Check Atmospheric profiles
author : Sylvie Dagoret-Campagne
affiliation : IJCLab/IN2P3/CNRS
date : November 27th 2018
Last update : November 7th 2023
These curves indicates the atmospheric model in libradtran is compatible with analytical physical formula.
[1]:
%load_ext autoreload
%autoreload 2
[2]:
import matplotlib.pyplot as plt
%matplotlib inline
import sys
import os
from astropy.io import fits
import numpy as np
import pandas as pd
[3]:
import matplotlib.colors as colors
import matplotlib.cm as cmx
[4]:
from astropy.io import fits
[5]:
# to enlarge the sizes
params = {'legend.fontsize': 'x-large',
'figure.figsize': (14, 10),
'axes.labelsize': 'x-large',
'axes.titlesize':'x-large',
'xtick.labelsize':'x-large',
'ytick.labelsize':'x-large'}
plt.rcParams.update(params)
[6]:
from analyticalmodels import libatmscattering as ray
#from analyticalmodels import N_A,R,M_air
N_A=6.0221409e23 # mol-1
M_air= 28.965338*1e-3 # u.g/u.mol (kg/mol)
m_air = M_air/N_A # single air molecule mass in grams
Configuration
[7]:
# from libradtran outout logs in verbose mode
# preselected sites
Dict_Of_sitesAltitudes = {'LSST':2.663,
'CTIO':2.207,
'OHP':0.65,
'PDM':2.8905,
'OMK':4.205,
'OSL':0.000,
}
# pressure calculated by libradtran at ground
Dict_Of_sitesPressures = {'LSST':731.50433,
'CTIO':774.6052,
'OHP':937.22595,
'PDM':710.90637,
'OMK':600.17224,
'OSL':1013.000,
}
# air molecule density at ground per cm3
Dict_Of_sitesMoleculesNDensity = {'LSST': 1.95597e19,
'CTIO': 2.04877e19,
'OHP': 2.39057e19,
'PDM': 1.91126e19,
'OMK': 1.66644e19,
'OSL': 2.54582e19,
}
# in grams per cm2 at ground
Dict_Of_sitesAirWeight = {'LSST':1.6e25,
'CTIO':1.6e25,
'OHP':2.0e25,
'PDM':1.5e25,
'OMK':1.3e25,
'OSL':2.2e25,
}
[8]:
dict2 = Dict_Of_sitesAltitudes
dict1 = Dict_Of_sitesPressures
dict3 = Dict_Of_sitesMoleculesNDensity
list_Of_sitesTags = [ 'OMK', 'PDM', 'LSST', 'CTIO', 'OHP', 'OSL']
list_Of_sitesNames = ['Mauna Kea', 'Pic du Midi','Rubin-LSST','CTIO','OHP','Sea Level']
List_Of_sitesPressures = [ dict1[sitetag] for sitetag in list_Of_sitesTags]
List_Of_sitesAltitudes = [ dict2[sitetag] for sitetag in list_Of_sitesTags]
List_Of_sitesAirNDensity = [ dict3[sitetag] for sitetag in list_Of_sitesTags] # number of molecules per cm3
NSites = len(list_Of_sitesTags)
[9]:
List_Of_siteMassDensity = [ n*m_air*1e6 for n in List_Of_sitesAirNDensity]
[10]:
lsst_altitude = 2.663 #km
lsst_pressure = 743 # hPa
airmass
[11]:
z_airmass=1.0
[12]:
all_airmass=np.array([1.0, 1.1, 1.2, 1.3, 1.4, 1.5,1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.5])
[13]:
NBSIMU=all_airmass.shape[0]
Seasons
[14]:
#YEARSTR="2018"
Winter
[15]:
pwv_w=10.2
o3_w=370.1
P_w=932.8
aer_w=0
Summer
[16]:
pwv_s=25.1
o3_s=318.1
P_s=938.9
aer_s=0
Altitude vs Pressure
[17]:
altitudes=np.linspace(0.,10000.,5000)
press=ray.Pressure_adiabatic(altitudes)
press=press/100.
[18]:
bbox = dict(boxstyle="round", fc="0.8")
bbox = dict(boxstyle="round", fc="0.8")
arrowprops = dict(facecolor='black',arrowstyle="->")
#connectionstyle="angle,angleA=0,angleB=90,rad=10")
plt.figure(figsize=(6,10))
plt.title("Adiabatic model for altitude vs pressure")
figname="AtmModel.png"
plt.plot(press,altitudes/1000,'k-')
plt.plot(List_Of_sitesPressures,List_Of_sitesAltitudes,'or',ms=5,label="Obs sites")
for idx in range(NSites):
pointed = (List_Of_sitesPressures[idx],List_Of_sitesAltitudes[idx])
point_text = (List_Of_sitesPressures[idx]-300,List_Of_sitesAltitudes[idx]-0.05)
#plt.annotate(list_Of_sitesNames[idx], xy=pointed, xytext=point_text,
# arrowprops=dict(facecolor='black', shrink=0.01),
# )
plt.annotate(list_Of_sitesNames[idx], xy=pointed, xytext=point_text,bbox=bbox, arrowprops=arrowprops)
plt.xlabel("Pressure (hPa)")
plt.ylabel("Altitude (km)")
#plt.axhline(lsst_altitude,color="red")
#plt.axvline(lsst_pressure,color="red")
plt.legend()
plt.grid()
plt.savefig(figname)
Altitude vs density
L’air humide est moins lourd (dense) que l’air sec
[19]:
T=273+20.
Hr=0.9
mass_density=ray.MassDensity_adiabatic(altitudes,T)
mass_density_h=ray.MassDensity_adiabatic_humid(altitudes,Hr,T)
kgm3_to_gl=1
[20]:
plt.figure(figsize=(6,10))
plt.title("Adiabatic model for altitude vs mass density")
figname="AirMassDensity.png"
plt.plot(mass_density,altitudes/1000,'k-',label='dry air')
plt.plot(mass_density_h,altitudes/1000,'r:',label='humid air')
# Not understood
#for idx in range(NSites):
# pointed = (List_Of_siteMassDensity[idx],List_Of_sitesAltitudes[idx])
# point_text = (List_Of_siteMassDensity[idx]-0.3,List_Of_sitesAltitudes[idx]-0.05)
# plt.annotate(list_Of_sitesNames[idx], xy=pointed, xytext=point_text,bbox=bbox, arrowprops=arrowprops)
plt.xlabel("Mass density (kg/m^3 or g/l)")
plt.ylabel("Altitude (km)")
plt.legend()
plt.grid()
plt.savefig(figname)
[ ]: