#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <math.h>

#ifndef _E_RADIUS
#define E_RADIUS 6371.0         /* average radius of the earth */
#endif

#ifndef _PI
#define PI 3.1415
#endif


double get_dist(double,double,double,double);

main(int argc, char *argv[]) {
/***********************************************************************
  calc_area.c           Keith Cherkauer             November 25, 1998

  This program computes the area being routed using the basin, or
  sub-basin flow direction file, and the drainage fraction file.  

***********************************************************************/

  FILE *fdirec, *ffract;

  int    i, row, col, nrows, ncols;
  int    direc;
  double cell_lat, cell_lng;
  double ll_lat, ll_lng;
  double dlat, tmplat;
  double cellsize;
  double fract;
  double tmpsum, areasum;
  double routarea, tmprout;
  double NODATA;
  double ddist;

  if(argc!=3) {
    fprintf(stderr,"Usage: %s <direc file> <fract file>\n",argv[0]);
    exit(0);
  }

  if((fdirec=fopen(argv[1],"r"))==NULL) {
    fprintf(stderr,"ERROR: Unable to open flow direction file %s\n",argv[1]);
    exit(0);
  }

  if((ffract=fopen(argv[2],"r"))==NULL) {
    fprintf(stderr,"ERROR: Unable to open drainage fraction file %s\n",argv[2]);
    exit(0);
  }

  /** Read Header **/
  fscanf(fdirec,"%*s %i", &ncols);
  fscanf(fdirec,"%*s %i", &nrows);
  fscanf(fdirec,"%*s %lf",&ll_lng);
  fscanf(fdirec,"%*s %lf",&ll_lat);
  fscanf(fdirec,"%*s %lf",&cellsize);
  fscanf(fdirec,"%*s %lf",&NODATA);

  /** Process Direction Number **/
  dlat  = cellsize / 10.;
  ddist = get_dist(0.,0.,dlat,0.);

  areasum  = 0.;
  routarea = 0;
  for(row=0;row<nrows;row++) {
    for(col=0;col<ncols;col++) {
      cell_lat = ((double)(nrows-row)-0.5)*cellsize+ll_lat;
      cell_lng = ((double)col+0.5)*cellsize+ll_lng;

      fscanf(fdirec,"%i",&direc);
      fscanf(ffract,"%lf",&fract);

      if(direc!=NODATA) {
	tmpsum = 0.;
	tmplat = cell_lat-0.5*cellsize + 0.5*dlat;
	for(i=0;i<10;i++) {
	  tmpsum += get_dist(tmplat,cell_lng-0.5*cellsize,
			     tmplat,cell_lng+0.5*cellsize) 
	    * ddist;
	  tmplat += dlat;
	}
	areasum += tmpsum * fract;

	routarea +=  E_RADIUS*E_RADIUS*fabs(cellsize)*PI/180*
	  fabs(sin((cell_lat-cellsize/2.0)*PI/180)-
	       sin((cell_lat+cellsize/2.0)*PI/180)) * fract;

      }
    }
  }

  fprintf(stdout,"Basin Area = %lf km^2.\n",areasum);
  fprintf(stdout,"Routing Model Basin Area = %lf km^2.\n",routarea);

  fclose(fdirec);
  fclose(ffract);

}

#undef E_RADIUS
#undef PI