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

Updates for molecular profile fit

parent a4555703
Pipeline #6894 failed with stages
in 0 seconds
......@@ -16,8 +16,9 @@ import numpy as np
import pandas as pd
from datetime import datetime
import matplotlib.pyplot as plt
from scipy import integrate
from scipy.optimize import curve_fit
from scipy.integrate import quad
from scipy.integrate import quad, quadpack
from scipy.interpolate import CubicSpline
from astropy.units.quantity import Quantity
from astropy import units as u
......@@ -95,7 +96,8 @@ class AtmosphericProfiles():
ABSMAX = 40 * u.km
def __init__(
self, models_data_filename: str, atmospheric_reader: AtmosphericDataReader
self, models_data_filename: str, atmospheric_reader: AtmosphericDataReader,
entry_data_filename: str = ''
) -> None:
"""
Instantiate attributes
......@@ -113,15 +115,68 @@ class AtmosphericProfiles():
atmospheric_reader: AtmosphericDataReader
Datatype to store all data from an atmospheric profile
"""
self.entry_data = {}
# Instantiate Models
self.models_class = models.Models()
self._get_models_data(models_data_filename)
self.atmospheric_reader = atmospheric_reader
self.calculations_done = False
if entry_data_filename:
self._read_entry_data(entry_data_filename)
self.calculations_done = True
else:
self.calculations_done = False
def _read_entry_data(self, filename: str) -> None:
"""Reads JSON file with the entry data"""
data_units = {
"height" : 1 * u.m,
"temperature" : 1 * u.K,
"numerical_density_relative_to_us_standard" : 1,
"log_numerical_density_relative_to_us_standard" : 1,
"drho" : 1,
"log_numerical_density_relative_to_us_standard_err" : 1,
"wind_speed" : 1 * u.m/u.s,
"wind_dir" : 1 * u.rad,
"rh" : 1,
"mass_overburden" : 1 * (u.g / (u.cm * u.cm))**2,
"log_mass_overburden" : 1,
}
self._get_models_data(models_data_filename)
with open(filename) as file:
entry_json = json.load(file)
try:
# Parser keys
# Change all the arrays to numpy arrays
for key, val in entry_json.items():
date = datetime.strptime(key,'%Y/%m/%d/%H')
self.entry_data[date] = val
for k, v in val.items():
if k == "models":
continue
if k in data_units.keys():
v = np.asarray(v)
if isinstance(data_units[k], u.Quantity):
v *= data_units[k].unit
self.entry_data[date][k] = v
##update models
self._append_models_data(date, self.entry_data[date]["height"])
print(self.entry_data)
except Exception as e:
raise AtmosphericProfilesError("JSON format error", e)
def _get_models_data(self, filename: str) -> None:
"""Reads JSON file with the models data"""
......@@ -582,12 +637,12 @@ class AtmosphericProfiles():
unit = x_data.unit
except:
unit = u.dimensionless_unscaled
ax.set_xlabel(plt_xlabel + '[{unit:FITS}]')
ax.set_xlabel(plt_xlabel + rf"[{unit:FITS}]")
try:
unit = y_data.unit
except:
unit = u.dimensionless_unscaled
ax.set_ylabel(plt_ylabel + '[{unit:FITS}]')
ax.set_ylabel(plt_ylabel + rf"[{unit:FITS}]")
if isinstance(x_data, Quantity):
x_data = x_data.to_value()
......@@ -750,7 +805,15 @@ class AtmosphericProfiles():
"""Integrate linear model params"""
# make sure limits are in the same unit
hi = hi.to(lo.unit)
try:
hi= hi.to(u.km)
lo= lo.to(u.km)
except:
raise AtmosphericProfilesError(
TypeError("hmin and hmax should be in meters")
)
chi_overlap = self.entry_data[date]["models"]["Linear Approx"]["chi2"]
popt1_overlap = self.entry_data[date]["models"]["Linear Approx"]["popt"][1].value
......@@ -773,7 +836,14 @@ class AtmosphericProfiles():
def integrate_full(self, date: datetime, lo: u.km, hi: u.km) -> float :
# make sure limits are in the same unit
hi = hi.to(lo.unit)
try:
hi= hi.to(u.km)
lo= lo.to(u.km)
except:
raise AtmosphericProfilesError(
TypeError("hmin and hmax should be in meters")
)
if hi > self.ABSMAX:
return -99
......@@ -804,15 +874,13 @@ class AtmosphericProfiles():
res, _ = quad(
self.models_class.stratosphere_exp_model,
max(lo, self.STRATOSPHERE_TRANSITION).value,
hi,
hi.value,
args=tuple(params + [self.STRATOSPHERE_TRANSITION.value]))
ans += res
return ans
def eval_overlap(self, date: datetime, h: u.km) -> float:
return (
self.entry_data[date]["models"]["Linear Approx"]["popt"][0].value
......@@ -851,6 +919,44 @@ 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: float, H0: u.km):
try:
h = h.to(u.km)
lo = lo.to(u.km)
H0 = H0.to(u.km)
except:
raise AtmosphericProfilesError(
TypeError("h, lo and H0 should be in km")
)
ans = np.zeros(h.shape)
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
* self.integrate_overlap(date, hlo_c, h_c)
)
elif h_c > self.STRATOSPHERE_TRANSITION:
ans[i] = (
self.eval_overlap(date, h_c) - 2 / costheta * Rayleigh.LIDAR_ratio_mol * beta_mol_0
* self.integrate_full(date, hlo_c, h_c)
)
else:
ans[i] = (
self.eval_sin(date, h_c) - 2 / costheta * Rayleigh.LIDAR_ratio_mol * beta_mol_0
* self.integrate_full(date, hlo_c,h_c)
)
return ans
def print_data(self, date: datetime) -> None:
"""Prints formatted data of all the information"""
......
import glob
import numpy as np
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 licel_reader import LICEL
from rayleigh import Rayleigh
from observatory import Observatory
from gdas_reader import GDASReader
......@@ -14,8 +15,8 @@ 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'
gdas_config = 'data_files/entry_data_gdas.json'
obs_config = 'data_files/ObservatoryData.json'
x_min = 0
x_max = 5000
......@@ -113,7 +114,7 @@ def show_file(filename):
return x, y, sr, ids
gdas_file = [ '/home/rgrau/Desktop/LIDAR/Analysis/gdas1.aug21.w4_magic1', ]
gdas_file = [ 'data_files/gdas1.aug21.w4_magic1', ]
obs_name = 'magic1'
obs = Observatory(obs_config, obs_name)
......@@ -122,7 +123,14 @@ 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
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)
......@@ -132,22 +140,17 @@ x_min //= 1000
x_max //= 1000
arr_y = np.zeros_like(x[ids])
x+=2200
#x+=2188.033
x/=1000
model = models.Models()
xs = np.zeros_like(x[ids])
for idx, i in enumerate(x[ids]):
arr_y[idx] = ap.eval_full(list(ap.entry_data.keys())[0], i * u.km)
xs[idx] = i
xx = x[ids]
ids2 = (xx>1) & (xx<5)
print(arr_y.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])
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(xs, arr_y)
plt.show()
plt.plot(x[ids], sr)
plt.plot(x[ids], y_fit)
plt.plot(xx[ids2], y_fit)
plt.xlim(0, 5)
plt.show()
......@@ -121,7 +121,7 @@ class Rayleigh():
self.calculate_sigma()
# calculate beta in units of km^-1 (Tomasi eq. 2)
self.beta = 1e5*N*self.sigma
self.beta = N.to(1/u.km**3) * self.sigma.to(u.km**2)
return self.beta
def calculate_n(self):
......
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