/*
** Copyright (2009) : Chris Velden, Tim Olander, Tony Wimmers,
**                    Jim Kossin, Howard Berger, and Dave Santek
**                    except where otherwise noted.
*/

/* include file containing all AODT library global variables */
#include "../inc/odtlib.h"
/* include file containing all AODT library variable definitions */
#include "../inc/odtlibdefs-x.h"

#define  RADIANSCONSTANT 0.017453292;
#define  ANGLE90DEGRADIANS 1.570797;

#define deg_to_rad(deg)  ((deg) * PI / 180.0)
#define rad_to_deg(rad)  ((rad) * 180.0 / PI)

double aodt_calctime( int, int );
float  aodt_slopecal( double,int );
float  aodt_getpwval( int,float,float,float );
float  aodt_PWlandsea( float );
float  aodt_ptovmax( float );
float  aodt_ptotno( float );
int    aodt_cmonth2julian( char * );
void   aodt_julian2cmonth( int, char * );
void   aodt_distance( float,float,float,float,int,float *,float * );
void   aodt_distance2( float, float, float, float, float *, float * );
float  aodt_atoif( char *, int, int );
void   aodt_calcskew( float *, int, float *, float *, float * );
int    aodt_idmyyd( int,int,int );
void   aodt_yddmy( int,int *,int *,int * );
int    aodt_oceanbasin( float,float );
int    aodt_sattypes( int,char * );
int    aodt_atcffilename( char *,int,char *,char *);
int    aodt_initcurrent(int);
int    aodt_stormspeedfromhistory(float *, float *);
int    aodt_checkIRdataHIE(void);
float  aodt_XLZA( float,float,float,float );

float aodt_slopecal(double SearchTimeInterval,
                    int SlopeInterceptFlag)
/* 
** Calculate slope or y-intercept of all points over SearchTimeInterval period
** Inputs  : SearchTimeInterval - time period to calculate slope or y-intercept
**           SlopeInterceptFlag - flag value indicating slope or y-intercept
**                                calculation and parameter to utilize
**                            1 = Final T# (will return slope)
**                            2 = Adjusted Raw T# (will return slope)
**                            3 = latitude (will return y-intercept)
**                            4 = longitude (will return y-intercept)
**           HistoryCurrentPointer_Global contains current analysis information
** Outputs : None
** Return  : Slope or Y-Intercept value of line over time period desired
*/
{
  int Counter=0;                        /* counter value */
  int CounterMinimum=4;                 /* counter minimum value */
  float SumX=0.0;                       /* X-axis sum value */
  float SumY=0.0;                       /* Y-axis sum value */
  float SumSquaresX=0.0;                /* X-axis sum squares value */
  float SumSquaresY=0.0;                /* Y-axis sum squares value */
  float SumXY=0.0;                      /* X*Y sum value */
  float XAxisValue;                     /* X-axis value */
  float YAxisValue;                     /* Y-axis value */
  float XMean;                          /* X-axis mean value */
  float YMean;                          /* Y-axis mean value */
  float XVariance;                      /* X-axis variance value */
  float YVariance;                      /* Y-axis variance value */
  float CovarianceXY;                   /* X*Y covariance value */
  float RValue;                         /* r (correlation?) value */
  float YIntercept;                     /* Y-intercept value */
  float Slope;                          /* slope value */
  float SlopeYIntReturnValue;           /* slope/y-intercept return value */
  float FirstLonValue;                  /* 1st lon value for extrap */
  double CurrentTime;                   /* current time */
  double HistoryRecTime;                /* history record time */
  double TimeThreshold;                 /* time threshold value */
  logical FoundValidRecordTF=FALSE;     /* found valid record logical */
  logical OverWaterTF;                  /* TC over water logical value */
  logical FirstLonValueTF=TRUE;         /* 1st lon value for extrap logical */
  struct HistoryFileRec_Structure *HistoryRecordPtr; /* history structure */

  HistoryRecordPtr=HistoryFirstRecordPointer_Global;

  CurrentTime=aodt_calctime(HistoryCurrentPointer_Global->IR.date,
                            HistoryCurrentPointer_Global->IR.time);
  TimeThreshold=CurrentTime-(SearchTimeInterval/24.0);

  while(HistoryRecordPtr!=0) {
     HistoryRecTime=aodt_calctime(HistoryRecordPtr->IR.date,
                                  HistoryRecordPtr->IR.time);
     if((HistoryRecTime<CurrentTime)&&(HistoryRecTime>=TimeThreshold)) {
        OverWaterTF=TRUE;
        if(SlopeInterceptFlag<3) {
           if(((LandFlagTF_Global)&&(HistoryRecordPtr->IR.land==1))||
              (HistoryRecordPtr->IR.Traw<1.0)) {
              OverWaterTF=FALSE;
           }
        }
/*
        if(((LandFlagTF_Global)&&(HistoryRecordPtr->IR.land==1))||
           (HistoryRecordPtr->IR.Traw<1.0)) { 
              OverWaterTF=FALSE;
        } 
*/
        if(OverWaterTF) {
           XAxisValue=(float)(CurrentTime-HistoryRecTime);
           if(SlopeInterceptFlag==1) {
              YAxisValue=HistoryRecordPtr->IR.Tfinal;
           }
           if(SlopeInterceptFlag==2) {
              YAxisValue=HistoryRecordPtr->IR.Traw;
           }
           if(SlopeInterceptFlag==3) {
              YAxisValue=HistoryRecordPtr->IR.latitude;
           }
           if(SlopeInterceptFlag==4) {
              YAxisValue=HistoryRecordPtr->IR.longitude;
              if(FirstLonValueTF) {
                 FirstLonValue=YAxisValue;
                 FirstLonValueTF=FALSE;
              } 
              else { 
                 if((FirstLonValue>100.0)&&(YAxisValue<-100.0)) {
                    /* dateline cross W to E */
                    YAxisValue=YAxisValue+360.0;
                 }
                 if((FirstLonValue<-100.0)&&(YAxisValue>100.0)) {
                    /* dateline cross E to W */
                    YAxisValue=YAxisValue-360.0;
                 }
              }
           }
           SumX=SumX+XAxisValue;
           SumY=SumY+YAxisValue;
           SumSquaresX=SumSquaresX+(XAxisValue*XAxisValue);
           SumSquaresY=SumSquaresY+(YAxisValue*YAxisValue);
           SumXY=SumXY+(YAxisValue*XAxisValue);
           Counter++;
           FoundValidRecordTF=TRUE;
        }
     }
     else {
        if(FoundValidRecordTF) {
           break;
        }
     }
     HistoryRecordPtr=HistoryRecordPtr->nextrec;
  }
  /* if calculating slope of Final T# values, add in current value */
  if(SlopeInterceptFlag<=2) {
     CounterMinimum=6;
     /* add current record to slope calculation */
     if(SlopeInterceptFlag==1) {
        YAxisValue=HistoryCurrentPointer_Global->IR.Tfinal;
     }
     if(SlopeInterceptFlag==2) {
        YAxisValue=HistoryCurrentPointer_Global->IR.Traw;
     }
     SumY=SumY+YAxisValue;
     SumSquaresY=SumSquaresY+(YAxisValue*YAxisValue);
     Counter++;
  }

  /* 
  ** compute least squares fit of line for data points
  ** using the equation  Y = Y* + r(varX/varY)(X - X*) = Mx + B
  ** X* = mean of X values (time values)  = XMean
  ** Y* = mean of Y values (T# values)    = YMean
  ** varX = variance of X                 = XVariance
  ** varY = variance of Y                 = YVariance
  ** r = covariance - Y*X*                = RValue
  ** M = slope of line (desired value)    = RValue*(sqrt(YVariance/XVariance))
  ** B = y-intercept                      = YMean-(slopecal*XMean)
  */

  /* must have more than CounterMinimum data values to calculate slope */
  if(Counter<CounterMinimum) {
     if((SlopeInterceptFlag==3)||(SlopeInterceptFlag==4)) {
        /* Slope=999.99; */
        Slope=MISSFP1;
     }
     else {
        Slope=0.0;
     }
     SlopeYIntReturnValue=Slope;
  }
  else {
     XMean=SumX/(float)Counter;
     YMean=SumY/(float)Counter;
     XVariance=(SumSquaresX/(float)Counter)-(XMean*XMean);
     YVariance=(SumSquaresY/(float)Counter)-(YMean*YMean);
     CovarianceXY=(SumXY/(float)Counter)-(XMean*YMean);
     RValue=CovarianceXY/A_SQRT(XVariance*YVariance);
     if((A_ABS(XVariance)<=0.0001)||(A_ABS(YVariance)<=0.0001)) {
        Slope=0.0;
     } 
     else {
        Slope=RValue*(A_SQRT(YVariance/XVariance));
     }
     Slope=(float)((int)(Slope*10.0))/10.0;
     if((SlopeInterceptFlag==3)||(SlopeInterceptFlag==4)) {
        YIntercept=YMean-(Slope*XMean);
        if(SlopeInterceptFlag==4) {
           if(YIntercept<-180.0) YIntercept=YIntercept+360.0;
           if(YIntercept>180.0) YIntercept=YIntercept-360.0;
        }
        /* y-intercept for latitude/longitude extrapolation */
        SlopeYIntReturnValue=YIntercept;
     }
     else {
        /* slope for Final T# slope calculation */
        SlopeYIntReturnValue=Slope;
     }
  }
  
  return SlopeYIntReturnValue;

}

