#!/usr/bin/env Python3
import numpy as np
from pdb import set_trace
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import sys
from AstroAtmosphere import *

'''
@ Author:		Joost van den Born 
@ Contact:		born@astron.nl
@ Description:	We plot the total refraction as a function of height.
				We find that most refraction happens in the lower atmosphere
				i.e. <10 km.
'''

def refractionIntegral2(n, R, z0, n0, R0):
	z0 = np.deg2rad(z0)
	return (n0 * R0 * np.sin(z0))/(n*np.sqrt((n*R)**2 - (n0 * R0 * np.sin(z0))**2 ))

def n_for_heights(n0, z0, R0=None, useStandardAtmosphere=True, heightData=None, rhoData=None):
	'''
		Using data that gives temperature, pressure and density as a function of height, 
		the gladstone-dale relation is invoked to calculate the path of a light ray 
		along the atmosphere. This is allows us to also include atmospheric activity in 
		different atmosphere layers.
		Recommended use with the US Standard Atmosphere from 1976.
		USSA1976 data generated from code supplied by http://www.pdas.com/atmos.html
	'''
	if useStandardAtmosphere:
		_h, _, _, _, _Temp, _Press, _Dens,_,_,_ = np.genfromtxt('si1py_modified.txt', unpack=True, skip_header=6)
	else:
		_Dens = rhoData
		_h 	  = heightData
	z0 		= np.deg2rad(z0)
	_Dens 	= np.asarray(_Dens)
	_h 	 	= np.asarray(_h) 
	_R 		= R0 + _h * 1000
	
	# Gladstone-Dale relation between refractive index and density
	_n_h = 1 + (n0 - 1) / _Dens[0] * _Dens
	_z_h = np.arcsin( (n0 * R0) / (_n_h * _R) * np.sin(z0) )

	return _h, _n_h, np.degrees(_z_h)


h, _, _, _, Temp, Press, Dens, vSound, _,_ = np.genfromtxt('si1py_modified.txt', unpack=True, skip_header=6)


# Conditions
T   = 288.15 	# k
P   = 101324  	# Pa
H   = 0.0
xc  = 314   	# ppm
lat = -24.5983 	# degrees
h0   = 0  	# m
l1  = 1.49 		# micron
l2  = 1.78 		# micron
zenith = 60 	# 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)

refr = refraction(lat, h0) 
refr.setReducedHeight(P, rho)
disp = dispersion(lat, h0)
disp.setReducedHeight(P, rho)


fig = plt.figure(figsize=(5,4))
ax = fig.add_subplot(111)
ax.grid(True, alpha=0.3)
ax.set_xlabel('Height (km)')
ax.set_ylabel('Normalized cumulative refraction')
ax.set_xlim([0, h[-1]])

for z in [30]:
	refList = []
	dispList = []
	for i in range(len(h)):
		refList.append(refr.refractionIntegral(n1, z, useStandardAtmosphere=False, heightData=h[:i+1], rhoData=Dens[:i+1])*3600)
		dispList.append(disp.refractionIntegral(n1, n2, z, useStandardAtmosphere=False, heightData=h[:i+1], rhoData=Dens[:i+1])*3.6e6)

	ax.plot(h, np.asarray(refList)/refList[-1], 'k-', label='$\lambda=1.49$ $\mu$m, z=%i$^{\\circ}$' %z)	
ax.set_ylim([0,1.1])

plt.tight_layout()
plt.show()