/* Program to compute area contributing to each pixel from pointers using  */
/* recursive algorithm  */

/* Created by David G Tarboton  */
                                                                           
/*    Input is read from file branch.in which has the following structure. */
/*    Record 1: File name of raw elevation data file (not used).           */
/*    Record 2: File name of pointer file input.                           */
/*    Record 3: File name of area file output.                             */
/*    Record 4: File name of adjusted elevation data file (not used).      */
/*    Subsequent records: (not used)                                       */
/*    Modified on 4/18/95 to output -1 when area in not unambigious        */
/*    due to edge effect.  */

#include <stdio.h> 
#define IGX 1401
#define IGY 1120       
#define MAXLN 80
   int nx,ny,igy;
   float dx,dy;  
   short int dir[IGX][IGY];
   int area[IGX][IGY],d1[9],d2[9];
main()
 {                
     
   char efile[MAXLN],pfile[MAXLN],afile[MAXLN]; 
   FILE   *fin;
   int i,j;      
/* 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("branch.in","r");
   fscanf(fin, "%s", efile);     
   fscanf(fin, "%s", pfile);     
   fscanf(fin, "%s", afile);     
/*  read pointers */
   igy=IGY;
   demread_(dir,pfile,&nx,&ny,&igy,&dx,&dy);
/* initialize area array to 0, -1 on boundary */
   for(i=0; i<ny; i++)
   {
      for(j=0; j<nx; j++)
        {
        if(i!=0 && i!=ny-1 && j!=0 && j!=nx-1 && dir[j][i]!= -1)
           area[j][i]=0;
        else area[j][i]= -1;
        }
   }  
/* call drainage area subroutine for each area */
/* work from middle outwards to avoid deep recursions  */
   for(i=ny/2; i<ny-1; i++)                    
   {  for(j=nx/2; j<nx-1; j++)
        darea(i,j);
      for(j=nx/2-1; j>=1; j--)
        darea(i,j);
   }  
   for(i=ny/2-1; i>=1; i--)                    
   {  for(j=nx/2; j<nx-1; j++)
        darea(i,j);
      for(j=nx/2-1; j>=1; j--)
        darea(i,j);
   }  
/* write out areas */
   demwrite_(area,afile,&nx,&ny,&igy,&dx,&dy);
 }     

/* function to compute area recursively */
darea(i,j)
  int i,j;
  {                                        
    int in,jn,k,con=0;
      /* con is a flag that signifies possible contaminatin of area
         due to edge effects  */
    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;                               
        for(k=1; k<=8; k++)
        {  in=i+d1[k];
           jn=j+d2[k];
/* test if neighbor drains towards cell excluding boundaryies */
           if(dir[jn][in]>=0 && (dir[jn][in]-k==4||dir[jn][in]-k==-4))
             {
                darea(in,jn);
                if(area[jn][in] < 0)con = -1;
                area[j][i]=area[j][i]+area[jn][in];
             }
           if(dir[jn][in] < 0)con = -1;
        }
        if(con == -1)area[j][i]= -1;
      }
    }
  }