double aodt_calctime(int InputJulianDate,
                     int InputHMSTime)
/* 
** Compute time in ADT xxxxx.yyy units, where xxxxx is the
** day and yyy is the percentage of the day.  This routine
** will also correct for Y2K problems.
** Inputs  : InputJulianDate - Julian date
**           InputHMSTime    - time in HHMMSS format
** Outputs : None 
** Return  : Date/Time value in xxxxx.yyy units 
*/
{
  int YearValue;                        /* year value */
  float SecondsValue;                   /* seconds value */
  float MinutesValue;                   /* minutes value */
  float HoursValue;                     /* hours value */
  float PartialDateValue;               /* part day value */
  double TimeReturnValue;               /* time return value */

  if(((InputJulianDate%1000)==0)||(InputHMSTime<0)) {
     TimeReturnValue=0.0;
  }
  else {
     YearValue=InputJulianDate/1000;  /* obtain year */
     /* check for millenium designation in the year.
        if it is not there, add it onto the beginning */
     if(YearValue<1900) {
        if(YearValue>70) {
           InputJulianDate=1900000+InputJulianDate;
        }
        else {
           InputJulianDate=2000000+InputJulianDate;
        }
     }
     SecondsValue=((float)(InputHMSTime%100))/3600.0;
     MinutesValue=((float)((InputHMSTime/100)%100))/60.0;
     HoursValue=(float)(InputHMSTime/10000);
     PartialDateValue=(HoursValue+MinutesValue+SecondsValue)/24.0;
     TimeReturnValue=(double)InputJulianDate+(double)PartialDateValue;
  }

  return TimeReturnValue;

}

float aodt_getpwval(int PressureWindIDValue,
                    float CIValue,
                    float Latitude_Input,
                    float Longitude_Input)
/* 
** Obtain pressure or wind value (for Atlantic or
** West Pacific storms) given the intensity estimate
** value.
** Inputs  : PressureWindIDValue  - flag for wind (1) or pressure (0) output
**           CIValue              - Current Intensity (CI) value
**           Latitude_Input       - Latitude value of analysis
**           Longitude_Input      - Longitude value of analysis
** Outputs : None
** Return  : Pressure or Wind Speed value
*/
{
  int XInc=2;                           /* X direction increment */
  int DomainID_Local;                   /* domain ID local value */
  int RetErr;                           /* fuction error return variable */
  int BasinID_Local;                    /* local basin ID value */
  float PWReturnValue;                  /* pressure/wind return value */
  float Vmax_Local;                     /* Max Wind Speed value from T# */
  float S;                              /* Normalized size parameter */
  float Vstorm_Local;                   /* Vmax adj for storm speed */
  float GaleRadius34_Climo;             /* Climatological Gale Radius */
  float StormSpeed;                     /* Storm Speed from history file */
  float StormDirection;                 /* Storm Direction from history file */
  float EXP_Local;
  float R34_Local;
  float ROCI_Local;
  float Vstorm_500;
  float Vstorm_500c;

  /* use traditional MSLP/Wind-T# conversion */
  /* determine correct pressure/wind array bin */
  while((CIValue>PW_TnoValues[XInc])&&(XInc<82)) { XInc++; }

  DomainID_Local=(DomainID_Global%2);
  /* convert CI value to wind/pressure value */
  if(PressureWindIDValue==1) {
     PWReturnValue=PW_WindValues[XInc];                      /* WIND */
  }
  else {
     PWReturnValue=PW_PressureValues[DomainID_Local][XInc];  /* PRESSURE */
     if(UseCKZTF_Global==1) {
        /* 
        ** use new Knaff/Zehr T#->MSLP conversion.
        ** Need Vmax calculation and two command line input values
        */
        BasinID_Local=aodt_oceanbasin(Latitude_Input,Longitude_Input);
        if(BasinID_Local==2) {
           GaleRadius34_Climo=82.0;
        }  
        else {
           GaleRadius34_Climo=107.0;
        }
        /* RetErr=aodt_stormspeedfromhistory(&StormSpeed,&StormDirection); */
        StormSpeed=11.0;  /* 11 knots -- climo value */
        /* printf("CKZGaleRadius34_Global=%f\n",CKZGaleRadius34_Global); */

        Vmax_Local=PW_WindValues[XInc];                    /* WIND */
        if(CKZGaleRadius34_Global<0.0) {
           ROCI_Local=A_ABS(CKZGaleRadius34_Global);
           R34_Local=(0.354*ROCI_Local)+13.3;
        }
        else {
           R34_Local=CKZGaleRadius34_Global;
        }
        /* printf("R34_local=%f\n",R34_Local); */
        /* printf("Latitude_Input=%f\n",Latitude_Input); */
        Vstorm_Local=Vmax_Local-(1.5*A_POW(StormSpeed,0.63));
        /* printf("Vstorm_local=%f\n",Vstorm_Local); */
        EXP_Local=0.1147+(0.0055*Vstorm_Local)-
                  (0.001*(A_ABS(Latitude_Input)-25.0));
        /* printf("EXP_local=%f\n",EXP_Local); */
        Vstorm_500=(R34_Local/9.0)-3.0;
        /* printf("Vstorm_500=%f\n",Vstorm_500); */
        Vstorm_500c=Vstorm_Local*A_POW(((66.785-(0.09102*Vstorm_Local)+
                    (1.0619*(A_ABS(Latitude_Input)-25.0)))/500.0),EXP_Local);
        /* printf("Vstorm_500c=%f\n",Vstorm_500c); */
        S=A_MAX((Vstorm_500/Vstorm_500c),0.4);
        /* printf("S=%f\n",S); */
        PWReturnValue=23.286-(0.483*Vstorm_Local)-
                      A_POW((Vstorm_Local/24.254),2)-
                      (12.587*S)-(0.483*A_ABS(Latitude_Input))+CKZPenv_Global;
        if(PWReturnValue>=(CKZPenv_Global-1.0)) {
           PWReturnValue=CKZPenv_Global-2.0;
        }
        /* printf("PWReturnValue=%f\n",PWReturnValue); */
        HistoryCurrentPointer_Global->IR.r34=(int)CKZGaleRadius34_Global;
        HistoryCurrentPointer_Global->IR.MSLPenv=(int)CKZPenv_Global;
     } 
  }
   
  return PWReturnValue;

}

