/* interpolation of satellite data to a regular grid */

/* this is the modified version created by S.Brodzik from the original used by
   U. Romatschke for interpolating PR TRMM data. 
   This is a faster and optimized code that use 2D matrices!!
*/



#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/types.h>

#define DEGREES 0.05
#define MAXSIZE 20
#define NODATA -999.0
#define PI 3.141592653589793
#define DEG2RAD (PI/180.0)
#define RAD2DEG 1/(PI/180.0)
#define MINDBZ 0.0

/* function prototypes. */

double degrees2distance(double gridLatitude, double latitude, 
			double gridLongitude, double longitude);
double cressman(double d2gp, float RADIUS);
void searcher(double distance[], int indices[], int indices2[], 
	      int *count, int *count2, float RADIUS);
void searcher2(double *latitude_2d, double *longitude_2d, 
	       double *gridLatitude, double *gridLongitude,  
	       double *reflectivity_2d, int indices2[], 
	       int *count2, int size2D, int i);

int main(int argv, char** argc)  {

  /* declaration of variables */
	
  int i, j;
  int count = 0, count2 = 0;
  FILE *fpLatitude, *fpLongitude;
  FILE *fpReflectivity;
  FILE *fpGridLatitude, *fpGridLongitude;
  FILE *fpNoSwath;
  FILE *fpNumSwaths, *fpGridRows, *fpGridColumns;
  FILE *fpGridReflectivity; 
  float RADIUS;
  int LEVELS, RAYS;
  int level, swath, ray;
  int size2D, size3D;
  int index, gridpoint;
  int *numSwathsI, *gridRowsI, *gridColumnsI;
  int numSwaths, gridRows, gridColumns;
  int *indices, *indices2;
  double *latitude, *latitude_2d;
  double *longitude, *longitude_2d;
  double *reflectivity, *reflectivity_2d;
  double *gridLatitude, *gridLongitude;
  double *noSwath;
  double *distance, *filteredDistance;
  double *weighting, *differences, *products; 
  double *interpolation;
  double interpolatedvalue, total, neighbour, proximity;
  double MAX;
  char datestring[50];
  char fname[128];

  /* check for correct number of input arguments */

  if (argv != 5) {
    fprintf(stdout, "Useage: %s [TRMM.XXXXX.YYYYMMDD.HHMMSS.] [RADIUS] [NUMLEVELS] [NUMRAYS]\n", argc[0]);
    exit(1);
  }

  /* get data from input parameters */

  strcpy(datestring,argc[1]);
  RADIUS = atof(argc[2]);
  LEVELS = atoi(argc[3]);
  RAYS = atoi(argc[4]);
  fprintf(stdout, "datestring=%s RADIUS=%f LEVELS=%d RAYS=%d\n", 
	  datestring, RADIUS, LEVELS, RAYS);

  /* open input files */

  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "latitude.txt");
  if ((fpLatitude = fopen(fname, "r")) == NULL)  {
    fprintf(stdout, "Problem opening input file %s\n", fname);
    exit(1);
  }
	
  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "longitude.txt");
  if ((fpLongitude = fopen(fname, "r")) == NULL)  {
    fprintf(stdout, "Problem opening input file %s\n", fname);
    exit(1);
  }

  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "reflectivity.txt");
  if ((fpReflectivity = fopen(fname, "r")) == NULL)  {
    fprintf(stdout, "Problem opening input file %s\n", fname);
    exit(1);
  }

  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "gridLatitude.txt");
  if ((fpGridLatitude = fopen(fname, "r")) == NULL)  {
    fprintf(stdout, "Problem opening input file %s\n", fname);
    exit(1);
  }

  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "gridLongitude.txt");
  if ((fpGridLongitude = fopen(fname, "r")) == NULL)  {
    fprintf(stdout, "Problem opening input file %s\n", fname);
    exit(1);
  }

  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "noSwath.txt");
  if ((fpNoSwath = fopen(fname, "r")) == NULL)  {
    fprintf(stdout, "Problem opening input file %s\n", fname);
    exit(1);
  }

  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "numSwaths.txt");
  if ((fpNumSwaths = fopen(fname, "r")) == NULL)  {
    fprintf(stdout, "Problem opening input file %s\n", fname);
    exit(1);
  }

  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "gridRows.txt");
  if ((fpGridRows = fopen(fname, "r")) == NULL)  {
    fprintf(stdout, "Problem opening input file %s\n", fname);
    exit(1);
  }

  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "gridColumns.txt");
  if ((fpGridColumns = fopen(fname, "r")) == NULL)  {
    fprintf(stdout, "Problem opening input file %s\n", fname);
    exit(1);
  }

  /* open output file */
  fname[0]=0;
  strcat(fname, datestring);
  strcat(fname, "gridReflectivity3D3L.txt");
  if ((fpGridReflectivity = fopen(fname, "w")) == NULL)  {
    fprintf(stdout, "Problem opening output file %s\n", fname);
    exit(1);
  }

  /* allocate space for input arrays */
  numSwathsI = (int *)malloc(sizeof(int));
  fscanf(fpNumSwaths, "%d", numSwathsI);
  gridRowsI = (int *)malloc(sizeof(int));
  fscanf(fpGridRows, "%d", gridRowsI);
  gridColumnsI = (int *)malloc(sizeof(int));
  fscanf(fpGridColumns, "%d", gridColumnsI);

  numSwaths = *numSwathsI;
  gridRows = *gridRowsI;
  gridColumns = *gridColumnsI;
  fprintf(stdout,"gridRows=%d and gridColumns=%d\n", gridRows,gridColumns);
  fprintf(stdout,"numRays =%d and numSwaths  =%d\n", RAYS,numSwaths);

  size3D = (LEVELS*numSwaths*RAYS);
  size2D = (numSwaths*RAYS);
  fprintf(stdout,"size3D=%ld and size2D=%ld\n",size3D,size2D);
 
  latitude = (double *)malloc(size3D*sizeof(double));
  latitude_2d = (double *)malloc(size2D*sizeof(double));
  longitude = (double *)malloc(size3D*sizeof(double));
  longitude_2d = (double *)malloc(size2D*sizeof(double));
  reflectivity = (double *)malloc(size3D*sizeof(double));
  reflectivity_2d = (double *)malloc(size2D*sizeof(double));
  gridLatitude = (double *)malloc(gridColumns*gridRows*sizeof(double));
  gridLongitude = (double *)malloc(gridColumns*gridRows*sizeof(double));
  noSwath = (double *)malloc(gridColumns*gridRows*sizeof(double));

  if (!latitude || !latitude_2d || !longitude || !longitude_2d || 
      !reflectivity || !reflectivity_2d || 
      !gridLatitude || !gridLongitude || !noSwath)  {
    fprintf(stdout, "Memory allocation error for one of input arrays.\n");
    exit(1);
  }

  /* allocate space for output array */
  interpolation = 
    (double *) malloc(LEVELS*gridColumns*gridRows*sizeof(double)); 

  if (!interpolation)  {
    fprintf(stdout, "Memory allocation error for output array.\n");
    exit(1);
  }

  /* read in input arrays */
  for (level=0; level<LEVELS; level++)  { 
    for (swath=0; swath<numSwaths; swath++)  {
      for (ray=0; ray<RAYS; ray++)  {
	fscanf(fpLatitude, "%lf", 
	       latitude + ((level*size2D) + (swath*RAYS) + ray));
	fscanf(fpLongitude, "%lf", 
	       longitude + ((level*size2D) + (swath*RAYS) + ray));
	fscanf(fpReflectivity, "%lf", 
	       reflectivity + ((level*size2D) + (swath*RAYS) + ray));
      }
    }
  }

  for (gridpoint=0; gridpoint<(gridColumns*gridRows); gridpoint++)  {
    fscanf(fpGridLatitude, "%lf", gridLatitude + (gridpoint));
    fscanf(fpGridLongitude, "%lf", gridLongitude + (gridpoint));
    fscanf(fpNoSwath, "%lf", noSwath + (gridpoint));
  }   	
    	
  /*****************************************************************************/
  /*****************************************************************************/

  /* calculate distance to each satellite footprint */
	
  for (level=0; level<LEVELS; level++)  { 
    if (level%10 == 0){ 
      fprintf(stdout, "level = %d\n", level);
    }

    /* grab reflectivity, latitude, and longitude for this level */
    for(j=0;j<size2D;j++)  {
      *(latitude_2d + j) = *(latitude + (level*size2D) + j);
      *(longitude_2d + j) = *(longitude + (level*size2D) + j);
      *(reflectivity_2d + j) = *(reflectivity + (level*size2D) + j);
    }

    /* FOR DEBUGGING GPM DATA ONLY */
    /*
    if (level==175) {
      MAX = NODATA;
      for (i=0; i<gridColumns*gridRows; i++)  {
	if (*(reflectivity_2d + j) > MAX) {
	  MAX = *(reflectivity_2d + j);
	}
      }
      if (MAX >= 0.0) {
	fprintf(stdout,"At 0.00km max refl = %f\n",10.0 * log10(MAX));
      }
      else {
	fprintf(stdout,"At 0.00km there is no positive value of zFactorCorrected\n");
      }
    }
    */

    for (i=0; i<gridColumns*gridRows; i++)  {

      if (noSwath[i]==1.0)  {
	interpolation[(level*gridRows*gridColumns + i)] = NODATA;
      } 

      else {

	/* allocate space for working arrays */
	weighting = (double *) malloc(size2D*sizeof(double)); 
	indices = (int *) malloc(size2D*sizeof(double));
	indices2 = (int *) malloc(size2D*sizeof(double));
	differences = (double *) malloc(size2D*sizeof(double)); 
	products = (double *) malloc(size2D*sizeof(double));
	distance = (double *) malloc(size2D*sizeof(double));
	filteredDistance = (double *) malloc(size2D*sizeof(double));

	if (!weighting || !indices || !indices2  || !differences || 
	    !products || !distance)  {
	  puts("Memory allocation error for one of working arrays.\n");
	  exit(1);
	} 

	/* initialize indices arrays */
	for(j=0;j<size2D;j++) {
	  indices2[j] = NODATA;
	  indices[j]  = NODATA;
	}

	/* find points within latitude/longitude box */
	searcher2(latitude_2d, longitude_2d, 
		  gridLatitude, gridLongitude, 
		  reflectivity_2d, indices2, 
		  &count2, size2D, i);

	for(j=0; j<count2; j++)  {
	  distance[indices2[j]] = 
	    degrees2distance(gridLatitude[i], latitude_2d[indices2[j]], 
			     gridLongitude[i], longitude_2d[indices2[j]]);
	}

	/* find points within RADIUS */
	searcher(distance, indices, indices2, &count, &count2, RADIUS);

	for(j=0; j<count; j++)  {
	  filteredDistance[j] = distance[indices[j]];
	}

	if (!count)  {

	  /* no points to interpolate */
	  interpolation[(level*gridRows*gridColumns + i)] = NODATA;

	} 

	else  { 

	  /* find weighting factors */
	  for(j=0;j<count;j++)  {
	    weighting[j] = cressman(filteredDistance[j],RADIUS);
	  }  

	  /* find nearest neighbor */
	  proximity = RADIUS;
	  for(j=0;j<count;j++)  {
	    if (filteredDistance[j] < proximity)  {
	      proximity = filteredDistance[j];
	      index = j;
	    }
	  }

	  /***********************************************************************/

	  /* this code is responsible for the interpolation */

	  neighbour = reflectivity_2d[indices[index]];
	  total = neighbour;

	  for(j=0;j<count;j++)  {
	    differences[indices[j]] = neighbour - reflectivity_2d[indices[j]];
	    products[j] = weighting[j]*differences[indices[j]];
	    interpolatedvalue = products[j] + total;
	    total = interpolatedvalue;
	    interpolation[(level*gridRows*gridColumns + i)] = interpolatedvalue;
	  }

	} /* end else (find weighting factors) */

	/* deallocate memory for working arrays */
	free(products);
	free(differences); 
	free(weighting);
	free(indices);
	free(indices2);
	free(distance);
	free(filteredDistance);

      } /* end else (allocate space for working arrays) */

    } /* for each gridpoint */

  } /* end for level... */

  /*****************************************************************************/
  /*****************************************************************************/

  /* For all non-NODATA values:
     1. convert linear values to log values and 
     2. filter out reflectivity values < MINDBZ
     3. Print all values to output file
  */

  for(j=0;j<LEVELS*gridRows*gridColumns;j++)  {
    if (interpolation[j] != NODATA)  {
      interpolation[j] = 10 * log10(interpolation[j]);
      if (interpolation[j] < MINDBZ )  {
	interpolation[j] = NODATA;
      }
    }
    fprintf(fpGridReflectivity, "%lf ", interpolation[j]);
  }
  fprintf(fpGridReflectivity, "\n");

  /*****************************************************************************/

  /* close files and deallocate memory */
  fclose(fpGridReflectivity);
  fclose(fpLatitude);
  fclose(fpLongitude);
  fclose(fpReflectivity);
  fclose(fpGridLatitude);
  fclose(fpGridLongitude);
  fclose(fpNoSwath);
  fclose(fpNumSwaths);
  fclose(fpGridRows);
  fclose(fpGridColumns);

  free(latitude);
  free(latitude_2d);
  free(longitude);
  free(longitude_2d);
  free(reflectivity);
  free(reflectivity_2d);
  free(gridLatitude);
  free(gridLongitude);
  free(noSwath);

  return(0);
		
} /* end of main /*

  /*******************************************************************************/

