import math

"""
The Purpose of this Program is to automate the task of Ventilation Surveying in Underground Mining Engineering.
This program can accept input of tabulated values (in csv file) for Leapfrogging and Roving-Base Altimeter
Indirect Method Surveys and provide data analysis and calculation output to a csv file.

Also provides many tools (Psychometric Properties of Air, Head Loss around Circuit, Specific Weight/Humidity Calc)
Analysis for Mine Ventilation Engineers.

pg 206, Mine Ventilation and Air Conditioning
Measured Values:
-Air Velocities (Vane Anemometer: v)
-Absolute and Differential Pressures or Heads (Manometer/Barometer/Altimeter: pb)
-Dry- and Wet-Bulb Temperatures (Thermometers: Tw, Td)
-Airway Dimensions (A)

Determine:
-Air Quantities
-Pressure Losses
-Air Specific Weights/Humidities
-Airway Resistance

Copyright Joaquin Roibal, August 2016, Latest Revision: 10/3/2016
All Rights Reserved
"""

def ManometerDirectMethod(TopList, BottomList, Manometer_Reading):
    """
    A Manometer is used to directly measure pressure difference (Direct Method).
    This function will perform the data analysis for Manometer Surveying given measured input values.
    Hl12 = (Hs1-Hs2)+(Hz1-Hz2)             #Head Loss (1->2) components
    Equations and Values from Example 6.3, page 208, Mine Ventilation and Air Conditioning Ramani and Wang
    """
    p3 = (BottomList[1] + (Manometer_Reading/27.69))  #Page 210, Pressure 3 = p2 + manometer reading (converted to psi)
    Ws = (TopList[-1]+BottomList[-1])/2               #Mean Specific Weight of Air in Hose
    Wh = (TopList[-1]+0.0705)/2                                  #Mean Specific Weight of Air in Shaft
    #Hl13 = (144*(Top[1]-p3)/5.2)+(Wh*(Top[0]-Bottom[0])/5.2)     #Page 209, Head Loss 1->3
    #Hl12 = (144*(Top[1]-Bottom[1])/5.2)+(Ws*(Top[0]-Bottom[0])/5.2) #Page 209, Head Loss 1 -> 2
    Hl12 = round((144/5.2)*(p3-BottomList[1])+(1/5.2)*(TopList[0]-BottomList[0])*(Ws-Wh), 3)
    return "Manometer Survey Results:\nWs: %s, Wh: %s, Head Loss in. Water: %s" %(Ws, Wh, Hl12)

def CalculateFrictionFactor(head_loss_f, length, diameter, quantity, spec_weight_air=0.075):
    """
    The Purpose of this function is to calculate the friction factor of an airway/ducting given parameters
    Utilizes Darcy-Weisbach equation and Atkinson Equation, 5.20 page 153 Mine Ventilation and Air Conditioning 3rd Edition
    """
    duct_perimeter = 2 * 3.14159 * (diameter / 2)
    area_opening = 3.14159 * (diameter / 2)**2
    rubbing_surface = length * duct_perimeter
    friction_factor_k = (spec_weight_air/0.075) * (head_loss_f*5.2*area_opening**3) / (duct_perimeter * length * quantity**2)
    return friction_factor_k
    
