#ifndef BAROMETRIC
#define BAROMETRIC

const static double mslp = 980.25; // mean sea level spressure
const static double int_lapse_rate = 0.0065; // a
const static double int_exponent = 5.255; // international barometric formular exponent calculated from (M * g) / (R * a)

//The height of the single floor levels.
const static double floor_height[5] = {0.0, 4.1, 3.4, 3.4, 3.4};


class BarometricFormular
{
    private:
        const double temperature; // T in Kelvin

        const double universal_gas_constant; // R
        const double molar_mass; // M
        const double gravitational_acceleration; // g
        const double lapse_rate; // a
        double _exponent;

    public:

        /** ctor */
        BarometricFormular(const double R, const double M, const double g, const double a, const double T):
        universal_gas_constant(R), molar_mass(M), gravitational_acceleration(g), lapse_rate(a), temperature(T){
            _exponent = (M * g) / (R * a);
        }

        /** ctor only with Temperature*/
        BarometricFormular(const double T) :
        universal_gas_constant(8.314), molar_mass(0.02896), gravitational_acceleration(9.80665), lapse_rate(0.0065), temperature(T){
            _exponent = (molar_mass * gravitational_acceleration) / (universal_gas_constant * lapse_rate);
        }

        /** Atmospheric Pressure Calculation */
        double getAtmosphericPressure(double p_0, double h_1) const{
            return p_0 * std::pow((1.0 - ((lapse_rate * h_1)/temperature)), _exponent);
        }

        /** Atmospheric Pressure Calculation above sea level*/
        double getAtmosphericPressure(double h_1) const{
            return mslp * std::pow((1.0 - ((lapse_rate * h_1)/temperature)), _exponent);
        }

        //TODO:: Height from pressure for the general formular

        //Static Functions

        /** International Barometric Formular*/
        static double s_getAtmosphericPressure(double p_0, double h_1, double kelvin){
            return p_0 * std::pow((1.0 - ((int_lapse_rate * h_1)/kelvin)), int_exponent);
        }

        /** International Barometric Formular above Sea Level*/
        static double s_getAtmosphericPressure(double h_1, double kelvin){
            return mslp * std::pow((1.0 - ((int_lapse_rate * h_1)/kelvin)), int_exponent);
        }

        /** International Barometric Formular above Sea Level at 15 degree*/
        static double s_getAtmosphericPressure(double height_above_sea_level){
            return mslp * std::pow((1.0 - ((int_lapse_rate * height_above_sea_level)/288.15)), int_exponent);
        }

        /** Get the height above sea level using a pressure measurment above sea level*/
        static double getHeightAboveSeaLevel(double p, double kelvin){
            // http://www.wolframalpha.com/input/?i=solve+for+h+++p+%3D+980.25*%281+-+0.0065+*+h%2FT%29^5.255
            return 41.4811 * ((3.70882 * kelvin) - (std::pow(p, 0.1902949571836346) * kelvin));
        }


        /** This is a helper Class only for gnupplot visualization for ipin2015*/
        static double getHeightForVisualizationOnly(double p, double z_0, double kelvin){

            // the height of the reference (first) pressure measurement
            double h_0 = 0.0;
            for(int i = 0; i <= z_0; i++){
                h_0 += floor_height[i];
            }

            // pressure value of h_0 above sea level
            // we define that the bottom of floor 0 is sea level ;).
            double p_0 = s_getAtmosphericPressure(h_0, kelvin);

            // pressure value of the current particle above floor 0 (sea level)
            double p_height = p_0 + p;

            // height of the particle above floor 0 (sea level)
            return getHeightAboveSeaLevel(p_height, kelvin);

        }

        static double s_getSeaLevelPressure(){
            return mslp;
        }

        static double getPressureOfFloorForVizualization(double z, double kelvin){

            int i = z + 0.5;
            double h_z = floor_height[i+1];
            double p_z = s_getAtmosphericPressure(h_z, kelvin);
            return std::abs(mslp - p_z);
        }

};

#endif // BAROMETRIC