/* util.c
 *          *
 *        Emilio            4/18/2001
 *          In SkipLines(), expanded the maximum line length to
 *          LONGLINELEN.
 *          *
 * Utilities file for all the C programs written to analyze binary and NetCDF
 * input and output for modelling using the Environmental Modeling Utilities,
 * or "EMU" package.
 ******************************************************************************/

#include "emu.h"


#define FREE_ARG char*
#define NR_END 0  /* this is supposed to provide a memory buffer - it makes allocation and
                     deallocation within the allocation functions a pain.  I've set it to 0.
                     If this is changed then memory allocation especially for VAR will be affected */


/* from function.c */
extern char * pszProgramName;

/* Local functions */
void ConvertToLowerCase(char * s);


/* ================================================================= */
/* find the month and year for the given time-step */

void getyrmonth(int tvar, int *year, int *month)
{
  int t;

  t = tvar + nFirstMonth - 1;
  *year = nFirstYear + (int)floor((float)t / 12.0);
  *month = t % 12;

}

/* ========================================================================= */
/*  increment yr and mo (where mo cycles, 0->11) */

void NEXTyrmo(int *yr, int *mo)
{
  if (*mo < 11)
    (*mo)++;
  else  /* if month is Dec, cycle *mo back to 0 (Jan), increment *yr */
  {
    *mo = 0;
    (*yr)++;
  }
}

/* ================================================================= */
/* return a cell index where the row-order is reversed from that used  *
 * in the cell index that's the first argument.  col and row are       *
 * base-0, and will be in the "right" order (left->right, top->bottom),*
 * the same implied by the first argument cell index                   */

long int getRowRevCol(long int i, int *col, int *row)
{
  long int irev;

  /* rows in NC files are arranged from S to N, in the
     OPPOSITE direction as rows in our binary image files */
  *col = i % nCols;
  *row = (int) ((float)(i - *col)/(float)nCols);
  irev = ((nRows - 1) - *row) * nCols + *col;

  return(irev);
}
/* ================================================================= */
/* calculate the are of cell i, using function area() */

double cellarea(int i)
{
  int row, col;
  double lat;

  (void) getRowRevCol(i, &col, &row);
  lat = fLatTop - (((double) row + 0.5) * fVScale);

  return (area(lat, fVScale, fHScale));
}

/* ================================================================= */
/*  library function revised by Matthew Thompson, 7/13/95 */

double area(double latitude, double vscale, double hscale)
/*
 * 'latitude' should be the latitude at the center of the cell (positive if
 * degrees N, negative if degrees S). 'vscale' is the number of vertical
 * degrees per cell - 1x1 map, 'vscale' would be '1'. 'hscale' is the number
 * of horizontal degrees per cell - 1x1 map, 'hscale' would be '1'. This
 * function assumes that the radius of the earth varies sinusoidally with
 * latitude - see PostScript documentation.
 */
/* returns area of a cell in m^2 */
{
  double          toplat, btmlat;
  /*
   * toplat = latitude at top of cell - btmlat = latitude at bottom of
   * cell (in degrees)
   */

  double          thetatop, thetabtm;
  /*
   * theta = latitude in radians from equator - thetatop = latitude at
   * top of cell in radians - thetabtm = latitude at bottom of cell in
   * radians
   */

  double          REQUATOR = 6378136;
  double          RPOLE = 6356751;
  double          D = REQUATOR - RPOLE;
  /* polar and equatorial radii of the earth in meters */
  /* from encyclopaedia brittanica online - the earth */

  double          NUMLAT = 180;
  /* number of vertical degrees in globe */

  double          NUMLONG = 360;
  /* number of horizontal degrees in globe */

  double          area = 0.0;
  /* calculated area */


  /* PERFORM CALCULATIONS */
  toplat = latitude + (vscale / 2);
  btmlat = latitude - (vscale / 2);

  thetatop = PI * (toplat / NUMLAT);
  thetabtm = PI * (btmlat / NUMLAT);

  area += ((RPOLE * RPOLE) + (2 * D * D / 3)) * sin(thetatop);
  area += D * RPOLE * cos(thetatop) * sin(thetatop);
  area += D * RPOLE * thetatop;
  area += (D * D / 3) * sin(thetatop) * cos(thetatop) * cos(thetatop);
  /* add first part of integral */

  area -= ((RPOLE * RPOLE) + (2 * D * D / 3)) * sin(thetabtm);
  area -= D * RPOLE * cos(thetabtm) * sin(thetabtm);
  area -= D * RPOLE * thetabtm;
  area -= (D * D / 3) * sin(thetabtm) * cos(thetabtm) * cos(thetabtm);
  /* subtract second part of integral */

  /*
   * now have area of entire latitudinal band defined by latitude and
   * vscale - next:  get area of individual cell, defined by hscale
   */
  return (2 * PI * area * (hscale / NUMLONG));
}
/* ================================================================= */
/* round off value to an integer */
int roundoff(double value)
{
  return((int)floor(value + 0.5));
}
/* ================================================================= */
float CorrCoef(float * pf1, float * pf2, int nOffset)
{
  int x,y;
  float fMean1, fMean2;
  float fDif1, fDif2;
  float f11 = 0.0F;
  float f22 = 0.0F;
  float f12 = 0.0F;
  float * * ppf1;
  float * * ppf2;

  ppf1 = VectorToMatrix(pf1);
  ppf2 = VectorToMatrix(pf2);

  fMean1 = Mean(ppf1, nOffset);
  fMean2 = Mean(ppf2, nOffset);

  for (y=nOffset; y<nRows-nOffset;y++)
  {
    for (x=nOffset; x<nCols-nOffset;x++)
    {
      fDif1 = ppf1[y][x] - fMean1;
      fDif2 = ppf2[y][x] - fMean2;
      f11 += fDif1 * fDif1;
      f22 += fDif2 * fDif2;
      f12 += fDif1 * fDif2;
    }
  }
  FreeVectorToMatrix(ppf1);
  FreeVectorToMatrix(ppf2);

  return(f12/(float) sqrt((double) (f11*f22)));
}

