///////////////////////////////////////////////////////////////////////////
// NAME: utility.cpp
// FUNCTION: Utility routines
// DESCRIPTION: Some useful utility routines
// REFERENCE: none
// CALLING SEQUENCE: none
// INPUTS: none
// OUTPUTS: none
// DEPENDENCIES: utility.h
// RESTRICTIONS: none
// HISTORY: none
//////////////////////////////////////////////////////////////////////////
#include <string>
#include "utility.h"
#include <cmath>
#include "ABI_const.h"
///////////////////////////////////////////////////////////////////////////
// NAME: void locate
// FUNCTION: Locates -x- in array -xt-.
// DESCRIPTION: Locates -x- in array -xt-, so that -x- is inside of an 
//              interval of two adjacent -xt- values: xt(i1) and xt(i2).
//              Returns i1, i2 and fraction frac = ( xt[i2] - x ) /
//              ( xt[i2] - xt[i1] ).
// REFERENCE: none
// CALLING SEQUENCE: locate(x,xt,nxt,i1,i2,frac)
// INPUTS: x, xt(array), nxt
// OUTPUTS: i1, i2, frac
// DEPENDENCIES: none
// RESTRICTIONS: none
// HISTORY: none
//////////////////////////////////////////////////////////////////////////
void locate(float32 x, float32 xt[], int32 nxt, int32& i1, int32& i2,
   float32& frac)
{
   if (x < xt[0] || nxt == 1) {
      i1 = 0;
      i2 = 0;
      frac = 1.0;
      return;
   }
   
   if (x > xt[nxt-1]) {
      i1 = nxt - 1;
      i2 = nxt - 1;
      frac = 1.0;
      return;
   }   
   
   for (int32 i=1; i<nxt; i++)
      if (x <= xt[i]) {
         i1 = i-1;
	 i2 = i;
	 if (abs(xt[i2] - xt[i1]) > 0.)
            frac = (xt[i2] - x)/(xt[i2] - xt[i1]);
         return;
      }      
}


//check the file as if for windows:  replace the '/'  with '\' and seach again
void replaceAll(char *pString, char originChar, char newchar)
{
   int len = strlen(pString);
   for (int i=0; i<len; i++) {
      if(pString[i] == originChar) pString[i] = newchar;

   }

}


//Set the cooresponding bit in Quality falg to 1
uint8 Set_QualityFlag_bit(uint8 Flag, uint8 Bit_No)
{
    Flag =Flag | uint8(pow(2.0,Bit_No-1));
 return Flag;
}


/******************************************************************************/
/*                                                                            */
/* NAME:                                                                      */
/*    RayRefl                                                                 */
/*                                                                            */
/* FUNCTION:                                                                  */
/*    Calculate Rayleigh reflectance                                          */
/*                                                                            */
/* DESCRIPTION:                                                               */
/*    This function calculate the molecular reflectance for a given           */
/*    molecular optical thickness, geometry and ABI channel.                  */
/*                                                                            */
/* REFERENCE                                                                  */
/*    Subroutine CHAND of 6S radiative transfer code                          */
/*    Vermote E.F. and Tanré, D., (1992)                                      */
/*      Analytical Expressions for Radiative Properties of Planar Rayleigh    */
/*      Scattering Media Including Polarization Contribution.                 */
/*      Journal of Quantitative Spectroscopy and Radiative Transfer ,         */
/*      47, (4): 305-314.                                                     */
/*                                                                            */
/* CALLING SEQUENCE                                                           */
/*    Called from PreProcess (ABIAOT_LndAlgo.cpp and ABIAOT_OcnAlgo.cpp)      */
/*                                                                            */
/*                                                                            */
/* INPUTS                                                                     */
/*    MolOD     Molecular optical thickness                                     */
/*    CosSza  Cosine of solar zenith angle                                    */
/*    CosVza  Cosine of viewing zenith angle                                  */
/*    RelAzi  Relative azimuth angle                                          */
/*    Chn     ABI channel index                                               */
/*                                                                            */
/* OUTPUTS                                                                    */
/*    Return Rayleigh reflectance                                             */
/*                                                                            */
/* DEPENDENCIES                                                               */
/*    None                                                                    */
/*                                                                            */
/* RESTRICTIONS                                                               */
/*    None                                                                    */
/*                                                                            */
/* HISTORY                                                                    */
/*                                                                            */
/*   DATE      AUTHOR              Description                                */
/*   --------  -----------    -----------------                               */
/*   02/10/07  NOAA/NESDIS    First Version                                   */
/*   03/20/08  NOAA/NESDIS    Second Version                                  */
/*                                                                            */
/******************************************************************************/

