#!/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.