#include <stdio.h>
#include <string.h>
#include <fcntl.h>
#include <math.h>
#include "nrutil.h"

/* arguments
  1) root name of the input bil file.  If omitted, the program will prompt.
  2) (optional) factor for scaling up 16-bit input files.  default 1000
     If the raw data are in feet, a factor = 304.8 yields meters.
  3) (optional) maximum value by which an elevation will be incremented.
     default = 900.  It should be smaller than the factor (2)
  4) (optional) factor to scale up area / base to retain integer precision.
  5) max number of downslope iterations.  Save time by avoiding valley bottoms.
/* The input image files are:
 * a) 16 or 32-bit integer filled elevations, the file name is the only program prompt
 * Harvey Greenberg
 * this should be adaptable to 4-pointers.
 * The raw elevation (eray) is scaled up by 1000 (4-byte)
 * Flat areas are sloped.
 * Then the flowaccumulation is run.
 * Output files include:
 * a) [root]sloped.bil - the filled DEM, scaled up by 1000 and sloped
 * b) [root]cum.bil - cumulative area
 * The output files can be converted to grids with the imagegrid command
 * at the Arc: prompt.  First copy the bilw and hdr files to [root]out.bilw,
 * etc., first editing the NBITS, BANDROWBYTES, and TOTALROWBYTES fields of the
 * out, sloped, and cum files to reflect their 32-bit format.
 *
 * 3/14/94 pm Instead of cummulative area, write out cummulative flux
 * 3/15/94 Use sqrt(2) factor on diagonals to compute flows
 * 4/25/94 optional second argument is factor
 * 6/21/94 removed the creation of hmgdir file
 * 10/24/94 removed a break which messed up sloping; added diags, and
 *          changed maxdif from 500 to 900
 * 8/21/95 changed all my close statements from afile to fp, etc., see blw file
 * 10/2/95 throw out some dead code; implement malloc; change name
 * 10/6/95 throw out last MAXROW reference, add an output slope grid
 * 10/10/95 No need to open filled.bil, because we have been assuming 
 *          that everything is filled.
 * 11/02/95 Various simplifications and fixes to make it work on the Queets.
 * 1/23/96 Fooled with maxdif and output file names.
 *	   Also extensive checks so the program will work well
 *	   on 32-bit input files, e.g. clayoquot. (see line 314.)
 * 6/27/96 Fooled with cell counter.  Make factor floating point.
 * 7/31/96 fool with matrix declarations to satify new compiler
 *	   Initialize inray for 32-bit input.
 * 8/28/96 Keep track of floating-point cell area, while using 100 internally.
 *	   added fourth parameter
 * 9/10/96 Fix the scaling of flux.
 * 11/7/97 1.081. no dimunition of flux on the diagonal
 * 2/24/98 1.082  simply diagonal stuff by ignoring flux, just doing area.  May help rounding errors.
 *		  drop the trueedge crap. the flux factor should be enough.
 * 10/13/98 1.083 Comment out floating point. Max iterations
 * 10/19/98 1.084 handling of real sinks which are coded -1
 * 11/03/98 1.085 juiced up algorithm for traversing array
 * 07/01/04 1.087 use the scaling factor with 4-bit input.
 */
FILE *fp;

short **inray,**dirray,**slopray;
int **eray,**lray,*downk;
float **fray;

/* these declarations define the 8-pointer method */
int nn = 8, nj[] = {1,1,0,-1,-1,-1,0,1},
	    ni[] = {0,1,1,1,0,-1,-1,-1},
	    nb[] = {1,2,4,8,16,32,64,128};
float dfactor[] = {1,1.4142,1,1.4142,1,1.4142,1,1.4142};

float ftweight,cum,edge,factor,area,fluxscale = 1.;

