/* Program to zero area of pixels that dont drain to a designated outlet */
/* pixel by calculating areas of pixels draining to the outlet pixel */
/*  using a recursive algorithm starting from the outlet pixel  */

/* Created by David G Tarboton  */
/* Modified by Shaleen Jain  3/28/95 */
/*  The modified code calculates the length of the longest stream.        */
/*  The code also computes the longest lengths of all sub basins in one   */
/*  go,so the area file and the ld file can be used to examine basin      */
/*  geomorphologic characteristics, such as, Hack's Law.                  */     

/*    Input is read from file ZAREA.IN which has the following structure. */
/*    Record 1: File name of pointer file input.                          */
/*    Record 2: File name of area file input (not used).                  */
/*    Record 3: File name of isolated area file output.                   */
/*    Subsequent Records: Row and Column of the outlet pixels isolate on. */
/*        More than one is possible to allow subnetworks to be computed   */
/*        first to prevent recursion depth or storage of two networks     */
/*        in the same file.                                               */


#include "lsm.h"
   char elfile[MAXLN];
   FILE *fout;
   float **ltp,**lengd,**gleng;
   double dist[9];
   int **area, **aa, tresh;
   void d2area();
void main()
 {
   char pfile[MAXLN],afile[MAXLN],lgfile[MAXLN],ltpfile[MAXLN],glfile[MAXLN];
   FILE  *fin;
   int i,j,row,col,err,ifp,jfp;
     double x=0., y=0.;
/* define directions */
   d1[1]=0; d1[2]= -1; d1[3]= -1; d1[4]= -1; d1[5]=0; d1[6]=1; d1[7]=1; d1[8]=1;
   d2[1]=1; d2[2]=1; d2[3]=0; d2[4]= -1; d2[5]= -1; d2[6]= -1; d2[7]=0; d2[8]=1;
/*  read in data  */
   fin=fopen("hack.in","r");
   fscanf(fin, "%s", pfile);
   fscanf(fin, "%s", afile); 
   fscanf(fin, "%s", lgfile);
   fscanf(fin, "%s", ltpfile);
   fscanf(fin, "%s", glfile);
   fscanf(fin, "%s", elfile);
   fscanf(fin, "%d", &tresh);   
   i=fscanf(fin,"%lf %lf",&x,&y);
   fclose(fin);
/*  read pointers */
/*    igy=IGY; 
    demread_(dir,pfile,&nx,&ny,&igy,&dx,&dy);*/

/*    read pointers  */
   err=gridread(pfile,(void ***)&dir,RPSHRDTYPE,&nx,&ny,&dx,&dy,
   bndbox,&csize,&ndv,&filetype);
/*   read areas  */
   err=gridread(afile,(void ***)&aa,RPINTDTYPE,&nx,&ny,&dx,&dy,
   bndbox,&csize,&ndv,&filetype);

  fprintf(stderr,"...now allocating memory...\n");
/* allocate memory to arrays */
   area= (int **)  matalloc(nx,ny,RPINTDTYPE);
   ltp= (float **) matalloc(nx,ny,RPFLTDTYPE);
   lengd= (float **) matalloc(nx,ny,RPFLTDTYPE);
   gleng= (float **) matalloc(nx,ny,RPFLTDTYPE);

  fprintf(stderr,"...allocation done...\n");
  fprintf(stderr,"%d ",tresh,"\n");

/*  Calculate Distances  */
    for(i=1; i<=8; i++){
     dist[i]=sqrt(d1[i]*d1[i]*dy*dy+d2[i]*d2[i]*dx*dx);
   }
/* open output file (ofile) */
/*   fout=fopen(ofile,"w");*/
/* get pixels to zero on */
/* There may be several with the first few judiciously chosen to */
/* restrict recursion depth   */
     /* initialize area array to 0 */
     for(i=0; i<ny; i++)
     {  for(j=0; j<nx; j++)
        {
          area[j][i]=0;
          ltp[j] [i]=0;
          lengd[j][i]=0;
          gleng[j][i]=0;
        }
     }
   if(x > 0. && y > 0.)   /*  Only compute area's for designated location  */
   {
		col= (int)floor((x-bndbox[0])/csize);
		row= (int)floor((bndbox[3]-y)/csize);
		if(row <0 || row > ny || col < 0 || col > nx)
		{
			printf("Given coords out of area\n");
			row=0; col=0;
		}
	}


     /* call drainage area subroutine for pixel to zero on  */
   fprintf(stderr,"..calcuating areas etc...\n");
  fout=fopen(elfile,"w");
     d2area(row,col,&ifp,&jfp);
/* write out to grid files */
/*   demwrite_(area,afile,&nx,&ny,&igy,&dx,&dy);   */
/*   demwrite_(lto,lfile,&nx,&ny,&igy,&dx,&dy);    */
/*   demwrite_(lengd,ldfile,&nx,&ny,&igy,&dx,&dy); */  
  fprintf(stderr,"..finished calcuating areas etc...\n");

err = gridwrite(lgfile,(void **)lengd,RPFLTDTYPE,nx,ny,dx,dy,
             bndbox,csize,-1.,filetype);
err = gridwrite(ltpfile,(void **)ltp,RPFLTDTYPE,nx,ny,dx,dy,
             bndbox,csize,-1.,filetype);
err = gridwrite(glfile,(void **)gleng,RPFLTDTYPE,nx,ny,dx,dy,
             bndbox,csize,-1.,filetype);

   fclose(fout);
   free(area[0]);
   free(area);
   free(aa[0]);
   free(aa);
   free(dir[0]);
   free(dir);
   free(ltp[0]);
   free(ltp);  
   free(lengd[0]);
   free(lengd);
   free(gleng[0]);
   free(gleng);
 }

/* function to compute area recursively */
void d2area(i,j,ifp,jfp)

int i,j,*ifp,*jfp;
{
    int in,jn,k,ifn,jfn;
    float ld;
    if(area[j][i]==0)
    {
      if(i!=0 && i!=ny-1 && j!=0 && j!=nx-1 && dir[j][i]!= -1)
      /* not on boundary  */
      {
        area[j][i]=1;
        *ifp = i;
        *jfp = j;
	
        for(k=1; k<=8; k++)
        {  in=i+d1[k];
           jn=j+d2[k];
	   /* test if neighbor drains towards cell excluding boundaries */
           if(dir[jn][in]>=0) 
            {
           if(aa[jn][in]>=tresh && (dir[jn][in]-k==4||dir[jn][in]-k==-4))
             {
                d2area(in,jn,&ifn,&jfn);
                area[j][i]=area[j][i]+area[jn][in];
				ld=lengd[jn][in]+dist[dir[jn][in]];
                ltp[j][i] += ltp[jn][in] + dist[dir[jn][in]];
		 if(ld > lengd[j][i])
		 {
                 lengd[j][i]=ld;
                 *ifp = ifn;
                 *jfp = jfn;
                 }
              }
           }

        }
	gleng[j][i] = sqrt((i-*ifp)*dy*(i-*ifp)*dy + 
                   (j-*jfp)*dx*(j-*jfp)*dx);
        fprintf(fout,"%d %d %d %d\n",i,j,*ifp,*jfp);
      }
    }
}