double degrees2distance(double gridLatitude, double latitude, 
			double gridLongitude, double longitude)  {
  double difflong = (DEG2RAD*longitude) - (DEG2RAD*gridLongitude);
  double absdifflong = fabs(difflong);

  return (111.24 * RAD2DEG * 
	  (acos(((sin(DEG2RAD*gridLatitude) * sin(DEG2RAD*latitude)) + 
		 (cos(DEG2RAD*gridLatitude) * cos(DEG2RAD*latitude) * 
		  cos(absdifflong))))));
}

/*******************************************************************************/

double cressman(double d2gp,float RADIUS)  {
  /* calculate weighting factor */
  return ((-1)*(((RADIUS*RADIUS)-(d2gp*d2gp))/((RADIUS*RADIUS)+(d2gp*d2gp))));
}

/*******************************************************************************/

void searcher(double *distance, int indices[], int indices2[], 
	      int *count, int *count2, float RADIUS)  {
  int i;
  int subi;

  *count = 0;
  subi = 0;

  for(i=0;i<*count2;i++)  {
    
    if (((*(distance + indices2[i])) < RADIUS) && 
	((*(distance + indices2[i])) >= 0.0))  {
      indices[subi] = (*(indices2 + i));
      (*count)++;
      subi++;
    }
  }
}

/*******************************************************************************/

void searcher2(double *latitude_2d, double *longitude_2d, 
	       double *gridLatitude, double *gridLongitude, 
	       double *reflectivity_2d, int indices2[], 
	       int *count2, int size2D, int i)  {

  int j;
  *count2 = 0;

  double temp1=(*(gridLatitude + i));
  double temp2=(*(gridLongitude + i));

  for(j=0;j<size2D;j++)  {
   		
    if (fabs(*(latitude_2d + j))<(fabs(temp1) + DEGREES) && 
	fabs(*(latitude_2d + j))>(fabs(temp1) - DEGREES) && 
	fabs(*(longitude_2d + j)<(fabs(temp2) + DEGREES)) && 
	fabs(*(longitude_2d + j))>(fabs(temp2) - DEGREES))  {
	 
      /* if (reflectivity_2d[j]>(-666))  { WHY WAS THIS HARDCODED IN? */
      if (reflectivity_2d[j] != NODATA)  {
	indices2[*count2] = j;
	(*count2)++;
      }						
    }
  }
}

/*******************************************************************************/
