import numpy as np
from AstroAtmosphere import *

'''
	Date: 02 May 2019
	Author: Joost van den Born
	Description: This script compares the uncertainty propagation results to a Monte Carlo method.
'''


# Parameters at Cerro Armazones
T   = 288.15 	# k
P   = 101325  	# Pa
H   = 0.0
xc  = 450   	# ppm
lat = -24.5983 	# degrees
h   = 3064  	# m
l1  = 1.49 		# micron
l2  = 1.78 		# micron
zenith = 45 	# degrees

dl1 = 0
dl2 = 0
dT  = 0
dP  = 2
dRH = 0
dCO2= 0
dz = 1/3600 # degrees

# Initializing dispersion model
at  = Observatory()

dPlist = [2, 20, 200, 0,    0,   0, 0,     0,    0,   0,  0, 0,   2,     20,   200]
dTlist = [0, 0,  0,   0.02, 0.2, 2, 0,     0,    0,   0,  0, 0,   0.02,  0.2,  2  ]
dRHlist= [0, 0,  0,   0,    0,   0, 0.002, 0.02, 0.2, 0,  0, 0,   0.002, 0.02, 0.2]
dCO2list=[0, 0,  0,   0,    0,   0, 0,     0,    0,   2, 20, 200, 2,     20,   200]

for i in range(len(dPlist)):
	nerr = at.dn_tph(l1, T, P, H, xc, dl=0, dT=dTlist[i], dP=dPlist[i], dRH=dRHlist[i], dCO2=dCO2list[i])
	print(nerr, dPlist[i], dTlist[i], dRHlist[i], dCO2list[i])






disp = dispersion(lat, h)
rho_obs = at.rho(P, T, H, xc)
disp.setReducedHeight(P, rho_obs)
zenith_true = 10

# We do a Monte Carlo error analysis on the dispersion.
print('MONTE-CARLO UNCERTAINTY ANALYSIS')
print('-'*35)

for j in range(len(dPlist)):
    sigmaP 	 = dPlist[j] # mbar
    sigmaT 	 = dTlist[j] # K
    sigmaH 	 = dRHlist[j]
    sigmaCO2 = dCO2list[j] # ppm

    N = 10000
    l1sample = np.random.normal(l1, 0, N)
    l2sample = np.random.normal(l2, 0, N)
    Tsample   = np.random.normal(T,   sigmaT,   N)
    Psample   = np.random.normal(P,   sigmaP,   N)
    Hsample   = np.random.normal(H,   sigmaH,   N)
    CO2sample = np.random.normal(xc, sigmaCO2, N)
    dzsample  = np.random.normal(zenith_true, dz, N)
    
    results = []
    nList = []

    for i in range(N):
        # Calculate the refraction at l1 and l2
        R1ca2 = at.n_tph(l=l1sample[i], T=Tsample[i], p=Psample[i], RH=Hsample[i], xc=CO2sample[i])
        R2ca2 = at.n_tph(l=l2sample[i], T=Tsample[i], p=Psample[i], RH=Hsample[i], xc=CO2sample[i])
        dispersionValue = disp.cassini(R1ca2, R2ca2, zenith_true)
        results.append(dispersionValue)
        nList.append(R1ca2)
    
    cassErrMC = np.std(results)*3.6e9
    nerrMC    = np.std(nList)
    print(nerrMC, dPlist[j], dTlist[j], dRHlist[j], dCO2list[j])