def NaturalVentilation(ShaftATop, ShaftABottom, ShaftBTop, ShaftBBottom):
    """
    The purpose of this function is to calculate the Natural Ventilation Head in Inches Water Gage.
    Inputs required: Lists in the following format: [DryBulbTemp, WetBulbTemp, Elevation, Pressure (in Hg)
    Method 2, page 297 in Ramani "Mine Ventilation And Air Conditioning" is used in commented example
    Equation used is from ME 440: Mine Ventilation with Dr. Bhattacharyya, ignoring vapor pressure
    :param ShaftATop:
    :param ShaftABottom:
    :param ShaftBTop:
    :param ShaftBBottom:
    :return:
    """
    """
    This Section is Commented Out Because NOT WORKING: Alternative Method Below
    spec_weight_air_shaft_a_top = psychrometricPropAir(ShaftATop[0], ShaftATop[1], ShaftATop[3])
    spec_weight_air_shaft_a_bottom = psychrometricPropAir(ShaftABottom[0], ShaftABottom[1], ShaftABottom[3])
    spec_weight_air_avg_upcast = (spec_weight_air_shaft_a_top[8] + spec_weight_air_shaft_a_bottom[8])/2
    spec_weight_air_shaft_b_top = psychrometricPropAir(ShaftBTop[0], ShaftBTop[1], ShaftBTop[3])
    spec_weight_air_shaft_b_bottom = psychrometricPropAir(ShaftBBottom[0], ShaftBBottom[1], ShaftBBottom[3])
    spec_weight_air_avg_downcast = (spec_weight_air_shaft_b_top[8] + spec_weight_air_shaft_b_bottom[8])/2
    L = ShaftBTop[2]-ShaftATop[2]
    print(L)
    print("Specific Weight Air Top A: ", spec_weight_air_shaft_a_top[9])
    print("Specific Weight Air Bottom A: ", spec_weight_air_shaft_a_bottom[9])
    print("Avg Spec Weight Upcast: ", spec_weight_air_avg_upcast)
    print("Avg Spec Weight Downcast: ", spec_weight_air_avg_downcast)
    inches_water_gage = (L/5.2)*(spec_weight_air_avg_downcast-spec_weight_air_avg_upcast)
    return inches_water_gage
    """
    
    #The Following Method Utilizes the equation from ME 440: Mine Ventilation by Dr. Bhattacharyya
    #NOTE: IGNORES VAPOR PRESSURE

    density_air_shaft_a_top = round((1.327/(460+ShaftATop[0]))*ShaftATop[-1], 6)
    print("Density Air Shaft A
 Top: ", density_air_shaft_a_top)
    density_air_shaft_a_bottom = round((1.327/(460+ShaftABottom[0])*ShaftABottom[-1]), 6)
    print("Density Air Shaft A Bottom: ", density_air_shaft_a_bottom)
    density_air_shaft_b_top = round((1.327
/(460+ShaftBTop[0])*ShaftBTop[-1]), 6)
    print("Density Air Shaft B Top: ", density_air_shaft_b_top)
    density_air_shaft_b_bottom = round((1.327/(460+ShaftBBottom[0])*ShaftBBottom[-1]), 6)
    print("Density Air Shaft B Bottom: ", density_air_shaft_b_bottom)
    density_avg_shaft_a = (density_air_shaft_a_bottom + density_air_shaft_a_top)/2
    density_avg_shaft_b = (density_air_shaft_b_bottom + density_air_shaft_b_top)/2

    pressure_diff = round(abs((density_avg_shaft_a - density_avg_shaft_b)), 6)
    elevation_diff = (ShaftBTop[-2]-ShaftABottom[-2])
    print("Pressure Difference: ", pressure_diff)
    print("Elevation Difference: ", elevation_diff)
    inches_water_gage = round((pressure_diff*elevation_diff)/5.2, 4)
    return inches_water_gage

def psychrometricPropAir(td, tw, pb):
    """
    The purpose of this function is to accept input of measured values (wet-bulb, dry-bulb temp, barometric pressure)
    to calculate the Psychrometric properties of Air (Spec Weight)and return a list of values calculated:
    Ps, Ps Prime, Pv, Phi, W, Ws, Mu, Pa, V, w, h . Will be used in other functions to calculate head loss, etc.

    Example Values and Equations from Page 13, Mine Ventilation And Air Conditioning Textbook by Ramani and Wang
    :param td: Dry Bulb Temperature
    :param tw: Wet Bulb Temperature
    :param pb: Pressure (in Hg)
    :return:
    """
    Td = (td + 459.67)       #Convert Temperature from Fahrenheit to Kelvin
    val_list = [td, Td, tw, pb]       #List of Values to be returned by function,
                        #Final Format for List: [td, Td, tw, pb, ps, ps_prime, pv, phi, W, Ws, mu, pa, v, w, h]

    m = 28.97           #Molecular Weight
    s = 1               #Specific Gravity
    R = 53.35           #ft*lb/lb mass*Degree Rankine, Gas Constant
    w = 0.0750          #lb/ft^3, Specific Weight at Standard Conditions
    standard_pb = 29.92 #in. Hg, Standard Barometric Pressure at Sea Level
    cp = 0.2403         #Btu/lb*degreeF, Specific Heat at Constant Pressure
    cv = 0.1714         #Btu/lb*degreeF, Specific Heat at Constant Volume
    gamma = 1.402       #Ratio of Spec heats at constant pressure and volume for diatomic gas


            #Calculate Saturation Vapor Pressures: (Page 15, Mine Ventilation and Air Conditioning)
    ps = 0.18079*math.exp((17.27*td-552.64)/(td+395.14))      #in. Hg, Saturation Vapor Pressure, Dry Bulb Temp, eq 2.2
    val_list.append(ps)
    ps_prime = 0.18079*math.exp((17.27*tw-552.64)/(tw+395.14))  #in. Hg, Saturation Vapor Pressure, Wet Bulb Temp
    val_list.append(ps_prime)

    pv = ps_prime - ((pb-ps_prime)*(td-tw))/(2800-1.3*tw)       #in. Hg, Partial Pressure of Water Vapor in Air, eq. 2.3
    val_list.append(pv)
    phi = pv/ps*100                         #Relative Humidity, eq. 2.4
    val_list.append(phi)

    W = 0.622*pv/(pb-pv)      #lb/lb dry air, Specific Humidity, Eq. 2.5
    val_list.append(W)
    W_grain = W*7000            #grains/lb dry air
    Ws = 0.622*ps/(pb-ps)       #lb/lb wet air, Specific Humidity, Eq. 2.5 (Wet