// airdensity.cc
//
// Michael Vera Oct 2005
//
// Given Dew point (dp), Temperature in Celsius (t), and
//  Measured Air Pressue (p) in millibars mb, return
//  Air Density in kg/m^3
//
// Determine Air Density given standard atmospheric variables
//  obtained from weahter stations of websites
//  (e.g. http://weatherunderground.com) such as
//  Humidity, Dew point, Temperature, and Pressure
//
// Density              = Density of Air (D)
// Pressure             = Total Pressure in Pascals (*100 mb to get Pascals) (P)
// DryAirPressure       = Dry Air Pressure (Pd)
// WaterVaporPressure   = Water Vapor Pressure (Pv)
// GasConstantDry       = Gas Constant for Dry Air, 287.05 J//(kg*degK) (Rd)
// GasConstantVapor     = Gas Constant for Water Vapor, 461.495 J/(kg*degK) (Rv)
// Temperature          = Temperature, degK = deg Celsius + 273.15 (T)
//
// 1. Density   = (DryAirPressure/(GasConstantDry*(Temperature+273.15))) +
//                  (WaterVaporPressure/(GasConstantVapor*(Temperature+273.15)))
//
// SaturationPressure   = Saturation Pressure of Water Vapor, mb or hPa (es)
// WaterVaporPressureAtDewPoint = Pv when Temperature is at Dew Point (Es)
// eso  = 6.1078
// p    = (c0+T*(c1+T*(c3+T*(c4+T*(c5+T*(c6+T*(c7+T*(c8+T(c9)))))))))
//REQUIRED EXTERNAL VARIABLE
// T    = Temperature in Celsius
// c0   = 0.99999683
// c1   = -0.90826951*(10**(-2))
// c2   = 0.78736169*(10**(-4))
// c3   = -0.61117958(10**(-6))
// c4   = 0.43884187(10**(-8))
// c5   = -0.29883885(10**(-10))
// c6   = 0.21874425*(10**(-12))
// c7   = -0.17892321*(10**(-14))
// c8   = 0.11112018*(10**(-16))
// c9   = -0.30994571*(10**(-19))
//

#include <iostream.h>

int main()
{

  double
    Temperature,
    PressureMeasured,
    DewPoint;

  // Get Temperature as input from the user
  cout << "\nPlease enter the current Temperature (Celsius): ";
  cin >> Temperature;

  // Get PressureMeasured as input from the user
  cout << "\nPlease enter the current Measured Pressure (hPa): ";
  cin >> PressureMeasured;

  // Get Dew Point as input from the user
  cout << "\nPlease enter the current Dew Point (Celsius): ";
  cin >> DewPoint;

  const double
    eso = 6.1078,
    c0   = 0.99999683,
    c1   = -0.90826951e-2,
    c2   = 0.78736169e-4,
    c3   = -0.61117958e-6,
    c4   = 0.43884187e-8,
    c5   = -0.29883885e-10,
    c6   = 0.21874425e-12,
    c7   = -0.17892321e-14,
    c8   = 0.11112018e-16,
    c9   = -0.30994571e-19;

  // 2.  Es     = eso*(p**(-8))
  //
  // Es = Pv when Temperature is at Dew Point
  // Es checked against 1958 Smithsonian Meterological Tables,
  //  Smithsonian Institute, Washington, D.C. for accuracy
  double p;
  p = c0+(DewPoint*(c1+(DewPoint*(c2+(DewPoint*(c3+(DewPoint*(c4+(DewPoint*(c5+(DewPoint*(c6+(DewPoint*(c7+(DewPoint*(c8+(DewPoint*c9)))))))))))))))));

  double
    SaturationPressure,
    DryAirPressure,
    WaterVaporPressure,
    Density;

  SaturationPressure = eso * pow(p, -8);
  WaterVaporPressure = SaturationPressure;

  DryAirPressure = PressureMeasured - WaterVaporPressure;

  const float
    GasConstantDry = 287.05,
    GasConstantVapor = 461.495;

  // Air Density Calculation
  Density = (DryAirPressure*100/(GasConstantDry*(Temperature+273.15))) +
              (WaterVaporPressure*100/(GasConstantVapor*(Temperature+273.15)));

  cout << "\nThe Air Density is " << Density << " kg/m^3.\n\n";

  return 0;
}