float aodt_PWlandsea( float PressureInputValue )
/* 
** obtain Landsea Vmax value given mslp value for a given 
** storm at a given latitude.
** lat < 25N : V = 12.016 * (1013 - P0)**(0.5337)
** 25-35N    : V = 14.172 * (1013 - P0)**(0.4778)
** 35-45N    : V = 16.086 * (1013 - P0)**(0.4333)
**
** Inputs  : PressureInputValue - input pressure value (in mb)
** Outputs : None
** Return  : Maximum Wind Speed (Vmax) in knots
*/
{
  int XInc=0;                           /* X direction increment */
  float MaxWindSpeedReturnValue;        /* max wind value return value */
  float LatitudeValue;                  /* latitude value */
  float MultFactorArray[3]={12.016,14.172,16.086}; /* mult factor array */
  float ExpFactorArray[3]={0.5337,0.4778,0.4333};  /* EXP factor array */

  LatitudeValue=HistoryCurrentPointer_Global->IR.latitude;
  if(LatitudeValue>=25.0) {
     XInc=1;
  }
  if(LatitudeValue>=35.0) {
     XInc=2;
  }

  MaxWindSpeedReturnValue=MultFactorArray[XInc]*
                          (A_POW((1013.0-PressureInputValue),
                                  ExpFactorArray[XInc]));

  return MaxWindSpeedReturnValue;
 
}

float aodt_ptovmax( float PressureInputValue )
/* 
** Derive wind speed (Vmax) value from pressure (MSLP) value
** Inputs  : Pressure value
** Outputs : None
** Return  : Maximum Wind Speed (Vmax) value 
*/
{
  int XInc=2;                           /* X direction increment */
  int DomainID_Local;                   /* domain ID local value */
  float MaxWindSpeedReturnValue;        /* max wind speed return value */
  float PartialValue;                   /* part day value */

  DomainID_Local=(DomainID_Global%2);
  while((PressureInputValue<PW_PressureValues[DomainID_Local][XInc])&&
        (XInc<82)) {
     XInc++;
  }

  if(XInc==2) {
     MaxWindSpeedReturnValue=PW_WindValues[2];
  }
  else {
     PartialValue=(PW_PressureValues[DomainID_Local][XInc-1]-
                   PressureInputValue)/
                  (PW_PressureValues[DomainID_Local][XInc-1]-
                   PW_PressureValues[DomainID_Local][XInc]);
     MaxWindSpeedReturnValue=PW_WindValues[XInc-1]+
                             ((PW_WindValues[XInc]-
                               PW_WindValues[XInc-1])*PartialValue);
  }

  return MaxWindSpeedReturnValue;

}

float aodt_ptotno( float PressureInputValue )
/* 
** derive T# value from pressure (MSLP) value
** Inputs  : pressure value
** Outputs : None
** Return  : T# value 
*/
{
  int XInc=2;                           /* X direction increment */
  int DomainID_Local;                   /* domain ID local value */
  float TNoReturnValue;                 /* T# return value */
  float PartialValue;                   /* partial T# valuee */

  DomainID_Local=(DomainID_Global%2);
  while((PressureInputValue<PW_PressureValues[DomainID_Local][XInc])&&
        (XInc<82)) {
     XInc++;
  }

  if(XInc==2) {
     TNoReturnValue=PW_TnoValues[2];
  } 
  else {
     PartialValue=(PW_PressureValues[DomainID_Local][XInc-1]-
                   PressureInputValue)/
                  (PW_PressureValues[DomainID_Local][XInc-1]-
                   PW_PressureValues[DomainID_Local][XInc]);
     TNoReturnValue=PW_TnoValues[XInc-1]+
                    ((PW_TnoValues[XInc]-
                      PW_TnoValues[XInc-1])*PartialValue);
  }

  return TNoReturnValue;

}

int aodt_cmonth2julian( char *InputDateCharString )
/* 
** Convert YYYYMMMDD format character string value to julian date Integer.
** Inputs  : InputDateCharString - character representation of date in 
**                                 YYYYmonDD format where DD and YYYY are 
**                                 integers and mon is a three character 
**                                 abbreviation for the month (e.g. 2000MAR13)
** Outputs : None
** Return  : Julian Date
*/
{
  int  DayValue;                        /* day value */
  int  YearValue;                       /* year value */
  int  MonthValue=0;                    /* month value */
  int  JulianDateReturnValue=0;         /* julian date return value */
  char MonthCharString[4]="\0";         /* month character string value */

  /* break down character string into various components */
  (void)sscanf(InputDateCharString,"%4d%3s%2d",
               &YearValue,MonthCharString,&DayValue);

  /* calculate integer month value */
  while((strncmp(MonthCharString,MonthID_String[MonthValue],3)!=0)&&
        (MonthValue<12)) {
     MonthValue++;
  }

  /* determine julian date from day/month/year */
  if(MonthValue<12) {
     JulianDateReturnValue=aodt_idmyyd(DayValue,MonthValue+1,YearValue);
  }
  else {
     JulianDateReturnValue=0.0;   /* Error */
  }

  return JulianDateReturnValue;

}

void aodt_julian2cmonth( int JulianDateInputValue,
                         char *ReturnDateCharString )
/* 
** convert julian date to YYYYMMMDD format for output.
** Inputs  : JulianDateInputValue - Julian date representation of date
** Outputs : ReturnDateCharString - character representation of date 
**                                  in YYYYmonDD format where DD and YYYY 
**                                  are integers and mon is a three
**           character abbreviation for the month (e.g. 2000MAR13)
** Return  : None
*/
{
  int DayValue;                         /* day value */
  int MonthValue;                       /* month value */
  int YearValue;                        /* year value */

  /* calculate date/month/year from julian date */
  (void)aodt_yddmy(JulianDateInputValue,&DayValue,&MonthValue,&YearValue);

  /* form character string representation from various components */
  (void)snprintf(ReturnDateCharString,10,"%4d%3s%2.2d",YearValue,
                MonthID_String[MonthValue-1],DayValue);

}

void aodt_distance(float StartLatitudeInput,
                   float StartLongitudeInput,
                   float EndLatitudeInput,
                   float EndLongitudeInput,
                   int UnitFlagID,
                   float *Distance_Return,
                   float *Angle_Return)
