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

First try of curve fit (under develpment)

parent 4c021168
Pipeline #6886 canceled with stages
......@@ -5,21 +5,22 @@ import pandas as pd
import json
import sys
from os.path import basename
from scipy.optimize import curve_fit
from licel_reader import LICEL
from observatory import Observatory
from gdas_reader import GDASReader
from atmospheric_profiles import AtmosphericProfiles
from astropy import units as u
import models
gdas_config = 'config/entry_data_gdas.json'
obs_config = 'config/ObservatoryData.json'
x_min = 0.
x_max = 5000.
y_min = 1e-4
y_max = 1e3
x_min = 0
x_max = 5000
y_min = 10
y_max = 16
is_groundline = False
sig_range = 140000
......@@ -48,8 +49,8 @@ def show_file(filename):
for (channel_id,channel) in enumerate(licel.channels):
#if (channel_id > 0):
# continue
if (channel_id != 0):
continue
channel.offset_range = offsets[nr] # + global_offset # correct for the different zero-offset
if (channel_id > 3):
......@@ -60,7 +61,11 @@ def show_file(filename):
(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)
#plt.plot(x,np.log(y*x**2),label=labels[channel_id])
ids = np.where(y>0)
sr = np.log(y[ids]*x[ids]**2)
plt.plot(x[ids], sr,label=labels[channel_id])
if (channel_id == 0):
x = x - 2 # move x by 2 meters to the left
......@@ -93,10 +98,10 @@ def show_file(filename):
ax.axhline(y=0., color='k')
ax.legend(loc='best')
else:
plt.plot(x,y,label=labels[channel_id])
#plt.plot(x,y,label=labels[channel_id])
plt.xlim(x_min,x_max)
plt.ylim(y_min,y_max)
plt.yscale('log')
#plt.yscale('log')
nr = nr + 1
......@@ -106,16 +111,40 @@ def show_file(filename):
plt.title(basename(filename))
plt.show()
gdas_file = [ '/Users/markusgaug/CTA/LIDAR/Measurements/LaPalma/GDAS/gdas1.aug21.w4_magic1', ]
return x, y, sr, ids
gdas_file = [ '/home/rgrau/Desktop/LIDAR/Analysis/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=True, save_to_file=True, metrics=True)
ap.calculate_results_for_all_profiles(plot=False, save_to_file=False, metrics=False)
for name in sorted(glob.glob(sys.argv[1])):
print(name)
x, y, sr, ids = show_file(name)
x_min //= 1000
x_max //= 1000
arr_y = np.zeros_like(x)
x+=2200
x/=1000
model = models.Models()
#for name in sorted(glob.glob(sys.argv[1])):
# print(name)
# show_file(name)
for idx, i in enumerate(x):
arr_y[idx] = ap.eval_full(list(ap.entry_data.keys())[0], i * u.km)
print(arr_y[ids].shape, sr.shape)
ans = curve_fit(model.quadratic_sin_model, x[ids], sr)
print(ans)
y_fit = model.quadratic_sin_model(x[ids], *ans[0])
print(*ans[0])
plt.plot(x[ids], sr)
plt.plot(x[ids], y_fit)
plt.show()
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