diff --git a/ADC.py b/ADC.py deleted file mode 100644 index 0436a7c..0000000 --- a/ADC.py +++ /dev/null @@ -1,178 +0,0 @@ -#!/usr/bin/env Python3 -import numpy as np -from scipy.optimize import fsolve - - -class ADC: - def __init__(self, wedge=(1.6554, 6.28927, -1.48306)): - self.p1 = prism('S-LAH71', 15, wedge[0], wedge[1]) - self.p2 = prism('S-FPL51', 15, wedge[1], wedge[2]) - self.p3 = prism('S-FPL51', 15, wedge[2], wedge[1]) - self.p4 = prism('S-LAH71', 15, wedge[1], wedge[0]) - self.theta = 90 - - def setOrientation(self, theta): - self.p1.setPrismOrientation(theta) - self.p2.setPrismOrientation(theta) - self.p3.setPrismOrientation(180-1*theta) - self.p4.setPrismOrientation(180-1*theta) - self.theta = theta - - def setFieldRotation(self, FR): - self.p1.setPrismOrientation(self.theta-FR) - self.p2.setPrismOrientation(self.theta-FR) - self.p3.setPrismOrientation(180-1*self.theta-FR) - self.p4.setPrismOrientation(180-1*self.theta-FR) - - def calculateDispersion(self, l1, l2, hx=0, hy=0): - ''' - Calculates the dispersion for the ADC using the current position, - for a normalized on sky field angle (hx, hy). Default field position - is in the center of the FoV. - ''' - _prismList = [self.p1, self.p2, self.p3, self.p4] - _d1x_list, _d1y_list = [], [] - _pupilRatio = 0.07943/38.542 - - hx = hx * 7.5e-3 - hy = hy * 7.5e-3 - ix = hx * 1/_pupilRatio - iy = hy * 1/_pupilRatio - - for l in [l1,l2]: - _deltaListX = np.zeros(len(_prismList)) - _deltaListY = np.zeros(len(_prismList)) - for p in range(4): - _n_g = _prismList[p].indexOfRefraction(l) - _b1x = _prismList[p].tilt1x - _b1y = _prismList[p].tilt1y - - _b2x = _prismList[p].tilt2x - _b2y = _prismList[p].tilt2y - - if p==0: - _t_iz_n_x = ix # Zero before first prism - _t_iz_n_y = iy # Zero before first prism - - # X direction - _t_i1_n_x = _t_iz_n_x #+ _b1x - _d1x = deviationCalculation2(_t_i1_n_x, _b1x, _b2x, _n_g) - _t_tz_n_x = _t_iz_n_x - _d1x - _t_iz_n_x = _t_tz_n_x # Set output to new input (for the next prism) - - # Y direction - _t_i1_n_y = _t_iz_n_y #+ _b1y - _d1y = deviationCalculation2(_t_i1_n_y, _b1y, _b2y, _n_g) - _t_tz_n_y = _t_iz_n_y - _d1y - _t_iz_n_y = _t_tz_n_y # Set output to new input (for the next prism) - - # Adding the deviations to the list with delta terms - _deltaListX[p] = _d1x - _deltaListY[p] = _d1y - _d1x_list.append(_deltaListX.sum()) - _d1y_list.append(_deltaListY.sum()) - - self.dispX = ( _d1x_list[1] - _d1x_list[0] ) * 3600 * 1000 * _pupilRatio#* 1.0326# * 0.9108 - self.dispY = ( _d1y_list[1] - _d1y_list[0] ) * 3600 * 1000 * _pupilRatio#* 1.0326# * 0.9108 - return self.dispX, self.dispY - - def optimize(self, inputDispersion, l1, l2, x0=50): - _solver = fsolve(optimizeFunction, x0=x0, args=(self, inputDispersion, l1, l2)) - return _solver[0] - - - - -class prism: - def __init__(self, glass, thickness, tilt1, tilt2=None): - self.glass = glass - self.thickness = thickness - self.tilt1x = tilt1 - self.tilt1y = 0 - self.tilt2x = tilt2 - self.tilt2y = 0 - - # Remember original tilts for reset: - self.tilt1x_OG = tilt1 - self.tilt1y_OG = 0 - self.tilt2x_OG = tilt2 - self.tilt2y_OG = 0 - - - def indexOfRefraction(self, l): - if self.glass == 'S-FPL51': - a = [1.170105050E+00, 4.757107830E-02, 7.638324450E-01] - b = [6.162039240E-03, 2.633728760E-02, 1.418826420E+02] - elif self.glass == 'S-LAH71': - a = [1.982800310E+00, 3.167584500E-01, 2.444726460E+00] - b = [1.189874560E-02, 5.271560010E-02, 2.132206970E+02] - return dispersionFormula(l, a, b) - - def setPrismOrientation(self, theta): - self.tilt1x, self.tilt1y = tiltRotation(self.tilt1x_OG, self.tilt1y_OG, theta) - self.tilt2x, self.tilt2y = tiltRotation(self.tilt2x_OG, self.tilt2y_OG, theta) - - - -def tiltRotation(tx,ty, theta): - theta = np.deg2rad(theta-90) - tx_new = tx * np.cos(theta) - ty * np.sin(theta) - ty_new = tx * np.sin(theta) + ty * np.cos(theta) - return tx_new, ty_new - - -def dispersionFormula(l, a, b): - n = np.sqrt( (a[0]* l**2)/(l**2 - b[0]) + (a[1] * l**2) / (l**2 - b[1]) + (a[2] * l**2) / (l**2 - b[2]) + 1) - return n - - -def deviationCalculation(t_i1, b1, b2, n): - # Returns deviation of the ray after the prism with apex angle alpha - # Units: degrees - alpha = b2 - b1 - delta = t_i1 + np.arcsin( np.sin(np.deg2rad(alpha)) * np.sqrt(n**2-(np.sin(np.deg2rad(t_i1)))**2) - np.sin(np.deg2rad(t_i1)) * np.cos(np.deg2rad(alpha)) )*180/np.pi - alpha - return delta - - -def deviationCalculation2(t_i1, b1, b2, n): - # Returns deviation of the ray after the prism with apex angle alpha - # only now using Snell's law only - # Units: degrees - t_i1 = np.deg2rad(t_i1) - b1 = np.deg2rad(b1) - b2 = np.deg2rad(b2) - - t_i = b1 - t_i1 - t_f = np.arcsin( 1 / n * np.sin( t_i ) ) - t_f_z = b1 - t_f - - t_i = b2 - t_f_z - t_f = np.arcsin( n / 1 * np.sin( t_i ) ) - t_f_z = b2 - t_f - - delta = t_i1 - t_f_z - delta = np.degrees(delta) - return delta - - -def optimizeFunction(theta, adcObject, inputDispersion, l1, l2, FR=None): - if FR: - FRdeg = np.deg2rad(FR) # 30 degrees field rotation - inputDispersionX = np.cos(FRdeg) * inputDispersion - inputDispersionY = np.sin(FRdeg) * inputDispersion - adcObject.setOrientation(theta) - ADCdisp = adcObject.calculateDispersion(l1,l2) - return np.sqrt((ADCdisp[1] - inputDispersionY)**2 + (ADCdisp[0] - inputDispersionX)**2) - else: - adcObject.setOrientation(theta) - ADCdisp = adcObject.calculateDispersion(l1,l2) - ADCdisp = ADCdisp[1] - return ADCdisp - inputDispersion - -def angleFunc(x, a): - ''' - Function to fit the ideal angle on - ''' - return np.arccos(a * np.tan(np.deg2rad(x))) * 180/np.pi - - diff --git a/[Tech note] Atmospheric dispersion models.pdf b/[Tech note] Atmospheric dispersion models.pdf deleted file mode 100644 index 05f55a9..0000000 Binary files a/[Tech note] Atmospheric dispersion models.pdf and /dev/null differ diff --git a/refractionModels.py b/refractionModels.py deleted file mode 100644 index 03722b2..0000000 --- a/refractionModels.py +++ /dev/null @@ -1,224 +0,0 @@ -#!/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)) - -