/* 
** Calculate distance and angle between two points 
** (StartLatitudeInput,StartLongitudeInput 
**  and EndLatitudeInput,EndLongitudeInput).
** Inputs  : StartLatitudeInput  - latitude of starting point
**           StartLongitudeInput - longitude of starting point
**           EndLatitudeInput    - latitude of ending point
**           EndLongitudeInput   - longitude of ending point
**           UnitFlagID          - flag for output unit type (1-km,2-mi,3-nmi)
** Outputs : Distance_Return     - distance between two points
**           Angle_Return        - angle between two points (0=north)
** Return  : None
*/
{
  float StartLatitude_Radians;          /* start point latitude in radians */
  float StartLongitude_Radians;         /* start point longitude in radians */
  float EndLatitude_Radians;            /* end point latitude in radians */
  float EndLongitude_Radians;           /* end point longitude in radians */
  float StartLatitude_Radians_ACOS;     /* start point lat ACOS value radians */
  float StartLongitude_Radians_ACOS;    /* start point lon ACOS value radians */
  float StartLatitude_Radians_ASIN;     /* start point lat ASIN value radians */
  float StartLongitude_Radians_ASIN;    /* start point lon ASIN value radians */
  float EndLatitude_Radians_ACOS;       /* end point lat ACOS value radians */
  float EndLongitude_Radians_ACOS;      /* end point lon ACOS value radians */
  float EndLatitude_Radians_ASIN;       /* end point lat ASIN value radians */
  float EndLongitude_Radians_ASIN;      /* end point lon ASIN value radians */
  float AdditiveValue;                  /* local additive value */
  float COSCOS_Difference;              /* end-start COS*COS values diff */
  float SINCOS_Difference;              /* end-start COS*SIN values diff */
  float LatitudeSIN_Difference;         /* latitude SIN value difference */
  float Distance_Intermediate;          /* intermediate distance value */
  float Distance_Final;                 /* final distance value */
  float Angle_Final;                    /* final angle value */

  Distance_Intermediate=0.0;
  StartLatitude_Radians=StartLatitudeInput*RADIANSCONSTANT;
  StartLongitude_Radians=StartLongitudeInput*RADIANSCONSTANT;
  EndLatitude_Radians=EndLatitudeInput*RADIANSCONSTANT;
  EndLongitude_Radians=EndLongitudeInput*RADIANSCONSTANT;
  StartLatitude_Radians_ACOS=A_COS(StartLatitude_Radians);
  StartLongitude_Radians_ACOS=A_COS(StartLongitude_Radians);
  StartLatitude_Radians_ASIN=A_SIN(StartLatitude_Radians);
  StartLongitude_Radians_ASIN=A_SIN(StartLongitude_Radians);
  EndLatitude_Radians_ACOS=A_COS(EndLatitude_Radians);
  EndLongitude_Radians_ACOS=A_COS(EndLongitude_Radians);
  EndLatitude_Radians_ASIN=A_SIN(EndLatitude_Radians);
  EndLongitude_Radians_ASIN=A_SIN(EndLongitude_Radians);
  COSCOS_Difference=(EndLatitude_Radians_ACOS*EndLongitude_Radians_ACOS)-
                    (StartLatitude_Radians_ACOS*StartLongitude_Radians_ACOS);
  SINCOS_Difference=(EndLatitude_Radians_ACOS*EndLongitude_Radians_ASIN)-
                    (StartLatitude_Radians_ACOS*StartLongitude_Radians_ASIN);
  LatitudeSIN_Difference=EndLatitude_Radians_ASIN-StartLatitude_Radians_ASIN;
  AdditiveValue=(COSCOS_Difference*COSCOS_Difference)+
                (SINCOS_Difference*SINCOS_Difference)+
                (LatitudeSIN_Difference*LatitudeSIN_Difference);
  Distance_Intermediate=A_SQRT(AdditiveValue);

  /* Distance_Final is distance in kilometers */
  Distance_Final=2.0*A_ASIN(Distance_Intermediate/2.0)*EARTHRADIUSKM;

  if(UnitFlagID==2) {
     /* Conversion to Miles */
     Distance_Final=((69.0*Distance_Final)+55)/111.0;
  }
  if(UnitFlagID==3) {
     /* Conversion to Nautical Miles */
     Distance_Final=((60.0*Distance_Final)+55)/111.0;
  }
  *Distance_Return=Distance_Final;

  /* Compute Final Angle */
  if(A_ABS(Distance_Final)>0.0001) { 
     Angle_Final=(A_SIN(StartLongitude_Radians-EndLongitude_Radians)*
                  A_SIN((PI/2.0)-EndLatitude_Radians))/
                  A_SIN(Distance_Intermediate);
  }
  else {
     Angle_Final=0.0;
  }
  if(A_ABS(Angle_Final)>1.0) {
     Angle_Final=A_SIGN(1.000,Angle_Final);
  }
  Angle_Final=A_ASIN(Angle_Final)/RADIANSCONSTANT;
  if(EndLatitude_Radians<StartLatitude_Radians) {
     Angle_Final=180.0-Angle_Final;
  }
  if(Angle_Final<0.0) {
     Angle_Final=360.0+Angle_Final;
  }
  *Angle_Return=Angle_Final;
}

void aodt_distance2(float StartLatitudeInput,
                    float StartLongitudeInput,
                    float DistanceInput,
                    float AngleInput,
                    float *EndLatitude_Return,
                    float *EndLongitude_Return)
/* 
** Calculate a latitude and longitude position from an 
** initial latitude/longitude and distance/angle values.
** Inputs  : StartLatitudeInput  - initial latitude
**           StartLongitudeInput - initial longitude
**           DistanceInput       - distance from initial position
**           AngleInput          - angle from initial position
** Outputs : EndLatitude_Return  - derived latitude
**           EndLongitude_Return - derived longitude
** Return  : None
*/
{
  int AngleInput_Integer;               /* input angle value */
  float StartLatitude_Radians;          /* start point latitude in radians */
  float StartLongitude_Radians;         /* start point longitude in radians */
  float StartLatitude_Radians_Flipped;  /* start point 90-X lat in radians */
  float Distance_Radians;               /* distance value in radians */
  float Angle_Radians;                  /* angle value in radians */
  float Latitude_AACOS;                 /* muliplicative latitude value */
  float Longitude_AASIN;                /* muliplicative longitude value */
  float ArgumentValue;                  /* local argument value */
  float ATANValue;                      /* Arc Tangent value */
  float Angle90Radians_Local=ANGLE90DEGRADIANS; /* 90degree angle in radians */
  float RadiansConstant_Local=RADIANSCONSTANT; /* radians constant value */

  StartLatitude_Radians=(90.0-StartLatitudeInput)*RADIANSCONSTANT;
  StartLatitude_Radians_Flipped=StartLatitude_Radians;
  StartLongitude_Radians=StartLongitudeInput*RADIANSCONSTANT;
  if(StartLatitudeInput<0.0) {
    StartLatitude_Radians_Flipped=-(90.0+StartLatitudeInput)*RADIANSCONSTANT;
    StartLongitude_Radians=(StartLongitudeInput-180.0)*RADIANSCONSTANT;
    AngleInput=360-AngleInput;
  }
  AngleInput_Integer=(int)AngleInput;
  Angle_Radians=-1.0*((float)((540-AngleInput_Integer)%360))*RADIANSCONSTANT;
  Distance_Radians=(DistanceInput/111.1)*RADIANSCONSTANT;
  Latitude_AACOS=A_ACOS((A_COS(StartLatitude_Radians)*
                         A_COS(Distance_Radians))+
                        (A_SIN(StartLatitude_Radians)*
                         A_SIN(Distance_Radians)*A_COS(Angle_Radians)));
  if(A_ABS(Latitude_AACOS)<0.0000001) {
    Longitude_AASIN=0.0;
  } 
  else {
     ArgumentValue=(A_SIN(Distance_Radians)*A_SIN(Angle_Radians))/
                    A_SIN(Latitude_AACOS);
     if(A_ABS(ArgumentValue)>1.0) {
        ArgumentValue=A_SIGN(1.0,ArgumentValue);
     }
     Longitude_AASIN=A_ASIN(ArgumentValue);
     ATANValue=(A_ATAN(A_SIN(Angle90Radians_Local-Angle_Radians)))/
                      (A_TAN(Angle90Radians_Local-Distance_Radians));
     if(ATANValue>StartLatitude_Radians_Flipped) {
        Longitude_AASIN=(2.0*Angle90Radians_Local)-Longitude_AASIN;
     }
  }
  Longitude_AASIN=StartLongitude_Radians-Longitude_AASIN;
  *EndLatitude_Return=90.0-(Latitude_AACOS/RadiansConstant_Local);
  *EndLongitude_Return=(float)
                      ((int)(10000*(Longitude_AASIN/
                                    RadiansConstant_Local))%3600000)/10000.0;
  if(*EndLongitude_Return<-180) {
     *EndLongitude_Return=*EndLongitude_Return+360;
  }

}

float aodt_atoif(char *StringPointer,
                 int BeginningPosition,
                 int EndingPosition)