/* ================================================================= */
float Mean(float * * ppfData, int nOffset)
{
  int y,x;
  float fSum = 0.0F;
  int nCnt = 0;

  for (y=nOffset; y<nRows-nOffset;y++)
    for (x=nOffset; x<nCols-nOffset;x++)
    {
      fSum +=  ppfData[y][x];
      nCnt++;
    }

  return(fSum/(float) nCnt);
}

/* ================================================================== */
void SwapBytes(BYTE * pbyData, int nSwapBytes, int nBytes)
{
  BYTE * pbyTemp;
  char s[81];
  int i;

  if (nSwapBytes == 2)
  {
    pbyTemp = alloc1d_uc(0,nBytes-1);

    memcpy(pbyTemp, pbyData, nBytes);
    swab(pbyTemp, pbyData, nBytes);

    free1d_uc(pbyTemp,0,nBytes-1);
  }
  else if (nSwapBytes == 4)
  {
    pbyTemp = alloc1d_uc(0,nBytes-1);
    memcpy(pbyTemp, pbyData, nBytes);

    for (i=0; i<nBytes; i++)
    {
      pbyData[3 + 4 * (i / 4) - (i % 4)] = pbyTemp[i];
    }

    free1d_uc(pbyTemp,0,nBytes-1);
  }
  else
  {
    sprintf(s, "Swap Bytes for %d bytes not supported\n", nSwapBytes);
    Warning(s, "");
  }
}
/* ================================================================== */
int IsBad(double dTest)
{
  char s[81];

#ifdef PC
  if (_isnan(dTest))
    return (1);

  sprintf(s, "%lf", dTest);
  if (strcmp("1.#INF00", s))
    return (0);
  else
    return (1);
#endif
#ifdef UNIX
  return (!finite(dTest));
/*
  if (isnan(dTest))
    return (1);

  sprintf(s, "%lf", dTest);
  if (strcmp("Inf", s))
    return (0);
  else
    return (1);
*/
#endif
}

/* ================================================================= */
/* returns 1 if var=type, 0 otherwise     */

int istype(char *var, char *type)
{
/*  ConvertToLowerCase(var);
  ConvertToLowerCase(type); */
  if (strncmp(var, type, strlen(type)) == 0)
    return(TRUE);
  else
    return(FALSE);
}

/* ================================================================= */
/* Allocate memory and assign string to pointer  */

void AssignString(char * * ppszTarget, char * pszSource)
{
	/* *ppszTarget = (char *) malloc( strlen(pszSource) + 1 ); */
  *ppszTarget = (char *) malloc((size_t)(strlen(pszSource)+1)*sizeof(char));
  strcpy(*ppszTarget, pszSource);
}

/* ================================================================== */

char * GetQuotedText(int * i, char * argv[])
{
  char s[LINELEN];
  char * pszText;


  if (argv[*i][0] == '"')
  {
    strcpy(s, &argv[*i][1]);
    (*i)++;
    while (argv[*i][strlen(argv[*i])-1] != '"')
    {
      strcat(s, argv[*i]);
      (*i)++;
    }
    strcat(s, argv[*i]);
    s[strlen(s)-1] = '\0';
  }
  else
    strcpy(s, argv[*i]);

  (*i)++;
  pszText = (char *) malloc(strlen(s)+1);
  strcpy(pszText, s);

  return(pszText);
}

/* ================================================================= */
/* drop trailing (surrounding) spaces, if present */
void DropTrailingSpaces(char *pszStr)
{
  char  s[LINELEN], *pCh;
  int   i;

  strcpy(s, pszStr);

  /* drop spaces at the beginning of the string */
  i = 0;
  while (s[i] == ' ')
    i++;
  pCh = &s[i];
  strcpy(s, pCh);

  /* drop spaces at the end of the string */
  pCh = &s[strlen(s)-1];
  while (*pCh == ' ')
    *(pCh--) = '\0';

  strcpy(pszStr, s);
}

/* ================================================================= */
/* drop surrounding quotes, if present */
void DropSurroundingQuotes(char *pszStr)
{
  char s[LINELEN], *pCh;

  strcpy(s, pszStr);

  /* drop quotes at the beginning of the string */
  pCh = &s[0];
  if (*pCh == '"')
  {
    *pCh = '\0';
    strcpy(s, ++pCh);
  }

  /* drop quotes at the end of the string */
  if (s[strlen(s)-1] == '"')
    s[strlen(s)-1] = '\0';

  strcpy(pszStr, s);
}

/* ================================================================= */
void ClearString(char *string, int StringLen)
{
  int c;

  for (c=0; c < StringLen; c++)
    string[c] = '\0';
/*
  for (c=0; c < StringLen-1; c++)
    string[c] = ' ';  * or fill it up with \0 ?? *

  string[StringLen-1] = '\0';
*/
}

/* ================================================================= */
int strcnt(char * psz, char c)
{
  int i = 0;
  int n = 0;

  while (psz[i])
  {
    if (psz[i++] == c)
      n++;
  }
  return(n);
}

/* ================================================================= */
/* NOTE: the space delimeter is a tricky case that must be handled with care;
   eg, multiple continuous spaces should be treated as a single space  */
/* NOTE: this function has memory. If a member of the array has been allocated in
   a previous call, b/c it's not cleared, the member will remain in future calls */

int GetLineElements(char **ppszArgs, char *pszLine, char delimeter)
{
  int   n = 0, i;
  char  s[LINELEN], *pCh, *pChSpace, *pChTmp;


  strcpy(s, pszLine);

  while ((pCh = strchr(s, delimeter)) != NULL)
  {
    /* space delimeter; treated in a special way */
    if (delimeter == ' ')
    {
      if (s[0] == delimeter)
      {
        /* drop all spaces at the beginning of the string, advance pointer */
        i = 0;
        while (s[i] == delimeter)
          i++;
        pChSpace = &s[i];
        strcpy(s, pChSpace);

        if ((pCh = strchr(s, delimeter)) == NULL) /* this is the last token */
          break;  /* exit while loop */
      }

      /* search for the quote character (") within this token;             *
       * extend token to the ending quote if token is surrounded by quotes */
      if (s[0] == '"')
      {
        if ((pChTmp = strchr(&s[1], '"')) != NULL) /* an enclosing quote is present */
          pCh = pChTmp;
      }

      /* split the string s */
      *pCh = '\0';
    }
    /* other delimeters */
    else
    {
      /* account for blank (empty) token */
      if (s[0] == delimeter)
      {
        AssignString(&ppszArgs[n++], "\0");
        *pCh = '\0';
        strcpy(s, ++pCh);
        continue;   /* continue to the next pass in the while loop */
      }

      /* search for the quote character (") within this token;             *
       * extend token to the ending quote if token is surrounded by quotes */
      if (s[0] == '"')
      {
        if ((pChTmp = strchr(&s[1], '"')) != NULL) /* an enclosing quote is present */
          pCh = pChTmp;
      }

      /* split the string s */
      *pCh = '\0';

      /* drop trailing spaces, if present */
      DropTrailingSpaces(s);
    }

    /* drop surrounding quotations, if present */
    DropSurroundingQuotes(s);

    /* assign token, then copy new split string */
    AssignString(&ppszArgs[n++], s);
    strcpy(s, ++pCh);
  }

  /* ----------- process the last entry ----------------- */

  /* Account for a blank (empty) last token, but only if the delimeter is NOT
     the space delimeter (that case is interpreted differently, because of
     complications with quoted elements and space delimeter, for the last entry.
     If there are no elements in the line (empty line), the line
     will be treated as having 0 elements.                                       */
  if (s[0] == '\n' || s[0] == '\0')
  {
    /* '\0' is returned if the element is empty */
    if (n > 0 && delimeter != ' ')
      AssignString(&ppszArgs[n++], "\0");
    return(n);
  }

  /* remove end-of-line character */
  if (s[strlen(s)-1] == '\n')
    s[strlen(s)-1] = '\0';

  /* drop trailing spaces and surrounding quotations, if present */
  DropTrailingSpaces(s);
  DropSurroundingQuotes(s);

  /* assign token, then split and copy new string */
  AssignString(&ppszArgs[n++], s);

  return(n);
}

