# -*- coding: utf-8 -*-
#TFY4165 Termisk fysikk 2024
#Øving 1
#Atmosfæren på Mars

import matplotlib.pyplot as plt
import numpy as np

#Konstanter
R = 8.314 #gasskonstanten (J/mol K)
g = 3.71  #tyngdens akselerasjon på Mars (m/s^2)
M = 0.04334 #midlere molare masse i atmosfæren på mars (kg/mol)
T0 = 234  #konstantleddet i T(z) (K)
p0 = 600  #midlere atmosfæretrykk på bakken (Pa)

H0 = R*T0/(g*M) #skalahøyde (m) med konst temp T0
z0 = 0  #bakkenivå
z1 = 16 #maks høyde over bakken (km)
N = 641 #antall z-verdier
h = (z1-z0)/(N-1) #steglengde (km)

#Nullstiller diverse arrays
integral = np.zeros(N)
Tinv = np.zeros(N) #invers temperatur (1/K)
height = np.zeros(N)
pressure = np.zeros(N)
pexp = np.zeros(N)
temperature = np.zeros(N)

for i in range(N):
    z=z0+(i-1)*h
    T=T0-2.25*z+14.0*np.exp(-2.0*z)
    Tinv[i]=1/T
    #Integrasjon med trapesmetoden og Simpsons metode
    #i=1 og i=2 behandles spesielt:
    if i==1:
        integral[i]=0.0
    elif i==2:
        integral[i]=(h/2)*(Tinv[1]+Tinv[2])
    elif np.mod(i,2)==0:     #i=2,4,6,...: trapesmetoden
        integral[i]=integral[i-2]+(h/2)*(Tinv[i-2] +
                    2*Tinv[i-1]+Tinv[i])
    else:                   #i=3,5,7,...: Simpsons metode
        integral[i]=integral[i-2]+(h/3)*(Tinv[i-2] + 
                    4*Tinv[i-1]+Tinv[i])

    height[i]=z
    #faktor 1000 for omregning fra km til m
    pressure[i]=p0*np.exp(-1000.0*(M*g/R)*integral[i])
    #med konstant H0 = RT/gM blir p(z) en ren exp-funksjon
    pexp[i]=p0*np.exp(-1000*z/H0)
    temperature[i]=T

    
#Plotting
plt.figure('Atmosfæren på Mars')
plt.plot(pressure,height,label='T(z) = 234 - 2.25z + 14 exp(-2z)')
plt.plot(pexp,height,label='T = 234 K')
plt.legend(loc='best')
plt.title('Trykket i atmosfæren på Mars',fontsize=20)
plt.ylabel('Høyde (km)',fontsize=24)
plt.xlabel('Trykk (Pa)',fontsize=24)
plt.xlim(0,600)
plt.ylim(0,16)
#plt.savefig("marstrykk.eps")
plt.grid()
plt.show()

#Plotting
plt.figure('Atmosfæren på Mars')
plt.plot(height,pressure,label='T(z) = 234 - 2.25z + 14 exp(-2z)')
plt.plot(height,pexp,label='T = 234 K')
plt.legend(loc='best')
plt.title('Trykket i atmosfæren på Mars',fontsize=20)
plt.xlabel('Høyde (km)',fontsize=24)
plt.ylabel('Trykk (Pa)',fontsize=24)
plt.ylim(0,600)
plt.xlim(0,16)
#plt.savefig("marstrykk1.eps")
plt.grid()
plt.show()