/*
** routine to convert from character to floating point 
** Inputs  : StringPointer     - pointer to beginning of character string
**           BeginningPosition - beginning character to convert
**           EndingPosition    - ending character to convert
** Outputs : None
** Return  : converted floating point value
*/
{
  int XInc;                             /* X direction increment */
  int CharASCIIValue;                   /* character ASCII integer value */
  int MultFactor=1;                     /* multiplier value */
  float FinalValue;                     /* local value */
  float MultFactorTens=1.0;             /* tens multiplier value */
  float ConvertValue=0.0;               /* conversion value */
  float Denominator=1.0;                /* denominator value */
  logical FoundValidTF=FALSE;           /* valid number found logical */

  for(XInc=EndingPosition-1;XInc>=BeginningPosition-1;XInc--) {
     if((StringPointer[XInc]=='S')||(StringPointer[XInc]=='E')) {
        MultFactor=-1;
     }
     CharASCIIValue=(int)(StringPointer[XInc])-48;
     if((CharASCIIValue>=0)&&(CharASCIIValue<=9)) {
        FoundValidTF=TRUE;
        ConvertValue=ConvertValue+((float)CharASCIIValue*MultFactorTens);
        MultFactorTens=MultFactorTens*10.0;
     } 
     else if(CharASCIIValue==-2) {
        Denominator=MultFactorTens;
     } 
     else {
        continue;
     }
  }

  if(FoundValidTF) {
    FinalValue=MultFactor*(ConvertValue/Denominator);
  }
  else {
    /* FinalValue=-999999.9; */
    FinalValue=MISSFN1;
  }

  return FinalValue;

}

void aodt_calcskew( float *DataArray_Input,
                    int DataArraySize,
                    float *ArrayAverage_Return,
                    float *ArrayStdDev_Return,
                    float *ArraySkew_Return )
/* Calculate average, standard deviation, and skew for a given data set.
    Inputs  : DataArray_Input  - data array
              DataArraySize - number of points in data array
    Outputs : ArrayAverage_Return  - average value in array
              ArrayStdDev_Return - standard deviation of data array
              ArraySkew_Return - skew of data array histogram
*/
{
  int XInc;                             /* X direction increment */
  float Average_Local;                  /* local average value */
  float StdDev_Local;                   /* local standard deviation value */
  float Skew_Local;                     /* local skew value */
  float DifferenceValue;                /* local difference value */
  float ArrayValuesSumSquared=0.0;      /* input array sum squared value */
  float ArrayValuesSumCubed=0.0;        /* input array sum cubed value */
  float ArrayValuesSum=0.0;             /* input array sum value */

  for(XInc=0;XInc<DataArraySize;XInc++) {
     ArrayValuesSum=ArrayValuesSum+DataArray_Input[XInc];
  }
  /* compute average value of data array */
  Average_Local=ArrayValuesSum/(float)DataArraySize;

  for(XInc=0;XInc<DataArraySize;XInc++) {
     DifferenceValue=DataArray_Input[XInc]-Average_Local;
     ArrayValuesSumSquared=ArrayValuesSumSquared+
                           (DifferenceValue*DifferenceValue);
     ArrayValuesSumCubed=ArrayValuesSumCubed+
                         (DifferenceValue*DifferenceValue*DifferenceValue);
  }
  if(DataArraySize<=1) {
     StdDev_Local=0.0;
     Skew_Local=0.0;
  } 
  else {
     /* calculate standard deviation of data array */
     StdDev_Local=sqrt((1.0/((float)DataArraySize-1.0))*ArrayValuesSumSquared);
     /* calculate skew of data array */
     Skew_Local=((1.0/((float)DataArraySize-1.0))*ArrayValuesSumCubed)/
                       (StdDev_Local*StdDev_Local*StdDev_Local);
  }

  /* return desired values */
  *ArrayAverage_Return=Average_Local;
  *ArrayStdDev_Return=StdDev_Local;
  *ArraySkew_Return=Skew_Local;

}

int aodt_idmyyd(int InputDay,
                int InputMonth,
                int InputYear)
/* 
** Convert dd/mm/yy to yyyyddd format.
** this routine was originally taken from McIDAS
** program idmyyd.for.
** Inputs  : InputDay   - day
**           InputMonth - month (integer)
**           InputYear  - year (YY or YYYY)
** Outputs : None
** Return  : Julian date or 0 (0=bad input data)
*/
{
  int JulianDate_Return;                /* julian date return value */
  int DayValue;                         /* local day value */
  int JulianDateMonthArray[12]={0,31,59,90,120,151,181,212,243,273,304,334};
                                        /* julian date array */

  /* perform a couple quality checks for day/month */
  if(((InputDay<1)||(InputDay>31))||((InputMonth<1)||(InputMonth>12))) {
     JulianDate_Return=0;
  } 
  else {
     DayValue=InputDay+JulianDateMonthArray[InputMonth-1];
     if(InputYear<1900) {
       if(InputYear>70) {
         InputYear=1900+InputYear;
       } 
       else {
         InputYear=2000+InputYear;
       }
     }
     /* Leap year check */
     if(((InputYear%4)==0)&&(InputMonth>2)) {
        DayValue=DayValue+1;
     }
     JulianDate_Return=(InputYear*1000)+DayValue;
  }

  return JulianDate_Return;

}

void aodt_yddmy(int JulianDateInput,
                int *Day_Return,
                int *Month_Return,
                int *Year_Return)
/* 
** Convert yyyyddd to dd/mm/yy format.
** Inputs  : JulianDateInput - Julian day (yyyyddd)
** Outputs : Day_Return      - date
**           Month_Return    - month
**           Year_Return     - year (yyyy)
** Return  : None
*/
{
  int DayNumberArray[2][13]={{0,31,59,90,120,151,181,212,243,273,304,334,365},
                             {0,31,60,91,121,152,182,213,244,274,305,335,366}};
                                        /* julian date (and leap year) values */
  int YearValue_Local;                  /* local year value */
  int DayValue_Local;                   /* local day value */
  int MonthValue_Local;                 /* local month value */
  int LeapYearFlag=0;                   /* leap year flag value */

  YearValue_Local=JulianDateInput/1000;
  if(YearValue_Local<1900) {
    if(YearValue_Local>70) {
      YearValue_Local=YearValue_Local+1900;
    } 
    else {
      YearValue_Local=YearValue_Local+2000;
    }
  }
  DayValue_Local=(JulianDateInput%1000);
  if((YearValue_Local%4)==0) {
     LeapYearFlag=1;
  }
  for(MonthValue_Local=0;MonthValue_Local<13;MonthValue_Local++) {
     if(DayValue_Local<=DayNumberArray[LeapYearFlag][MonthValue_Local]) {
        break;
     }
  }
  *Month_Return=MonthValue_Local;
  *Day_Return=DayValue_Local-DayNumberArray[LeapYearFlag][MonthValue_Local-1];
  *Year_Return=YearValue_Local;

}

int aodt_sattypes(int SatelliteIDInput,
                  char *SatelliteName_Return)
