Commit 31f4f71a authored by Markus Gaug's avatar Markus Gaug
Browse files

fixes after intens tests of atmospheric profiles

parent a4555703
......@@ -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,81 @@ 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
self.popt_units = {
"c0": u.dimensionless_unscaled,
"c1": 1 / u.km,
"c2": 1 / (u.km ** 2),
"c3": 1 / (u.km ** 3),
"c4": 1 / (u.km ** 4),
"c5": 1 / (u.km ** 5),
"A": u.dimensionless_unscaled,
"omega": 1 / u.km,
"phi": u.km,
}
self._get_models_data(models_data_filename)
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,
}
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 in data_units.keys():
v = np.array(v)
if isinstance(data_units[k], u.Quantity):
v *= data_units[k].unit
self.entry_data[date][k] = v
for model_name, model_data in self.entry_data[date]["models"].items():
self.entry_data[date]["models"][model_name]["function"] = self.models_class.models[model_name]["function"]
for i, constant in enumerate(self.entry_data[date]["models"][model_name]["p0_names"]):
self.entry_data[date]["models"][model_name]["popt"][i] = u.Quantity(
self.entry_data[date]["models"][model_name]["popt"][i],
self.popt_units[constant])
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"""
......@@ -325,18 +393,6 @@ class AtmosphericProfiles():
parameters. By default, this array will contain information about the height
"""
self.popt_units = {
"c0": u.dimensionless_unscaled,
"c1": 1 / u.km,
"c2": 1 / (u.km ** 2),
"c3": 1 / (u.km ** 3),
"c4": 1 / (u.km ** 4),
"c5": 1 / (u.km ** 5),
"A": u.dimensionless_unscaled,
"omega": 1 / u.km,
"phi": u.km,
}
self.entry_data[date]["models"] = {
model_name: {
"function": self.models_class.models[model_name]["function"],
......@@ -582,12 +638,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,10 +806,20 @@ 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
popt0_overlap = self.entry_data[date]["models"]["Linear Approx"]["popt"][0].value
print("popt",popt0_overlap,popt1_overlap)
if chi_overlap > chi2limit:
return -1
......@@ -765,15 +831,22 @@ class AtmosphericProfiles():
return 0
return (
np.exp(chi_overlap) / popt1_overlap
np.exp(popt0_overlap) / popt1_overlap
* (np.exp(popt1_overlap * hi.value) - np.exp(popt1_overlap * lo.value))
)
) * u.km
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
......@@ -787,32 +860,35 @@ class AtmosphericProfiles():
params = np.array([i.value for i in self.entry_data[date]["models"]["Quadratic Sine Approx"]["popt"]])
res, _ = quad(
self.models_class.quadratic_sin_model,
self.models_class.quadratic_sin_exp_model,
max(lo, self.OVERLAP_TRANSITION).value,
min(hi, self.STRATOSPHERE_TRANSITION).value,
args=tuple(params))
ans += res
ans += res * u.km
if hi > self.STRATOSPHERE_TRANSITION:
#TODO:
params = [
self.entry_data[date]["models"]["Quadratic Sine Approx"]["popt"][i].value
for i in range(self.models_class.models["Stratosphere model"]["nargs"])
]
params[0] = self.eval_sin(date, self.STRATOSPHERE_TRANSITION)
params[1] = 0.059055 # Average for La Palma winter
params[2] = 6.63877e-4 # Average for La Palma winter
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
ans += res * u.km
return ans
def eval_overlap(self, date: datetime, h: u.km) -> float:
return (
self.entry_data[date]["models"]["Linear Approx"]["popt"][0].value
......@@ -832,6 +908,10 @@ class AtmosphericProfiles():
]
params[0] = self.eval_sin(date, self.STRATOSPHERE_TRANSITION)
params[1] = 0.059055 # Average for La Palma winter
params[2] = 6.63877e-4 # Average for La Palma winter
print("parameters", params)
res = self.models_class.stratosphere_model(
h.value, *params,
self.STRATOSPHERE_TRANSITION.value)
......@@ -851,6 +931,45 @@ 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):
try:
h = h.to(u.km)
lo = lo.to(u.km)
H0 = H0.to(u.km)
beta_mol_0 = beta_mol_0.to(1/u.km)
except:
raise AtmosphericProfilesError(
TypeError("h, lo, betal_mol_0 and H0 should have units")
)
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_strato(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"""
......@@ -894,4 +1013,4 @@ class AtmosphericProfiles():
print(data.to_markdown())
print("Model info")
print(self.entry_data[date]["models"])
np.set_printoptions(**oldoptions)
np.set_printoptions(**oldoptions)
\ No newline at end of file
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
......@@ -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*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)
......@@ -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()
......@@ -88,7 +88,7 @@ def show_file(filename):
if (bg_range == 0):
ax.set_ylim(-0.1,0.5)
else:
ax.set_ylim(-1.,2.5)
ax.set_ylim(-0.01,0.05)
ax.axhline(y=0., color='k')
ax.legend(loc='best')
else:
......
{
"Linear Approx": [
null,
[
5000,
2000
],
[
-Infinity,
Infinity
]
],
"Quadratic Approx": [
[
0,
-0.95,
-0.17
],
[
Infinity,
3000
],
[
-Infinity,
Infinity
]
],
"Cubic Approx": [
[],
[
Infinity,
3000
],
[
-Infinity,
Infinity
]
],
"Quartic Approx": [
[],
[
Infinity,
3000
],
[
-Infinity,
Infinity
]
],
"Quintic Approx": [
[],
[
Infinity,
3000
],
[
-Infinity,
Infinity
]
],
"Quadratic Sine Approx": [
[
0.2,
0.3,
1.0
],
[
21000,
3000
],
[
[
-Infinity,
-Infinity,
-0.45,
-25,
0.2,
-30
],
[
Infinity,
Infinity,
0.40,
35,
2.2,
15
]
]
]
}
\ No newline at end of file
......@@ -38,7 +38,7 @@ class Models():
}
def linear_model(self, x: np.array, c0: float, c1: float) -> float:
return c0 + c1 * x / 10
return c0 + c1 * x
def quadratic_model(self, x: np.array, c0: float, c1: float, c2: float) -> float:
return c0 + c1 * x / 10 + c2 * x * x / 100
......
{
"magic1": {
"Latitude": [
"28,45,42.462",
"+"
],
"Longitude": [
"17,53,26.525",
"-"
],
"Height": 2199.4,
"ObservatoryName": "Observatorio del Roque de los Muchachos (Magic1)",
"MagneticFieldX": 29.28,
"MagneticFieldY": 3.6,
"MagneticFieldZ": -23.0,
"GeoidOffset": 43.89
},
"WuerzburgCity": {
"Latitude": [
"51,38,48.0",
"+"
],
"Longitude": [
"9,56,36.0",
"+"
],
"Height": 300,
"ObservatoryName": "Wuerzburg City",
"MagneticFieldX": -1,
"MagneticFieldZ": -1,
"GeoidOffset": 45.7
},
"Hess": {
"Latitude": [
"23,16,16.5",
"-"
],
"Longitude": [
"16,30,02.1",
"+"
],
"Height": 1820,
"ObservatoryName": "H.E.S.S. location, Khomas Highland, Namibia",
"MagneticFieldX": -1,
"MagneticFieldZ": -1,
"GeoidOffset": 32.7
},
"Veritas": {
"Latitude": [
"31,40,30.21",
"+"
],
"Longitude": [
"110,57,07.77",
"-"
],
"Height": 1268,
"ObservatoryName": "VERITAS location, Tucson, Arizona",
"MagneticFieldX": -1,
"MagneticFieldZ": -1,
"GeoidOffset": -29.6
},
"AtmoscopeNamibia": {
"Latitude": [
"23,14,44.8",
"-"
],
"Longitude": [
"16,30,59.2",
"+"
],
"Height": 1746,
"ObservatoryName": "Location of Atmoscope in Namibia, close to H.E.S.S. site",
"MagneticFieldX": -1,
"MagneticFieldZ": -1,
"GeoidOffset": 32.7
},
"AtmoscopeTenerife1": {
"Latitude": [
"28,16,45.0",
"+"
],
"Longitude": [
"16,32,0.6",
"-"
],
"Height": 2275,
"ObservatoryName": "Location of Atmoscope in Tenerife, center of CTA candidate site",
"MagneticFieldX": -1,
"MagneticFieldZ": -1,
"GeoidOffset": 50.5
},
"AtmoscopeTenerife2": {
"Latitude": [
"28,16,23.9",
"+"
],
"Longitude": [
"16,32,33.0",
"-"
],
"Height": 2277,
"ObservatoryName": "Location of Atmoscope in Tenerife, edge of CTA candidate site, view to OT",
"MagneticFieldX": -1,
"MagneticFieldZ": -1,
"GeoidOffset": 50.5
}
}
\ No newline at end of file
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