/* ma.c */ 
/* Program to compute area contributing to each Pixel in DEM */
/* for multiple direction cell outflow */

/* Program modification of D.G. Tarboton single-flow-path  */
/* algorithm by T.H. Jackson, Aug. 1991 @ UWRL, Logan, Utah */

/* Input from work_files.in with following structure:           */
/* Record 1 - File name of adjusted elevation data          */
/* Record 2 - File name of pointer matrix                   */
/* Record 3 - File name of area outlet file                 */
/* Record 4 - File name  not used                           */

#include <stdio.h>
#include <math.h>
#define IGX 801
#define IGY 501
#define MAXLN 80
#define Nocells IGX*IGY 
   int nx,ny,igy;
   float dx,dy;
   short int elev[IGX][IGY],dir[IGX][IGY];
   int d1[9],d2[9];
   float slope[Nocells],prop[IGX][IGY][9],w[9],distance[9], area[IGX][IGY];
main()
{
   char efile[MAXLN],pfile[MAXLN],afile[MAXLN];
   FILE   *fin;
   int i,j,k,in,jn;
   float S_pos_tot;
/* 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;

/* Multiple Outlet Weighting Factors from Quinn et al (1991).  */
w[1]=0.5; w[2]=.354; w[3]=0.5; w[4]=.354;
w[5]=0.5; w[6]=.354; w[7]=0.5; w[8]=.354;

/* Read in data */
fin=fopen("ma.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);
demread_(elev,efile,&nx,&ny,&igy,&dx,&dy);
/* Calculate distances between cells */                                      
for(k=1; k<=8; k++)
   distance[k]=sqrt(d1[k]*d1[k]*dy*dy + d2[k]*d2[k]*dx*dx);

/* Compute proportions of every cell's flow to each of its neighbors */

for(i=0; i<ny; i++)
    {for(j=0; j<nx; j++)
       {for(k=1; k<=8; k++)       /* initialize prop[j][i][k] */
           prop[j][i][k] = 0.0;
        }
    }
for(i=1; i< ny-1; i++)
    {for(j=1; j< nx-1; j++)
       {S_pos_tot = 0.0;
        for(k=1; k<=8; k++)
           {in = i + d1[k];
            jn = j + d2[k];  
            slope[k]=(elev[j][i] - elev[jn][in])/distance[k];
            if(slope[k]>=0.0)
               S_pos_tot = S_pos_tot + slope[k]*w[k];
           }
        if(S_pos_tot>0.0)   
             {for(k=1; k<=8; k++)
                 {if(slope[k]>=0.0)
                    prop[j][i][k] = w[k]*slope[k]/S_pos_tot;
                  else
                    prop[j][i][k] = 0.0;
                  }
              }
        else 
        /* The cell is in a flat area so use the single direction */
        /* pointer read calculated in a previous program */
         prop[j][i][dir[j][i]] = 1.0;       
       }
    }
/* initialize area array to 0 */
    for(i=0; i<ny; i++)
       {for(j=0; j<nx; j++)
          area[j][i] = 0.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 */

dparea(i,j)
int i,j;
{
   int in,jn,k,kback,con=0;
   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 && dir[j][i] != -1)
       /* i.e. the cell isn't on the cell grid boundary */
         {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;
             if(prop[jn][in][kback]>0.)
                {dparea(in,jn);
                 if(area[jn][in] < -0.5) con = -1;
                 area[j][i] = area[j][i] + area[jn][in]*prop[jn][in][kback];
                }          
             if(dir[jn][in] < 0)con = -1;
            }
         if(con == -1)area[j][i]= -1;
         }
     }
}