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

changes on the linear model

parent 3d1ad752
Pipeline #6971 failed with stages
in 0 seconds
......@@ -87,6 +87,7 @@ class AtmosphericProfiles():
# (m a.s.l): GPS: 2231.917 m - 43.883 m (EGM2008)
H0_WS = 2188.033 * u.m
H_ns = 10300 # (m^-1) approximate scale height of particle density at ORM
H0_LIDAR = 2188.24 * u.m
#Temperature_mean = 8.89 # mean temperature at the MAGIC site
......@@ -161,8 +162,6 @@ class AtmosphericProfiles():
"log_mass_overburden" : 1,
}
with open(filename) as file:
entry_json = json.load(file)
......@@ -631,7 +630,6 @@ class AtmosphericProfiles():
One PNG Matplotlib plot saved to local directory
"""
fig, ax = plt.subplots()
ax.set_title(plt_title)
......@@ -655,6 +653,9 @@ class AtmosphericProfiles():
ax.scatter(x_data, y_data, label="Real Data")
#plot vertical line
ax.axvline(x=self.H0_LIDAR, color="black", linestyle="--")
for model_name in models:
#Drop astropy units for all the elements in list
params = np.array([e.to_value() if isinstance(e, Quantity) else e for e in self.entry_data[date]["models"][model_name]["popt"]])
......@@ -673,8 +674,15 @@ class AtmosphericProfiles():
try:
plt.savefig("images/" + self.atmospheric_reader.name + "/" + date.strftime("%Y_%m_%d_%H") + plt_title.replace(" ", "").replace("-", "_") + '.png')
except FileNotFoundError:
os.mkdir(os.listdir(os.path.join(os.getcwd(), "images", self.atmospheric_reader.name)))
plt.savefig("images/" + self.atmospheric_reader.name + "/" + date.strftime("%Y_%m_%d_%H") + plt_title.replace(" ", "").replace("-", "_") + '.png')
try:
os.mkdir(os.path.join(os.getcwd(), "images"))
plt.savefig("images/" + self.atmospheric_reader.name + "/" + date.strftime("%Y_%m_%d_%H") + plt_title.replace(" ", "").replace("-", "_") + '.png')
except Exception as e:
try:
os.mkdir(os.path.join(os.getcwd(), "images", self.atmospheric_reader.name))
plt.savefig("images/" + self.atmospheric_reader.name + "/" + date.strftime("%Y_%m_%d_%H") + plt_title.replace(" ", "").replace("-", "_") + '.png')
except Exception as e:
AtmosphericProfilesError(e)
plt.close(fig)
......@@ -819,9 +827,10 @@ class AtmosphericProfiles():
chi_overlap = self.entry_data[date]["models"]["Linear Approx"]["chi2"]
popt1_overlap = self.entry_data[date]["models"]["Linear Approx"]["popt"][1].value
popt1_overlap = self.entry_data[date]["models"]["Linear Approx"]["popt"][1].value /10
popt0_overlap = self.entry_data[date]["models"]["Linear Approx"]["popt"][0].value
print("popt",popt0_overlap,popt1_overlap)
#print("popt",popt0_overlap,popt1_overlap)
if chi_overlap > chi2limit:
return -1
......@@ -894,7 +903,7 @@ class AtmosphericProfiles():
def eval_overlap(self, date: datetime, h: u.km) -> float:
return (
self.entry_data[date]["models"]["Linear Approx"]["popt"][0].value
+ self.entry_data[date]["models"]["Linear Approx"]["popt"][1].value * h.value
+ self.entry_data[date]["models"]["Linear Approx"]["popt"][1].value/10 * h.value
)
def eval_sin(self, date: datetime, h: u.km) -> float:
......@@ -951,6 +960,7 @@ class AtmosphericProfiles():
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
......@@ -1015,3 +1025,20 @@ class AtmosphericProfiles():
print("Model info")
print(self.entry_data[date]["models"])
np.set_printoptions(**oldoptions)
if __name__ == "__main__":
from gdas_reader import GDASReader
OBS = "data_files/ObservatoryData.json"
OBSNAME = "magic1"
obs = Observatory(OBS, OBSNAME)
# Instantiate Gdas
FILENAME = [ "data_files/gdas1.feb19.w1_magic1"]
gdas = GDASReader(FILENAME, obs)
# Instantiate AtmosphericProfiles
ap = AtmosphericProfiles("data_files/entry_data_gdas.json", gdas)
ap.calculate_results_for_all_profiles(plot=True, save_to_file=False, metrics=False)
date = list(ap.entry_data.keys())[0]
......@@ -114,7 +114,7 @@ def show_file(filename):
return x, y, sr, ids
gdas_file = [ 'data_files/gdas1.aug21.w4_magic1', ]
gdas_file = [ '../gdas1.aug21.w4_magic1', ]
obs_name = 'magic1'
obs = Observatory(obs_config, obs_name)
......
......@@ -38,7 +38,7 @@ class Models():
}
def linear_model(self, x: np.array, c0: float, c1: float) -> float:
return c0 + c1 * x
return c0 + c1 * x / 10
def quadratic_model(self, x: np.array, c0: float, c1: float, c2: float) -> float:
return c0 + c1 * x / 10 + c2 * x * x / 100
......
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