mirror of
https://gitlab.astro-wise.org/micado/atmosphericmodel.git
synced 2026-10-29 20:51:05 +00:00
353 lines
13 KiB
Python
353 lines
13 KiB
Python
#!/usr/bin/env Python
|
|
from __future__ import division
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from mpl_toolkits.axes_grid1 import host_subplot
|
|
from AstroAtmosphere import *
|
|
from tabulate import tabulate
|
|
|
|
|
|
|
|
'''
|
|
@ Author: Joost van den Born
|
|
@ Contact: born@astron.nl
|
|
@ Description: An exploration of the parameter space for the Ciddor-Cassini model.
|
|
The results here are given in table 3 and Appendix A of the paper.
|
|
|
|
'''
|
|
|
|
|
|
if __name__ == '__main__':
|
|
|
|
plt.rc('font', family='serif')
|
|
plt.rc('axes', labelsize='large')
|
|
|
|
at = Observatory()
|
|
|
|
pdrFigs = False
|
|
errorSpace = False
|
|
generate_images = True
|
|
|
|
|
|
T = 279.65 # k
|
|
P = 71200 # Pa
|
|
H = 0.22
|
|
xc = 450 # ppm
|
|
lat = -24.5983 # degrees
|
|
h = 3064 # m
|
|
l1 = 1.49 # micron
|
|
l2 = 1.78 # micron
|
|
zenith = 30 # degrees
|
|
|
|
rho = at.rho(p=P, T=T, RH=H, xc=xc)
|
|
disp = dispersion(lat, h)
|
|
disp.setReducedHeight(P, rho)
|
|
|
|
# Change dl, dT, dP, dRH, dCO2 to the values of interest.
|
|
n1 = at.n_tph(l=1.49, T=279.65, p=71200, RH=0.22, xc=450)
|
|
n2 = at.n_tph(l=1.78, T=279.65, p=71200, RH=0.22, xc=450)
|
|
R1ca = at.dn_tph(l=1.49, T=279.65, p=71200, RH=0.22, xc=450,
|
|
dl=0.001, dT=0.4, dP=10, dRH=0.01, dCO2=20)
|
|
R2ca = at.dn_tph(l=1.78, T=279.65, p=71200, RH=0.22, xc=450,
|
|
dl=0.001, dT=0.4, dP=10, dRH=0.01, dCO2=20)
|
|
|
|
|
|
# print (n1, at.dn_tph(l=1.49, T=288.0, p=101325, RH=0.0, xc=450, dl=0.0, dT=0.5, dP=0.0, dRH=0.0, dCO2=0) )
|
|
# print(disp.cassiniError( zenith, l1, l2, T, P, H, xc, dl1=0.001, dl2=0.001, dT=0.4, dP=10, dRH=0.01, dCO2=20 )*3.6e9)
|
|
|
|
zenith_true = 30 # degrees
|
|
dispValue = disp.cassini(n1, n2, zenith_true)
|
|
defaultDispCA = disp.cassini(n1, n2, zenith_true)*3.6e6
|
|
errorDispValue = disp.cassiniError(zenith, l1, l2, T, P, H, xc, dl1=0.001, dl2=0.001, dT=0.4, dP=10, dRH=0.01, dCO2=20 )
|
|
print ('Dispersion between 1.49 and 1.78 micron at z = %i deg: %.3f +/- %.3f mas' %(zenith_true, dispValue*3.6e6, errorDispValue*3.6e6))
|
|
|
|
|
|
R1s = at.n_tph(l=1.49, T=288.15, p=101325, RH=0.0, xc=450)
|
|
R2s = at.n_tph(l=1.78, T=288.15, p=101325, RH=0.0, xc=450)
|
|
|
|
# Varying the parameters in the non-errror formula:
|
|
if pdrFigs:
|
|
Trange = np.linspace(262, 298, 1000)
|
|
Prange = np.linspace(69000, 73000, 1000)
|
|
Hrange = np.linspace(0.1, 0.3, 1000)
|
|
|
|
fig = plt.figure()
|
|
|
|
# Variation in T
|
|
dispListS, dispListCA = [], []
|
|
for t in Trange:
|
|
R1ca = at.n_tph(l=1.49, T=t, p=71200, RH=0.22, xc=450)
|
|
R2ca = at.n_tph(l=1.78, T=t, p=71200, RH=0.22, xc=450)
|
|
dispValueS = disp.cassini(R1s, R2s, zenith_true) * 3.6e6
|
|
dispValueCA = disp.cassini(R1ca, R2ca, zenith_true) * 3.6e6
|
|
dispListS.append(dispValueS)
|
|
dispListCA.append(dispValueCA)
|
|
ax1 = fig.add_subplot(131)
|
|
ax1.grid(True, alpha=0.3)
|
|
ax1.plot(Trange, np.abs(dispListCA-defaultDispCA), 'r-', label='Cerro Armazones (p=712 mbar, T=6.5C, H=22%)')
|
|
ax1.set_xlabel('Temperature (K)')
|
|
ax1.set_ylabel('Dispersion (mas)')
|
|
ax1.tick_params(axis ='both', direction='in')
|
|
|
|
# Variation in P
|
|
dispListS, dispListCA = [], []
|
|
for p in Prange:
|
|
R1ca = at.n_tph(l=1.49, T=279.65, p=p, RH=0.22, xc=450)
|
|
R2ca = at.n_tph(l=1.78, T=279.65, p=p, RH=0.22, xc=450)
|
|
dispValueS = disp.cassini(R1s, R2s, zenith_true) * 3.6e6
|
|
dispValueCA = disp.cassini(R1ca, R2ca, zenith_true) * 3.6e6
|
|
dispListS.append(dispValueS)
|
|
dispListCA.append(dispValueCA)
|
|
ax2 = fig.add_subplot(132)
|
|
ax2.grid(True, alpha=0.3)
|
|
ax2.plot(Prange/100, np.abs(dispListCA-defaultDispCA), 'r-', label='Cerro Armazones (p=712 mbar, T=6.5C, H=22%)')
|
|
ax2.set_xlabel('Pressure (mbar)')
|
|
ax2.set_ylabel('Dispersion (mas)')
|
|
ax2.tick_params(axis ='both', direction='in')
|
|
|
|
# Variation in T
|
|
dispListS, dispListCA = [], []
|
|
for h in Hrange:
|
|
R1ca = at.n_tph(l=1.49, T=279.65, p=71200, RH=h, xc=450)
|
|
R2ca = at.n_tph(l=1.78, T=279.65, p=71200, RH=h, xc=450)
|
|
dispValueS = disp.cassini(R1s, R2s, zenith_true) * 3.6e6
|
|
dispValueCA = disp.cassini(R1ca, R2ca, zenith_true) * 3.6e6
|
|
dispListS.append(dispValueS)
|
|
dispListCA.append(dispValueCA)
|
|
ax3 = fig.add_subplot(133)
|
|
ax3.grid(True, alpha=0.3)
|
|
ax3.plot(Hrange*100, np.abs(dispListCA-defaultDispCA), 'r-', label='Cerro Armazones (p=712 mbar, T=6.5C, H=22%)')
|
|
ax3.set_xlabel('Humidity (%)')
|
|
ax3.set_ylabel('Dispersion (mas)')
|
|
ax3.tick_params(axis ='both', direction='in')
|
|
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
|
|
# Plotting variations in parameters with corresponding errors.
|
|
if errorSpace:
|
|
Trange = np.linspace(-50, 50, 100)
|
|
Prange = np.linspace(-50, 50, 10)
|
|
Hrange = np.linspace(-0.025, 0.025, 10)
|
|
|
|
fig = plt.figure()
|
|
# Variation in T
|
|
errorRange = [0.1, 0.2, 0.3, 1.0]
|
|
errorRange = errorRange[::-1]
|
|
|
|
ax1 = fig.add_subplot(111)
|
|
ax1.grid(True, alpha=0.3)
|
|
|
|
for e in range(len(errorRange)):
|
|
dispList = []
|
|
errorList = []
|
|
for t in Trange:
|
|
R1n = at.n_tph(l=1.49, T=279.65+t, p=71200, RH=0.22, xc=450)
|
|
R2n = at.n_tph(l=1.78, T=279.65+t, p=71200, RH=0.22, xc=450)
|
|
dispValue = disp.cassini(R1n, R2n, zenith_true) * 3.6e6
|
|
errorValue = disp.cassiniError(zenith_true, l1, l2, T, P, H, xc, dl1=0.0, dl2=0.0, dT=errorRange[e], dP=0, dRH=0, dCO2=0 ) * 3.6e6
|
|
dispList.append(dispValue)
|
|
errorList.append(errorValue)
|
|
dispList = np.asarray(dispList)
|
|
errorList = np.asarray(errorList)
|
|
ax1.fill_between(Trange+279.65, dispList-errorList, dispList+errorList, label='$\delta p$ = %s $^{\circ}$C' %(errorRange[e]))
|
|
|
|
ax1.plot(Trange+279.65, dispList, 'k-')
|
|
ax1.set_xlabel('Temperature (K)')
|
|
ax1.set_ylabel('Dispersion (mas)')
|
|
ax1.set_xlim([np.min(Trange)+279.65,np.max(Trange)+279.65])
|
|
ax1.set_ylim([19,32])
|
|
handles, labels = ax1.get_legend_handles_labels()
|
|
ax1.legend(handles[::-1], labels[::-1], title='Error in measurements')
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
|
|
# Varying the parameters in the uncertainty formula:
|
|
if generate_images:
|
|
################################################
|
|
# Listing some values and the results in a table
|
|
pList = [101325, 101325, 71200, 71200, 71200]
|
|
TList = [288.15, 288.15, 279.65, 279.65, 279.65]
|
|
HList = [0, 0, 0.22, 0.22, 0.6]
|
|
CO2List = [450] * len(pList)
|
|
lList = [1.35] * len(pList)
|
|
dpList = [0, 20, 20, 10, 20]
|
|
dTlist = [0, 0.2, 0.2, 0.1, 0.2]
|
|
dRHlist = [0, 0.02, 0.02, 0.01, 0.02]
|
|
dCO2list = [0, 20, 20, 10, 20]
|
|
dlList = [0, 0.001, 0.001, 0.001, 0.001]
|
|
|
|
refractionList, sigmaList = [], []
|
|
for i in range(len(pList)):
|
|
n = at.n_tph(l=lList[i], T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i])
|
|
s = at.dn_tph(l=lList[i], T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i],
|
|
dl=dlList[i], dT=dTlist[i], dP=dpList[i], dRH=dRHlist[i], dCO2=dCO2list[i])
|
|
refractionList.append((n-1)*1e8)
|
|
sigmaList.append(s*1e8)
|
|
|
|
tableData = []
|
|
for i in range(len(pList)):
|
|
tableData.append([pList[i], TList[i], HList[i], CO2List[i], lList[i], dpList[i], dTlist[i], dRHlist[i], dCO2list[i], dlList[i], refractionList[i], sigmaList[i]])
|
|
print('REFRACTION: ')
|
|
print (tabulate(tableData, headers = ['p', 'T', 'H', 'x_c', 'l', 'dp', 'dT', 'dH', 'dx_c', 'dl', 'n', 's']))
|
|
|
|
#################################################
|
|
print ('\n\n\n')
|
|
|
|
################################################
|
|
# Listing some values and the results in a table
|
|
pList = [101325, 101325, 101325, 101325, 101325]
|
|
TList = [288.15, 288.15, 288.15, 288.15, 288.15]
|
|
HList = [0.0, 0.0, 0.0, 0.0, 0.0]
|
|
CO2List = [450] * len(pList)
|
|
# lList = [1.35] * len(pList)
|
|
dpList = [0, 10, 10, 20, 100]
|
|
dTlist = [0, 0.1, 0.2, 0.2, 1.0]
|
|
dRHlist = [0, 0.01, 0.01, 0.02, 0.05]
|
|
dCO2list = [0, 20, 20, 20, 20]
|
|
# dlList = [0, 0.001, 0.001, 0.001, 0.001] # Uncomment this and comment the next
|
|
# line to check the impact with uncertainty in wavelength
|
|
dlList = [0, 0, 0, 0, 0]
|
|
dzList = [0] * len(pList)
|
|
|
|
dispersionList, errorList = [], []
|
|
for i in range(len(pList)):
|
|
zenith_true = 30 # degrees
|
|
n1 = at.n_tph(l=1.49, T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i])
|
|
n2 = at.n_tph(l=1.78, T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i])
|
|
dn1 = at.dn_tph(l=1.49, T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i],
|
|
dl=dlList[i], dT=dTlist[i], dP=dpList[i], dRH=dRHlist[i], dCO2=dCO2list[i])
|
|
dn2 = at.dn_tph(l=1.78, T=TList[i], p=pList[i], RH=HList[i], xc=CO2List[i],
|
|
dl=dlList[i], dT=dTlist[i], dP=dpList[i], dRH=dRHlist[i], dCO2=dCO2list[i])
|
|
dispValue = disp.cassini(n1, n2, zenith_true) * 3.6e6
|
|
errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, TList[i], pList[i], HList[i], CO2List[i],
|
|
dl1=dlList[i], dl2=dlList[i], dT=dTlist[i], dP=dpList[i], dRH=dRHlist[i], dCO2=dCO2list[i], dz=dzList[i] ) * 3.6e9
|
|
dispersionList.append(dispValue)
|
|
errorList.append(errorDispValue)
|
|
|
|
tableData = []
|
|
for i in range(len(pList)):
|
|
tableData.append([pList[i], TList[i], HList[i], CO2List[i], dpList[i], dTlist[i], dRHlist[i], dCO2list[i], dlList[i], dispersionList[i], errorList[i]])
|
|
print('DISPERSION: ')
|
|
print (tabulate(tableData, headers = ['p', 'T', 'H', 'x_c', 'dp', 'dT', 'dH', 'dx_c', 'dl', 'dispersion', 'error']))
|
|
|
|
# Variation of the uncertainty in P
|
|
pRange = np.linspace(0, 100, 1000)
|
|
errorList = []
|
|
for p in pRange:
|
|
errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450,
|
|
dl1=0.001, dl2=0.001, dT=0.2, dP=p, dRH=0.01, dCO2=20, dz=0 ) * 3.6e9
|
|
errorList.append(errorDispValue)
|
|
|
|
fig = plt.figure(figsize=(4,3))
|
|
ax = fig.add_subplot(111)
|
|
ax.grid(True, alpha=0.3)
|
|
ax.plot(pRange, errorList, 'k-')
|
|
ax.set_xlabel('$\sigma_p$ (Pa)')
|
|
ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
|
|
ax.set_xlim([np.min(pRange), np.max(pRange)])
|
|
ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
|
|
ax.tick_params(axis ='both', direction='in')
|
|
plt.tight_layout()
|
|
|
|
|
|
# Variation of the uncertainty in T
|
|
TRange = np.linspace(0.0, 1, 1000)
|
|
errorList = []
|
|
for T in TRange:
|
|
errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450,
|
|
dl1=0.001, dl2=0.001, dT=T, dP=10, dRH=0.01, dCO2=20, dz=0 ) * 3.6e9
|
|
errorList.append(errorDispValue)
|
|
|
|
fig = plt.figure(figsize=(4,3))
|
|
ax = fig.add_subplot(111)
|
|
ax.grid(True, alpha=0.3)
|
|
ax.plot(TRange, errorList, 'k-')
|
|
ax.set_xlabel('$\sigma_T$ (K)')
|
|
ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
|
|
ax.set_xlim([np.min(TRange), np.max(TRange)])
|
|
ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
|
|
ax.tick_params(axis ='both', direction='in')
|
|
plt.tight_layout()
|
|
|
|
# Variation of the uncertainty in H
|
|
HRange = np.linspace(0.000, 0.05, 1000)
|
|
errorList = []
|
|
for H in HRange:
|
|
errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450,
|
|
dl1=0.001, dl2=0.001, dT=0.2, dP=10, dRH=H, dCO2=20, dz=0 ) * 3.6e9
|
|
errorList.append(errorDispValue)
|
|
|
|
fig = plt.figure(figsize=(4,3))
|
|
ax = fig.add_subplot(111)
|
|
ax.grid(True, alpha=0.3)
|
|
ax.plot(HRange, errorList, 'k-')
|
|
ax.set_xlabel('$\sigma_H$')
|
|
ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
|
|
ax.set_xlim([np.min(HRange), np.max(HRange)])
|
|
ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
|
|
ax.tick_params(axis ='both', direction='in')
|
|
plt.tight_layout()
|
|
|
|
|
|
# Variation of the uncertainty in xc
|
|
CO2Range = np.linspace(0, 200, 1000)
|
|
errorList = []
|
|
for CO2 in CO2Range:
|
|
errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450,
|
|
dl1=0.001, dl2=0.001, dT=0.2, dP=10, dRH=0.01, dCO2=CO2, dz=0 ) * 3.6e9
|
|
errorList.append(errorDispValue)
|
|
|
|
fig = plt.figure(figsize=(4,3))
|
|
ax = fig.add_subplot(111)
|
|
ax.grid(True, alpha=0.3)
|
|
ax.plot(CO2Range, errorList, 'k-')
|
|
ax.set_xlabel('$\sigma_{x_c}$ (ppm)')
|
|
ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
|
|
ax.set_xlim([np.min(CO2Range), np.max(CO2Range)])
|
|
ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
|
|
ax.tick_params(axis ='both', direction='in')
|
|
plt.tight_layout()
|
|
|
|
# Variation of the uncertainty in wavelength
|
|
lRange = np.linspace(0.0, 0.001, 1000)
|
|
errorList = []
|
|
for l in lRange:
|
|
errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450,
|
|
dl1=l, dl2=l, dT=0.2, dP=10, dRH=0.01, dCO2=20, dz=0 ) * 3.6e9
|
|
errorList.append(errorDispValue)
|
|
|
|
fig = plt.figure(figsize=(4,3))
|
|
ax = fig.add_subplot(111)
|
|
ax.grid(True, alpha=0.3)
|
|
ax.plot(lRange*1000, errorList, 'k-')
|
|
ax.set_xlabel('$\sigma_{\lambda}$ (nm)')
|
|
ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
|
|
ax.set_xlim([np.min(lRange)*1000, np.max(lRange)*1000])
|
|
ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
|
|
ax.tick_params(axis ='both', direction='in')
|
|
plt.tight_layout()
|
|
|
|
|
|
# Variation of the uncertainty in zenith
|
|
zRange = np.linspace(0, 60/3600, 1000)
|
|
errorList = []
|
|
for z in zRange:
|
|
errorDispValue = disp.cassiniError( zenith_true, 1.49, 1.78, T=279.65, p=71200, RH=0.22, xc=450,
|
|
dl1=0.001, dl2=0.001, dT=0.2, dP=10, dRH=0.01, dCO2=20, dz=z ) * 3.6e9
|
|
errorList.append(errorDispValue)
|
|
|
|
fig = plt.figure(figsize=(4,3))
|
|
ax = fig.add_subplot(111)
|
|
ax.grid(True, alpha=0.3)
|
|
ax.plot(zRange*3600, errorList, 'k-')
|
|
ax.set_xlim([np.min(zRange)*3600, np.max(zRange)*3600])
|
|
ax.set_ylim([np.min(errorList)-0.05*(np.max(errorList) - np.min(errorList)), np.max(errorList)])
|
|
ax.set_xlabel('$\sigma_z$ (arcseconds)')
|
|
ax.set_ylabel('$\sigma_{\Delta R}$ ($\mu$as)')
|
|
ax.tick_params(axis ='both', direction='in')
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|