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

new features

parent 1ebfeb66
Pipeline #7077 failed with stages
in 0 seconds
......@@ -9,6 +9,7 @@ import glob
# LIDAR analysis modules
from licel_reader import LICEL
from lidar_analysis import glue_lange
from rayleigh import Rayleigh
from observatory import Observatory
from gdas_reader import GDASReader
......@@ -34,15 +35,26 @@ model = models.Models() # intialize the fitting mode
# initialize the axis ranges
x_min = 0 # in km
x_max = 5. # in km
y_min = 10
y_max = 16
y_min = 11
y_max = 17.5
# molecular fit ranges
x_mol0 = 2.0 # in km
x_mol1 = 5. # in km
max_ranges_visualized = [ 6. , 0.5, # 355 nm in km
6. , 6., # 389 nm in km
6., 0.5, # 532 nm in km
0.8, 0.8 ] # near range in km
# 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)' ]
offsets = [ 2 , -7 , 10, -5 , -7, -9, -3, -5 ] # internal Licel electronic time offsets
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
......@@ -69,6 +81,32 @@ def show_file(filename):
costheta = np.cos(licel.zenith_ang * np.pi / 180.) # needed for LIDAR inversion
integrate_mol_from = 0 * u.km # start point of extinction integration
x_max_used = x_max * costheta
is_glue = True
id_amp = 2
id_pc = 3
testplot = True # plot a test plot to control the gluing
rebin_glue = 0 # do not rebin for gluing
deltaR = 1000 # initial gluing window
deltae = 1.5
fmin = 0.01
fmax = 180.
id0 = 30
licel.channels[id_amp].offset_range = offsets[id_amp] # correct for the different zero-offset in channel_id 1
licel.channels[id_pc].offset_range = offsets[id_pc] # correct for the different zero-offset in channel_id 1
if (is_glue):
(a_Raman, b_Raman) = glue_lange(licel,id_amp,id_pc,rebin_glue,
bg_range=bg_range+(int)(licel.channels[id_amp].num_data/2),bg_range2=(int)(licel.channels[id_amp].num_data/2),
deltaR=deltaR,
deltaa=deltae,fmin=fmin,fmax=fmax,plot=testplot,id0=id0) #2000.,7500.,deadtime,testplot)
if (testplot):
plt.show()
for (channel_id,channel) in enumerate(licel.channels): # loop over all channels
......@@ -90,7 +128,14 @@ def show_file(filename):
if (channel_id == 0):
x_range = x_range - 2 # move x by 2 meters to the left for fine-adjustment
ids = np.where((y_sig>0) & (x_range>0)) # discard any negative data for the logarithm in plot
if (is_glue and channel.analog): # apply the gluing here!!
print ('apply glueing here: ', a_Raman, b_Raman)
y_sig = y_sig*a_Raman + b_Raman
# discard any negative data for the logarithm in plot
# discard range with too high noise
ids = np.where((y_sig>0) & (x_range>0) & (x_range < max_ranges_visualized[channel_id]*1000))
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]
......@@ -106,23 +151,23 @@ def show_file(filename):
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
ids2 = (x_height > 2) & (x_height < 5)
ids2 = (x_height > x_mol0*costheta) & (x_height < x_mol1*costheta)
# 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])
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],p0=[ 15. ])
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,
plt.text(x_min+(x_max_used-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,
plt.text(x_min+(x_max_used-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.xlim(x_min,x_max_used)
plt.ylim(y_min,y_max)
#plt.yscale('log')
......
......@@ -128,7 +128,7 @@ def impose_rebinning(y,b_idx,offset,offseterr,verbose=False):
def glue_lange(lic_reader,channel_id1=0, channel_id2=1, rebin=0, bg_range=500,deltaR=100,deltaa=1., fmin=8., fmax=180.,plot=False, id0=100):
def glue_lange(lic_reader,channel_id1=0, channel_id2=1, rebin=0, bg_range=500, bg_range2=2, deltaR=100,deltaa=1., fmin=8., fmax=180.,plot=False, id0=100):
"""
Apply the method of Lange D., Kumar D., Rocadenbosch F., Sicard M. and Comeron A., Rev. Boliviana de Fisica 20s (2011) 4-6.
......@@ -154,8 +154,8 @@ def glue_lange(lic_reader,channel_id1=0, channel_id2=1, rebin=0, bg_range=500,de
"""
#lic_reader.channels[channel_id2].apply_deadtime_correction(deadtime)
(x0,y0,yerr0,offset,offset_err) = lic_reader.channels[channel_id1].get_data(rebin,bg_range) # Analog
(x1,y1,yerr1,offset,offset_err) = lic_reader.channels[channel_id2].get_data(rebin,bg_range) # PC
(x0,y0,yerr0,offset,offset_err) = lic_reader.channels[channel_id1].get_data(rebin=rebin,bg_range=bg_range,bg_range2=bg_range2) # Analog
(x1,y1,yerr1,offset,offset_err) = lic_reader.channels[channel_id2].get_data(rebin=rebin,bg_range=bg_range) # PC
print ('\n\nenter in glueing here\n')
......@@ -203,7 +203,7 @@ def glue_lange(lic_reader,channel_id1=0, channel_id2=1, rebin=0, bg_range=500,de
e2 = lambda x: (np.square(vpc[ids] - x[0]*va[ids] - x[1])/vpc[ids]).sum() # divide here by vpc for better stability! (not done in Lange et al.)
res = minimize(e2,(0.1,1.0e-4),tol=1e-8,bounds=bnds)
#print ('id0=',id0,'id1=',id1, 'a=',res.x[0], 'b=', res.x[1], 'e2=',res.fun)
print ('id0=',id0,'id1=',id1, 'x[id0]=',x0[id0],'x[id1]=',x0[id1],'a=',res.x[0], 'b=', res.x[1], 'e2=',res.fun, 'flast=', flast, 'fmin=',fmin)
id_mean = id0 + int((id1-id0)/2)
#if (res.x[0] > 0):
......
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