/* zma2.c */ 
/* Program to compute area contributing to each Pixel in DEM */
/* for cell outflow based on angles*/

/* Input from stdin with following structure:           */
/* Record 1 - File name of angle matrix                   */
/* Record 2 - File name of area file to output                 */
/* Subsequent records are pixel co-ordinates to initiate calculation */
/* 0 0 calculates the whole region */

#include <stdio.h>
#include <math.h>
#define IGX 1401
#define IGY 1120
#define MAXLN 80
#define PI 3.14159265359
   int nx,ny,igy;
   float dx,dy;
   int dd1[8],dd2[8],*d1,*d2;
   float ang[IGX][IGY],aref[10], area[IGX][IGY];

/* Declare program functions */

        float 
        prop();
        void dparea(),demwrite_();
 
main()
{
   char pfile[MAXLN],afile[MAXLN];
   FILE   *fin;
   int i,j,k,in,jn,irz,icz;
/* Define directions */
d1=dd1-1; d2=dd2-1;
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 */
printf("Input direction (angle) file name (must be .a4 or .r4)\n");
scanf("%s", pfile);
printf("Input file name for output areas (must be .a4 or .r4)\n");
scanf("%s", afile);

/* read angles */
igy = IGY;
demread_(ang,pfile,&nx,&ny,&igy,&dx,&dy);
/* Calculate angles for interpolation */                                      
aref[0]= -atan2(dy,dx);
aref[1]=0.;
aref[2]= -aref[0];
aref[3]=0.5*PI;
aref[4]=PI-aref[2];
aref[5]=PI;
aref[6]=PI+aref[2];
aref[7]=1.5*PI;
aref[8]=2.*PI-aref[2];
aref[9]=2.*PI;
/* initialize area array to 0 */
    for(i=0; i<ny; i++)
       {
       for(j=0; j<nx; j++)
         {
         if(i!=0 && i!=ny-1 && j!=0 && j!=nx-1 && ang[j][i]>= -0.5)
            area[j][i]=0;
         else area[j][i]= -1;
         }
       }
   printf("Input row and column to initiate calculation (outlet)\n");
   printf("0 0 calculates the whole file\n");
   i=scanf("%d %d",&irz,&icz);
   while(i == 2 && irz > 0 && icz > 0)
   {
/* call drainage area subroutine for pixel to zero on  */
     irz=irz-1;   /* decrease index for c indexing starting from 0 */
     icz=icz-1;
     dparea(irz,icz);
     printf("Input another row and column to initiate calculation (outlet)\n");
     printf("[return] to end\n");
     i=scanf("%d %d",&irz,&icz);
   }
  if(irz == 0){

/* Call drainage area subroutine for each cell           */
/* working from the middle out to avoid deep recursions  */ 

/* "lower block" */
for(i=ny/2; i<ny-1; i++)
   {for(j=nx/2; j<nx-1; j++)
      dparea(i,j);
    for(j=nx/2-1; j>=1; j--)
      dparea(i,j);
   }
 /* "upper block" */
for(i=ny/2-1; i>=1; i--)
   {for(j=nx/2; j<nx-1; j++)
      dparea(i,j);
    for(j=nx/2-1; j>=1; j--)
      dparea(i,j);
   }
}
/*  Write out areas  */
   demwrite_(area,afile,&nx,&ny,&igy,&dx,&dy);
}  

/* function to compute partial areas recursively */

void dparea(i,j)
int i,j;
{
   int in,jn,k,kback,con=0;
    /* con is a flag that signifies possible contaminatin of area
         due to edge effects  */
   float p;
   if(area[j][i]==0.0)    /* i.e., if the cell hasn't been looked at yet */
      {
        if(i!=0 && i!=ny-1 && j!=0 && j!=nx-1 && ang[j][i] >= -0.5)
       /* i.e. the cell isn't outside domain */
         {
           area[j][i] = 1.0;
           for(k=1; k<=8; k++)
           {
             in = i + d1[k];
             jn = j + d2[k];
/* for each neighboring cell, find the proportion of of  outflow */
/* draining back to the cell in question */
             kback = k + 4;
             if(kback>8)kback = kback - 8;
             p=prop(ang[jn][in],kback);
             if(p>0.)
                {
                 dparea(in,jn);
                 if(area[jn][in] < -0.5)con = -1;
                 area[j][i] = area[j][i] + area[jn][in]*p;
                }
             if(ang[jn][in] < -0.5)con = -1;         
            }
           if(con == -1)area[j][i]= -1.;
         }
     }
} 

float prop(a,k)
float a;
int k;
{
float p=0.;
if(k == 1 && a > PI)a=a-2.*PI;
if(a > aref[k-1] && a < aref[k+1]){
  if( a > aref[k])
       p=(aref[k+1]-a)/(aref[k+1]-aref[k]);
  else 
      p=(a-aref[k-1])/(aref[k]-aref[k-1]);
  }
return(p);
}