using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;

namespace blahgeopotpressuretest
{
    class Program
    {
        protected static double VaporPressLiquid(double T)
        {
            //convert the input temperature from deg C to kelvin
            T += 273.15;

            double Psat = Math.Exp(-0.58002206E4 / T + 0.13914993E1 - 0.48640239E-1 * T + 0.41764768E-4 * Math.Pow(T, 2.0)
                - 0.14452093E-7 * Math.Pow(T, 3.0) + 0.65459673E1 * Math.Log(T)) / 100.0;

            return Psat;
        }

        protected static double VirtualTemp(double p, double t, double u)
        {
            double Es,         // saturated vapor pressure (mb)
                   Tv;         // calculated virtual temperature (K)

            if ((u > 0.0) && (u < 100.0) && (p > 0.0))
            {                            // validate P & U data
                Es = VaporPressLiquid(t); //use the Hyland and Wexler vapor pressure equation (previously used the Buck 1981 equation)
                Tv = (t + 273.15) * (1 + 0.37821 * ((u * Es) / (100.0 * p)));
            }
            else
            {
                Tv = t + 273.15;
            }
            return Tv;
        }

        protected static double CalculateGeoPotPressure(double currentTemperature, double currentHumidity, double currentAltitude, double previousTemperature, double previousHumidity, double previousAltitude, double previousGeoPotPressure)
        {
            //virtual temperature (K)
            double virtualTemperature;
            //calculated pressure (mb)
            double currentGeoPotPressure;

            currentGeoPotPressure = previousGeoPotPressure;
            if (previousGeoPotPressure > 0.0)
            {
                if (previousAltitude != currentAltitude)
                {
                    //average virtual temperature using known pressure
                    virtualTemperature = (VirtualTemp(previousGeoPotPressure, previousTemperature, previousHumidity) + VirtualTemp(previousGeoPotPressure, currentTemperature, currentHumidity)) / 2.0;
                    currentGeoPotPressure = previousGeoPotPressure * Math.Exp((-9.80665 * (currentAltitude - previousAltitude)) / (287.04 * virtualTemperature));
                }
                else
                {
                    currentGeoPotPressure = previousGeoPotPressure;
                }
            }

            return currentGeoPotPressure;
        }

        static void Main(string[] args)
        {
            char[] splitter = new char[] { ' ', '\t' };
            StreamReader reader = new StreamReader(new FileStream("GeoPotPressTest.txt", FileMode.Open));
            StreamWriter writer = new StreamWriter(new FileStream("output.csv", FileMode.Create));
            writer.WriteLine("RH [%], Temperature [deg C], Altitude [km], Geopot Pressure [mb]");

            double previousTemperature = double.NaN;
            double previousHumidity = double.NaN;
            double previousAltitude = double.NaN;
            double previousGeopotPressure = double.NaN;

            reader.ReadLine();
            while (!reader.EndOfStream)
            {
                string[] lineArr = reader.ReadLine().Split(splitter, StringSplitOptions.RemoveEmptyEntries);
                if (lineArr.Length == 3)
                {
                    double temperature = double.Parse(lineArr[1]);
                    double humidity = double.Parse(lineArr[0]);
                    double altitude = double.Parse(lineArr[2]) * 1000;

                    if (double.IsNaN(previousTemperature))
                    {
                        previousTemperature = temperature;
                        previousHumidity = humidity;
                        previousAltitude = altitude;
                        previousGeopotPressure = 973.72;
                    }

                    double geopotPressure = CalculateGeoPotPressure(temperature, humidity, altitude, previousTemperature, previousHumidity, previousAltitude, previousGeopotPressure);

                    writer.WriteLine("{0:0.0}, {1:0.00}, {2:0.000}, {3:0.000}", humidity, temperature, altitude/1000, geopotPressure);

                    previousTemperature = temperature;
                    previousHumidity = humidity;
                    previousAltitude = altitude;
                    previousGeopotPressure = geopotPressure;
                }
            }

            reader.Close();
            writer.Close();
        }
    }
}