/* mktransect.c
 *                *
 *              Emilio          4/1/99
 *                Updated program.
 *              Emilio          4/21/98
 *                *
 * COMPILE AS: gccemu mktransect
 * RUN AS: mktransect <lon_o>   
 * eg, mktransect -70 1 -55 -10  OR  mktransect -70 4 -70 -10
 ******************************************************************************/

#include "emu.h"

const char * cpszUsage = {"usage: mktransect \n\
Create a binary or NC file with containing a transect mask.\n\
\n\
Command-line switches specific to this program are:\n\
\t-xy lon_o lat_o lon_f lat_f \n\
\n\
"};

/* local functions */
void Setup(void);
void Process(void);
void Output(void);

int lonlatToidx(float lon, float lat);
void idxTolonlat(int i, float *lon, float *lat);
int getNormalidx(char *dir, float normalangle, float delDegdist, float lon, float lat);
float distlonlatTokm(float lon1, float lat1, float lon2, float lat2);

/* local variables */
float lon_o, lat_o, lon_f, lat_f;
WORD *transect;
float scale_factor = 0.1;
int errnc;


int main(int argc, char *argv[])
{
  AssignString(&pszNCtitle, "Transect Mask");
  AssignString(&pszNChist, "transect mask created by mktransect");

  Initialize(&argc, argv, NOMASK);

  Setup();
  Process();
  Output();

  CleanUp();

  return(0);
}

/* ================================================================= */
void Setup(void)
{
  int i;

  transect = alloc1d_us(0,nRows*nCols-1);
  for (i = 0; i < nRows*nCols; i++)
    transect[i] = -999;
}