/* ================================================================== */
/* Function to parse a line of fixed-width columns  */
/* indices are base-0; maximum line length is LONGLINELEN characters long */

void GetFixedWidthElements(char *ppszArgs[], int TotElements, char *pszLine, int indices[][2])
{
  char s[LONGLINELEN], *pCh;
  int token, len;

  for (token=0; token < TotElements; token++)
  {
    pCh = &pszLine[indices[token][0]];
    strcpy(s, pCh);
    len = indices[token][1] - indices[token][0] + 1;
    s[len] = '\0';

    AssignString(&ppszArgs[token], s);
  }
}

/* ================================================================= */
void FreeLineElements(char **StringArray, int TotElements)
{
  int i;

  for (i=0; i < TotElements; i++)
    free((char *)StringArray[i]);
}

/* ================================================================= */

void DoesFileExist(char * pszFile)
{
  FILE * fp;

  if ((fp = fopen(pszFile, "r")) == NULL)
    Warning("File %s does not exist\n", pszFile);
  else
    fclose(fp);
}

/* ========================================================================= */
/* read a line from fin, write line to char pointer argument "line"
 * with \n removed;  return 0 if EOF was reached, otherwise line length */

int fgetline(char *line, int lim, FILE *fin)
{
  register int c;
  register char *cs;

  cs = line;
  while((--lim > 0) && ((c = fgetc(fin)) != EOF) && (c != '\n'))
    *cs++ = c;

  *cs = '\0';

  return ((c == EOF && cs == line) ? 0 : strlen(line));
}
/* ================================================================= */
/* Count total number of lines in file, in one sweep through file;   *
 * when done rewind file pointer locator to the beginning            */
long GetTotalLines(FILE *fid)
{
  int TotalLines = 0;

  while (SkipLines(fid, 1))
    TotalLines++;
  rewind(fid);

  return(TotalLines);
}
/* ================================================================= */
/* Skip a number of lines equal to SkipSize, in the file       *
 * maximum length of char array ("width" of line): LONGLINELEN */
BOOL SkipLines(FILE *fid, long SkipSize)
{
  long n;
  char line[LONGLINELEN];

  for (n=0; n < SkipSize; n++)
    if (fgetline(line, LONGLINELEN, fid) == 0)  /* end of file has been reached */
      return(FALSE);

  return(TRUE);
}

/* ========================================================================= */

void PrintTime(void)
{
  /* THESE VARIABLES ARE USED TO PROVIDE CALENDAR TIMES TO THE USER */
  time_t          TIME;   /* time.h defined type - used to store local
                           * time for the function localtime() - this
                           * type is defined in time.h */
  struct tm      *DATE;   /* structure containing time information */
  char           *sTime;  /* string used to print time */
  char           *spTime; /* pointer to sTime */


  TIME = time(NULL);
  if ((DATE = localtime(&TIME)) == NULL)  /* no internal clock */
  {
    fprintf(stderr, "Local time not available.\n");
    fprintf(stdout, "\n");
  }
  else
  {
    sTime = asctime(DATE);  /* asctime() - provides allocated
                             * memory to this string */

    /* replace \n string provided by asctime() with \0 */
    spTime = sTime;
    while (*spTime != '\n')
      spTime++;
    *spTime = '\0';

    fprintf(stdout, "%s\n", sTime);
  }
}
/* ================================================================== */

void Warning(char * pszWarningString, char * pszArg)
{
  if (!bQuiet)
  {
    printf("%s :\n", pszProgramName);
    printf(pszWarningString, pszArg);
  }
}
/* ================================================================= */
void ExitError(char * pszErrorString, char * arg)
{
  printf("%s :\n", pszProgramName);
  printf(pszErrorString, arg);
  if (bDebug)
    printf("\nuse %s -h to get a help screen, continuing in debug mode\n", pszProgramName);
  else
  {
    printf("\nuse %s -h to get a help screen, exiting to system\n", pszProgramName);
    CleanUp();
    exit(1);
  }
}
/* ================================================================== */
/* standard error handler
  (obsolescent; eventually will be replaced by ExitError)             */

void nrerror(char error_text[], char *errorvar)
{
  fprintf(stderr,"Run-time error:  ");
  fprintf(stderr, error_text, errorvar);
  fprintf(stderr,"...now exiting to system...\n");
  exit(1);
}


/* ================================================================= */
/* =========== MEMORY ALLOCATION/DEALLOCATION FUNCTIONS ============ */

