#Processing WVR tipping curve data
#BEGIN-------------------------------------------------------------------------
#------------------------------------------------------------------------------
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy import stats
#------------------------------------------------------------------------------
#Retrieving the antenna temperature using eq 5 from the technical documentation
def RetrieveAntTemp(V1, Vh, Vc, Th, Tc):
    g = (Th-Tc)/(Vh-Vc)
    o = Tc - (g*Vc)
    return(g*V1) + o
#------------------------------------------------------------------------------
def DataExtract(data, Th_23, Th_31): 
    #Function to extract the angle and and the brightness tempeatures for each angle
    angles = [0]*8
    temp_23 = [0]*8
    temp_31 = [0]*8
    for i in range(0,8):
         angles[i] = float(data[i][2])
         Tc = data[i][10]
         temp_23[i] = RetrieveAntTemp(data[i][8], data[i][7], data[i][6], Th_23, Tc)
         Tc = data[i][12]
         temp_31[i] = RetrieveAntTemp(data[i][5], data[i][4], data[i][3], Th_31, Tc)
    return angles, temp_23, temp_31   
#------------------------------------------------------------------------------
#Get the line characteristics
def LineFit(x_data, y_data):
    slope, intercept, r, p, std_err = stats.linregress(x_data, y_data)
    return slope, intercept
#------------------------------------------------------------------------------
def PlotTipCurves(angles, temp_23, temp_31, title):
    #Function to plot the tip curve
    plt.rcParams['font.family'] = 'Times New Roman'
    plt.rcParams['font.size'] = 12
    
    airmass_2 = [0]*9
    airmass = [0]*8
    for i in range(8):
        airmass[i] = Airmass(angles[i])
        airmass_2[i+1] = Airmass(angles[i])
    slope_23, intercept_23 = LineFit(airmass, temp_23)
    slope_31, intercept_31 = LineFit(airmass, temp_31)
    line_23 = [0]*9
    line_31 = [0]*9
    for i in range(len(airmass_2)):
        line_23[i] = (airmass_2[i] * slope_23) + intercept_23
        line_31[i] = (airmass_2[i] * slope_31) + intercept_31
    print("23 GHz 0 airmass intercept of ", str(intercept_23), " K")
    print("31 GHz 0 airmass intercept of ", str(intercept_31), " K")
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 6))
    #fig.suptitle(title)
    ax1.plot(airmass_2, line_23, "--", color = "black")
    ax1.plot(airmass, temp_23, 'x', color="blue")
    ax1.set_title("23 GHz channel", font="Times New Roman")
    ax1.set_xlabel("Airmass")
    ax1.set_ylabel("Brightness temperature (K)")
    ax1.grid()
    ax2.plot(airmass_2, line_31, "--", color = "black")
    ax2.plot(airmass, temp_31, 'x', color="blue")
    ax2.set_title("31 GHz channel", font="Times New Roman")
    ax2.set_xlabel("Airmass")
    ax2.grid()
    plt.show()
    
    return
#------------------------------------------------------------------------------
#Calculate airmass
def Airmass(A):
   return 1/np.sin(np.pi*A/180) 
#------------------------------------------------------------------------------
#Reading in the data
data = []
with open("T240312.csv") as data_file:
    i = 0
    for line in data_file:
        data_1 = line.split(",")
        data_2 = [0]*len(data_1)
        for i in range(len(data_1)):
            if (i == 1) or (i == 2):
                data_2[i] = data_1[i]
            else:
                data_2[i] = (float(data_1[i]))
        data.append(data_2)  
TippCurve = []
for i in range(8):
    TippCurve.append(data[i])
print(TippCurve) 
#------------------------------------------------------------------------------
Th_nom = 700 # nominal or starting guess for hot load temperature
angles, temp_23, temp_31 = DataExtract(TippCurve, Th_nom, Th_nom)

#Creating an array of arrmass form the angles
airmass = [0]*8
for i in range(8):
    airmass[i] = Airmass(angles[i])
    
#We can plot the uncalibrated tip curve
print("Raw tip curve data")
print("Expected hot load temperature: 700 K")
PlotTipCurves(angles, temp_23, temp_31, "Uncalibrated Tip curves")
print("")