Commit 1ebfeb66 authored by Markus Gaug's avatar Markus Gaug
Browse files

added lot of doc, fixed mol. fits

parent 3893000c
Pipeline #7015 failed with stages
in 0 seconds
import glob
from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import json
import sys
from os.path import basename
from scipy.optimize import curve_fit
from astropy import units as u
import glob
# LIDAR analysis modules
from licel_reader import LICEL
from rayleigh import Rayleigh
from observatory import Observatory
from gdas_reader import GDASReader
from atmospheric_profiles import AtmosphericProfiles
from astropy import units as u
import models
# required configuration files
gdas_config = 'config/entry_data_gdas.json'
obs_config = 'config/ObservatoryData.json'
x_min = 0
x_max = 10000
# initialize the molecular profile from GDAS
# location of the corresponding GDAS file, needs to be adapted!!!
gdas_file = [ '../../Measurements/LaPalma/GDAS/gdas1.aug21.w4_magic1', ]
obs_name = 'magic1'
obs = Observatory(obs_config, obs_name) # get the observatory coordinates and altitude
gdas = GDASReader(gdas_file, obs) # read the GDAS file
ap = AtmosphericProfiles(gdas_config, gdas) # initialize the atmospheric profiles of the GDAS file
H0 = ap.H0_WS # altitude of the observatory
ap.calculate_results_for_all_profiles(plot=False, save_to_file=False, metrics=False) # calculuate the atmospheric profiles of the GDAS file
dates = list(ap.entry_data.keys())[0] # available dates in GDAS file
model = models.Models() # intialize the fitting models
# initialize the axis ranges
x_min = 0 # in km
x_max = 5. # in km
y_min = 10
y_max = 16
is_groundline = False
sig_range = 140000
# standard channel configuration
labels = [ '355 nm (analog)', '355 nm (photon counting)',
'389 nm (analog)', '389 nm (photon counting)',
'532 nm (analog)', '532 nm (photon counting)',
'near range (analog)', 'near range (photon counting)' ]
rebin = -1
rebin_thres = 0
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)
offsets = [ 2 , -7 , 10, -5 , -7, -9, -3, -5 ] # internal Licel electronic time offsets
bg_range = 5000 # Background subtraction (both atmospheric and electronic)
used_elastic_channels = [0,4] # Used elastic channels for molecular profile fits and Klett inversion
dx_range = 2. # current uncertainty of range in m
elast = [0,4]
# apply a constant re-binning (-1 means no re-binning)
rebin = -1
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)
# apply a dynamic re-binning threshold in photo-electrons
rebin_thres = -1 # 0.1 # 4 # in p.e.
rebin_seq = 50 # increase re-binning every x slices
rebin_fac = 2 # increase re-binning by x after every rebin_seq slices
# function to evaluate the molecular profile at altitude h
def eval_molecular_profile(h: np.array, param: float, hlo: float, costheta: float, beta_mol_0: float):
#print (h,hlo,h,beta_mol_0,H0)
return param - np.log(Rayleigh.ptMAGIC) + ap.eval_f(dates, hlo, h*u.km, costheta, beta_mol_0, H0)
def show_file(filename):
licel = LICEL(filename)
licel.print_header()
nr = 0
plt.figure(figsize=[15,7])
if (is_groundline):
plt.subplot(421)
plt.figure(figsize=[15,7]) # initialize a new figure for each data file
for (channel_id,channel) in enumerate(licel.channels):
licel = LICEL(filename) # initialize Licel reader for filename
licel.print_header() # print header information
costheta = np.cos(licel.zenith_ang * np.pi / 180.) # needed for LIDAR inversion
integrate_mol_from = 0 * u.km # start point of extinction integration
for (channel_id,channel) in enumerate(licel.channels): # loop over all channels
if ((channel_id != 0) and (channel_id != 4)):
continue
channel.offset_range = offsets[nr] # + global_offset # correct for the different zero-offset
if (channel_id > 3):
channel.offset_range = offsets[channel_id] # correct for the different zero-offset
if (channel_id > 3): # work-around, until David has set correct bin widths for the new Licels
channel.bin_width = 3.75
channel.print_header() # print channel header
channel.print_header() # print channel header
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)
# get the data to plot
if (channel_id == 2): # channel 2 has a different ground line offset at about half of the signal range, and need separate background estimation
(x_range,y_sig,y_err,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=rebin_thres,rebin_seq=rebin_seq,rebin_fac=rebin_fac,
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=4,rebin_seq=30,rebin_fac=2,bg_range=bg_range,verbose=True)
(x_range,y_sig,y_err,offset,offseterr) = channel.get_data(rebin=rebin,rebin_thres=rebin_thres,rebin_seq=rebin_seq,rebin_fac=rebin_fac,
bg_range=bg_range,
verbose=True)
if (channel_id == 0):
x_range = x_range - 2 # move x by 2 meters to the left for fine-adjustment
ids = np.where(y>0)
sr = np.log(y[ids]*x[ids]**2)
plt.plot(x[ids], sr,label=labels[channel_id], alpha = 0.5)
ids = np.where((y_sig>0) & (x_range>0)) # discard any negative data for the logarithm in plot
x_range = x_range[ids]
x_height= x_range/1000. * costheta # convert x-axis from range to height for MAGIC analysis
y_sig = y_sig[ids]
y_err = y_err[ids]
sr = np.log(y_sig * x_range**2) # standard logarithm of range-corrected signal
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
y_err = 1/sr * (y_err * x_range**2 + y_sig * 2 * x_range * dx_range )
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))
plt.plot(x_height, sr,label=labels[channel_id], alpha = 0.5)
if channel_id in used_elastic_channels: # molecular fit and Klett only for selected elastic channels
if (channel_id == 0):
x = x - 2 # move x by 2 meters to the left
if (channel.analog == False):
y = y * 0.008
if (is_groundline):
ax = plt.subplot(4,2,channel_id+1)
ax.plot(x,y,label=labels[channel_id])
ax.set_xlim(0.,sig_range)
if (channel_id > 3):
ax.set_xlim(0.,sig_range/2.)
ax.set_xlabel('distance (m)')
if (channel.analog):
ax.set_ylabel('raw signal (mV)')
if (bg_range == 0):
ax.set_ylim(4.,6.)
else:
if (channel_id == 6):
ax.set_ylim(-0.01,0.01)
else:
ax.set_ylim(-0.03,0.03)
else:
ax.set_ylabel('raw signal (MHz)')
if (bg_range == 0):
ax.set_ylim(-0.1,0.5)
else:
ax.set_ylim(-1.,2.5)
ax.axhline(y=0., color='k')
ax.legend(loc='best')
else:
#plt.plot(x,y,label=labels[channel_id])
plt.xlim(x_min,x_max)
plt.ylim(y_min,y_max)
#plt.yscale('log')
ray = Rayleigh(channel.wavelength_nm * u.nm) # initialize Rayleigh backscatter coefficient for used wavelength
beta_mol_0 = ray.dbeta_domega((180 * u.deg).to(u.rad)) # mol. backscatter coefficient at altitude 0 from LIDAR
nr = nr + 1
ids2 = (x_height > 2) & (x_height < 5)
# fit molecular profile to data and obtain the LIDAR constant c0
ans = curve_fit(lambda xl, c0: eval_molecular_profile(xl, c0, integrate_mol_from, costheta, beta_mol_0), x_height[ids2] + H0.to(u.km).value, sr[ids2])
y_fit = eval_molecular_profile(x_height[ids2] + H0.to(u.km).value, *ans[0], integrate_mol_from, costheta, beta_mol_0)
chi2_red = np.sum(((sr[ids2] - y_fit)/y_err[ids2])**2) / (y_fit.size - 1)
print(*ans[0], chi2_red, y_err[ids2])
plt.plot(x_height[ids2], y_fit, label = r"molecular fit {:.0f} nm, $\chi^{{2}}$/NDF={:.1f}".format(channel.wavelength_nm,chi2_red))
plt.text(x_min+(x_max-x_min)*0.45, y_min+(y_max-y_min)*0.92,
u'zenith={:.0f}\u00B0'.format(licel.zenith_ang),
fontsize='x-large')
plt.text(x_min+(x_max-x_min)*0.45, y_min+(y_max-y_min)*0.85,
'{:.0f} shots'.format(licel.laser1['shots']),
fontsize='x-large')
plt.xlim(x_min,x_max)
plt.ylim(y_min,y_max)
#plt.yscale('log')
plt.legend(loc='best')
plt.grid(axis='x')
plt.title(basename(filename))
plt.xlabel('Altitude above LIDAR (km)')
plt.ylabel(r'Logarithm range-corrected signal (log(p.e. $\cdot m^{{2}}) )$')
plt.show()
if __name__ == "__main__":
for name in sorted(glob.glob(sys.argv[1])):
print(name)
show_file(name)
for name in sorted(glob.glob(sys.argv[1])):
print(name)
show_file(name)
x_min //= 1000
x_max //= 1000
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