LONG l1d_VAR = 0;

VAR * alloc1d_VAR(long nl, long nh)
/* allocate a VAR _1d with subscript range v[nl..nh] */
{
  VAR * v;
  int i;

  v=(VAR *) malloc((size_t) (nh-nl+1+NR_END)*sizeof(VAR) );
  if (!v)
    nrerror("allocation failure in alloc1d_VAR()","");

  for(i=nl+NR_END;i<=nh+NR_END;i++)
  {
    v[i].nParams    = 0;
    v[i].nStrParams = 0;

    v[i].nDyn = 0;
    v[i].pstDyn = NULL;
  }

  l1d_VAR++;
  return v-nl+NR_END;
}
/* ================================================================== */

void free1d_VAR(VAR * pstVarTmp, long nl, long nh)
/* free a VAR _1d allocated with alloc1d_VAR() */
{
  int i, j, v;
  DYNFILE * pstDyn;

  l1d_VAR--;
  for (v=nl+NR_END; v<=nh+NR_END; v++)
  {
    if (pstVarTmp[v].nParams > 0)
      free1d_f(pstVarTmp[v].P, 0,pstVarTmp[v].nParams-1);

    if (pstVarTmp[v].nStrParams > 0)
      for (i=0; i < pstVarTmp[v].nStrParams; i++)
        free(pstVarTmp[v].szP[i]);

    if (pstVarTmp[v].nDyn)
    {
      for (j=0; j < pstVarTmp[v].nDyn; j++)
      {
        pstDyn = &(pstVarTmp[v].pstDyn[j]);
        if (pstDyn->nElements)
        {
          for (i=0; i < pstDyn->nElements; i++)
            free(pstDyn->ppszElement[i]);
          free(pstDyn->ppszElement);
        }
      }
      free(pstVarTmp[v].pstDyn);
    }
  }

  free((FREE_ARG) (pstVarTmp+nl-NR_END));
}
/* ================================================================== */

MASK *alloc1d_MASK(long nl, long nh)
/* allocate a MASK _1d with subscript range v[nl..nh] */
{
  MASK *v;

  v=(MASK *) malloc((size_t) (nh-nl+1+NR_END)*sizeof(MASK) );
  if (!v) nrerror("allocation failure in alloc1d_MASK()","");

  return v-nl+NR_END;
}

/* ================================================================= */
int *alloc1d_i(long nl, long nh)
/* allocate an int _1d with subscript range v[nl..nh] */
{
  int *v;

  v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
  if (!v) nrerror("allocation failure in alloc1d_i()","");
  return v-nl+NR_END;
}

/* ================================================================= */
void free1d_i(int *v, long nl, long nh)
/* free an int _1d allocated with alloc1d_i() */
{
  free((FREE_ARG) (v+nl-NR_END));
}

/* ================================================================= */
LONG lBytesAllocated;
LONG l1d_s = 0;
short *alloc1d_s(long nl, long nh)
/* allocate a short _1d with subscript range v[nl..nh] */
{
  short *v;
  size_t lBytes;

  lBytes = (size_t) (nh-nl+1+NR_END)*sizeof(short);

  v=(short *) malloc(lBytes);
  if (!v)
    nrerror("allocation failure in alloc1d_s()","");

  lBytesAllocated += lBytes;

  l1d_s++;
  return v-nl+NR_END;
}

/* ================================================================== */
void free1d_s(short *v, long nl, long nh)
/* free a short _1d allocated with alloc1d_s() */
{
  free((FREE_ARG) (v+nl-NR_END));
  l1d_s--;
}

/* ================================================================= */
LONG l1d_l = 0;
long *alloc1d_l(long nl, long nh)
/* allocate a long _1d with subscript range v[nl..nh] */
{
  long *v;
  size_t lBytes;

  lBytes = (size_t) (nh-nl+1+NR_END)*sizeof(long);

  v=(long *) malloc(lBytes);
  if (!v)
    nrerror("allocation failure in alloc1d_l()","");

  lBytesAllocated += lBytes;

  l1d_l++;
  return v-nl+NR_END;
}

/* ================================================================== */
void free1d_l(long *v, long nl, long nh)
/* free a long _1d allocated with alloc1d_l() */
{
  free((FREE_ARG) (v+nl-NR_END));
  l1d_l--;
}

/* ================================================================== */
double *alloc1d_d(long nl, long nh)
/* allocate a double _1d with subscript range v[nl..nh] */
{
  double *v;

  v=(double *) malloc((size_t) (nh-nl+1+NR_END)*sizeof(double) );
  if (!v) nrerror("allocation failure in alloc1d_d()","");

  return v-nl+NR_END;
}