int n,i,j,k,l,fd,fe,ff,nrows,ncols,dbits, length, thiscount,totalcount,
    maxe = -1,maxdrop, flag, w[8],zmin,zmax,mindif,dif,tweight,
    bigits,sinkcells,downz,fixes,delv,sinkc,ii,jj,ir,jr,itcounter,thiscountold,
    sinkk,sinkmax,pourk,pourmax,sinki[99],sinkj[99],pouri[99],pourj[99],
    sinkelv,pourelv,pourc,dlength,checksum,maxdif,dx,dy,dist,elv,dz,maxits=9999;
char root[80],afile[80],bfile[80],ffile[80],sfile[80],progname[80];
main (int argc, char *argv[]) {
 strcpy(progname,argv[0]);
 printf("%s Version 1.087: 07/01/2004\n",progname);
 fprintf(stderr,"%s Version 1.087: 07/01/2004\n",progname);
if(strcmp(progname,"sloparea") == 0) printf("%s is sloparea\n",progname);
if(strcmp(progname,"sloparea_nocheck") == 0) printf("%s is sloparea_nocheck\n",progname);
if(strcmp(progname,"whatever") == 0) printf("%s is whatever\n",progname);
 if(sizeof(int) != 4 || sizeof(short) != 2 || sizeof(float) != 4){
   fprintf(stderr,"Int, short and float must be 4,2, and 4, not %d, %d, and %d\n",
	sizeof(int),sizeof(short),sizeof(float));
   return 1;
 }
 if(argc == 1){
    fprintf(stderr,"Enter file root name: ");
    scanf("%s",root);
 }
 else
    strcpy(root,argv[1]);
 factor = 1000;
 if(argc >= 3){
   sscanf(argv[2],"%f",&factor);
   printf("Using vertical scale factor of %f\n",factor);
 }
 maxdif = factor - 2;
 if(argc >= 4){
   sscanf(argv[3],"%d",&maxdif);
   printf("Using maximum increment of %d\n",maxdif);
 }
 if(argc >= 5){
   sscanf(argv[4],"%f",&fluxscale);
   printf("Output cum file will be scaled up by %f to reduce integer truncation errors.\n",fluxscale);
 }
 if(argc >= 6){
   sscanf(argv[5],"%d",&maxits);
   printf("Processing will stop after %d iterations.\n",maxits);
 }
    strcpy(afile,root);
    strcat(afile,".hdr");
    if ((fp = fopen(afile,"r")) != NULL) {
      if(fscanf(fp,"%*s%*s%*s%*s%*s%d%*s%d%*s%*s%*s%d%",
            &nrows,&ncols,&dbits)) {
	/* keywords would be safer */
      printf("Elevation rows, cols & bits: %d, %d and %d\n",nrows,ncols,dbits);
      fprintf(stderr,"Elevation rows, cols & bits: %d, %d and %d\n",nrows,ncols,dbits);
      fclose(fp);
      }
    }
    else {
      fprintf(stderr,"Unable to open header file %s\n",afile);
      return 1; /* or query for info */
    }
    eray = imatrix(0,nrows,0,ncols);
    fray = fmatrix(0,nrows,0,ncols);
    inray = smatrix(0,nrows,0,ncols);
    dirray = smatrix(0,nrows,0,ncols);
    slopray = smatrix(0,nrows,0,ncols); /* slope to lowest neighbor. not really used */
    lray = imatrix(0,nrows,0,ncols);
    downk = ivector(1,nrows);
    
    if(dbits == 32){
	fprintf(stderr,"The input file is 32-bits, and is assumed to have no sinks, except for legitimate sinks, which are coded -1.\n");
    }
    else if(dbits != 16){
        fprintf(stderr,"The input file must be 16 or 32 bits; it is %d bits\n",dbits);
        printf("The input file must be 16 or 32 bits; it is %d bits\n",dbits);
        return 1;
    }

    strcpy(afile,root);
    strcat(afile,".bilw");
    area = 900;
    edge = 30.;
    if ((fp = fopen(afile,"r")) != NULL) {
      if(fscanf(fp,"%f",&edge)){
	area = (edge * edge);
      }
      else {
        fprintf(stderr,"Unable to read cell size, using 900\n");
      }
    fclose(fp); 
    }
    else{
      strcpy(afile,root);
      strcat(afile,".blw");	/* new 9660 naming convention */
      if ((fp = fopen(afile,"r")) != NULL) {
        if(fscanf(fp,"%f",&edge)){
	  area = (edge * edge);
          }  
          else {
            printf("Unable to read cell size, using 900\n");
          }  
        fclose(fp);
      }
    }

    printf("Cell area is %f square units\n",area);
     
    strcpy(afile,root);
    strcat(afile,".bil");
   if(dbits == 16){
    if ((fd = open(afile,O_RDONLY,0)) > 0) { 
      length = dbits * ncols / 8;
      dlength = ncols * 4;
      for(i=0;i < nrows;i++){
	if((n = read(fd,inray[i],length)) != length){
	  fprintf(stderr,"Failed to read line %i - %i\n",i,n);
	  printf("Failed to read line %i - %i\n",i,n);
	  return 1;
	}
      }
      close(fd);
    }
    else {
      fprintf(stderr,"Failed to open %s\n",afile);
      printf("Failed to open %s\n",afile);
      return 1;
    }

    for(i=0;i < nrows;i++){
      for(j=0;j < ncols;j++)
	if(inray[i][j] == -1) /* added 10/19/98 for legitimate sinks */
	  eray[i][j] = -1;
	else
	  eray[i][j] = inray[i][j] * factor;
    }
  }
  else if (dbits == 32) {
    if ((fd = open(afile,O_RDONLY,0)) > 0) { 
      length = ncols * 2;
      dlength = ncols * 4;
      for(i=0;i < nrows;i++){
	if((n = read(fd,eray[i],dlength)) != dlength){
	  fprintf(stderr,"Failed to read line %i - %i\n",i,n);
	  printf("Failed to read line %i - %i\n",i,n);
	  return 1;
	}
      }
      if(argc >= 3){  /* 7/1/2004: if scale factor is an argument,
		         use it on 32-bit input as well.   */
        for(i=0;i < nrows;i++){
          for(j=0;j < ncols;j++)
	    if(eray[i][j] != -1) /* for legitimate sinks */
		eray[i][j] = eray[i][j] * factor;
	}
      }

      for(i=0;i < nrows;i++){
	for(j=0;j < ncols;j++)
          inray[i][j] = 1; /* we use for a flag later */
      }
      close(fd);
    }
    else {
      fprintf(stderr,"Failed to open %s\n",afile);
      printf("Failed to open %s\n",afile);
      return 1;
    }

  } /* got our 32-bit eray, one way or another */
    strcpy(bfile,root);
    strcat(bfile,"sloped.bil");
    if ((fe = creat(bfile,0664)) == -1) { 
      fprintf(stderr,"Failed to creat %s\n",bfile);
      printf("Failed to creat %s\n",bfile);
      return 1;
    }
/* code any zero elevation on edge as -1, so we can spot 0-meter sinks  */
    for(j=0;j < ncols;j++){
      dirray[0][j] = nb[6]; /* pours off upwards*/
      if(eray[0][j] == 0) eray[0][j]--;
      dirray[nrows-1][j] = nb[2]; /* pours off down */
      if(eray[nrows-1][j] ==0) eray[nrows-1][j]--;
    }
    for(i=1;i < nrows-1;i++){
      dirray[i][0] = nb[4]; /* pours off left */
      if(eray[i][0] == 0) eray[i][0]--;
      dirray[i][ncols-1] = nb[0]; /* pours off right */
      if(eray[i][ncols-1] == 0) eray[i][ncols-1]--;
    }
 if (strcmp(progname,"sloparea") == 0){
    fprintf(stderr,"Propogate any -1 elevations in from the edge\n");
    bigits = 0;
    fixes = 999;
    while(fixes > 0){
      fixes = 0;
      bigits++;
      for(i=1;i < nrows-1;i++){
        for(j=1;j < ncols-1;j++){
          if (eray[i][j] == 0){
                for(k=0;k < nn;k++){
                  if(eray[i + ni[k]][j + nj[k]] == -1){
                    eray[i][j]--;
                    fixes++;
                    break;
                  }
                }   
          }       
        }
      }   
    }      /* propogating elevation change 0 -> -1 from edge */
    fprintf(stderr,"Done in %d iterations\n",bigits);
 } else if (strcmp(progname,"sloparea_nocheck") == 0){
    fprintf(stderr,"We are assuming all elevations <= 0 are sinks\n");
    for(i=0;i < nrows;i++){
      for(j=0;j < ncols;j++){
        if(eray[i][j] <= 0) eray[i][j] = -1;
      }
    }
 } else{
    fprintf(stderr,"The program must be called sloparea or sloparea_nocheck\n");
    return 1;
 }

/* CALCULATE ESRI-style directions.   */
    fprintf(stderr,"Calculating directions on filled DEM\n");
    for(i=1;i < nrows-1;i++){
      for(j=1;j < ncols-1;j++){
        elv = eray[i][j];
   if(elv == -1)
     dirray[i][j] = -1;
   else {
        maxdrop = -1;
        flag = 0;
        for(k=0;k < nn;k++){
          dz = elv - eray[i + ni[k]][j + nj[k]];

          if(dz > maxdrop) {  		/* was >= Was that right? */
            maxdrop = dz;
	    pourk = k;
            /* flag = nb[k]; LAST ONE FOUND HAS PRIORITY  */
          }
        }
        for(k=0;k < nn;k++){
          dz = elv - eray[i + ni[k]][j + nj[k]];
          if(dz == maxdrop) flag = nb[k];
           /* flag = flag | nb[k]; would allow multi-pointing, but it's bad */ 
        }
        if(maxdrop == -1){
	  dirray[i][j] = /* - */   1;
	  fprintf(stderr,"Danger Will Robinson, sink at %d %d %d\n",i,j,eray[i][j]);
	  printf("Danger Will Robinson, sink at %d %d %d\n",i,j,eray[i][j]);
	  if(elv == 0){
	    printf("Fixing a zero point\n");
	    k = 0;
	    if(eray[i-1][j-1] > 0){
		elv += eray[i-1][j-1];
		k++;
	    }
	    if(eray[i-1][j] > 0){
		elv += eray[i-1][j];
		k++;
	    }
	    if(eray[i-1][j+1] > 0){
		elv += eray[i-1][j+1];
		k++;
	    }
	    if(eray[i][j-1] > 0){
		elv += eray[i][j-1];
		k++;
	    }
	    if(eray[i][j+1] > 0){
		elv += eray[i][j+1];
		k++;
	    }
	    if(eray[i+1][j-1] > 0){
		elv += eray[i+1][j-1];
		k++;
	    }
	    if(eray[i+1][j] > 0){
		elv += eray[i+1][j];
		k++;
	    }
	    if(eray[i+1][j+1] > 0){
		elv += eray[i+1][j+1];
		k++;
	    }
	    if(k != 8)
	      fprintf(stderr,"Hey Will, all %d neighbors are non-zero\n",k);
	      printf("Hey Will, all %d neighbors are non-zero\n",k);
	    eray[i][j] = (int)( (float)elv / (float)k + 1 );
	/* ignore the remote posibility that the dirray should need redoing */
	  }
	  else{
	    fprintf(stderr,"Sink elevation is %d; bailing out\n",elv);
	    printf("Sink elevation is %d; bailing out\n",elv);
	    return 1;
	  }
	}
        else if(maxdrop == 0) dirray[i][j] = 0;
          else {
	    dirray[i][j] = flag;
	    slopray[i][j] = maxdrop / ( edge * dfactor[pourk]);
          }
      }
   } /* elv != 0 -1 */
    }
/* At this point, a good cell has a dirray value to indicate one or more pour points.
 * A  sink-point is -1.
 * A flat spot, possibly a sink, is 0
 * Cells pouring into a sink are not differentiated
 */
/* BEGIN HEART OF SLOPER */
    fprintf(stderr,"sloping on the filled array\n");
    bigits = 0;
    fixes = 999;
    while(fixes > 0){
      fixes = 0;
      bigits++;
      for(i=1;i < nrows-1;i++){
        for(j=1;j < ncols-1;j++){
         if(eray[i][j] > 0){ /* was != before 11/4/95       */  
          if(dirray[i][j] == 0){ /* same elev as a neighbor; doesn't pour */
            maxdrop = bigits + 1; /* was maxdif */
            for(k=0;k < nn;k++){ /* look for neighbor which pours, but the
*		desired neighbor may not be designated a pourer until a later
*		iteration.  First, make sure that the neighbor is
*		not pouring back this way.  Also, do not allow maxdrop
*		greater than bigits.  This slows things done a little,
*		but is safer. */
/*
printf("pours? at %d %d on k = %d, dir = %d\n",i + ni[k],j + nj[k],k,dirray[i +
ni[k]][j + nj[k]]);
*/
              if(dirray[i + ni[k]][j + nj[k]] > 0){
	       if((nb[k] != (16 * dirray[i + ni[k]][j + nj[k]])) &&
	       ((nb[k] * 16) != dirray[i + ni[k]][j + nj[k]])){ /* not backandforth*/
                delv = eray[i + ni[k]][j + nj[k]] - eray[i][j];
                if(delv < maxdrop){ /* if neighbors are equal, keep the last */
                  maxdrop = delv;
                  pourk = k;
                }
               }   
              }   
            }   
            if(maxdrop <= bigits){
                  dirray[i][j] = nb[pourk];
                  eray[i][j] = eray[i + ni[pourk]][j + nj[pourk]] + 1;
                  slopray[i][j] = 1.0; /* i.e. .001, that's close enough */
                  fixes++;
/* but what if you have raised this up to the level another cell pouring in?
*  This is a problem only with our 32-bit files. */
		  for(k=0;k < nn;k++){ /* for each neighbor */
		    if(eray[i + ni[k]][j + nj[k]] <= eray[i][j]){ /* lower neigbor */
		      if((nb[k] == (16 * dirray[i + ni[k]][j + nj[k]])) ||
		      ((nb[k] * 16) == dirray[i + ni[k]][j + nj[k]])){
			fprintf(stderr,"Special adjustment at %d %d from %d to %d\n",
			  i,j,eray[i + ni[k]][j + nj[k]],eray[i][j]);
			printf("Special adjustment at %d %d from %d to %d\n",
			  i,j,eray[i + ni[k]][j + nj[k]],eray[i][j]);
		/*eray[i + ni[k]][j + nj[k]] = eray[i][j] + 1; */
/* don't try to fix the other cell; just flag it as not done */
			dirray[i + ni[k]][j + nj[k]] = 0;
		      }
		    }
		  }
              /*  break;  commented out 10/24/94   */
            }
          } /* for each flat spot */
         }      /* elevation > 0 */
        } 
      }  
        fprintf(stderr,"finishes iteration %d, fixing %d\n",bigits,fixes);
    }
 
/* END HEART OF SLOPER */

/*  Is this right?  if(bigits != 1 && fixes != 0){  */
/* write out fixed DEM (scaled by 1000)  */

    if ((fe = open(bfile,O_WRONLY,0)) <= 0) {
          printf("Failed to open filled bil file to write %d",fe);
          fprintf(stderr,"Failed to open filled bil file to write %d",fe);
          return 1;
    }
      for(i=0;i < nrows;i++){
        if((n = write(fe,eray[i],dlength)) != dlength){
          printf("Failed to write line %i\n",n);
          fprintf(stderr,"Failed to write line %i\n",n);
          return 1;
       }
      } 
      close(fe);
/*  } */
    printf("Wrote out fixed elevation file: %s\n",bfile);

/* write out slope grid (slope is to lowest neighbor)  */
/* not now
    strcpy(sfile,root);
    strcat(sfile,"slop.bil");
    if ((fe = creat(sfile,0664)) == -1) { 
      printf("Failed to creat %s\n",sfile);
      return 1;
    }
    if ((fe = open(sfile,O_WRONLY,0)) <= 0) {
          printf("Failed to open slope bil file %s to write %d",sfile,fe);
          return 1;
    }
      for(i=0;i < nrows;i++){
        if((n = write(fe,slopray[i],length)) != length){
          printf("Failed to write line %i\n",n);
          return 1;
       }
      } 
      close(fe);
*/

/* water always flows off the edge - not calculating to edge for now */
/* seed each cell with its own area 
 * We are throwing out the separate flux tally.  The B unit length is always
* the edge length as test show.  The only in which this might disturb
 * people is that flux is calculated exiting the cell.
 * use fray, the floating array.
*/
      for(i=0;i < nrows;i++){
        for(j=0;j < ncols;j++){
          fray[i][j] = area;
        }
      }  
/* reuse inray to flag processed cells          */
/* flag edges as done */
  for(i=0;i < nrows;i++){
    inray[i][0] = -999;
    inray[i][ncols -1] = -999;
  }  
  for(j=0;j < ncols;j++){
     inray[0][j] = -999;
     inray[nrows -1][j] = -999;
  }

    totalcount = (nrows - 2) * (ncols - 2);
    fprintf(stderr,"Calculating flow for %d cells\n",totalcount);
    bigits = 0;
    itcounter = 0;
    thiscountold = 9999999;

    thiscount = 99;
    while(thiscount > 0){
      bigits++;
      if (((float)thiscount / (float)thiscountold) < 0.9) itcounter++;
/*fprintf(stderr,"bigits=%i,itcounter=%i,thiscount=%i thiscountold=%i : %f v. 0.9",bigits,itcounter,thiscount,thiscountold,(float)thiscount/(float)thiscountold);
*/
      thiscountold = thiscount;
      thiscount = 0;

      for(ir=1;ir < nrows-1;ir++){
	if(itcounter % 2)
	  i = ir;
	else
	  i = nrows-1-ir;
/* if(ir < 4 || ir > 297) printf("ir=%i, i = %i\n",ir,i);
*/
	for(jr=1;jr < ncols-1;jr++){
	  if(itcounter % 4 < 2)
	    j = jr;
	  else
	    j = ncols-1-jr;
/*if((ir < 4 || ir > 297) && (jr < 4 || jr > 297))printf(" jr=%i, j = %i\n",jr,j);
*/
	  elv = eray[i][j];
          for(k=0;k < nn;k++){
	    if(inray[i][j] < 0)
		break; /* this cell done  or edge of quad */ 
	    if((eray[i + ni[k]][j + nj[k]] > elv) /*adj cells can have same elv*/
		&& (inray[i + ni[k]][j + nj[k]] > 0)){
		break; /* not ready for this cell */
	    }
	  }
	  if(k == 8){  /* no unaccounted-for uphill flows */
	    thiscount++;
	    tweight = 0;
	    cum = fray[i][j];
	    maxdrop = -9999888;
	    pourk = -1;
            for(k=0;k < nn;k++){
	      w[k] = elv - eray[i+ni[k]][j+nj[k]];
	      if(w[k] >0){ /* altering Berkeley algorithm re diagonals */
		if(((k&1) == 1) && (w[k] != 1)) /* and nn == 8 */
		  w[k] *= .7071;
		tweight += w[k];
		if(w[k] > maxdrop){
		  maxdrop = w[k];
		  pourk = k;
		}
	      }
	    }
	    ftweight = (float)tweight;
	    for(k=0;k < nn;k++){
	      if(w[k] >0){
		fray[i+ni[k]][j+nj[k]] +=  (float)w[k] * cum / ftweight;
	      }
	    }
	    if(pourk == -1){
	      if(elv >= 0)printf("WARNING: no pour from cell %d %d\n",i,j);
	      if(elv >= 0)fprintf(stderr,"WARNING: no pour from cell %d %d\n",i,j);
	      dirray[i][j] = 0;
	    }
	    else
	      dirray[i][j] = nb[pourk];
	    inray[i][j] = -bigits; /* flag this cell done */
	  } /* if */
	}   /* j  */
      }     /* i  */
      totalcount = totalcount - thiscount;
      fprintf(stderr," iteration %d: %d more done, %d to go\n",bigits,thiscount,totalcount);
      if(bigits > maxits){
	printf("Max iterations exceeded; leaving with %d unprocessed cells\n",totalcount);
	fprintf(stderr,"Max iterations exceeded; leaving with %d unprocessed cells\n",totalcount);
	break;
      }
    }       /* while */
/* New innovation: write the results as a float file
    strcpy(bfile,root);
    strcat(bfile,"area");
    printf("Writing out unscaled floating-point cummulative area to %s\n",bfile);
    if ((fe = creat(bfile,0664)) == -1) { 
      printf("Failed to creat %s\n",bfile);
      fprintf(stderr,"Failed to creat %s\n",bfile);
      return 1;
    }

    if ((fe = open(bfile,O_WRONLY,0)) <= 0) { 
	  printf("Failed to open cum bil file to write %d",fd);
	  fprintf(stderr,"Failed to open cum bil file to write %d",fd);
	  return 1;
    }
      for(i=0;i < nrows;i++){
	if((n = write(fe,fray[i],dlength)) != dlength){
	  printf("Failed to write line %i\n",n);
	  fprintf(stderr,"Failed to write line %i\n",n);
	  return 1;
       }
      }
      close(fe);
*/
/* The traditional area/B as long integers */
      for(i=0;i < nrows;i++){
        for(j=0;j < ncols;j++)
          lray[i][j] = (fray[i][j] / edge) * fluxscale + 0.5 ;
      }

    strcpy(bfile,root);
    strcat(bfile,"cum.bil");
    printf("Writing out cummulative flux to %s\n",bfile);
    if(fluxscale != 1.0)
	printf("Output cum file will be scaled up by %f\n",fluxscale);
    if ((fe = creat(bfile,0664)) == -1) { 
      printf("Failed to creat %s\n",bfile);
      fprintf(stderr,"Failed to creat %s\n",bfile);
      return 1;
    }

    if ((fe = open(bfile,O_WRONLY,0)) <= 0) { 
	  printf("Failed to open cum bil file to write %d",fd);
	  fprintf(stderr,"Failed to open cum bil file to write %d",fd);
	  return 1;
    }
      for(i=0;i < nrows;i++){
	if((n = write(fe,lray[i],dlength)) != dlength){
	  printf("Failed to write line %i\n",n);
	  fprintf(stderr,"Failed to write line %i\n",n);
	  return 1;
       }
      }
      close(fe);

/* write out the direction array
    strcpy(bfile,root);
    strcat(bfile,"hmgdir.bil");
    if ((fe = creat(bfile,0664)) == -1) {
      printf("Failed to creat %s\n",bfile);
      return 1;
    }
      
    if ((fe = open(bfile,O_WRONLY,0)) <= 0) {
          printf("Failed to open hmgdir bil file to write %d",fd);
          return 1;
    }
      for(i=0;i < nrows;i++){
        if((n = write(fe,dirray[i],length)) != length){
          printf("Failed to write line %i\n",n);
          return 1;
       }
      }   
      close(fe);
*/
return 0;
}