#!/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()