Files
AtmosphericModel/scripts/ADCsteps.py
T

87 lines
2.6 KiB
Python

#!/usr/bin/env Python3
import numpy as np
from AstroAtmosphere import *
import matplotlib.pyplot as plt
'''
@ Author: Joost van den Born
@ Contact: born@astron.nl
@ Description: This script uses the equations of the paper to calculate the
residual dispersion due to the finite step size of the ADC.
'''
def angleFunc(x, a):
'''
Function to fit the ideal angle on
'''
return np.arccos(a * np.tan(np.deg2rad(x))) * 180/np.pi
def angleFuncDeriv(x, a):
return (a * (1+np.tan(np.deg2rad(x))**2)) / np.sqrt(1 - a**2 * np.tan(np.deg2rad(x))**2)
def cassiniDispersionDeriv(n1, n2, zenith, alpha):
a1 = n1 * alpha * np.cos(np.deg2rad(zenith)) / np.sqrt(1 - n1**2 * alpha**2 * np.sin(np.deg2rad(zenith))**2)
a2 = n2 * alpha * np.cos(np.deg2rad(zenith)) / np.sqrt(1 - n2**2 * alpha**2 * np.sin(np.deg2rad(zenith))**2)
return a1 - a2
# Parameters at Cerro Armazones
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
at = dispersion(lat, h)
at.setReducedHeight(P, at.rho(P, T, H, xc))
n1 = at.n_tph(l1, T, P, H, xc)
n2 = at.n_tph(l2, T, P, H, xc)
alpha = at.rc / (at.rc+at.h)
# c_fH = 0.499413 # From Zemax
c_fH = 0.46816061517198737 # Completely analytical cf
# f_fH = 0.46462883308144837 # Best fit to analytical ADC model
dispersion = at.cassini(n1, n2, zenith)*3.6e6
dispersionDeriv = cassiniDispersionDeriv(n1, n2, zenith, alpha)*3.6e6
optAngle = angleFunc(zenith, c_fH)
optAngleDeriv = angleFuncDeriv(zenith, c_fH)
maxDev = optAngleDeriv / dispersionDeriv
print('Zenith = %s degrees' %zenith)
print('Dispersion in mas: ', dispersion)
print('Dispersion/degree zenith at z=%s: %s mas/degree z' %(zenith, dispersionDeriv))
print('Optimum angle ADC: %s' %(optAngle))
print('Optimum angle ADC/degree zenith at z=%s: %s' %(zenith, optAngleDeriv))
print('Max deviation per mas dispersion: %s degrees' %(maxDev))
print('At designed stepsize at most %s mas of residual dispersion' %(1/maxDev*(360/1262.5)))
# Plotting as a function of zenith angle
dataList = []
for zenith in np.linspace(5,65,100):
dispersion = at.cassini(n1, n2, zenith)*3.6e6
dispersionDeriv = cassiniDispersionDeriv(n1, n2, zenith, alpha)*3.6e6
optAngle = angleFunc(zenith, c_fH)
optAngleDeriv = angleFuncDeriv(zenith, c_fH)
maxDev = optAngleDeriv / dispersionDeriv
dataList.append(1/maxDev*(360/1262.5))
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(np.linspace(5,65,100), dataList)
ax.grid(True, alpha=0.2)
ax.set_xlabel('Zenith (degrees)')
ax.set_ylabel('Residual dispersion (mas)')
plt.show()