Atmosphere.h

Go to the documentation of this file.

// $Id: Atmosphere.h,v 1.2 2014/06/04 17:09:38 mpound Exp $
#include "carma/services/Units.h"
#include "carma/services/Angle.h"
#include "carma/services/Pressure.h"
#include "carma/services/Temperature.h"
#include "carma/services/Velocity.h"

#ifndef CARMA_SERVICES_ATMOSPHERE_H
#define CARMA_SERVICES_ATMOSPHERE_H
namespace carma  {
  namespace services {
  class Atmosphere {
      public:

      //------------------------------------------------------
      //             PUBLIC STATIC CONSTANTS 
      //------------------------------------------------------
      

      static const double DEFAULT_ATM_PRESSURE = 780.0;

      static const double STANDARD_ATM_PRESSURE = 1013.25;

      static const double MIN_ATM_PRESSURE = 680; // 0.9*DEFAULT for CARMA

      static const double MAX_ATM_PRESSURE = 1115.0; //1.1atm at sealevel

      //"room temperature"
      static const double DEFAULT_AIR_TEMP = 300.0; //80 Fahrenheit

      //??need reasonable limits
      static const double MIN_AIR_TEMP = 242.0; // -25 Fahrenheit

      //??need reasonable limits
      static const double MAX_AIR_TEMP = 316.5; // 110 Fahrenheit 

      //??need reasonable limits
      static const double DEFAULT_RH = 25.0; 

      //??need reasonable limits
      static const double MIN_RH = 5.0; 

      static const double MAX_RH = 100.0; 

      // This is the dewpoint for the DEFAULT_AIRTEMP (300K)
      // and DEFAULT_RH (50%)
      //  carma/environment/Test/tAtmosphere RH=50 Ta=26.85
      // 25% = 277.388
      //  carma/environment/Test/tAtmosphere RH=25 Ta=26.85
      // 50% = 288.761; 
      static const double DEFAULT_DEW_POINT = 278.388;

      //For MIN_AIRTEMP, MIN_RH
      static const double MIN_DEW_POINT = 221.9;

      //dewpoint at 100% RH is the airtemp
          //value is initialized with the storage declaration in the .cc file
      static const double MAX_DEW_POINT;

      static const double MIN_ELEVATION = 0.0174533; // radians

      static const double DEFAULT_WIND_SPEED = 0.0;

      static const double MIN_WIND_SPEED = 0.0;

      static const double MAX_WIND_SPEED = 200.0;

      typedef enum dewPointMethodEnum { 
          CLAUSIUS_CLAPEYRON, 
          MAGNUS_TETENS,
          HOFFMAN_WELCH 
      } dewPointMethodType ;

      //------------------------------------------------------
      //            PUBLIC METHODS
      //------------------------------------------------------
      
      Atmosphere();

      virtual ~Atmosphere();

    double computeSaturatedPressure(double airTemp) const;

    double waterColumn(double airTemp, double relHumid) const;

    double waterPartialPressure(double airTemp, double relHumid) const;

    double waterVaporDensity(double airTemp, double relHumid) const;


    double computeDewPoint(double airTemp, 
        double relHumid, 
        dewPointMethodType method = CLAUSIUS_CLAPEYRON ) const;


    double computeHumidity(double airTemp, 
        double dewTemp, 
        dewPointMethodType method = CLAUSIUS_CLAPEYRON ) const;

    double computeZenithRefractivity(double airTemp, 
                     double atmPressure,
                     double relHumid,
                     double frequency) const;

    double computeRefractionCorrection(double airTemp, 
                       double atmPressure,
                       double relHumid,
                       double elevation,
                       double frequency,
                       double altitude) const;
    double ulichRefractionCorrection(double airTemp, 
                     double atmPressure, 
                     double relHumid,
                     double elevation,
                     double frequency) const;
    double computePathlength(double airTemp,
                 double atmPressure,
                 double relHumid,
                 double elevation,
                 double frequency,
                 double altitude = 0.0) const;


    double computeHorizonPathlength(double airTemp, 
                    double atmPressure,
                    double relHumid,
                    double frequency) const;

    static double safeAirTemperature(double airTemp);

    static bool isSafeAirTemperature(double airTemp);

    static carma::services::Temperature safeAirTemperature(
        carma::services::Temperature airTemp);
    static bool isSafeAirTemperature( carma::services::Temperature airTemp);
    
    carma::services::Temperature lastSafeAirTemperature( double airTemp );

    carma::services::Temperature lastSafeAirTemperature( 
       carma::services::Temperature airTemp );

    carma::services::Temperature lastSafeAirTemperature( ) const;

