Commit f0abd979 authored by Markus Gaug's avatar Markus Gaug
Browse files

Merge branch 'master' of gitlab.cta-observatory.org:cta-array-elements/ccf/LIDAR_Analysis

because it's necessary
parents 7cf1a6ec af024050
Pipeline #7012 failed with stages
in 0 seconds
......@@ -87,6 +87,7 @@ class AtmosphericProfiles():
# (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
H0_LIDAR = 2188.24 * u.m
#Temperature_mean = 8.89 # mean temperature at the MAGIC site
......@@ -161,8 +162,6 @@ class AtmosphericProfiles():
"log_mass_overburden" : 1,
}
with open(filename) as file:
entry_json = json.load(file)
......@@ -251,9 +250,13 @@ class AtmosphericProfiles():
try:
plt.savefig("images/" + self.atmospheric_reader.name + "/" + "metrics.png")
except FileNotFoundError:
os.mkdir(os.listdir(os.path.join(os.getcwd(), "images", self.atmospheric_reader.name)))
plt.savefig("images/" + self.atmospheric_reader.name + "/" + "metrics.png")
try:
d = os.path.join(os.getcwd(), "images", self.atmospheric_reader.name)
if not os.path.exists(d):
os.makedirs(d)
plt.savefig("images/" + self.atmospheric_reader.name + "/" + "metrics.png")
except Exception as e:
AtmosphericProfilesError(e)
def _plot_metrics(self) -> None:
"""Save plot to images/atmospheric_reader.name/chi2*.png of the chi2 histogram"""
......@@ -276,8 +279,13 @@ class AtmosphericProfiles():
try:
plt.savefig("images/" + self.atmospheric_reader.name + "/" + ("chi2_hist_of_" + filename + ".png").replace("/", ""))
except FileNotFoundError:
os.mkdir(os.listdir(os.path.join(os.getcwd(), "images", self.atmospheric_reader.name)))
plt.savefig("images/" + self.atmospheric_reader.name + "/" + ("chi2_hist_of_" + self.atmospheric_reader.filename + ".png").replace("/", ""))
try:
d = os.path.join(os.getcwd(), "images", self.atmospheric_reader.name)
if not os.path.exists(d):
os.makedirs(d)
plt.savefig("images/" + self.atmospheric_reader.name + "/" + ("chi2_hist_of_" + filename + ".png").replace("/", ""))
except Exception as e:
AtmosphericProfilesError(e)
plt.close()
......@@ -631,7 +639,6 @@ class AtmosphericProfiles():
One PNG Matplotlib plot saved to local directory
"""
fig, ax = plt.subplots()
ax.set_title(plt_title)
......@@ -655,6 +662,9 @@ class AtmosphericProfiles():
ax.scatter(x_data, y_data, label="Real Data")
#plot vertical line
ax.axvline(x=self.H0_LIDAR.value, color="black", linestyle="--")
for model_name in models:
#Drop astropy units for all the elements in list
params = np.array([e.to_value() if isinstance(e, Quantity) else e for e in self.entry_data[date]["models"][model_name]["popt"]])
......@@ -673,9 +683,15 @@ class AtmosphericProfiles():
try:
plt.savefig("images/" + self.atmospheric_reader.name + "/" + date.strftime("%Y_%m_%d_%H") + plt_title.replace(" ", "").replace("-", "_") + '.png')
except FileNotFoundError:
os.mkdir(os.listdir(os.path.join(os.getcwd(), "images", self.atmospheric_reader.name)))
plt.savefig("images/" + self.atmospheric_reader.name + "/" + date.strftime("%Y_%m_%d_%H") + plt_title.replace(" ", "").replace("-", "_") + '.png')
try:
d = os.path.join(os.getcwd(), "images", self.atmospheric_reader.name)
if not os.path.exists(d):
os.makedirs(d)
plt.savefig("images/" + self.atmospheric_reader.name + "/" + date.strftime("%Y_%m_%d_%H") + plt_title.replace(" ", "").replace("-", "_") + '.png')
except Exception as e:
AtmosphericProfilesError(e)
plt.close(fig)
def get_block_elements_at_plevel(
......@@ -819,9 +835,10 @@ class AtmosphericProfiles():
chi_overlap = self.entry_data[date]["models"]["Linear Approx"]["chi2"]
popt1_overlap = self.entry_data[date]["models"]["Linear Approx"]["popt"][1].value
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)
#print("popt",popt0_overlap,popt1_overlap)
if chi_overlap > chi2limit:
return -1
......@@ -894,7 +911,7 @@ class AtmosphericProfiles():
def eval_overlap(self, date: datetime, h: u.km) -> float:
return (
self.entry_data[date]["models"]["Linear Approx"]["popt"][0].value
+ self.entry_data[date]["models"]["Linear Approx"]["popt"][1].value * h.value
+ self.entry_data[date]["models"]["Linear Approx"]["popt"][1].value/10 * h.value
)
def eval_sin(self, date: datetime, h: u.km) -> float:
......@@ -934,7 +951,7 @@ class AtmosphericProfiles():
return self.eval_strato(date, h)
##NOT tested
def eval_f(self, date: datetime, lo: u.km, h: u.km, costheta: float, beta_mol_0: 1/u.km, H0: u.km):
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)
lo = lo.to(u.km)
......@@ -951,6 +968,7 @@ class AtmosphericProfiles():
for i, h_i in enumerate(h):
hlo_c = lo + H0
h_c = h_i + H0
if h_c < self.OVERLAP_TRANSITION:
ans[i] = (
self.eval_overlap(date, h_c) - 2 / costheta * Rayleigh.LIDAR_ratio_mol * beta_mol_0
......
......@@ -19,7 +19,7 @@ gdas_config = 'config/entry_data_gdas.json'
obs_config = 'config/ObservatoryData.json'
x_min = 0
x_max = 5000
x_max = 10000
y_min = 10
y_max = 16
......@@ -37,7 +37,25 @@ offsets = [ 2 , -7 , 10, -5 , -7, -9, -3, -5 ]
bg_range = 5000 # 5000 # Background subtraction (both atmospheric and electronic)
if (is_groundline):
bg_range = 5000
gdas_file = [ '../gdas1.aug21.w4_magic1', ]
obs_name = 'magic1'
obs = Observatory(obs_config, obs_name)
gdas = GDASReader(gdas_file, obs)
ap = AtmosphericProfiles(gdas_config, gdas)
ap.calculate_results_for_all_profiles(plot=False, save_to_file=False, metrics=False)
elast = [0,4]
H0 = ap.H0_WS
dates = list(ap.entry_data.keys())[0]
model = models.Models()
def eval_GDAS(x: u.km, param: float, hlo: float, costheta: float, beta_mol_0: float):
return param - np.log(Rayleigh.ptMAGIC) + ap.eval_f(dates, hlo, x*u.km, costheta, beta_mol_0, H0)
def show_file(filename):
licel = LICEL(filename)
......@@ -50,7 +68,7 @@ def show_file(filename):
for (channel_id,channel) in enumerate(licel.channels):
if (channel_id != 0):
if ((channel_id != 0) and (channel_id != 4)):
continue
channel.offset_range = offsets[nr] # + global_offset # correct for the different zero-offset
......@@ -61,12 +79,29 @@ def show_file(filename):
if (channel_id == 2):
(x,y,yerr,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=0,rebin_seq=50,rebin_fac=0,bg_range=bg_range+(int)(channel.num_data/2),bg_range2=(int)(channel.num_data/2),verbose=False)
else:
(x,y,yerr,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=0,rebin_seq=50,rebin_fac=0,bg_range=bg_range,verbose=False)
(x,y,yerr,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=4,rebin_seq=30,rebin_fac=2,bg_range=bg_range,verbose=True)
ids = np.where(y>0)
sr = np.log(y[ids]*x[ids]**2)
plt.plot(x[ids], sr,label=labels[channel_id])
plt.plot(x[ids], sr,label=labels[channel_id], alpha = 0.5)
if channel_id in elast:
hlo = 0 * u.km
costheta = 1
ray = Rayleigh(channel.wavelength_nm * u.nm)
beta_mol_0 = ray.dbeta_domega((180 * u.deg).to(u.rad)) #8.242519245726348e-06*1/u.m
arr_y = np.zeros_like(x[ids])
xx = x[ids]/1000
ids2 = (xx>2) & (xx<8)
print(arr_y.shape, sr.shape)
ans = curve_fit(lambda x, c0: eval_GDAS(x, c0, hlo, costheta, beta_mol_0), xx[ids2] + H0.to(u.km).value, sr[ids2])
y_fit = eval_GDAS(xx[ids2] + H0.to(u.km).value, *ans[0], hlo, costheta, beta_mol_0)
print(*ans[0])
plt.plot(xx[ids2]*1000, y_fit, label = "molecular fit {:.0f} nm".format(channel.wavelength_nm))
if (channel_id == 0):
x = x - 2 # move x by 2 meters to the left
......@@ -112,45 +147,11 @@ def show_file(filename):
plt.title(basename(filename))
plt.show()
return x, y, sr, ids
gdas_file = [ 'data_files/gdas1.aug21.w4_magic1', ]
obs_name = 'magic1'
obs = Observatory(obs_config, obs_name)
gdas = GDASReader(gdas_file, obs)
ap = AtmosphericProfiles(gdas_config, gdas)
ap.calculate_results_for_all_profiles(plot=False, save_to_file=False, metrics=False)
H0 = ap.H0_WS
dates = list(ap.entry_data.keys())[0]
hlo = 0 * u.km
costheta = 1
beta_mol_0 = 8.242519245726348e-06*1/u.m
def eval_GDAS(x: u.km, param: float):
return param - np.log(Rayleigh.ptMAGIC) + ap.eval_f(dates, hlo, x*u.km, costheta, beta_mol_0, H0)
for name in sorted(glob.glob(sys.argv[1])):
print(name)
x, y, sr, ids = show_file(name)
show_file(name)
x_min //= 1000
x_max //= 1000
arr_y = np.zeros_like(x[ids])
#x+=2188.033
x/=1000
model = models.Models()
xx = x[ids]
ids2 = (xx>1) & (xx<5)
print(arr_y.shape, sr.shape)
ans = curve_fit(eval_GDAS, xx[ids2]+2.188, sr[ids2])
y_fit = eval_GDAS(xx[ids2]+2.188, *ans[0])
print(*ans[0])
plt.plot(x[ids], sr)
plt.plot(xx[ids2], y_fit)
plt.xlim(0, 5)
plt.show()
......@@ -454,7 +454,7 @@ class LICEL:
print ('new offset=',offset,' +- ',offset_err)
# optional additional general rebinning
rebin_start = 0 # bin number
rebin_start = 10 # bin number
#rebin_thres = 4 # threshold in p.e.
#rebin_seq = 10 # increase re-binning every 50 bins
#rebin_fac = 2 # multiply binsize by 3 eveny rebin_seq
......
......@@ -38,7 +38,7 @@ class Models():
}
def linear_model(self, x: np.array, c0: float, c1: float) -> float:
return c0 + c1 * x
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
......
......@@ -28,7 +28,6 @@ Adapted from MRayleigh, written by Markus Gaug <markus.gaug@uab.cat>, 04/2013
"""
import matplotlib.pyplot as plt
from math import pi, cos, log
import numpy as np
from astropy import units as u
from astropy.units import cds
......@@ -49,7 +48,7 @@ class Rayleigh():
#TODO: Add to conf file
ptMAGIC = p_mean / pressure_us_standard * temperature_us_standard / T_mean # mean air temperature at MAGIC site [K]
logptMAGIC = log(ptMAGIC) # mean air temperature at MAGIC site [K]
logptMAGIC = np.log(ptMAGIC) # mean air temperature at MAGIC site [K]
def __init__(self, wavelength : u.nm, p : u.hPa = pressure_us_standard,
T : u.K = temperature_us_standard, RH : float = relative_humidity_us_standard,
......@@ -102,7 +101,7 @@ class Rayleigh():
N = self.numerical_density_us_standard * self.p/self.pressure_us_standard * self.temperature_us_standard/self.T
N2 = N*N
pi3 = pow(pi, 3)
pi3 = pow(np.pi, 3)
F = self.king()
self.sigma = 24*pi3*(n2 - 1)**2 / (l4*N2*(n2 + 2)**2) * F # units of cm^-2
return self.sigma
......@@ -206,7 +205,7 @@ class Rayleigh():
F = (c1*F1 + c2*F2 + c3*F3 + c4*F4 + c5*F5)/(c1 + c2 + c3 + c4 + c5) # Tomasi eq. 22
return F
def phase_function(self, angle):
def phase_function(self, angle : u.rad):
"""
Calculates the Chandrasekhar phase function.
......@@ -216,12 +215,18 @@ class Rayleigh():
Returns:
Chandrasekhar phase function for scattering of natural light
"""
try:
angle.to(u.rad)
except Exception as e:
raise Exception("Input angle should be in radians or equivalent", e)
rho = self.depolarization()
# need to solve Chandrasekhar eq. 254 for gamma as a function of rho
f1 = (2 + 2*rho)/(2 + rho)
f2 = (1 - rho)/(1 + rho)
return 0.75*f1*(1 + f2*pow(cos(angle), 2)) # Chandrasekhar eq. 255
return 0.75*f1*(1 + f2*pow(np.cos(angle), 2)) # Chandrasekhar eq. 255
def depolarization(self):
"""
......@@ -234,7 +239,7 @@ class Rayleigh():
return 6*(F - 1)/(3 + 7*F) # Tomasi eq. 5, solved for rho
# TODO: where does this come from? Needs unit test
def dbeta_domega(self, angle):
def dbeta_domega(self, angle : u.rad):
"""
Get the back-scattering coefficient for a given scattering angle.
......@@ -244,13 +249,19 @@ class Rayleigh():
Returns:
back-scattering coefficient in units of km^-1
"""
try:
angle.to(u.rad)
except Exception as e:
raise Exception("Input angle should be in radians or equivalent", e)
try: # check that beta has been calculated
self.beta
except NameError:
self.calculate_beta()
phase_func = self.phase_function(angle)
return phase_func*self.beta/(4*pi)
return phase_func*self.beta/(4*np.pi)
def print_params(self):
"""Prints Rayleigh scattering parameters."""
......@@ -319,7 +330,7 @@ class Rayleigh():
for angle in angles:
f = self.phase_function(np.deg2rad(angle))
y_sigma.append(f)
y_beta.append(f*self.beta/(4*pi))
y_beta.append(f*self.beta/(4*np.pi))
l_str = r'Wavelength $\lambda$={:.1f} nm'.format(self.wavelength)
p_str = r'Pressure p={:.2f} hPa'.format(self.p)
......
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