mirror of
https://gitlab.astro-wise.org/micado/atmosphericmodel.git
synced 2026-10-29 20:51:05 +00:00
Removed some unnessary files
This commit is contained in:
Joost van den Born committed
1 parent
14482ef5f5
commit
51ce3213c9
3 files changed
-402
No files matched your search
@@ -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
|
||||
|
||||
|
||||
Binary file not shown.
@@ -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))
|
||||
|
||||
|
||||
Reference in new issue
Block a user