/* ================================================================== */
void free1d_MASK(MASK *v, long nl, long nh)
/* free a MASK _1d allocated with alloc1d_MASK() */
{
  free((FREE_ARG) (v+nl-NR_END));
}

/* ================================================================== */
void free1d_d(double *v, long nl, long nh)
/* free a double _1d allocated with alloc1d_d() */
{
  free((FREE_ARG) (v+nl-NR_END));
}

/* ================================================================== */
double **alloc2d_d(long nrl, long nrh, long ncl, long nch)
/* allocate a double _2d with subscript range m[nrl..nrh][ncl..nch] */
{
  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  double **m;

  /* allocate pointers to rows */
  m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
  if (!m) nrerror("allocation failure 1 in alloc2d_d()","");
  m += NR_END;
  m -= nrl;

  /* allocate rows and set pointers to them */
  m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
  if (!m[nrl]) nrerror("allocation failure 2 in alloc2d_d()","");
  m[nrl] += NR_END;
  m[nrl] -= ncl;

  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

  /* return pointer to array of pointers to rows */
  return m;
}

/* ================================================================== */
short ***alloc3d_s(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a short _3d with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
  long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
  short ***t;

  /* allocate pointers to pointers to rows */
  t=(short ***) malloc((size_t)((nrow+NR_END)*sizeof(short**)));
  if (!t) nrerror("allocation failure 1 in alloc3d_f()","");
  t += NR_END;
  t -= nrl;

  /* allocate pointers to rows and set pointers to them */
  t[nrl]=(short **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(short*)));
  if (!t[nrl]) nrerror("allocation failure 2 in alloc3d_s()","");
  t[nrl] += NR_END;
  t[nrl] -= ncl;

  /* allocate rows and set pointers to them */
  t[nrl][ncl]=(short *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(short)));
  if (!t[nrl][ncl]) nrerror("allocation failure 3 in alloc3d_s()","");
  t[nrl][ncl] += NR_END;
  t[nrl][ncl] -= ndl;

  for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
  for(i=nrl+1;i<=nrh;i++) {
    t[i]=t[i-1]+ncol;
    t[i][ncl]=t[i-1][ncl]+ncol*ndep;
    for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
  }

  /* return pointer to array of pointers to rows */
  return t;
}