    static double safeAtmPressure(double atmPressure);

    static bool isSafeAtmPressure(double atmPressure);

    static carma::services::Pressure safeAtmPressure(
        carma::services::Pressure atmPressure );
    static bool isSafeAtmPressure( carma::services::Pressure atmPressure );
    
    carma::services::Pressure lastSafeAtmPressure( double atmPressure );

    carma::services::Pressure lastSafeAtmPressure( 
       carma::services::Pressure atmPressure );

    carma::services::Pressure lastSafeAtmPressure( ) const;

    static double safeRelativeHumidity(double relHumid);

    static bool  isSafeRelativeHumidity(double relHumid);

    double lastSafeRelativeHumidity( double relHumid );

    double lastSafeRelativeHumidity( ) const;
    
    static double safeDewPoint(double dewPoint);

    static bool  isSafeDewPoint(double dewPoint);

    static carma::services::Temperature safeDewPoint(
        carma::services::Temperature dewPoint );

    static bool  isSafeDewPoint( carma::services::Temperature dewPoint );

    bool isConsistentDewPoint(services::Temperature dewPoint,
                  services::Temperature airTemp, 
                  double relHumid,
        dewPointMethodType method = CLAUSIUS_CLAPEYRON,
            double consistencyPercent = 10.0 );

    carma::services::Temperature lastSafeDewPoint( double dewPoint );

    carma::services::Temperature lastSafeDewPoint( 
       carma::services::Temperature dewPoint );

    carma::services::Temperature lastSafeDewPoint( ) const;
    
    static double safeElevation(double elevation);

    static bool isSafeElevation(double elevation);

    static carma::services::Angle safeElevation(
        carma::services::Angle elevation );
    static bool isSafeElevation( carma::services::Angle elevation );

    static carma::services::Velocity safeWindSpeed(
        carma::services::Velocity windSpeed );

    /*
     * @return true is wind speed 
     * is within reasonable bounds for the CARMA site. 
     * @see #safeWindSpeed(services::Velocity windSpeed)
     */
    static bool isSafeWindSpeed( carma::services::Velocity windSpeed );
    
    static double safeWindSpeed(double windSpeed);

    /*
     * @return true is wind speed 
     * is within reasonable bounds for the CARMA site. 
     * @see #safeWindSpeed(double windSpeed)
     */
    static bool isSafeWindSpeed( double windSpeed );

      private:
      //------------------------------------------------------
      //            PRIVATE STATIC CONSTANTS 
      //------------------------------------------------------
      static const double CC_A = 6.105; // leading coeff.
      static const double CC_B = 25.22; // T coeff.
      static const double CC_C = 5.31;  // ln(T) coeff.

      static const double AWIPS_C1  = 0.0091379024;
      static const double AWIPS_C2  = 6106.396;
      static const double AWIPS_C3  = 223.1986;
      static const double AWIPS_C4  = 0.0182758048;
      static const double AWIPS_C15 = 26.66082;


      static const double MT_A = 17.27;
      static const double MT_B = 237.7;

      static const double HW_A = 1.598E9;
      static const double HW_B = 5370;

      static const double SW_DRY_AIR = 77.6;

      static const double BW_DRY_AIR = 77.6;

      static const double SW_INDUCED_DIPOLE = 0.128;
      static const double BW_INDUCED_DIPOLE = 5.6;

      static const double SW_PERM_DIPOLE = 3.776E3;
      static const double BW_PERM_DIPOLE = 3.75E3;

      static const double REFRAC_COEFF_A = -2.8354E-5;

      static const double REFRAC_COEFF_B = 5.4012E-7;

      static const double MIN_REFRAC_FREQ = 50.0; //50 GHz

      // Equations 13. (radio)
       static const double A1_RAD_1     =  0.5753868; // constant
       static const double A1_RAD_P     =  0.5291E-4; // * P
       static const double A1_RAD_PPW   = -0.2819E-4; // * Pwater
       static const double A1_RAD_PPWSQ = -0.9381E-6; // * Pwater^2
       static const double A1_RAD_T     = -0.5958E-3; // * T
       static const double A1_RAD_TSQ   =  0.2657E-5; // * T^2

       static const double A2_RAD_1     =  1.301211;  // constant
       static const double A2_RAD_P     =  0.2003E-4; // * P
       static const double A2_RAD_PPW   = -0.7285E-4; // * Pwater
       static const double A2_RAD_PPWSQ =  0.2579E-5; // * Pwater^2
       static const double A2_RAD_T     = -0.2595E-2; // * T
       static const double A2_RAD_TSQ   =  0.8509E-5; // * T^2