/* 
** obtain satellite name given ID number
** Inputs  : SatelliteIDInput - Internal ADT satellite ID number
** Outputs : csat - Satellite character string ID
** Return  : None
*/
{
  int Ret_ID;                      /* return ID variable */

  if (SatelliteIDInput<0) {
     snprintf(SatelliteName_Return,7,"%s","MSNG");
  } else {
     snprintf(SatelliteName_Return,7,"%s"," OTHER");
     if(SatelliteIDInput==21) snprintf(SatelliteName_Return,7,"%s"," GOES1");
     if(SatelliteIDInput==23) snprintf(SatelliteName_Return,7,"%s"," GOES2");
     if(SatelliteIDInput==25) snprintf(SatelliteName_Return,7,"%s"," GOES3");
     if(SatelliteIDInput==27) snprintf(SatelliteName_Return,7,"%s"," GOES4");
     if(SatelliteIDInput==29) snprintf(SatelliteName_Return,7,"%s"," GOES5");
     if(SatelliteIDInput==31) snprintf(SatelliteName_Return,7,"%s"," GOES6");
     if(SatelliteIDInput==33) snprintf(SatelliteName_Return,7,"%s"," GOES7");
     if(SatelliteIDInput==34) snprintf(SatelliteName_Return,7,"%s","  FY2B");
     if(SatelliteIDInput==35) snprintf(SatelliteName_Return,7,"%s","  FY2C");
     if(SatelliteIDInput==36) snprintf(SatelliteName_Return,7,"%s","  FY2D");
     if(SatelliteIDInput==37) snprintf(SatelliteName_Return,7,"%s","  FY2E");
     if(SatelliteIDInput==38) snprintf(SatelliteName_Return,7,"%s","  FY2F");
     if(SatelliteIDInput==39) snprintf(SatelliteName_Return,7,"%s","  FY2G");
     if(SatelliteIDInput==40) snprintf(SatelliteName_Return,7,"%s","  FY2H");
     if((SatelliteIDInput>=42)&&(SatelliteIDInput<=45)) {
        snprintf(SatelliteName_Return,7,"%s","  NOAA");
     }
     if(SatelliteIDInput==51) snprintf(SatelliteName_Return,7,"%s","  MSG1");
     if(SatelliteIDInput==52) snprintf(SatelliteName_Return,7,"%s","  MSG2");
     if(SatelliteIDInput==53) snprintf(SatelliteName_Return,7,"%s","  MSG3");
     if(SatelliteIDInput==54) snprintf(SatelliteName_Return,7,"%s","  MET3");
     if(SatelliteIDInput==55) snprintf(SatelliteName_Return,7,"%s","  MET4");
     if(SatelliteIDInput==56) snprintf(SatelliteName_Return,7,"%s","  MET5");
     if(SatelliteIDInput==57) snprintf(SatelliteName_Return,7,"%s","  MET6");
     if(SatelliteIDInput==58) snprintf(SatelliteName_Return,7,"%s","  MET7");
     if((SatelliteIDInput>=60)&&(SatelliteIDInput<=69)) {
        snprintf(SatelliteName_Return,7,"%s","  NOAA");
     }
     if(SatelliteIDInput==70) snprintf(SatelliteName_Return,7,"%s"," GOES8");
     if(SatelliteIDInput==72) snprintf(SatelliteName_Return,7,"%s"," GOES9");
     if(SatelliteIDInput==74) snprintf(SatelliteName_Return,7,"%s","GOES10");
     if(SatelliteIDInput==76) snprintf(SatelliteName_Return,7,"%s","GOES11");
     if(SatelliteIDInput==78) snprintf(SatelliteName_Return,7,"%s","GOES12");
     if(SatelliteIDInput==83) snprintf(SatelliteName_Return,7,"%s","  GMS4");
     if(SatelliteIDInput==83) snprintf(SatelliteName_Return,7,"%s","  GMS5");
     if(SatelliteIDInput==84) snprintf(SatelliteName_Return,7,"%s","MTSAT1");
     if(SatelliteIDInput==85) snprintf(SatelliteName_Return,7,"%s","MTSAT2");
     if(SatelliteIDInput==86) snprintf(SatelliteName_Return,7,"%s"," HIM-8");
     if((SatelliteIDInput>=87)&&(SatelliteIDInput<=95)) {
        snprintf(SatelliteName_Return,7,"%s","  DMSP");
     }
     if(SatelliteIDInput==95) snprintf(SatelliteName_Return,7,"%s","  FY1B");
     if(SatelliteIDInput==96) snprintf(SatelliteName_Return,7,"%s","  FY1C");
     if(SatelliteIDInput==97) snprintf(SatelliteName_Return,7,"%s","  FY1D");
     if((SatelliteIDInput>=101)&&(SatelliteIDInput<=171)) {
        snprintf(SatelliteName_Return,7,"%s"," MODIS");
     }
     if(SatelliteIDInput==180) snprintf(SatelliteName_Return,7,"%s","GOES13");
     if(SatelliteIDInput==182) snprintf(SatelliteName_Return,7,"%s","GOES14");
     if(SatelliteIDInput==184) snprintf(SatelliteName_Return,7,"%s","GOES15");
     if(SatelliteIDInput==186) snprintf(SatelliteName_Return,7,"%s","GOES16");
     if(SatelliteIDInput==188) snprintf(SatelliteName_Return,7,"%s","GOES17");
     if((SatelliteIDInput>=195)&&(SatelliteIDInput<=196)) {
        snprintf(SatelliteName_Return,7,"%s","  DMSP");
     }
     if(SatelliteIDInput==230) snprintf(SatelliteName_Return,7,"%s","KLPNA1");
     if(SatelliteIDInput==240) snprintf(SatelliteName_Return,7,"%s","MetOpA");
     if(SatelliteIDInput==241) snprintf(SatelliteName_Return,7,"%s","MetOpB");
     if(SatelliteIDInput==242) snprintf(SatelliteName_Return,7,"%s","MetOpC");
  }

  Ret_ID=0;

  return Ret_ID;

}


int aodt_oceanbasin(float LatitudeInput,
                    float LongitudeInput)
/* 
** determine ocean basin given latitude and longitude position of storm
** Inputs  : LatitudeInput - latitude
**           LongitudeInput - longitude
** Outputs : None
** Return  : basin type (0=atlantic,1=west pacific,2,east pacific,3=Indian)
*/
{
  int   BasinID_Local;                  /* local basin ID value */
  float AtlEPacDivision;                /* atlantic/east pac division value */

  if(LatitudeInput>=0.0) {
     /* northern hemisphere */
     if(LongitudeInput<=-100.0) {
        BasinID_Local=1;      /* West Pacific */
     } 
     else if((LongitudeInput>-100.0)&&(LongitudeInput<=-20.0)) {
        BasinID_Local=3;      /* Indian */
     }
     else if(LongitudeInput>=100.0) {
        BasinID_Local=2;      /* East Pacific */
     }
     else {
        /* -20 to ~+100 */
        if(LatitudeInput>20.0) {
           BasinID_Local=0;      /* Atlantic */
        }
        else if(LatitudeInput<10.0) {
           if(LongitudeInput<80.0) {
              BasinID_Local=0;      /* Atlantic */
           }
           else {
              BasinID_Local=2;      /* East Pacific */
           }
        }
        else {
           /* latitude between 10 and 20 north */
           /* 
           ** slope of line between (100W,20N) and (80W,10N) is 2
           ** if slope of new point and (100,20)  > 2, storm is in atlantic
           ** if slope of new point and (100,20) <= 2, storm is in pacific
           */
           AtlEPacDivision=(100.0-LongitudeInput)/(20-LatitudeInput);
           if(AtlEPacDivision>2.0) {
              BasinID_Local=0;      /* Atlantic */
           }
           else {
              BasinID_Local=2;      /* East Pacific */
           }
        }
     }
  }
  else {
     /* southern hemisphere */
     if(LongitudeInput<=-135.0) {
        BasinID_Local=1;      /* West Pacific */
     }
     else if((LongitudeInput>-135.0)&&(LongitudeInput<=-20.0)) {
        BasinID_Local=3;      /* Indian */
     }
     else if((LongitudeInput>-20.0)&&(LongitudeInput<=67.0)) {
        BasinID_Local=0;      /* Atlantic */
     }
     else {
        BasinID_Local=2;      /* East Pacific */
     }
  }

  return BasinID_Local;

}

int aodt_atcffilename(char *ATCFDirectoryPath_Input,
                      int InputValue,
                      char *StormName_Input, 
                      char *ATCFFileName_Return)
