mirror of
https://gitlab.astro-wise.org/micado/atmosphericmodel.git
synced 2026-10-30 05:01:02 +00:00
242 lines
8.1 KiB
Python
242 lines
8.1 KiB
Python
#!/usr/bin/env Python3
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from matplotlib.lines import Line2D
|
|
from matplotlib import cm
|
|
from matplotlib.colors import LogNorm
|
|
import sys
|
|
from AstroAtmosphere import *
|
|
from scipy import interpolate
|
|
|
|
'''
|
|
@ Author: Joost van den Born
|
|
@ Contact: born@astron.nl
|
|
@ Description: This script uses the US Standard Atmosphere and perturbs it
|
|
at several points. Then full integration is used to calculate
|
|
the atmospheric dispersion. This is then compared to when we
|
|
would have only used the atmospheric conditions directly at
|
|
the observatory. This is done to assess whether we should
|
|
monitor the weather at every height along the line of sight of
|
|
the telescope.
|
|
'''
|
|
|
|
def perturbAtmosphere(h, Temp_perturbed, Press_perturbed):
|
|
for _i in range(len(h)):
|
|
if h[_i] >= 1.0 and h[_i] <= 2.0:
|
|
Tnoise = np.random.normal(0,5,1)
|
|
Temp_perturbed[_i] = Temp_perturbed[_i] + Tnoise
|
|
Pnoise = np.random.normal(0,1000,1)
|
|
Press_perturbed[_i] = Press_perturbed[_i] + Pnoise
|
|
if h[_i] >= 7.0 and h[_i] <= 10.0:
|
|
Tnoise = np.random.normal(0,5,1)
|
|
Temp_perturbed[_i] = Temp_perturbed[_i] + Tnoise
|
|
Pnoise = np.random.normal(0,500,1)
|
|
Press_perturbed[_i] = Press_perturbed[_i] + Pnoise
|
|
if h[_i] >= 25.0 and h[_i] <= 30.0:
|
|
Tnoise = np.random.normal(0,10,1)
|
|
Temp_perturbed[_i] = Temp_perturbed[_i] + Tnoise
|
|
Pnoise = np.random.normal(0,250,1)
|
|
Press_perturbed[_i] = Press_perturbed[_i] + Pnoise
|
|
if h[_i] >= 50.0 and h[_i] <= 60.0:
|
|
Tnoise = np.random.normal(0,15,1)
|
|
Temp_perturbed[_i] = Temp_perturbed[_i] + Tnoise
|
|
Pnoise = np.random.normal(0,4,1)
|
|
Press_perturbed[_i] = Press_perturbed[_i] + Pnoise
|
|
return Temp_perturbed, Press_perturbed
|
|
|
|
plt.rc('font', family='serif')
|
|
plt.rc('axes', labelsize='large')
|
|
|
|
h, _, _, _, Temp, Press, Dens, vSound, _,_ = np.genfromtxt('si1py_modified.txt', unpack=True, skip_header=6)
|
|
|
|
# Conditions at Cerro Armazones
|
|
T = 288.15 # k
|
|
P = 101324 # Pa
|
|
H = 0.0
|
|
xc = 314 # ppm
|
|
lat = -24.5983 # degrees
|
|
h0 = 0 # m
|
|
l1 = 1.49 # micron
|
|
l2 = 1.78 # micron
|
|
zenith = 60 # degrees
|
|
|
|
# Calculate the index of refraction
|
|
at = Observatory()
|
|
n1 = at.n_tph(l=l1, T=T, p=P, RH=H, xc=xc)
|
|
n2 = at.n_tph(l=l2, T=T, p=P, RH=H, xc=xc)
|
|
rho = at.rho(p=P, T=T, RH=H, xc=xc)
|
|
|
|
refr = refraction(lat, h0)
|
|
refr.setReducedHeight(P, rho)
|
|
disp = dispersion(lat, h0)
|
|
disp.setReducedHeight(P, rho)
|
|
|
|
|
|
|
|
# We must create a cube of weather data:
|
|
# We assume a flat atmosphere for now, so that h does not need to be adjusted
|
|
# x_max = 200 km, gives zenith_max = ~66.7 degrees
|
|
# Every 20 km in horizontal direction a new perturbation and interpolation
|
|
# for every 500 m.
|
|
|
|
x_max = 200 # km
|
|
|
|
Temp_cube = []
|
|
Press_cube = []
|
|
Temp_cube_n = []
|
|
Press_cube_n = []
|
|
Temp_cube.append(Temp)
|
|
Press_cube.append(Press)
|
|
Temp_cube_n.append(Temp)
|
|
Press_cube_n.append(Press)
|
|
|
|
|
|
|
|
for x in range(int(x_max/(2))):
|
|
# Adding noise
|
|
# 1-2 km: 2 degree noise, 100 Pa noise
|
|
# 7-10 km: 3 degree noise, 150 Pa noise
|
|
# 25-30 km: 5 degree noise, 250 Pa noise
|
|
Temp_perturbed = Temp.copy()
|
|
Press_perturbed = Press.copy()
|
|
|
|
Temp_perturbed, Press_perturbed = perturbAtmosphere(h, Temp_perturbed, Press_perturbed)
|
|
Temp_cube.append(Temp_perturbed)
|
|
Press_cube.append(Press_perturbed)
|
|
Temp_cube_n.append(Temp)
|
|
Press_cube_n.append(Press)
|
|
|
|
|
|
f_Temp_n = interpolate.interp1d(np.linspace(0,x_max, int(x_max/2)+1), Temp_cube_n, axis=0)
|
|
f_Press_n = interpolate.interp1d(np.linspace(0,x_max, int(x_max/2)+1), Press_cube_n, axis=0)
|
|
f_Temp = interpolate.interp1d(np.linspace(0,x_max, int(x_max/2)+1), Temp_cube, axis=0)
|
|
f_Press = interpolate.interp1d(np.linspace(0,x_max, int(x_max/2)+1), Press_cube, axis=0)
|
|
|
|
Temp_interpolated_n = f_Temp_n(np.linspace(0,x_max, int(x_max/0.5)+1) )
|
|
Press_interpolated_n = f_Press_n(np.linspace(0,x_max, int(x_max/0.5)+1) )
|
|
Temp_interpolated = f_Temp(np.linspace(0,x_max, int(x_max/0.5)+1) )
|
|
Press_interpolated = f_Press(np.linspace(0,x_max, int(x_max/0.5)+1) )
|
|
|
|
x, y = np.mgrid[0:200:(int(x_max/0.5)+1)*1j, 0:np.max(h):len(h)*1j]
|
|
|
|
fig = plt.figure()
|
|
ax1 = fig.add_subplot(211)
|
|
im = ax1.pcolor(x,y, Temp_interpolated, cmap=cm.coolwarm)
|
|
fig.colorbar(im, ax=ax1, label='Temperature (K)')
|
|
ax1.set_aspect(1)
|
|
ax1.set_ylabel('Height (km)')
|
|
|
|
ax2 = fig.add_subplot(212)
|
|
im2 = ax2.pcolor(x,y, Press_interpolated, cmap=cm.coolwarm, norm=LogNorm())
|
|
ax2.set_aspect(1)
|
|
ax2.set_xlabel('Horizontal distance (km)')
|
|
ax2.set_ylabel('Height (km)')
|
|
fig.colorbar(im2, ax=ax2, label='Pressure (Pa)')
|
|
plt.tight_layout()
|
|
plt.savefig('2d_atmosphere.png', bbox_inches='tight')
|
|
plt.savefig('2d_atmosphere.pdf', bbox_inches='tight')
|
|
plt.show(block=False)
|
|
|
|
|
|
# Now for a given zenith angle we need to find the temp and press along the line of sight
|
|
x = np.linspace(0,200,(int(x_max/0.5)+1))
|
|
y = h
|
|
|
|
f_Temp_2d_n = interpolate.interp2d(y,x,Temp_interpolated_n, kind='linear')
|
|
f_Press_2d_n = interpolate.interp2d(y,x,Press_interpolated_n, kind='linear')
|
|
f_Temp_2d = interpolate.interp2d(y,x,Temp_interpolated, kind='linear')
|
|
f_Press_2d = interpolate.interp2d(y,x,Press_interpolated, kind='linear')
|
|
|
|
fig2 = plt.figure()
|
|
ax = fig2.add_subplot(121)
|
|
ax.grid(True, alpha=0.3)
|
|
ax.plot(Temp, h, 'k-')
|
|
|
|
ax2 = fig2.add_subplot(122)
|
|
ax2.grid(True, alpha=0.3)
|
|
ax2.plot(Press, h, 'k-')
|
|
|
|
|
|
zenithRange = np.linspace(30,65,100)
|
|
refrDiff, dispDiff = [], []
|
|
|
|
for zenith in zenithRange:
|
|
z = np.deg2rad(zenith)
|
|
y_los = np.linspace(0, 86, len(h))
|
|
x_los = np.linspace(0, 86*np.tan(z), len(h))
|
|
|
|
Temp_los, Press_los = [], []
|
|
Temp_los_n, Press_los_n = [], []
|
|
for i in range(len(h)):
|
|
Temp_los.append(f_Temp_2d(y_los[i], x_los[i])[0])
|
|
Press_los.append(f_Press_2d(y_los[i], x_los[i])[0])
|
|
Temp_los_n.append(f_Temp_2d_n(y_los[i], x_los[i])[0])
|
|
Press_los_n.append(f_Press_2d_n(y_los[i], x_los[i])[0])
|
|
|
|
Temp_los = np.asarray(Temp_los)
|
|
Press_los = np.asarray(Press_los)
|
|
Temp_los_n = np.asarray(Temp_los_n)
|
|
Press_los_n = np.asarray(Press_los_n)
|
|
rho_los = at.rho(p=Press_los, T=Temp_los, RH=H, xc=xc)
|
|
rho_los_n = at.rho(p=Press_los_n, T=Temp_los_n, RH=H, xc=xc)
|
|
|
|
ax.plot(Temp_los, h, label='z=%i' %zenith)
|
|
ax2.plot(Press_los, h, label='z=%i' %zenith)
|
|
|
|
refrInt = refr.refractionIntegral(n1, z, useStandardAtmosphere=False, heightData=h, rhoData=rho_los_n)*3600*1000
|
|
refrInt_pert = refr.refractionIntegral(n1, z, useStandardAtmosphere=False, heightData=h, rhoData=rho_los)*3600*1000
|
|
refrDiff.append(refrInt_pert - refrInt)
|
|
|
|
dispInt = disp.refractionIntegral(n1, n2, z, useStandardAtmosphere=False, heightData=h, rhoData=rho_los_n)*3.6e6
|
|
dispInt_pert = disp.refractionIntegral(n1, n2, z, useStandardAtmosphere=False, heightData=h, rhoData=rho_los)*3.6e6
|
|
dispDiff.append(dispInt_pert - dispInt)
|
|
|
|
ax.legend()
|
|
ax2.legend()
|
|
|
|
|
|
|
|
# fig3 = plt.figure()
|
|
# ax1 = fig3.add_subplot(211)
|
|
# ax2 = fig3.add_subplot(212)
|
|
# ax1.tick_params(axis ='both', direction='in')
|
|
# ax2.tick_params(axis ='both', direction='in')
|
|
# ax1.grid(True, alpha=0.3)
|
|
# ax2.grid(True, alpha=0.3)
|
|
# ax1.semilogy(zenithRange, np.abs(refrDiff), 'k-')
|
|
# ax2.semilogy(zenithRange, np.abs(dispDiff), 'k-')
|
|
# ax1.set_xlabel('Zenith (degrees)')
|
|
# ax2.set_xlabel('Zenith (degrees)')
|
|
# ax1.set_ylabel('|$R_{p} - R_{SA}$| (mas)')
|
|
# ax2.set_ylabel('|$\Delta R_{p} - \Delta R_{SA}$| (mas)')
|
|
# plt.tight_layout()
|
|
# plt.savefig('differential_USSA_for_different_z.png', bbox_inches='tight')
|
|
# plt.savefig('differential_USSA_for_different_z.pdf', bbox_inches='tight')
|
|
# plt.show()
|
|
|
|
fig3 = plt.figure(figsize=(9,7))
|
|
ax1 = fig3.add_subplot(211)
|
|
ax2 = fig3.add_subplot(212)
|
|
ax1.tick_params(axis ='both', direction='in', which='both', labelsize=16)
|
|
ax2.tick_params(axis ='both', direction='in', which='both', labelsize=16)
|
|
ax1.grid(True, alpha=0.3)
|
|
ax2.grid(True, alpha=0.3)
|
|
|
|
|
|
ax1.plot(zenithRange, np.log10(np.abs(refrDiff)), 'k-', linewidth=2)
|
|
ax2.plot(zenithRange, np.log10(np.abs(dispDiff)), 'k-', linewidth=2)
|
|
ax1.set_xlabel('Zenith (degrees)', fontsize=18)
|
|
ax2.set_xlabel('Zenith (degrees)', fontsize=18)
|
|
ax1.set_ylabel('$\log_{10}($|$R - R_{o}$|$)$ (mas)', fontsize=18)
|
|
ax2.set_ylabel('$\log_{10}($|$\Delta R - \Delta R_{o}$|$)$ (mas)', fontsize=18)
|
|
plt.tight_layout()
|
|
# plt.savefig('differential_USSA_for_different_z.png', bbox_inches='tight')
|
|
# plt.savefig('differential_USSA_for_different_z.pdf', bbox_inches='tight')
|
|
plt.show()
|
|
|
|
|
|
|
|
|
|
# Comparison of different air columns.
|
|
|
|
|