/* ================================================================== */
void free3d_s(short ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* free a short _3d allocated by alloc3d_s() */
{
  free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
  free((FREE_ARG) (t[nrl]+ncl-NR_END));
  free((FREE_ARG) (t+nrl-NR_END));
}

/* ================================================================== */
float ***alloc3d_f(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a float _3d with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
  long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
  float ***t;

  /* allocate pointers to pointers to rows */
  t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
  if (!t) nrerror("allocation failure 1 in alloc3d_f()","");
  t += NR_END;
  t -= nrl;

  /* allocate pointers to rows and set pointers to them */
  t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
  if (!t[nrl]) nrerror("allocation failure 2 in alloc3d_f()","");
  t[nrl] += NR_END;
  t[nrl] -= ncl;

  /* allocate rows and set pointers to them */
  t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
  if (!t[nrl][ncl]) nrerror("allocation failure 3 in alloc3d_f()","");
  t[nrl][ncl] += NR_END;
  t[nrl][ncl] -= ndl;

  for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
  for(i=nrl+1;i<=nrh;i++) {
    t[i]=t[i-1]+ncol;
    t[i][ncl]=t[i-1][ncl]+ncol*ndep;
    for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
  }

  /* return pointer to array of pointers to rows */
  return t;
}

/* ================================================================== */
void free3d_f(float ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* free a float _3d allocated by alloc3d_f() */
{
  free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
  free((FREE_ARG) (t[nrl]+ncl-NR_END));
  free((FREE_ARG) (t+nrl-NR_END));
}

/* ================================================================== */
double ***alloc3d_d(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a double _3d with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
  long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
  double ***t;

  /* allocate pointers to pointers to rows */
  t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
  if (!t) nrerror("allocation failure 1 in alloc3d_d()","");
  t += NR_END;
  t -= nrl;

  /* allocate pointers to rows and set pointers to them */
  t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
  if (!t[nrl]) nrerror("allocation failure 2 in alloc3d_d()","");
  t[nrl] += NR_END;
  t[nrl] -= ncl;

  /* allocate rows and set pointers to them */
  t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
  if (!t[nrl][ncl]) nrerror("allocation failure 3 in alloc3d_d()","");
  t[nrl][ncl] += NR_END;
  t[nrl][ncl] -= ndl;

  for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
  for(i=nrl+1;i<=nrh;i++) {
    t[i]=t[i-1]+ncol;
    t[i][ncl]=t[i-1][ncl]+ncol*ndep;
    for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
  }

  /* return pointer to array of pointers to rows */
  return t;
}

/* ================================================================== */
void free3d_d(double ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* free a double _3d allocated by alloc3d_d() */
{
  free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
  free((FREE_ARG) (t[nrl]+ncl-NR_END));
  free((FREE_ARG) (t+nrl-NR_END));
}

/* ================================================================== */
LONG l2d_f = 0;
float **alloc2d_f(long nrl, long nrh, long ncl, long nch)
/* allocate a float _2d with subscript range m[nrl..nrh][ncl..nch] */
{
  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  float **m;

  /* allocate pointers to rows */
  m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
  if (!m) nrerror("allocation failure 1 in alloc2d_f()","");
  m += NR_END;
  m -= nrl;

  /* allocate rows and set pointers to them */
  m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
  if (!m[nrl]) nrerror("allocation failure 2 in alloc2d_f()","");
  m[nrl] += NR_END;
  m[nrl] -= ncl;

  for(i=nrl+1;i<=nrh;i++)
    m[i]=m[i-1]+ncol;

  /* return pointer to array of pointers to rows */
  l2d_f++;
  return m;
}

/* ================================================================== */
void free2d_f(float **m, long nrl, long nrh, long ncl, long nch)
/* free a float _2d allocated by alloc2d_f() */
{
  l2d_f--;
  free((FREE_ARG) (m[nrl]+ncl-NR_END));
  free((FREE_ARG) (m+nrl-NR_END));
}

/* ================================================================== */
LONG l2d_l = 0;
LONG **alloc2d_l(long nrl, long nrh, long ncl, long nch)
/* allocate a LONG _2d with subscript range m[nrl..nrh][ncl..nch] */
{
  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  LONG **m;

  /* allocate pointers to rows */
  m=(LONG **) malloc((size_t)((nrow+NR_END)*sizeof(LONG*)));
  if (!m) nrerror("allocation failure 1 in alloc2d_l()","");
  m += NR_END;
  m -= nrl;

  /* allocate rows and set pointers to them */
  m[nrl]=(LONG *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(LONG)));
  if (!m[nrl]) nrerror("allocation failure 2 in alloc2d_l()","");
  m[nrl] += NR_END;
  m[nrl] -= ncl;

  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

  /* return pointer to array of pointers to rows */
  l2d_l++;
  return m;
}

/* ================================================================= */
void free2d_l(LONG **m, long nrl, long nrh, long ncl, long nch)
/* free a LONG _2d allocated by alloc2d_l() */
{
  l2d_l--;
  free((FREE_ARG) (m[nrl]+ncl-NR_END));
  free((FREE_ARG) (m+nrl-NR_END));
}

/* ================================================================== */
LONG l2d_s = 0;
short **alloc2d_s(long nrl, long nrh, long ncl, long nch)
/* allocate a short _2d with subscript range m[nrl..nrh][ncl..nch] */
{
  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  short **m;

  /* allocate pointers to rows */
  m=(short **) malloc((size_t)((nrow+NR_END)*sizeof(short*)));
  if (!m) nrerror("allocation failure 1 in alloc2d_s()","");
  m += NR_END;
  m -= nrl;

  /* allocate rows and set pointers to them */
  m[nrl]=(short *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(short)));
  if (!m[nrl]) nrerror("allocation failure 2 in alloc2d_s()","");
  m[nrl] += NR_END;
  m[nrl] -= ncl;

  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

  /* return pointer to array of pointers to rows */
  l2d_s++;
  return m;
}

/* ================================================================= */
void free2d_s(short **m, long nrl, long nrh, long ncl, long nch)
/* free a short _2d allocated by alloc2d_s() */
{
  l2d_s--;
  free((FREE_ARG) (m[nrl]+ncl-NR_END));
  free((FREE_ARG) (m+nrl-NR_END));
}

/* ================================================================== */
void free2d_d(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double _2d allocated by alloc2d_d() */
{
  free((FREE_ARG) (m[nrl]+ncl-NR_END));
  free((FREE_ARG) (m+nrl-NR_END));
}

/* ================================================================= */
LONG l1d_f = 0;
float *alloc1d_f(long nl, long nh)
/* allocate a float _1d with subscript range v[nl..nh] */
{
  float *v;

  v=(float *) malloc((size_t) (nh-nl+1+NR_END)*sizeof(float) );
  if (!v)
    nrerror("allocation failure in alloc1d_f()","");

  l1d_f++;
  return v-nl+NR_END;
}
/* ================================================================== */
void free1d_f(float *v, long nl, long nh)
/* free a float _1d allocated with alloc1d_f() */
{
  l1d_f--;
  free((FREE_ARG) (v+nl-NR_END));
}

/* ================================================================== */
LONG l1d_uc = 0;
unsigned char *alloc1d_uc(long nl, long nh)
/* allocate an unsigned char _1d with subscript range v[nl..nh] */
{
  unsigned char *v;

  v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
  if (!v)
    nrerror("allocation failure in alloc1d_uc()","");
  l1d_uc++;
  return v-nl+NR_END;
}

/* ================================================================== */
void free1d_uc(unsigned char *v, long nl, long nh)
/* free an unsigned char _1d allocated with alloc1d_uc() */
{
  l1d_uc--;
  free((FREE_ARG) (v+nl-NR_END));
}

/* ================================================================== */
LONG l1d_us = 0;
unsigned short *alloc1d_us(long nl, long nh)
/* allocate an unsigned short _1d with subscript range v[nl..nh] */
{
  unsigned short *v;

  v=(unsigned short *) malloc((size_t) (nh-nl+1+NR_END)*sizeof(unsigned short) );
  if (!v) nrerror("allocation failure in alloc1d_us()","");

  l1d_us++;
  return v-nl+NR_END;
}

/* ================================================================= */
void free1d_us(unsigned short *v, long nl, long nh)
/* free an unsigned short _1d allocated with alloc1d_us() */
{
  l1d_us--;
  free((FREE_ARG) (v+nl-NR_END));
}

/* ================================================================== */
LONG l2d_by = 0;
BYTE **alloc2d_by(long nrl, long nrh, long ncl, long nch)
/* allocate a short _2d with subscript range m[nrl..nrh][ncl..nch] */
{
  long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  BYTE **m;

  /* allocate pointers to rows */
  m=(BYTE **) malloc((size_t)((nrow+NR_END)*sizeof(BYTE*)));
  if (!m) nrerror("allocation failure 1 in alloc2d_by()","");
  m += NR_END;
  m -= nrl;

  /* allocate rows and set pointers to them */
  m[nrl]=(BYTE *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(BYTE)));
  if (!m[nrl]) nrerror("allocation failure 2 in alloc2d_s()","");
  m[nrl] += NR_END;
  m[nrl] -= ncl;

  for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

  /* return pointer to array of pointers to rows */
  l2d_by++;
  return m;
}

/* ================================================================= */
void free2d_by(BYTE **m, long nrl, long nrh, long ncl, long nch)
/* free a short _2d allocated by alloc2d_by() */
{
  l2d_by--;
  free((FREE_ARG) (m[nrl]+ncl-NR_END));
  free((FREE_ARG) (m+nrl-NR_END));
}
/* ================================================================= */

void CheckDeallocation(void)
{
  /* missing l2d_d check */
  if (l2d_f)
    printf("alloc2d_f = %ld\n", l2d_f);
  if (l2d_l)
    printf("alloc2d_l = %ld\n", l2d_l);
  if (l2d_s)
    printf("alloc2d_s = %ld\n", l2d_s);
  if (l2d_by)
    printf("alloc2d_by = %ld\n", l2d_by);
  if (l1d_VAR)
    printf("alloc1d_VAR = %ld\n", l1d_VAR);
  /* missing l1d_d check */
  if (l1d_f)
    printf("alloc1d_f = %ld\n", l1d_f);
  /* missing l1d_i check */
  if (l1d_l)
    printf("alloc1d_l = %ld\n", l1d_l);
  if (l1d_s)
    printf("alloc1d_s = %ld\n", l1d_s);
  if (l1d_us)
    printf("alloc1d_us = %ld\n", l1d_us);
  if (l1d_uc)
    printf("alloc1d_uc = %ld\n", l1d_uc);
}
/* ================================================================= */
float * * VectorToMatrix (float * pfVector)
{
  int i;
  float * * ppfMatrix;

  ppfMatrix = (float * *) malloc (sizeof(float *) * nRows);
  for (i=0; i< nRows; i++)
    ppfMatrix[i] = &pfVector[i*nCols];

  return(ppfMatrix);
}
/* ================================================================= */
void FreeVectorToMatrix (float * * ppfMatrix)
{
  int i;

  for (i=0; i< nRows; i++)
    free(ppfMatrix[i]);
  free(ppfMatrix);
}
/* ================================================================= */

void freeArgsStr(float **m, long nrl, long nrh, long ncl, long nch)
{
  free((FREE_ARG) (m[nrl]+ncl-NR_END));
}

/* ================================================================= */
/* ====================== LOCAL FUNCTIONS ========================== */
/* doesn't work on the sun with gcc, don't use for now */

void ConvertToLowerCase(char * s)
{
  int i = 0;
  while (s[i])
  {
    s[i] = (char) tolower(s[i]);
    i++;
  }
}