/*
** Derived ATCF file name using various input parameters
** Inputs  : ATCFDirectoryPath_Input - directory path for output file
**           InputValue  - Not currently used
**           StormName_Input - Storm ID character string value (e.g. 12L)
** Outputs : ATCFFileName_Return - Complete file name and path for ATCF file
** Return  : 0
*/
{
  int Ret_ID;                           /* return ID variable */
  int BasinID_Local;                    /* local basin ID value */
  int YearValue;                        /* local year value */
  int DayValue;                         /* local day value */
  int MonthValue;                       /* local month value */
  int HHMMTime;                         /* hours/minutes time value */
  float LatitudeValue_Local;            /* local latitude value */
  float LongitudeValue_Local;           /* local longitude value */
  char *Local_String;                   /* local allocated character string */

  Local_String=calloc((size_t)50,sizeof(char));

  Ret_ID=0;
  LatitudeValue_Local=HistoryCurrentPointer_Global->IR.latitude;
  LongitudeValue_Local=HistoryCurrentPointer_Global->IR.longitude;
  BasinID_Local=aodt_oceanbasin(LatitudeValue_Local,LongitudeValue_Local);
  /* rawflag=InputValue/1000; */
  /* finalflag=(InputValue-(rawflag*1000))/100; */
  (void)aodt_yddmy(HistoryCurrentPointer_Global->IR.date,
                   &DayValue,&MonthValue,&YearValue);
  HHMMTime=HistoryCurrentPointer_Global->IR.time/100;
  snprintf(Local_String,50,"%10s_%4d%2.2d%2.2d%4.4d_%3s_%3s",
         "NESDIS_DVTO",YearValue,MonthValue,DayValue,HHMMTime,
         StormName_Input,"FIX");
  strcpy(ATCFFileName_Return,ATCFDirectoryPath_Input);
  strncat(ATCFFileName_Return,Local_String,strlen(Local_String));

  free(Local_String);

  return Ret_ID;;

}

int aodt_initcurrent( int InitializeResetIDFlag )
/* 
** initialize HistoryCurrentPointer_Global array or 
** reset values for land interaction situations 
** Inputs  : InitializeResetIDFlag - Reset value in Current History structure
** Outputs : None
** Return  : Success/Failure Flag
*/
{
  int Ret_ID;                           /* return ID variable */
  int SizeDouble;                       /* double size */
  int SizeFloat;                        /* float size */
  char CommentString_Local[50]="\0";    /* local comment string value */

  Ret_ID=0;
  if(InitializeResetIDFlag==0) {
     HistoryCurrentPointer_Global=(struct HistoryFileRec_Structure *)
                   malloc(sizeof(struct HistoryFileRec_Structure));
     HistoryCurrentPointer_Global->IR.latitude=MISSFP2;
     HistoryCurrentPointer_Global->IR.longitude=MISSFP2;
     HistoryCurrentPointer_Global->IR.land=0;
     HistoryCurrentPointer_Global->IR.autopos=0;
     HistoryCurrentPointer_Global->IR.sattype=0;
     strcpy(HistoryCurrentPointer_Global->IR.comment,CommentString_Local);
     DiagnosticString_Global=calloc((size_t)50000,sizeof(char));
     HistoryFileName_Global=calloc((size_t)200,sizeof(char));
     ForecastFileName_Global=calloc((size_t)200,sizeof(char));
     HistoryFilePath_Global=calloc((size_t)200,sizeof(char));
     TopographyFilePath_Global=calloc((size_t)200,sizeof(char));
     ListFilePath_Global=calloc((size_t)200,sizeof(char));
     ForecastFilePath_Global=calloc((size_t)200,sizeof(char));
     ImageDataFilePath_Global=calloc((size_t)200,sizeof(char));
     PMWHDFFilePath_Global=calloc((size_t)200,sizeof(char));
     
     SizeFloat=sizeof(float);
     SizeDouble=sizeof(double);
     ForecastLatitudes_Global=calloc((size_t)10,SizeFloat);
     ForecastLongitudes_Global=calloc((size_t)10,SizeFloat);
     ForecastTimes_Global=calloc((size_t)10,SizeDouble);
  }

  HistoryCurrentPointer_Global->IR.Traw=0.0;
  HistoryCurrentPointer_Global->IR.TrawO=0.0;
  HistoryCurrentPointer_Global->IR.Tfinal=0.0;
  HistoryCurrentPointer_Global->IR.CI=0.0;
  HistoryCurrentPointer_Global->IR.eyet=MISSFP1;
  HistoryCurrentPointer_Global->IR.warmt=MISSFP1;
  HistoryCurrentPointer_Global->IR.cloudt=MISSFP1;
  HistoryCurrentPointer_Global->IR.cloudt2=MISSFP1;
  HistoryCurrentPointer_Global->IR.cwcloudt=MISSFP1;
  HistoryCurrentPointer_Global->IR.warmlatitude=MISSFP2;
  HistoryCurrentPointer_Global->IR.warmlongitude=MISSFP2;
  HistoryCurrentPointer_Global->IR.eyecdosize=0.0;
  HistoryCurrentPointer_Global->IR.eyestdv=0.0;
  HistoryCurrentPointer_Global->IR.cloudsymave=0.0;
  HistoryCurrentPointer_Global->IR.eyescene=0;
  HistoryCurrentPointer_Global->IR.cloudscene=0;
  HistoryCurrentPointer_Global->IR.eyesceneold=-1;
  HistoryCurrentPointer_Global->IR.cloudsceneold=-1;
  HistoryCurrentPointer_Global->IR.rule9=0;
  HistoryCurrentPointer_Global->IR.rule8=0;
  HistoryCurrentPointer_Global->IR.LBflag=0;
  HistoryCurrentPointer_Global->IR.rapiddiss=0;
  HistoryCurrentPointer_Global->IR.eyefft=0;
  HistoryCurrentPointer_Global->IR.cloudfft=0;
  HistoryCurrentPointer_Global->IR.cwring=0;
  HistoryCurrentPointer_Global->IR.ringcb=0;
  HistoryCurrentPointer_Global->IR.ringcbval=0;
  HistoryCurrentPointer_Global->IR.ringcbvalmax=0;
  HistoryCurrentPointer_Global->IR.CIadjp=0.0;
  HistoryCurrentPointer_Global->IR.rmw=MISSFN1;
  HistoryCurrentPointer_Global->IR.mwscore=MISSFN1;
  HistoryCurrentPointer_Global->IR.mwdate=1900001;
  HistoryCurrentPointer_Global->IR.mwtime=0;
  HistoryCurrentPointer_Global->IR.r34=MISSIN1;
  HistoryCurrentPointer_Global->IR.MSLPenv=MISSIN2;
  if(InitializeResetIDFlag==0) {
     HistoryCurrentPointer_Global->nextrec=NULL;
  }
  ErrorBit1=-1;
  ErrorBit0=-1;
  DiagBit1=-1;
  DiagBit0=-1;

  return Ret_ID;

}

aodt_stormspeedfromhistory(float *StormMotionSpeed, 
                           float *StormMotionDirection)