      // Equations 17. (optical)
       static const double A1_OPT_1      =  0.5787089; // constant
       static const double A1_OPT_P      =  0.5609E-4; // * P
       static const double A1_OPT_T      = -0.6229E-3; // * T
       static const double A1_OPT_TSQ    =  0.2824E-5; // * T^2
       static const double A1_OPT_PPW    =  0.5177E-3; // * Pwater
       static const double A1_OPT_PPWSQ  =  0.2900E-6; // * Pwater^2
       static const double A1_OPT_WAVE   = -0.1644E-1; // * lambda
       static const double A1_OPT_WAVESQ =  0.4910E-1; // * lambda^2

       static const double A2_OPT_1      =  1.302474;  // constant
       static const double A2_OPT_P      =  0.2142E-4; // * P
       static const double A2_OPT_T      =  0.1287E-2; // * Pwater
       static const double A2_OPT_TSQ    =  0.6500E-6; // * Pwater^2
       static const double A2_OPT_PPW    = -0.2623E-2; // * T
       static const double A2_OPT_PPWSQ  =  0.8776E-5; // * T^2
       static const double A2_OPT_WAVE   = -0.6298E-2; // * lambda
       static const double A2_OPT_WAVESQ =  0.1890E-1; // * lambda^2

       // Zero point offsets used by Mangum (2001)
       // Same as Yan (1996) except temperature has been converted
       // to Kelvin. These are all to be ADDED to the measured
       // values.
       //
       // pressure offset, mbar
       static const double P_OFFSET    = -1013.25;
       // temperature offset, K
       static const double T_OFFSET    = -258.15; 
       // wavelength offset, micron
       static const double WAVE_OFFSET = -0.5320; 

       // coefficients for normalized effective zenith
       // Yan equation 12, Mangum equations 24 and 27.
       // R = 8314.14 J/(kmol K) = kg/(m^2 s^2 kmol K)
       // M = 28.970 kg/kmol
       // g = 9.7840 m/s^2 (in vertical column of air at tropopause)
       // R/(M*g) = 29.33895 meters
       // From units: 
       // You have: 1 K*gasconstant/(28.970 kg kmol^-1*9.7840 m s^-2)
       // You want: 
       //         Definition: 29.333895 m
       static const double RMG = 29.333895;

       // some constants used in the mapping function
       // see computeMappingFunction(...)
       static const double MAP1 = 13.24969;
       static const double MAP2 = 173.4233;

      //------------------------------------------------------
      //          PRIVATE NON-STATIC NON-CONSTANT DATA 
      //------------------------------------------------------
      carma::services::Temperature lastSafeAirTemperature_;
      carma::services::Pressure lastSafeAtmPressure_;
      carma::services::Temperature lastSafeDewPoint_;
      double lastSafeRelativeHumidity_;
       
      //------------------------------------------------------
      //            STRUCTS USED BY THIS CLASS
      //------------------------------------------------------
      
      typedef struct {
          double ir;

          double dry;

          double wet;

          double opacity;

      } scaleHeightType;

      static const scaleHeightType scaleHeight;

      const carma::services::Units units;

    //------------------------------------------------------
    //            PRIVATE METHODS
    //------------------------------------------------------
    double freqDependence(double frequency);

    double computeMappingFunction(
                   double airTemp, 
                   double atmPressure, 
                   double relHumid,
                   double elevation,
                   double frequency, 
                   double altitude) const;

    /*
     * This is to test against the Ulich (1981)
     * mapping function from NRAO.
     * @param elevation The <em>in vacuo</em> 
     *        (uncorrected) elevation in radians
     */
    double ulichMappingFunction(double elevation) const;

    /*
     * Calculates the normalized effective zenith parameter, I,
     * which is used in calculating the refraction mapping function.
     * This calculates I<sup>2</sup>csc(elevation).
     * @param elevation The <em>in vacuo</em> 
     *        (uncorrected) elevation in radians
     * @param  airTemp The ambient air temperature in Kelvin
     * @param  the altitude (station height) of the antenna above 
     *         the array reference plane, meters.
     * @return the normalized effective zenith parameter (unitless).
     * @see Mangum (2001) equation 27
     */
    double normalizedEffectiveZenith(double airTemp, 
                     double elevation, 
                     double altitude) const;

    bool isOptical(double frequency) const;

    double a1Radio(double p0, double t0, double ppw) const;

    double a2Radio(double p0, double t0, double ppw) const;

    double a1Optical(double p0, double t0, double ppw, double w0) const;

    double a2Optical(double p0, double t0, double ppw, double w0) const;

  };

  }
}
#endif //CARMA_SERVICES_ATMOSPHERE_H