Commit 6415e528 authored by Roger.Grau.Haro's avatar Roger.Grau.Haro
Browse files

klett updates

parent 7c3c1466
Pipeline #7949 failed with stages
in 0 seconds
......@@ -85,8 +85,9 @@ class AtcaWizard():
ValueError(
"Lidar ratio must be bigger than 0; cannot initialize LIDAR... "
))
self.slope_guess = 1/AtmosphericProfiles.H_ns + self.beta_mol_0 * Rayleigh.ptMAGIC * 16/3 * np.pi / self.coszd
def exchange_point_for_klett(self, height: np.ndarray, sr: np.ndarray, href: float, param: float) -> np.ndarray:
"""Exchange point at height href for the Klett alpha"""
......
......@@ -86,9 +86,9 @@ class AtmosphericProfiles():
# Data measurement from private company for the MAGIC Collaboration 28/07/15 11:25
# (m a.s.l): GPS: 2231.917 m - 43.883 m (EGM2008)
H0_WS = 2188.033 * u.m
H_ns = 10300 # (m^-1) approximate scale height of particle density at ORM
H_ns = 10300 * u.m # (m^-1) approximate scale height of particle density at ORM
H0_LIDAR = 2188.24 * u.m
#Temperature_mean = 8.89 # mean temperature at the MAGIC site
Temperature_mean = 8.89 * u.deg_C # mean temperature at the MAGIC site
OVERLAP_TRANSITION = 4.4 * u.km
......@@ -837,17 +837,15 @@ class AtmosphericProfiles():
chi_overlap = self.entry_data[date]["models"]["Linear Approx"]["chi2"]
popt1_overlap = self.entry_data[date]["models"]["Linear Approx"]["popt"][1].value /10
popt0_overlap = self.entry_data[date]["models"]["Linear Approx"]["popt"][0].value
#print("popt",popt0_overlap,popt1_overlap)
if chi_overlap > chi2limit:
return -1
return -1 * u.km
if chi_overlap < 0:
return -99
return -99 * u.km
if popt1_overlap == 0:
return 0
return 0 * u.km
return (
np.exp(popt0_overlap) / popt1_overlap
......@@ -868,9 +866,9 @@ class AtmosphericProfiles():
)
if hi > self.ABSMAX:
return -99
return -99 * u.km
ans = 0
ans = 0 * u.km
if lo < self.OVERLAP_TRANSITION:
ans += self.integrate_overlap(date, lo, min(hi, self.OVERLAP_TRANSITION))
......@@ -887,7 +885,6 @@ class AtmosphericProfiles():
ans += res * u.km
if hi > self.STRATOSPHERE_TRANSITION:
#TODO:
params = [
self.entry_data[date]["models"]["Quadratic Sine Approx"]["popt"][i].value
for i in range(self.models_class.models["Stratosphere model"]["nargs"])
......@@ -932,7 +929,6 @@ class AtmosphericProfiles():
params[1] = 0.059055 # Average for La Palma winter
params[2] = 6.63877e-4 # Average for La Palma winter
print("parameters", params)
res = self.models_class.stratosphere_model(
h.value, *params,
self.STRATOSPHERE_TRANSITION.value)
......@@ -952,7 +948,6 @@ class AtmosphericProfiles():
return self.eval_strato(date, h)
##NOT tested
def eval_f(self, date: datetime, lo: u.km, h: np.ndarray, costheta: float, beta_mol_0: 1/u.km, H0: u.km):
try:
h = h.to(u.km)
......
......@@ -395,8 +395,12 @@ def show_file(filename):
alpha_fitted_green = alpha_fitted
aod_err_green = aod_err
aod_green = aod_fitted
offset = 0
offset = 0
Klett1 = lidar_klett(costheta, x_height*u.m, sr, b_widths*u.m, 33, x_mol0*u.m, x_mol1*u.m, ans1[0])
Klett2 = lidar_klett(costheta, x_height*u.m, sr, b_widths*u.m, 33, x_mol2*u.m, x_mol3*u.m, ans2[0])
alpha_err = np.sqrt(aod_err**2/(x_mol2 - x_mol1)**2 + aod_fitted**2/(x_mol2 - x_mol1)**2 /12)
print ('AOD fitted: ', aod_fitted, ' ans2: ', ans2[0], ' ans1: ', ans1[0], 'alpha_fitted:', alpha_fitted, ' alpha_err:', alpha_err)
......
......@@ -133,10 +133,11 @@ def impose_rebinning(y,b_idx,offset,offseterr,verbose=False):
return (np.resize(yy,cnt),np.resize(yyerr,cnt))
def lidar_klett(coszd: float, height: np.ndarray, sr: np.ndarray, bin_width: np.ndarray, lidar_ratio: float, href_l: u.Qunatity,
href_h: u.Quantity, klett_param: float, atmospheric_profiles: AtmosphericProfiles = None):
def lidar_klett(coszd: float, height: np.ndarray, sr: np.ndarray, bin_width: np.ndarray, lidar_ratio: float, href_l: u.Quantity,
href_h: u.Quantity, klett_param: float, date: datetime = None, filename: str = None, atmospheric_profiles: AtmosphericProfiles = None):
atca = AtcaWizard(coszd, atmospheric_profiles)
#atca.get_weather_from_file(date, filename)
sr_exchanged = atca.exchange_point_for_klett(height, sr, href=href_l, param=klett_param)
......
import numpy as np
from astropy import units as u
class Models():
"""Class containing all the approximations models
......@@ -38,7 +39,7 @@ class Models():
}
def linear_model(self, x: np.array, c0: float, c1: float) -> float:
return c0 + c1 * x / 10
return c0 + c1 * x/10
def quadratic_model(self, x: np.array, c0: float, c1: float, c2: float) -> float:
return c0 + c1 * x / 10 + c2 * x * x / 100
......@@ -125,4 +126,11 @@ class Models():
stratosphere_transition_const: float
) -> float:
return np.exp(self.stratosphere_model(x, c0, c1, c2, stratosphere_transition_const))
\ No newline at end of file
return np.exp(self.stratosphere_model(x, c0, c1, c2, stratosphere_transition_const))
def exp_model(
self, x: np.ndarray, c0: float, c1: float
) -> float:
return c0 - c1*x
......@@ -152,7 +152,7 @@ class Rayleigh():
Mw = 0.018015 * (u.kg / u.mol) # molar mass of water vapor [kg/mol]
Xw = humidity.MolarFractionWaterVapor(self.p, self.T, self.RH) # molar fraction of water vapor in moist air
Za = humidity.Compressibility(self.pressure_us_standard, self.temperature_us_standard, 0) # compressibility of dry air
Zw = humidity.Compressibility(13.33 * u.hPa, 293.15 * u.K, 1) # compressibility of pure water vapor
Zw = humidity.Compressibility(13.33 * u.hPa, 293.15 * u.K, 1) # compressibility of pure water vapor
Zm = humidity.Compressibility(self.p, self.T, Xw) # compressibility of moist air
# density of dry air at standard p and T
......@@ -490,4 +490,4 @@ class Rayleigh():
axes[1].ticklabel_format(axis='y', style='sci', useOffset=False, scilimits=(0, 0))
axes[1].minorticks_on()
fig.subplots_adjust(left=0.07, right=0.96, bottom=0.09, top=0.93)
fig.subplots_adjust(left=0.07, right=0.96, bottom=0.09, top=0.93)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment