#!/usr/bin/env Python import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import host_subplot from atmosphericModel import * from pdb import set_trace import pandas as pd ''' @ Author: Joost van den Born @ Contact: born@astron.nl @ Description: Modeling the atmospheric refraction using multiple models of the atmosphere. Models included: - Plane-parallel atmosphere - Cassini's model - Tan5 model (A*tan(z) + B*tan^3(z) + C*tan^5(z)) - Oriani's theorem (A*tan(z) + B*tan^3(z)) - Erf model (from Corbard et al., 2019) - I do not yet consider full integration of the refraction integral. ''' if __name__ == '__main__': plt.rc('font', family='serif') # plt.rc('xtick', labelsize='x-small') # plt.rc('ytick', labelsize='x-small') plt.rc('axes', labelsize='large') # Conditions at Cerro Armazones # T = 279.65 # k # P = 71200 # Pa # H = 0.22 # xc = 450 # ppm T = 288.15 # k P = 101325 # Pa H = 0.0 xc = 314 # ppm lat = -24.5983 # degrees h = 3064 # m l1 = 1.49 # micron l2 = 1.78 # micron zenith = 30 # 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) refractionAnalysis = True dispersionAnalysis = True tableCompare = True if dispersionAnalysis: # Initiate the refraction class ref = dispersion(lat, h) ref.setReducedHeight(P, rho) # Compare the refraction models to check if the output is okay print ('Plane - parallel atmosphere dispersion: %s mas' %(ref.planeParallel(n1, n2, zenith) * 3.6e6)) print ('Cassini dispersion: %s mas' %(ref.cassini(n1, n2, zenith) * 3.6e6)) print ('Tan5 dispersion: %s mas' %(ref.tan5(n1, n2, zenith) * 3.6e6)) print ('Oriani dispersion: %s mas' %(ref.oriani(n1, n2, zenith) * 3.6e6)) print ('Corbard dispersion: %s mas' %(ref.corbard(n1, n2, zenith)*3.6e6)) print ('Full integration: %s mas' %(ref.refractionIntegral(n1, n2, zenith)*3.6e6)) # Now we calculate the dispersion over a wide range of zenith angles # for ideal gas atmosphere zenithRange = np.linspace(5, 85, 1000) ppList = ref.planeParallel(n1, n2, zenithRange) * 3.6e6 cassiniList = ref.cassini(n1, n2, zenithRange) * 3.6e6 tan5List = ref.tan5(n1, n2, zenithRange) * 3.6e6 orianiList = ref.oriani(n1, n2, zenithRange) * 3.6e6 corbardList = [ref.corbard(n1, n2, z) * 3.6e6 for z in zenithRange] fullIntegrationList = [ref.refractionIntegral(n1, n2, z) * 3.6e6 for z in zenithRange] # for isothermal atmosphere cassiniList_isotherm = ref.cassini(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6 tan5List_isotherm = ref.tan5(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6 orianiList_isotherm = ref.oriani(n1, n2, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e6 corbardList_isotherm = [ref.corbard(n1, n2, z, reduced_height=ref.H_isotherm(T)) * 3.6e6 for z in zenithRange] # check when the variation between the different models pass different values dispersionWidth = [] widthList = [] for z in range(len(zenithRange)): compList = [cassiniList[z], tan5List[z], orianiList[z], corbardList[z], cassiniList_isotherm[z], tan5List_isotherm[z], orianiList_isotherm[z], corbardList_isotherm[z]]#, fullIntegrationList[z]] width = np.max(compList) - np.min(compList) widthList.append(width) dispersionWidth.append((zenithRange[z], width)) masBorder = np.argmax(np.asarray(widthList)>1) p1masBorder = np.argmax(np.asarray(widthList)>0.1) p01masBorder = np.argmax(np.asarray(widthList)>0.01) p001masBorder = np.argmax(np.asarray(widthList)>0.001) # Plotting it all, with vertical lines where the variation between the models exceed 1, 0.1, 0.01 and 0.001 mas. fig = plt.figure(figsize=(9,5)) ax = fig.add_subplot(111) ax.set_xlim([30,85]) ax.set_xlabel('Zenith (degrees)') ax.set_ylabel('Dispersion (mas)') ax.grid(True, alpha=0.3) ax.axvline(dispersionWidth[masBorder][0], linestyle='--', color='k', linewidth=0.5) ax.text(dispersionWidth[masBorder][0]-0.5, 402, '1 mas', ha='center', va='bottom', rotation='vertical') ax.axvline(dispersionWidth[p1masBorder][0], linestyle='--', color='k', linewidth=0.5) ax.text(dispersionWidth[p1masBorder][0]-0.5, 402, '0.1 mas', ha='center', va='bottom', rotation='vertical') ax.axvline(dispersionWidth[p01masBorder][0], linestyle='--', color='k', linewidth=0.5) ax.text(dispersionWidth[p01masBorder][0]-0.5, 402, '10 $\mu$as', ha='center', va='bottom', rotation='vertical') ax.axvline(dispersionWidth[p001masBorder][0], linestyle='--', color='k', linewidth=0.5) ax.text(dispersionWidth[p001masBorder][0]-0.5, 402, '1 $\mu$as', ha='center', va='bottom', rotation='vertical') ax.plot(zenithRange, ppList, 'k-', label='Plane-parallel', alpha=0.2) ax.plot(zenithRange, cassiniList, 'g-', label='Cassini') ax.plot(zenithRange, tan5List, 'b-', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$') ax.plot(zenithRange, orianiList, 'y-', label='$Atan(z) - Btan^3(z)$') ax.plot(zenithRange, corbardList, 'r-', label='Error function') ax.plot(zenithRange, fullIntegrationList, 'c-', label='Full integration') ax.plot(zenithRange, cassiniList_isotherm, 'g--', label='Cassini (isotherm)') ax.plot(zenithRange, tan5List_isotherm, 'b--', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$ (isotherm)') ax.plot(zenithRange, orianiList_isotherm, 'y--', label='$Atan(z) - Btan^3(z)$ (isotherm)') ax.plot(zenithRange, corbardList_isotherm, 'r--', label='Error function (isotherm)') ax.tick_params(axis ='both', direction='in') ax.legend() plt.tight_layout() plt.show() if refractionAnalysis: # Initiate the refraction class ref = refraction(lat, h, P, rho) # Compare the refraction models to check if the output is okay print ('Plane - parallel atmosphere refraction: %s arcsec' %(ref.planeParallel(n1, zenith) * 3.6e3)) print ('Cassini refraction: %s arcsec' %(ref.cassini(n1, zenith) * 3.6e3)) print ('Tan5 refraction: %s arcsec' %(ref.tan5(n1, zenith) * 3.6e3)) print ('Oriani refraction: %s arcsec' %(ref.oriani(n1, zenith) * 3.6e3)) print ('Corbard refraction: %s arcsec' %(ref.corbard(n1, zenith)*3.6e3)) print ('Full integration: %s mas' %(ref.refractionIntegral(n1, zenith)*3.6e3)) # Now we calculate the dispersion over a wide range of zenith angles # for ideal gas atmosphere zenithRange = np.linspace(0, 88, 1000) ppList = ref.planeParallel(n1, zenithRange) * 3.6e3 cassiniList = ref.cassini(n1, zenithRange) * 3.6e3 tan5List = ref.tan5(n1, zenithRange) * 3.6e3 orianiList = ref.oriani(n1, zenithRange) * 3.6e3 corbardList = [ref.corbard(n1, z) * 3.6e3 for z in zenithRange] fullIntegrationList = [ref.refractionIntegral(n1, z) * 3.6e3 for z in zenithRange] # for isothermal atmosphere cassiniList_isotherm = ref.cassini(n1, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e3 tan5List_isotherm = ref.tan5(n1, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e3 orianiList_isotherm = ref.oriani(n1, zenithRange, reduced_height=ref.H_isotherm(T)) * 3.6e3 corbardList_isotherm = [ref.corbard(n1, z, reduced_height=ref.H_isotherm(T)) * 3.6e3 for z in zenithRange] # Plotting it all, with vertical lines where the variation between the models exceed 1, 0.1, 0.01 and 0.001 arcsec. fig = plt.figure(figsize=(9,5)) ax = fig.add_subplot(111) ax.set_xlim([80,88]) ax.set_xlabel('Zenith (degrees)') ax.set_ylabel('Refraction (arcsec)') ax.grid(True, alpha=0.3) ax.plot(zenithRange, ppList, 'k-', label='Plane-parallel', alpha=0.2) ax.plot(zenithRange, cassiniList, 'g-', label='Cassini') ax.plot(zenithRange, tan5List, 'b-', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$') ax.plot(zenithRange, orianiList, 'y-', label='$Atan(z) - Btan^3(z)$') ax.plot(zenithRange, corbardList, 'r-', label='Error function') ax.plot(zenithRange, fullIntegrationList, 'c-', label='Full integration') ax.plot(zenithRange, cassiniList_isotherm, 'g--', label='Cassini (isotherm)') ax.plot(zenithRange, tan5List_isotherm, 'b--', label='$Atan(z) - Btan^3(z) + Ctan^5(z)$ (isotherm)') ax.plot(zenithRange, orianiList_isotherm, 'y--', label='$Atan(z) - Btan^3(z)$ (isotherm)') ax.plot(zenithRange, corbardList_isotherm, 'r--', label='Error function (isotherm)') ax.tick_params(axis ='both', direction='in') ax.legend() plt.tight_layout() plt.show() if tableCompare: refTableData = [] dispTableData = [] ref = refraction(lat, h, P, rho) disp = dispersion(lat, h) disp.setReducedHeight(P, rho) zenithRange = [10, 30, 50, 60, 70, 80, 85] headerList = np.array(['z', 'Plane-parallel', 'Cassini', 'tan3', 'tan5', 'Error function', 'Full integration']) refTableData.append(headerList) dispTableData.append(headerList) for z in zenithRange: ppValue = round(ref.planeParallel(n1, z) * 3.6e3, 3) cassiniValue = round(ref.cassini(n1, z) * 3.6e3, 3) tan3Value = round(ref.oriani(n1, z) * 3.6e3, 3) tan5Value = round(ref.tan5(n1, z) * 3.6e3, 3) erfValue = round(ref.corbard(n1, z) * 3.6e3, 3) fullIntValue = round(ref.refractionIntegral(n1, z) * 3.6e3, 3) refTableData.append([z, ppValue, cassiniValue, tan3Value, tan5Value, erfValue, fullIntValue]) ppValue = round(disp.planeParallel(n1, n2, z) * 3.6e6, 3) cassiniValue = round(disp.cassini(n1, n2, z) * 3.6e6, 3) tan3Value = round(disp.oriani(n1, n2, z) * 3.6e6, 3) tan5Value = round(disp.tan5(n1, n2, z) * 3.6e6, 3) erfValue = round(disp.corbard(n1, n2, z) * 3.6e6, 3) fullIntValue = round(disp.refractionIntegral(n1, n2, z) * 3.6e6, 3) dispTableData.append([z, ppValue, cassiniValue, tan3Value, tan5Value, erfValue, fullIntValue]) dfr = pd.DataFrame(np.asarray(refTableData)) dfd = pd.DataFrame(np.asarray(dispTableData)) print(dfr.to_string(index=False)) print(dfd.to_string(index=False))