/* 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  */

/*    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.                                               */

/* This program uses Input/Output routines from CUTIL.C  */

#include <stdio.h> 
#define IGX 1401
#define IGY 1401       
#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,irz,icz;      
/* 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("zarea.in","r");
   fscanf(fin, "%s", pfile);     
   fscanf(fin, "%s", efile);  /* this file name not used */   
   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;
        }
   }  
/* get pixels to zero on */
/* There may be several with the first few judicuously chosen to */
/* restrict recursion depth   */
   i=fscanf(fin,"%d %d",&irz,&icz); 
   while(i == 2)
   {
/* call drainage area subroutine for pixel to zero on  */
     irz=irz-1;   /* decrease index for c indexing starting from 0 */
     icz=icz-1;   
     darea(irz,icz);
     i=fscanf(fin,"%d %d",&irz,&icz); 
   }
/* 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;
      }
    }
  }