#ifndef __ATMOSPHERE_H__
#define __ATMOSPHERE_H__

#include <cmath>

/**
 * 標準大気モデル
 *
 */
class InternationalStrandardAtmosphere {
  public:
    static const double p0; ///< Pressure at sea level
    static const double rho0; ///< Density at sea level
    static const double T0; ///< Temperature at sea level
    static const double a0; ///< Speed of sound at sea level
    static const double g0; ///< Acceleration of gravity at sea level
    static const double R; ///< Real gas constant for air
    static const double gamma; // adiabatic index @see http://en.wikipedia.org/wiki/Speed_of_sound
    static const double L; //

    /**
     * ある高度における気温を求める
     *
     * @param h 高度(m) TODO: 高度の拡張
     * @param deltaT 温度のずれ(K)
     * @return 気温(K)
     */
    static double T(const double &h, const double &deltaT = 0) {
      return T0 + (L * h) + deltaT;
    }
    /**
     * ある高度における気圧を求める
     *
     * @param h 高度(m)
     * @return 気圧(N/m^2)
     */
    static double p(const double &h){
      return p0 * std::pow((1.0 + L * h / T0), -g0 / (L * R));
    }
    /**
     * ある高度における密度を求める
     *
     * @param h 高度(m)
     * @param deltaT 温度のずれ(K)
     * @return 密度(kg/m^3)
     */
    static double rho(const double &h, const double &deltaT = 0) {
      return p(h) / (R * T(h, deltaT));
    }

    /**
     * ある高度における音速を求める
     *
     * @param h 高度(m)
     * @param deltaT 温度のずれ(K)
     * @return 音速(m/s)
     */
    static double mach1(const double &h, const double &deltaT = 0){
      return std::sqrt(gamma * R * T(h, deltaT));
    }

    // TODO: 要検証 {
    static double tas(const double &staticT, const double &p_static, const double &p_total){
      return std::sqrt((2.0 * gamma * R * staticT / (gamma - 1))
          * (std::pow((p_total / p_static), ((gamma - 1) / gamma)) - 1) );
    }

    static double eas(const double &p_static, const double &p_total){
      return tas(T0, p_static, p_total);
    }

    static double cas(const double &p_dynamic){
      return tas(T0, p0, p_dynamic + p0);
    }

    static double tas2eas(const double &tas, const double &staticT){
      return tas * std::sqrt(T0 / staticT);
    }

    static double eas2cas(const double &eas, const double &p_static){
      double p_total(p_static * std::pow((eas * eas / ((2.0 * gamma * R * T0) / (gamma - 1)) + 1), (gamma / (gamma - 1))) );
      return cas(p_total - p_static);
    }

    static double cas2eas(const double &cas, const double &p_static){
      double p_dynamic(p0 * (std::pow((cas * cas / ((2.0 * gamma * R * T0) / (gamma - 1)) + 1), (gamma / (gamma - 1))) - 1));
      return eas(p_static, p_static + p_dynamic);
    }

    static double eas2tas(const double &eas, const double &staticT){
      return eas * std::sqrt(staticT / T0);
    }
    // TODO: 要検証 }
};

const double InternationalStrandardAtmosphere::p0 = 101325; // N/m^2
const double InternationalStrandardAtmosphere::rho0 = 1.225; // kg/m^3
const double InternationalStrandardAtmosphere::T0 = 288.15; // K
const double InternationalStrandardAtmosphere::a0 = 340.294; // m/s
const double InternationalStrandardAtmosphere::g0 = 9.80665; // m/s^2
const double InternationalStrandardAtmosphere::R = 287.04; //m^2 / K s^2
const double InternationalStrandardAtmosphere::gamma = 1.4; //
const double InternationalStrandardAtmosphere::L = -0.0065; //

#endif /* #ifndef __ATMOSPHERE_H__ */