/* ================================================================= */
void Process(void)
{
  int i, i_o, i_f, col, row;        int bug=0,flag=0;
  WORD distval, Totpts=0;
  int col_o, row_o, col_f, row_f;
  float lon, lat, del_lon;
  float angle, normalangle, theta, delDegdist, dist, Totdist;


  i_o = lonlatToidx(lon_o, lat_o);
  idxTolonlat(i_o, &lon_o, &lat_o);
  i_f = lonlatToidx(lon_f, lat_f);
  idxTolonlat(i_f, &lon_f, &lat_f);

  (void) getRowRevCol(i_o, &col_o, &row_o);
  (void) getRowRevCol(i_f, &col_f, &row_f);
  if (bDebug)  {
    printf("HERE:%d, lon_o=%.3f lat_o=%.3f lon_f=%.3f lat_f=%.3f \n", bug++, lon_o,lat_o,lon_f,lat_f);
    printf("HERE:%d, col_o=%d row_o=%d col_f=%d row_f=%d \n", bug++, col_o,row_o,col_f,row_f);
    }

  /* angle is angle from 1st endpoint (_o) to last endpoint (_f), with respect to the
     "x-axis" (>0 in the increasing-longitude direction); range: -179.999... to 180  */
  angle = RAD2DEG * atan2((lat_f - lat_o) , (lon_f - lon_o));

  /* normalangle is the angle perpendicular to the vector, in "right-handed" direction  */
  normalangle = (angle >= -90.0 && angle < -180.0) ? (angle + 270.0) : (angle - 90.0);

  /* delDegdist: distance step along vector direction; range: fHScale*(1 -> sqrt(2)) */
  theta = fabs(fabs(angle) - 90.0);
  theta = (theta > 45) ? (90.0 - theta) : theta;   /* theta: 0->45 deg */
  delDegdist = fHScale/cos(theta/RAD2DEG);

  Totdist = distlonlatTokm(lon_f, lat_f, lon_o, lat_o);

  /* Put origin value into transect array */
  dist    = 10.0;        /* The origin will be 1 km */
  distval = (WORD) floor((double)dist + 0.5);

  lon = lon_o;
  lat = lat_o;
  i = i_o;

  transect[i] = distval;
  i = getNormalidx("right", normalangle, delDegdist, lon, lat);
  transect[i] = distval;
  i = getNormalidx("left", normalangle, delDegdist, lon, lat);
  transect[i] = distval;
  Totpts+=3;

  if (bDebug)  printf("HERE:%d, Totdist=%.3f (0.1km), angle=%.4f, normalangle=%.4f,
         delDegdist=%.4f,  i=%d\n", bug++, Totdist,angle,normalangle,delDegdist,i);

  while (dist <= (Totdist + 10.0))
  {
    dist   += (DEG2KM / scale_factor) * delDegdist;
    distval = (WORD) floor((double)dist + 0.5);

    lon += delDegdist * cos(angle/RAD2DEG);
    lat += delDegdist * sin(angle/RAD2DEG);
    i = lonlatToidx(lon, lat);

    (void) getRowRevCol(i, &col, &row);  flag++;
    if (!bQuiet && ((flag % 20) == 0 || flag==1))
      printf("  i=%d,  lon=%.3f,lat=%.3f, col=%d,row=%d,  dist=%.3f \n",
             i,lon,lat,col,row,dist);

    transect[i] = distval;
    i = getNormalidx("right", normalangle, delDegdist, lon, lat);
    transect[i] = distval;
    i = getNormalidx("left", normalangle, delDegdist, lon, lat);
    transect[i] = distval;
    Totpts+=3;
  }

  if (!bQuiet)  printf("\n   Total Number of Points Selected: %d \n", Totpts);
}

/* ================================================================= */
/* ==================== WRITE OUT BINARY FILE ====================== */
void Output(void)
{
  FILE *fout;
  int ncID;               /* netCDF file ID */
  int tID;                /* time axis ID   */
  int varID;
  char varName[NAMELEN];      /* Name of NC variable(s) to be written      */
  char varLongName[NAMELEN];  /* Long Name of NC variable(s) to be written */


  if (bOutBinary)    /* write output to binary files */
  {
    fout = OUTPUTfile(pszOutputFile, ".bil");
    fwrite(transect, sizeof(WORD), nRows*nCols, fout);
    fclose(fout);
  }
  else
  {
    /* setup new NC file from NC headers */
    SetupDefaultAxes(FALSE, FALSE);
    ncID = newNCsetup(&tID, pszOutputFile, pszNCtitle, pszNChist);

    strcpy(varName, "transect");
    sprintf(varLongName,"%s mask", varName);
    /* the NC lib doesn't support unsigned short, that MAY cause a problem */
    varID = DefNCvar(ncID, varName, EMU_SHORT, "xy", varLongName, "km");
    AddNCCompressionAttributes(ncID, varID, scale_factor, 0.0);
    errnc = nc_put_var_short(ncID, varID, transect);
    errnc = nc_close(ncID);
  }
}

/* ================================================================= */
void ProcessCommandLineArgs(int * pArgc, char * argv[])
{
  int i;

  for (i=1; i<*pArgc; i++)
    if (strcmp(argv[i], "-xy") == 0)
    {
      lon_o = atof(argv[++i]);
      lat_o = atof(argv[++i]);
      lon_f = atof(argv[++i]);
      lat_f = atof(argv[++i]);
    }
}

/* ================================================================= */
void LocalCleanUp(void)
{
  free(pszOutputFile);
  free1d_us(transect,0,nRows*nCols-1);
}


/* ================================================================= */
/* get the 1-d location index from (lon,lat)         */
int lonlatToidx(float lon, float lat)
{
  int row, col;

  col = (int) floor((double)(lon - fLonLeft)/fHScale);
  row = (int) floor((double)(fLatTop - lat)/fVScale);

  return (row*nCols + col);
}

/* reshift (lon,lat) point to the center of the cell */
void idxTolonlat(int i, float *lon, float *lat)
{
  int row, col;

  (void) getRowRevCol(i, &col, &row);
  *lon = fLonLeft + fHScale*((float)col + 0.5);
  *lat = fLatTop  - fVScale*((float)row + 0.5);
}

int getNormalidx(char *dir, float normalangle, float delDegdist, float lon, float lat)
{
  float lonNorm, latNorm, angle;

  if (istype(dir, "left"))
     angle = (normalangle >= 0.0) ? (normalangle-180.0) : (normalangle+180.0);
  else
     angle = normalangle;

  lonNorm = lon + delDegdist * cos(angle/RAD2DEG);
  latNorm = lat + delDegdist * sin(angle/RAD2DEG);

  return (lonlatToidx(lonNorm, latNorm));
}

float distlonlatTokm(float lon1, float lat1, float lon2, float lat2)
{
  return ((float) sqrt((double)(SQR(lon2 - lon1) + SQR(lat2 - lat1))) * DEG2KM * 10);
}