mirror of
https://gitlab.astro-wise.org/micado/atmosphericmodel.git
synced 2026-10-30 05:01:02 +00:00
86 lines
2.3 KiB
Python
86 lines
2.3 KiB
Python
import numpy as np
|
|
from AstroAtmosphere import *
|
|
|
|
'''
|
|
Date: 02 May 2019
|
|
Author: Joost van den Born
|
|
Description: This script compares the uncertainty propagation results to a Monte Carlo method.
|
|
'''
|
|
|
|
|
|
# Parameters at Cerro Armazones
|
|
T = 288.15 # k
|
|
P = 101325 # Pa
|
|
H = 0.0
|
|
xc = 450 # ppm
|
|
lat = -24.5983 # degrees
|
|
h = 3064 # m
|
|
l1 = 1.49 # micron
|
|
l2 = 1.78 # micron
|
|
zenith = 45 # degrees
|
|
|
|
dl1 = 0
|
|
dl2 = 0
|
|
dT = 0
|
|
dP = 2
|
|
dRH = 0
|
|
dCO2= 0
|
|
dz = 1/3600 # degrees
|
|
|
|
# Initializing dispersion model
|
|
at = Observatory()
|
|
|
|
dPlist = [2, 20, 200, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 20, 200]
|
|
dTlist = [0, 0, 0, 0.02, 0.2, 2, 0, 0, 0, 0, 0, 0, 0.02, 0.2, 2 ]
|
|
dRHlist= [0, 0, 0, 0, 0, 0, 0.002, 0.02, 0.2, 0, 0, 0, 0.002, 0.02, 0.2]
|
|
dCO2list=[0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 20, 200, 2, 20, 200]
|
|
|
|
for i in range(len(dPlist)):
|
|
nerr = at.dn_tph(l1, T, P, H, xc, dl=0, dT=dTlist[i], dP=dPlist[i], dRH=dRHlist[i], dCO2=dCO2list[i])
|
|
print(nerr, dPlist[i], dTlist[i], dRHlist[i], dCO2list[i])
|
|
|
|
|
|
|
|
|
|
|
|
|
|
disp = dispersion(lat, h)
|
|
rho_obs = at.rho(P, T, H, xc)
|
|
disp.setReducedHeight(P, rho_obs)
|
|
zenith_true = 10
|
|
|
|
# We do a Monte Carlo error analysis on the dispersion.
|
|
print('MONTE-CARLO UNCERTAINTY ANALYSIS')
|
|
print('-'*35)
|
|
|
|
for j in range(len(dPlist)):
|
|
sigmaP = dPlist[j] # mbar
|
|
sigmaT = dTlist[j] # K
|
|
sigmaH = dRHlist[j]
|
|
sigmaCO2 = dCO2list[j] # ppm
|
|
|
|
N = 10000
|
|
l1sample = np.random.normal(l1, 0, N)
|
|
l2sample = np.random.normal(l2, 0, N)
|
|
Tsample = np.random.normal(T, sigmaT, N)
|
|
Psample = np.random.normal(P, sigmaP, N)
|
|
Hsample = np.random.normal(H, sigmaH, N)
|
|
CO2sample = np.random.normal(xc, sigmaCO2, N)
|
|
dzsample = np.random.normal(zenith_true, dz, N)
|
|
|
|
results = []
|
|
nList = []
|
|
|
|
for i in range(N):
|
|
# Calculate the refraction at l1 and l2
|
|
R1ca2 = at.n_tph(l=l1sample[i], T=Tsample[i], p=Psample[i], RH=Hsample[i], xc=CO2sample[i])
|
|
R2ca2 = at.n_tph(l=l2sample[i], T=Tsample[i], p=Psample[i], RH=Hsample[i], xc=CO2sample[i])
|
|
dispersionValue = disp.cassini(R1ca2, R2ca2, zenith_true)
|
|
results.append(dispersionValue)
|
|
nList.append(R1ca2)
|
|
|
|
cassErrMC = np.std(results)*3.6e9
|
|
nerrMC = np.std(nList)
|
|
print(nerrMC, dPlist[j], dTlist[j], dRHlist[j], dCO2list[j])
|
|
|