#include <stdio.h>
#include <stdlib.h>

namespace AtmosphereNS {

  const float c_cmperns = 29.9792458; 
  const float hGroundRef= 1452.e2;

  class AtmModel {
  public:
    float ai[5],bi[5],ci[5];
    float hi[5]; // cm
    int atmflag; // 1...12 Montly models  10,...,50 Seasonal models
  };

  class Atmosphere {    
  private:

    AtmModel amonth[12];
    AtmModel aseason[4];
    AtmModel acurrent;
    bool ok; //Initialized to false, only set to true when acurrent is set

  public:
    Atmosphere(bool UseGDAS=false);
    ~Atmosphere();
    AtmModel *GetCurrentAtmosphere(){return &acurrent;};
    int GetCurrentAtmosphereMonth(){ return acurrent.atmflag; };
    void SetCurrentAtmosphere(AtmModel theAtm);
    void SetCurrentAtmosphere(int imonth);  // imonth=1,...,12
    void SetCurrentAtmosphere(float *ai, float *bi, float *ci, float *hi); //
    void printCurrentAtmosphere(FILE *fp);
    void ClearCurrentAtmosphere();

    //Using current atmosphere
    float VertDepthToHeight(float vdepth);
    float SlantDepthToHeight(float sdepth, float theta );
    float HeightToVertDepth(float h);
    float HeightToSlantDepth(float h, float theta);
    float GetDensity_i(float h);
    float GetDensity(float h);
    float GetDensityAboveGround(float hground, float theta, float nRad);
    float GetDensityAboveGround_i(float hground, float theta, float nRad);
    float GetSlant0(){return acurrent.ai[4]*acurrent.ci[4]/acurrent.bi[4];};
    float GetC0(){return acurrent.ci[0];};

    float GetDX_rectilinear(float r, float psi, float SlantDepth,     float theta, float hground );
    float GetDX_rectilinear_i(float r, float psi, float SlantDepth,     float theta, float hground );
    float GetDX_diffusive(  float r, float psi, float SlantDepth,     float theta, float hground );
    float GetDX_diffusive_i(  float r, float psi, float SlantDepth,     float theta, float hground );

    float GetTimePF_delay(float x, float y, float z,  float theta, float azi);
    float GetTimePF(float r, float psi, float SlantDepth, float theta, float hground);
    float GetTimePF_h(float r, float psi, float hStart, float theta, float hground);

    float GetTimeCF_h(float r, float psi, float h1, float theta, float hground);
    float GetTimeCF(float r, float psi, float SlantDepth, float theta, float hground);
    
  };
}