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