/* 
** calculate storm speed and direction from history file.  Programs
** from aodtwinds.c routines.
*/
{
  int RetErr;                           /* fuction error return variable */
  float SpeedReturn;                    /* Local speed value */
  float DirectionReturn;                /* Local direction value */
  float TimeDifference;                 /* time difference */
  float Latitide_Current;               /* current storm latitude */
  float Longitude_Current;              /* current storm longitude */
  float Latitude_Previous;              /* previous storm latitude */
  float Longitude_Previous;             /* previous storm longitude */
  float LatitudeMinus6hrs;              /* latitude position @ T-6hrs */
  float LongitudeMinus6hrs;             /* longitude position @ T-6hrs */
  float LatitudeMinus6hrs_interp;       /* interp latitude position @ T-6hrs */
  float LongitudeMinus6hrs_interp;      /* interp longitude position @ T-6hrs */
  double CurrentTime;                   /* current time */
  double CurrentTimeMinus6hrs;          /* current time @ T-6hrs */
  double HistoryRecTime;                /* history record time */
  double HistoryRecTimePrevious;        /* previous history record time */
  double HistoryRecTimeMinus6hrs;       /* previous history time @ T-6hrs */
  struct HistoryFileRec_Structure *HistoryRecordPtr; /* history structure */
  logical FoundMinus6hrRecordTF=FALSE;  /* found record @ T-6hrs logical */

  /* get curreint info */
  Latitide_Current=HistoryCurrentPointer_Global->IR.latitude;
  Longitude_Current=HistoryCurrentPointer_Global->IR.longitude;

  CurrentTime=aodt_calctime(HistoryCurrentPointer_Global->IR.date,
                            HistoryCurrentPointer_Global->IR.time);
  CurrentTimeMinus6hrs=aodt_calctime(HistoryCurrentPointer_Global->IR.date,
                            HistoryCurrentPointer_Global->IR.time)-0.25;

  /* get previous storm position for storm motion calculation */
  Latitude_Previous=Latitide_Current;
  Longitude_Previous=Longitude_Current;
  if((HistoryFirstRecordPointer_Global==0)||
     (strlen(HistoryFileName_Global)==0)) {
     /* dont calculate storm motion... no previous record */
     SpeedReturn=0.0;
     DirectionReturn=0.0;
     RetErr=0;
  }
  else {
     HistoryRecordPtr=HistoryFirstRecordPointer_Global;
     HistoryRecTimePrevious=CurrentTime;
     while(HistoryRecordPtr!=0) {
        HistoryRecTime=aodt_calctime(HistoryRecordPtr->IR.date,
                                     HistoryRecordPtr->IR.time);
        if(HistoryRecTime>=CurrentTime) {
           break;
        }
        Latitude_Previous=HistoryRecordPtr->IR.latitude;
        Longitude_Previous=HistoryRecordPtr->IR.longitude;
        HistoryRecTimePrevious=HistoryRecTime;
        /* determine six hour previous point */
        if((HistoryRecTime>=CurrentTimeMinus6hrs)&&
           (!FoundMinus6hrRecordTF)) {
          LatitudeMinus6hrs=HistoryRecordPtr->IR.latitude;
          LongitudeMinus6hrs=HistoryRecordPtr->IR.longitude;
          HistoryRecTimeMinus6hrs=HistoryRecTime;
          LatitudeMinus6hrs_interp=aodt_slopecal((double)6.0,3);
          LongitudeMinus6hrs_interp=aodt_slopecal((double)6.0,4);
          if((Latitide_Current<999.0)&&
             (Longitude_Current<999.0)) {
             FoundMinus6hrRecordTF=TRUE;
          }
        }
        HistoryRecordPtr=HistoryRecordPtr->nextrec;
     }
     SpeedReturn=0.0;
     DirectionReturn=0.0;
     if(FoundMinus6hrRecordTF) {
        /* printf("USING 6 HOUR AVERAGE MOTION \n"); */
        TimeDifference=(float)(CurrentTime-HistoryRecTimeMinus6hrs);
        if(LatitudeMinus6hrs_interp<999.) {
           RetErr=aodt_motion(LatitudeMinus6hrs,LongitudeMinus6hrs,
                              LatitudeMinus6hrs_interp,LongitudeMinus6hrs_interp,
                              TimeDifference,
                              &SpeedReturn,&DirectionReturn);
        }
     }
     else {
        /* printf("USING 1 HOUR MOTION \n"); */
        TimeDifference=(float)(CurrentTime-HistoryRecTimePrevious);
        if(TimeDifference>0.0) {
           RetErr=aodt_motion(Latitude_Previous,Longitude_Previous,
                              Latitide_Current,Longitude_Current,
                              TimeDifference,
                              &SpeedReturn,&DirectionReturn);
        }
     }
     *StormMotionSpeed=SpeedReturn;
     *StormMotionDirection=DirectionReturn;
  }

  return RetErr;
}

int aodt_checkIRdataHIE(void)
{
/* Perform Line/Element temperatue data check for input image array 
** to make sure values are w/in gross error check temperature thresholds.
** This is only done for HIE version of the ADT code and is called from
** the routine odtloadIRimageHIE.c.
** Inputs  : None
** Outputs : None
** Return  : -17 : data read off edge of image
**            12 : good read of image
*/

  int RetErr;                           /* fuction error return variable */
  int XInc;                             /* X direction increment */
  int YInc;                             /* Y direction increment */
  int BadPointCounter;                  /* counter for bad pixel values */
  int BadLineCounter;                   /* counter for bad lines values */
  float TemperatureValue;               /* pixel temperatue value */
  float MinimumThresholdValue=120.0;    /* minimum temperatue threshold value */
  float MaximumThresholdValue=320.0;    /* maximum temperatue threshold value */
  char *OutStr;                         /* local output string */

  OutStr=calloc((size_t)5000,sizeof(char));

  RetErr=12;
  BadLineCounter=0;
  for(YInc=0;YInc<IRDataGrid_Global->numy;YInc++) {
     BadPointCounter=0;
     for(XInc=0;XInc<IRDataGrid_Global->numx;XInc++) {
        TemperatureValue = IRDataGrid_Global->temp[YInc][XInc];
        if((TemperatureValue>=MaximumThresholdValue)||
           (TemperatureValue<=MinimumThresholdValue)) {
           /* data at point is bad (>47C or <-153C).  Count number
              of bad points in line... if there are more than 10 bad 
              points, replace entire line with previous line */
           BadPointCounter++;
           if(XInc>0) {
              IRDataGrid_Global->temp[YInc][XInc]=
                          IRDataGrid_Global->temp[YInc][XInc-1];
           }
        }
     }
     if((BadPointCounter>10)&&(YInc>0)) {
        /* data line in image is bad... replace with previous */
        snprintf(OutStr,5000,"BAD DATA POINTS DETECTED AT LINE %d, REPLACING WITH PREVIOUS LINE\n",YInc);
        strcat(DiagnosticString_Global,OutStr);
        for(YInc=0;YInc<IRDataGrid_Global->numy;YInc++) {
           IRDataGrid_Global->temp[YInc][XInc]=
                      IRDataGrid_Global->temp[YInc-1][XInc];
        }
        BadLineCounter++;
        if(BadLineCounter>=10) {
           snprintf(OutStr,5000,"BAD IMAGE... EXITING\n");
           strcat(DiagnosticString_Global,OutStr);
           YInc=IRDataGrid_Global->numy;
           RetErr=-17;
        }
     }
  }

  free(OutStr);
  return RetErr;

}

float aodt_XLZA(float rlat,float rlon,float plat,float plon)
/* approximate local zenith angle (angle from vertical to the satellite)
   Inputs  : rlat - latitude of target position on earth
             rlon - longitude of target position on earth
             plat - latitude of sub-satellite point
             plon - longitude of sub-satellite point
   Outputs : ang  - local zenith angle
*/
{
  double dif, r, sd, cosal, yy, zz, dp, alpha, zaz;

  r = 6371.229;           /*  mean earth radius   */
  sd = r + 35790.0;       /*  mean distance from satellite to earth center  */

  rlat = fabs(rlat);
   /*  For the latitude in East of Greenwich  */
  if (rlon < 0) rlon = rlon + 360.0;
  rlon = fabs(rlon);
  dif = plon - rlon;
  dif = fabs(dif);

  cosal = cos(deg_to_rad(rlat)) * cos(deg_to_rad(dif));
  dp = sqrt(pow(r, 2) + pow(sd, 2) - (2 * r * sd * cosal));
  yy = sqrt(1-pow(cosal, 2)) / cosal;
  alpha = atan(yy);
  zz = sd * sin(alpha) / dp;
  yy = zz / (sqrt(1 - pow(zz, 2)));
  zaz = atan(yy);
  zaz = rad_to_deg(zaz);
  return zaz;

}