# -*- coding: utf-8 -*-
#   Atmosfæærenà jorda
#   Termisk Fysikk 2024

import matplotlib.pyplot as plt
import numpy as np

R=8.314    # gasskonstanten (J/mol K)
g=9.81     # tyngdens akselerasjon jorda (m/s^2)
M=0.029    # midlere molar masse i atmosfaeren paa jorda (kg/mol)
T0=288     # konstant verdi brukt for temperatur (K)

H0=R*T0/(g*M)   # skalahoyde (m) med konstant temperatur T0

z0=0.001      # velger hoyde z0 = 0.001 like over bakken (km),
              # ikke z0=0, for åaa unngaa aa dele med null nedenfor
z1=350.0      # velger maksimal hoyde over bakken z1 = 350 km
N=100         # antall integrasjonspunkter
h=(z1-z0)/(N-1)   # intervall-lengde (km)

#p350 = publisert verdi for trykket i hoyde 350 km
p350 = 2.5e-6  #Pa
#p0 = kjent verdi for trykket paa bakkenivaa
p0 = 1.0e5   #Pa
#initialisering av verdi for integralet
F=0  
#tilpasningsparametre for T(z)
a1 = 275
a2 = 8.4
a3 = 0.290
a4 = 0.002565
a5 = 0.000
a6 = 90
a7 = 5
a8 = 1000
a9 = 288
# nullstiller for beregning av integralet \int_0^z dz'/T(z')
integral=zeros(N)
# nullstiller for Tinv = 1/T
Tinv=zeros(N)
#for plotting
height=zeros(N)
pressure=zeros(N)
pdown=zeros(N)
pup=zeros(N)
temperature=zeros(N)
for i in range(1,N,1):
     z=z1-(i-1)*h
     #Tilpasning til maalt T(z)
     Thigh=(a8*(a1/z)^2*exp(-a1/z)/(-exp(-a1/z)+1)^2)/(1+exp((a6-z)/a7))
     Tlow=(a9-a2*z+a3*z^2-a4*z^3)/(1+exp((z-a6)/a7))
     T = Tlow + Thigh
     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
     pressure[i]=p350*np.exp(1000.0*(M*g/R)*integral[i])
     pdown[i]=p350*np.exp((z1-z)*1000/H0)
     pup[i]=p0*np.exp((z0-z)*1000/H0)
     temperature[i]=T
 
#plotter trykk mot høyde i semilog-plott
#semilogy(height,pressure,height,pdown,height,pup);
#legend('tilpasset T','T = 288 K, kjent p(350km)','T = 288 K, kjent p(0)');
#ylabel('Trykk (Pa)');
#xlabel('Høyde (km)');
#grid on;
#figure;
#%plotter temperatur mot høyde i lineært plott (evt loglog plott)
#%loglog(height,temperature);
#plot(height,temperature);
#ylabel('Temperatur (K)');
#xlabel('Høyde (km)');
#grid on;