import numpy as np
from matplotlib import pyplot as plt

alpha = 149597870.691  #km to AU conversion
mPs_to_auPday = 86.4/alpha  #m/s to au/day conversion
kg_to_g = 1000      #kg to g conversion
Ndays = 80
TempValues = 40 #1000
mN2 = 28/kg_to_g #Mass of N2 gas in kg/mole
R = 8.314*mPs_to_auPday**2 #Joules/mole*Kelvin = kg*m^2/(s^2*mol*K) -> kg*au^2/(day^2*mol*K)
tsteps = 80

day = np.zeros(Ndays,dtype=float)
V = np.zeros(Ndays,dtype=float)
VO = np.zeros(Ndays,dtype=float)
diffV = np.zeros(tsteps, dtype=np.float64)
MassRatio = np.zeros(tsteps*TempValues, dtype=np.float64).reshape(TempValues, tsteps)

f = open('Velocities2.csv','r')
f.readline()
for i in range(0,Ndays):
    d, vtemp, votemp = f.readline().split(',')
    #print(d,vtemp,votemp)
    day[i] = float(d)
    V[i] = float(vtemp)
    VO[i] = float(votemp)
f.close

#Defining function to calculate the velocity of nitrogen gas\
#Uses Boltzman-Weighted Gas Velocity and D&J's Efficiency Factor tau
def CalcN2V(MassRatio,TempValues,R,mN2,diffV):
    from numpy import exp
    from numpy import pi
    tau = 0.45
    for Temp in range(0,TempValues):
        #ve will be in au per day
        ve = tau*np.sqrt(8/pi*R*(Temp+1)/mN2)
        #print(ve)
        for t in range(0, tsteps):
            MassRatio[Temp, t] = exp(-diffV[t]/ve)

CalcN2V(MassRatio,TempValues,R,mN2,diffV)