float32 RayRefl(float32 MolOD, float32 CosSza, float32 CosVza,
                float32 RelAzi, int32 Chn)
{
   /* Cosine of relative azimuth angle used for Fourier expansion           */

   float32 Phi = PIc - RelAzi * DEG2RAD;    /* PI - (SolAzi - SatAzi)        */
   float32 CosPhi0 = 1.;                   /* Cosine of N*Phi  (N=0,1,2)    */
   float32 CosPhi1 = cos(Phi);
   float32 CosPhi2 = cos(2.*Phi);


   /* Fourier expansion of Rayleigh Scattering phase function (N=0,1,2)     */

   float32 DepFac = 0.0279;                    /* Depolarization factor   */
   float32 DepCorr = DepFac / (2.0 - DepFac);  /* Depolarization correction */
   DepCorr = (1.0 - DepCorr)/(1.0 + 2.* DepCorr);
   float32 PhaseFun0 = 1.0 + (3.0 * CosSza * CosSza - 1.0) *
                             (3.0 * CosVza * CosVza - 1.0) * DepCorr / 8.0;
   float32 PhaseFun1 = -1.0 * CosSza * CosVza * sqrt(1.0 - CosSza * CosSza) *
                       sqrt(1.0 - CosVza * CosVza) * DepCorr * 0.5 * 1.5;
   float32 PhaseFun2 = (1.0 - CosSza * CosSza) * (1.0 - CosVza * CosVza) *
                       DepCorr * 0.5 * 0.375;


   /* Fourier expansion of single scattering contribution (N=0,1,2)         */

//   float32 Tmp = (1.0 - exp(-MolOD * (1.0 / CosSza + 1.0 / CosVza))) /
//                 (4.0 * (CosSza + CosVza));
   float32 Tmp = 0.0;
   if (abs(CosSza) > 0.0 && abs(CosVza) > 0.0)
      Tmp = (1.0 - exp(-MolOD * (1.0 / CosSza + 1.0 / CosVza))) /
            (4.0 * (CosSza + CosVza));
   float32 SingleScat0 = PhaseFun0 * Tmp;
   float32 SingleScat1 = PhaseFun1 * Tmp;
   float32 SingleScat2 = PhaseFun2 * Tmp;

   /* Calculate DeltaMolOD   (N=0,1,2)                                        */

   float32 DeltaMolODCoef[10];   /* Coefficients used for estimating DeltaMolOD */
   DeltaMolODCoef[0] = 1.0;
   DeltaMolODCoef[1] = log(MolOD);
   DeltaMolODCoef[2] = CosSza + CosVza;
   DeltaMolODCoef[3] = DeltaMolODCoef[1] * DeltaMolODCoef[2];
   DeltaMolODCoef[4] = CosSza * CosVza;
   DeltaMolODCoef[5] = DeltaMolODCoef[1] * DeltaMolODCoef[4];
   DeltaMolODCoef[6] = CosSza * CosSza + CosVza * CosVza;
   DeltaMolODCoef[7] = DeltaMolODCoef[1] * DeltaMolODCoef[6];
   DeltaMolODCoef[8] = CosSza * CosSza * CosVza * CosVza;
   DeltaMolODCoef[9] = DeltaMolODCoef[1] * DeltaMolODCoef[8];

   float32 DeltaMolOD0 = 0.0;


   for (int32 i = 0; i < 10; i++)
   {
      DeltaMolOD0 += DeltaMolODCoef[i] * RAYCOEFF0[i];
   }
   float32 DeltaMolOD1 = DeltaMolODCoef[0] * RAYCOEFF1[0] +
                       DeltaMolODCoef[1] * RAYCOEFF1[1];
   float32 DeltaMolOD2 = DeltaMolODCoef[0] * RAYCOEFF2[0] +
                       DeltaMolODCoef[1] * RAYCOEFF2[1];


   /* Calculate Rayleigh Reflection
       Summation of first 3 (N=0,1,2) Fourier expansion terms               */

   if (abs(CosSza) > 0.0 && abs(CosVza) > 0.0)
       Tmp = (1.0 - exp(-MolOD / CosSza)) * (1.0 - exp(-MolOD/ CosVza));

   float32 RayRef0 = (SingleScat0 + PhaseFun0 * DeltaMolOD0 * Tmp) * CosPhi0;
   float32 RayRef1 = (SingleScat1 + PhaseFun1 * DeltaMolOD1 * Tmp) * CosPhi1
                      * 2.0;
   float32 RayRef2 = (SingleScat2 + PhaseFun2 * DeltaMolOD2 * Tmp) * CosPhi2
                      * 2.0;

   return (RayRef0 + RayRef1 + RayRef2);
}


//Calculating sunglint angle 
float32 SunGlintAng(float32 SOLZA, float32 SENZA, float32 SOLAZ, float32 SENAZ)
{
    float32 SunGlintA = Missing_Value;
    float32 DegreeToRadiance = 3.1415926/180.0;
    float32 RadianceToDegree = 180.0/3.1415926;
    float32 Temp1 = cos(SOLZA*DegreeToRadiance)*cos(SENZA*DegreeToRadiance);
    float32 Temp2 = sin(SOLZA*DegreeToRadiance)*sin(SENZA*DegreeToRadiance);
    float32 Temp3 = cos((180.0-(SOLAZ-SENAZ))*DegreeToRadiance);
    SunGlintA = acos(Temp1+Temp2*Temp3)*RadianceToDegree;

 return SunGlintA;
}
float32 ConfidenceLevel_U(float32 test, float32 Threshold)
{
    float32 offset1=abs(Threshold)*0.01;
    float32 offset2=abs(Threshold)*0.02;
    if (Threshold ==0.0 ){
    offset1=0.01;
    offset2=0.02;
    }
    float32 Con_Lev = 0.0;
    if ((test >=(Threshold+offset1)) && (test <(Threshold+offset2))) {
       Con_Lev=0.5;
     }
    if (test >=(Threshold+offset2)) {
       Con_Lev=1.0;
     }
 return Con_Lev;
}

float32 ConfidenceLevel_L(float32 test, float32 Threshold)
{
    float32 offset1=abs(Threshold)*0.01;
    float32 offset2=abs(Threshold)*0.02;
    if (Threshold ==0.0 ){
    offset1=0.01;
    offset2=0.02;
    }

    float32 Con_Lev = 0.0;
    if ((test <=(Threshold-offset1)) && (test >(Threshold-offset2))) {
       Con_Lev=0.5;
     }
    if (test <=(Threshold-offset2)) {
       Con_Lev=1.0;
     }
 return Con_Lev;
}

float32 ConfidenceLevel_M(float32 test, float32 Threshold1, float32 Threshold2)
{
    float32 Con_Lev = 0.5;
    if ((test >=(Threshold1+abs((Threshold2-Threshold1))/3.0)) && (test <=(Threshold2-abs((Threshold2-Threshold1))/3.0))) {
       Con_Lev=1.0;
     }
    if ((test < (Threshold1+abs((Threshold2-Threshold1))/3.0)) || (test >(Threshold2-abs((Threshold2-Threshold1))/3.0))) {
       Con_Lev=0.0;
     }
 return Con_Lev;
}

/* Function: find files by pattern, and copy the filename of the first file
   to filename1 & return 1 if file is found.
   Limitation: there is file saperator "/" or "\" in the string "pattern"
*/
int findFile(const char *pattern, char *filename1)
{
   char syscommand[STRMAX];           /* Glob pattern. */
   char tempFile[STRMAX]  ;

   sprintf(tempFile,"temp.txt");
   /* get file names based on year and day, and put them in a file */

   char *pos= strstr(pattern, "/");
   sprintf(syscommand, "ls     %s > %s", pattern, tempFile);

   if(pos == NULL)
   sprintf(syscommand, "dir /b %s > %s", pattern, tempFile);

   //printf("   Search command= '%s'\n", syscommand);
   system(syscommand);

   /* reaf file names from this file */
   FILE *pFile = fopen (tempFile, "rt");
   if (pFile==NULL) return 0;

   fscanf(pFile, "%s", filename1);
   fclose (pFile);

   remove(tempFile);
   return 1;
}