/* VModis_interpol.cpp */

/*
remap
Copyright (C) 2006 Fabrice Ducos, fabrice.ducos@icare.univ-lille1.fr

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

// interpolation code for Modis, J. Descloitres (adapted by F. Ducos)

#include <math.h>
#include "VModis_interpol.h"

void get_step_and_offset(int n1, int n2, float *offset, float *step, 
			 int rowsperscan, int georowsamplingstep, unsigned char crosstrack)
{
  // new implementation by J. Descloitres and F. Ducos

  // Map i1 into i2:  i2 = offset + step * i1
  if (n1 >= n2) *step = 1. / (int) ((float)n1 / n2 + 0.5);         
  // Must add 0.5 to handle 1354/271 
  else  *step = (int) ((n2 + n1 - 1) / n1);

  *offset = ( *step - 1 ) / 2;

  if (crosstrack && georowsamplingstep > 1) {
    // this is to be validated (not tested)
    if (rowsperscan >= 20) {
      *offset += 0.5*10/(20 * georowsamplingstep);
    }
    if (rowsperscan >= 40) {
      *offset += 0.5*10/(40 * georowsamplingstep);
    }
    assert(rowsperscan < 80);
  }
}

void get_virtual_row(const coord_type *lon1, const coord_type *lat1, const coord_type *lon2, const coord_type *lat2, double t, int ncols, double *lon3, double *lat3)
{
  int geocol;
  const coord_type PROJUNDEF = DEFAULT_COORD_FILL_VALUE; // F. Ducos
  
  for (geocol=0; geocol<ncols; geocol++) {
    
    /* Draft version - UNDEF value must not be hard-coded eventually */
    if ( lon1[geocol] == 0  &&  lat1[geocol] == 0  ||   /* Filter (0,0)  
							   because of MOD09CRS problem */
	 lon2[geocol] == 0  &&  lat2[geocol] == 0  ||
	 lon1[geocol] == -999  ||
	 lon2[geocol] == -999  ||
	 lat1[geocol] == -999  ||
	 lat2[geocol] == -999 ) {
      lon3[geocol] = PROJUNDEF;                /* It is enough to mark lon3 as bad */
      continue;
    }
    lon3[geocol] = (1 - t) * lon1[geocol] + t * lon2[geocol];
    lat3[geocol] = (1 - t) * lat1[geocol] + t * lat2[geocol];
    if (lon1[geocol] - lon2[geocol] > 180) lon3[geocol] += t * 360;      
    /* Interpolation across the dateline */
    if (lon2[geocol] - lon1[geocol] > 180) lon3[geocol] -= t * 360;      
    /* Change lon3 but not lon1 and lon2 */
    if (lon3[geocol] >  180) lon3[geocol] -= 360;
    if (lon3[geocol] < -180) lon3[geocol] += 360;
    
  } 
}

void get_georows(double fractrow, PRODUCT *input, int iscan, int *georow1, int *georow2, double *t)
{
  int georowsperscan;
  
  georowsperscan = (int)(10 * input->Np_geo / 1354. + 0.5);     /* Number of geolocation rows per scan */
  
  *georow1 = (int) floor(input->rowoffset + fractrow * input->rowstep);
  *georow2 = *georow1 + 1;
  if (*georow1 < iscan * georowsperscan) {
    *georow1 = iscan * georowsperscan + 1;
    *georow2 = iscan * georowsperscan;
  }
  if (*georow2 > (iscan + 1) * georowsperscan - 1) {
    *georow1 = (iscan + 1) * georowsperscan - 1;
    *georow2 = (iscan + 1) * georowsperscan - 2;
  }
  
  *t = (input->rowoffset + fractrow * input->rowstep - *georow1) /  
    (*georow2 - *georow1);
  
}

    
void resample_row_geolocation(INTERP_OUTPUT output, double *oldlon, 
			      double *oldlat, PRODUCT *input, coord_type *newlon, coord_type *newlat)
{
  int jcol, geocol1, geocol2, npoints;
  double u, coloffset;
  const coord_type PROJUNDEF = DEFAULT_COORD_FILL_VALUE;
  
  switch (output) {
  case INTERP_BOUNDS:
    npoints = input->sds.Np + 1;
    coloffset = -0.5;
    break;
  case INTERP_CENTERS:
    npoints = input->sds.Np;
    coloffset = 0;
    break;
  default :
    assert(! APPNAME ": rsample_row_geolocation: unexpected value for parameter 'output'");
    break;
  }
  for (jcol=0; jcol<npoints; jcol++) {
    geocol1 = (int) floor(input->coloffset + (jcol + coloffset) * input->colstep);
    geocol2 = geocol1 + 1;
    if (geocol1 < 0) {
      geocol1 = 1;
      geocol2 = 0;
    }
    if (geocol2 > input->Np_geo - 1) {
      geocol1 = input->Np_geo - 1;
      geocol2 = input->Np_geo - 2;
    }
      
    /* F. Ducos, we do not compare oldlon with PROJUNDEF for it could be a 'Not a Number' */

    /* old code
      if ( oldlon[geocol1] == PROJUNDEF  ||
      oldlon[geocol2] == PROJUNDEF ) {
      newlon[jcol] = PROJUNDEF; // It is enough to mark newlon as bad
      continue;   // Skip to next pixel in current row
      }
    */
 
    if ( ((-180. <= oldlon[geocol1] && oldlon[geocol1] <= 180.) == false) 
	 || ((-180. <= oldlon[geocol2] && oldlon[geocol2] <= 180.) == false) ) {
      newlon[jcol] = PROJUNDEF;
      continue;
    }

    u = (input->coloffset + (jcol + coloffset) * input->colstep - 
	 geocol1) / (geocol2 - geocol1);
    newlon[jcol] = (1 - u) * oldlon[geocol1] + u * oldlon[geocol2];
    newlat[jcol] = (1 - u) * oldlat[geocol1] + u * oldlat[geocol2];
    if (oldlon[geocol1] - oldlon[geocol2] >= 180)       /* 
							   Interpolation across the dateline */
      newlon[jcol] += u * 360;
    if (oldlon[geocol2] - oldlon[geocol1] >= 180)
      newlon[jcol] -= u * 360;
    if (newlon[jcol] >  180) newlon[jcol] -= 360;
    if (newlon[jcol] < -180) newlon[jcol] += 360;
  }
  
}