Skip to content

GitLab

Commit b296cd78 authored by Scott Griffiths's avatar Scott Griffiths
Browse files

Incremental progress on gdas_entry and msis_entry. Added preliminary msis_reader.

parent 11ecd72b
Expand all Show whitespace changes
Inline Side-by-side
atca_wizard.py
View file @ b296cd78
...@@ -47,6 +47,7 @@ class AtcaWizard: ...@@ -47,6 +47,7 @@ class AtcaWizard:
self.fMSIS = None self.fMSIS = None
self.fGDAS = None self.fGDAS = None
# TODO: This function needs revision
def LoadWeather(self, filename="WS2013.root"): def LoadWeather(self, filename="WS2013.root"):
if self.fWeather: if self.fWeather:
delete self.fWeather delete self.fWeather
...@@ -87,6 +88,7 @@ class AtcaWizard: ...@@ -87,6 +88,7 @@ class AtcaWizard:
print(args[0], " ", args[1], " ", args[2]) print(args[0], " ", args[1], " ", args[2])
return (args) return (args)
# TODO: This function needs revision
def LoadMSIS(self, filename="msis2013.root", verbose=False): def LoadMSIS(self, filename="msis2013.root", verbose=False):
if self.fMSIS: if self.fMSIS:
delete self.fMSIS delete self.fMSIS
...@@ -121,6 +123,7 @@ class AtcaWizard: ...@@ -121,6 +123,7 @@ class AtcaWizard:
if verbose: if verbose:
self.fMSIS.Print() self.fMSIS.Print()
# TODO: This function needs revision
def CreateGDAS(self, files="gdas1.*13.w*_magic1", outfile="GDAS1.root"): def CreateGDAS(self, files="gdas1.*13.w*_magic1", outfile="GDAS1.root"):
tlist = [] tlist = []
MGDASFileRead r(files) MGDASFileRead r(files)
...@@ -132,6 +135,7 @@ class AtcaWizard: ...@@ -132,6 +135,7 @@ class AtcaWizard:
tlist.append(w) tlist.append(w)
w.Close() w.Close()
# TODO: This function needs revision
def LoadGDAS(self, filename="gdas2013_magic1.root", verbose=False): def LoadGDAS(self, filename="gdas2013_magic1.root", verbose=False):
if self.fGDAS: if self.fGDAS:
delete self.fGDAS delete self.fGDAS
...@@ -170,6 +174,7 @@ class AtcaWizard: ...@@ -170,6 +174,7 @@ class AtcaWizard:
b = self.fGDAS.GetBranch("MGDASEntry.") b = self.fGDAS.GetBranch("MGDASEntry.")
b.SetAutoDelete(False) b.SetAutoDelete(False)
# TODO: This function needs revision
def GetWeather(self, mjd, verbose=False): def GetWeather(self, mjd, verbose=False):
if not self.fWeather: if not self.fWeather:
raise RuntimeError("No weather file 'WS*.root' has been loaded, please use LoadWeather first.") raise RuntimeError("No weather file 'WS*.root' has been loaded, please use LoadWeather first.")
...@@ -224,6 +229,7 @@ class AtcaWizard: ...@@ -224,6 +229,7 @@ class AtcaWizard:
dust = args[10] dust = args[10]
return temperature, pressure, humidity, dust return temperature, pressure, humidity, dust
# TODO: This function needs revision
def GetMSIS(self, mjd, verbose=False): def GetMSIS(self, mjd, verbose=False):
if not self.fMSIS: if not self.fMSIS:
raise RuntimeError("No MSIS file 'MSIS*.root' has been loaded, please use LoadMSIS first.") raise RuntimeError("No MSIS file 'MSIS*.root' has been loaded, please use LoadMSIS first.")
...@@ -268,6 +274,7 @@ class AtcaWizard: ...@@ -268,6 +274,7 @@ class AtcaWizard:
break break
return msis_entry return msis_entry
# TODO: This function needs revision
def GetGDAS(self, mjd, verbose=False): def GetGDAS(self, mjd, verbose=False):
if not self.fGDAS: if not self.fGDAS:
raise RuntimeError("No GDAS file 'GDAS*.root' has been loaded, please use LoadGDAS first.") raise RuntimeError("No GDAS file 'GDAS*.root' has been loaded, please use LoadGDAS first.")
......
gdas_entry.py
View file @ b296cd78
This diff is collapsed.
msis_entry.py
View file @ b296cd78
...@@ -5,6 +5,7 @@ Copyright: MAGIC Software Development, 2000-2013 ...@@ -5,6 +5,7 @@ Copyright: MAGIC Software Development, 2000-2013
from __future__ import division, print_function from __future__ import division, print_function
from math import pi, sqrt, atan2, log, exp from math import pi, sqrt, atan2, log, exp
from collections import namedtuple
from datetime import datetime from datetime import datetime
from astropy.time import Time from astropy.time import Time
import numpy as np import numpy as np
...@@ -17,16 +18,6 @@ class MSISEntry: ...@@ -17,16 +18,6 @@ class MSISEntry:
gkLoschmidt = 2.079153E25 # Na divided by molar mass of air: 0.0289644 gkLoschmidt = 2.079153E25 # Na divided by molar mass of air: 0.0289644
def __init__(self): def __init__(self):
self.fLatitude(-99)
self.fLongitude(-99)
self.fDebug(False)
self.fHeight = []
self.fDensity = []
self.fLogDensity = []
self.fMassOverburden = []
self.fLogMassOverburden = []
self.fTemperature = []
self.clear() self.clear()
def fill_date(self, line): def fill_date(self, line):
...@@ -35,12 +26,12 @@ class MSISEntry: ...@@ -35,12 +26,12 @@ class MSISEntry:
where YYYY is the year, MM is the month, DD is the day, and HH is the hour. where YYYY is the year, MM is the month, DD is the day, and HH is the hour.
""" """
dt = datetime.strptime(line, "%Y %m %d %H") dt = datetime.strptime(line, "%Y %m %d %H")
apt = Time(dt) # use astropy for mjd conversion astropy_time = Time(dt) # use astropy for mjd conversion
self.fMJD = apt.mjd self.mjd = astropy_time.mjd
self.fHour = dt.hour self.hour = dt.hour
self.fDay = dt.day self.day = dt.day
self.fMonth = dt.month self.month = dt.month
self.fYear = dt.year self.year = dt.year
def fill_observatory(self, line): def fill_observatory(self, line):
n = sscanf(line.Data(), "latitude=%f, longitude=%f, %n", lat, lon, len) n = sscanf(line.Data(), "latitude=%f, longitude=%f, %n", lat, lon, len)
...@@ -49,8 +40,8 @@ class MSISEntry: ...@@ -49,8 +40,8 @@ class MSISEntry:
print(lat, "-", lon, "-") print(lat, "-", lon, "-")
print(" n is:", n, " instead of 2\n") print(" n is:", n, " instead of 2\n")
return False return False
self.fLatitude = lat self.latitude = lat
self.fLongitude = lon self.longitude = lon
def fill_profile(self, line, obs, ver=0): def fill_profile(self, line, obs, ver=0):
n = sscanf(line.Data(), "%f %g %f %n", height, rho, T, len) n = sscanf(line.Data(), "%f %g %f %n", height, rho, T, len)
...@@ -62,188 +53,180 @@ class MSISEntry: ...@@ -62,188 +53,180 @@ class MSISEntry:
rho *= self.gkLoschmidt*1e3 # multiply by Na, and divide by molecular weight of dry air, transform from g/cm^3 to kg/m^3 rho *= self.gkLoschmidt*1e3 # multiply by Na, and divide by molecular weight of dry air, transform from g/cm^3 to kg/m^3
rho /= self.gkNs # transform in units of n_s rho /= self.gkNs # transform in units of n_s
height = self.get_altitude_from_geopotential_height(height)
size = len(self.fHeight) self.height.append(height)
self.density.append(rho)
self.log_density.append(log(rho) if rho > 0 else -1)
self.temperature.append(T)
self.fHeight = np.zeros(size + 1) def get_altitude_from_geopotential_height(self, geop_height):
self.fDensity = np.zeros(size + 1) """
self.fLogDensity = np.zeros(size + 1) Convert from geopotential height to geometric altitude.
self.fTemperature = np.zeros(size + 1)
height = obs.GetAltitudeFromGeopotentialHeight(height*1000)/1000 # use internally km, for precision of fit geop_height : height in km
"""
# TODO: FIXME: Using magic1_latitude for testing.
magic1_latitude = 5.01988021533642415e-01 # obtained from MARS
self.fHeight[size] = height cos2lat = cos(2*magic1_latitude)
self.fDensity[size] = rho Z = (1 + 0.002644*cos2lat)*geop_height + (1 + 0.0089*cos2lat)*geop_height**2/6245
self.fLogDensity[size] = log(rho) if rho > 0 else -1 return Z
self.fTemperature[size] = T
def clear(self): def clear(self):
self.fHeight = [] self.height = []
self.fDensity = [] self.density = []
self.fLogDensity = [] self.log_density = []
self.fMassOverburden = [] self.mass_overburden = []
self.fLogMassOverburden = [] self.log_mass_overburden = []
self.fTemperature = [] self.temperature = []
self.fMJD = -1 self.params_total = {}
self.fHour = -1 self.params_0_to_4_km = {}
self.fDay = -1 self.params_4_to_10_km = {}
self.fMonth = -1 self.params_10_to_40_km = {}
self.fYear = -1 self.params_40_to_100_km = {}
self.fP0Tot = -99 def fit(self, lower=None, upper=None, param_guess=None):
self.fP1Tot = -99 def quadratic(x, c0, c1, c2):
self.fP2Tot = -99 return c0 + c1*x/10 + c2*x*x/100
self.fChi2Tot = -99 xdata, ydata = self.select_in_xrange(lower, upper)
popt, pcov = curve_fit(quadratic, xdata, ydata, p0=param_guess)
self.fP00to4km = -99 chi2 = np.sum((quadratic(xdata, *popt) - ydata)**2)
self.fP10to4km = -99 ndf = len(xdata) - len(popt)
self.fP20to4km = -99 popt[1] /= 10
self.fChi20to4km = -99 popt[2] /= 100
return tuple(popt) + (chi2/ndf,)
self.fP04to10km = -99
self.fP14to10km = -99 def select_in_xrange(lower, upper):
self.fP24to10km = -99 """
self.fChi24to10km = -99 Select a subset of xdata, ydata for a given range of x values.
self.fP010to40km = -99 Args:
self.fP110to40km = -99 lower, upper : bounds for x (inclusive)
self.fP210to40km = -99
self.fChi210to40km = -99 Returns:
x, y data over the interval lower <= x <= upper
self.fP040to100km = -99 """
self.fP140to100km = -99 if lower is None:
self.fP240to100km = -99 lower = -np.inf
self.fChi240to100km = -99 if upper is None:
upper = np.inf
def fit(self, g, low, up, o="RWMQ0+", go="goff"): xdata = self.height
g.Fit(fQuad, o, go, low, up) ydata = self.log_density
p0 = self.fQuad.GetParameter(0) idx = np.where((xdata >= lower) & (xdata <= upper))
p1 = self.fQuad.GetParameter(1)/10. return xdata[idx], ydata[idx]
p2 = self.fQuad.GetParameter(2)/100.
chi2 = self.fQuad.GetChisquare() / self.fQuad.GetNDF()
return p0, p1, p2, chi2
def fit_profiles(self): def fit_profiles(self):
grho = new TGraph(len(self.fHeight) - 1, self.fHeight, self.fLogDensity) param_names = ['p0', 'p1', 'p2', 'chi2']
self.fit(grho, 0, 25, self.fP0Tot, self.fP1Tot, self.fP2Tot, self.fChi2Tot) # grho = new TGraph(len(self.height) - 1, self.height, self.log_density)
self.fit(grho, 0, 4, self.fP00to4km, self.fP10to4km, self.fP20to4km, self.fChi20to4km)
self.fit(grho, 4, 10, self.fP04to10km, self.fP14to10km, self.fP24to10km, self.fChi24to10km)
self.fit(grho, 10, 40, self.fP010to40km, self.fP110to40km, self.fP210to40km, self.fChi210to40km)
self.fit(grho, 40, 100, self.fP040to100km, self.fP140to100km, self.fP240to100km, self.fChi240to100km)
TSpline3 sst("spline", self.fHeight, self.fTemperature, len(self.fHeight)) self.params_total = zip(param_names, self.fit(lower=0, upper=25))
self.fTemperatureMAGIC = sst.Eval(AtcaWizard.H0_WS/1000) self.params_0_to_4_km = zip(param_names, self.fit(lower=0, upper=4))
self.params_4_to_10_km = zip(param_names, self.fit(lower=4, upper=10))
self.params_10_to_40_km = zip(param_names, self.fit(lower=10, upper=40))
self.params_40_to_100_km = zip(param_names, self.fit(lower=40, upper=100))
TSpline3 ss("spline", self.fHeight, self.fLogDensity, len(self.fHeight)) self.temperature_MAGIC = get_temperature(AtcaWizard.H0_WS/1000) # TODO: do we need this?
overburden = 0
size = len(self.height)
self.mass_overburden = np.zeros(size)
self.mass_overburden[size-1] = 0
self.log_mass_overburden = np.zeros(size)
self.log_mass_overburden[size-1] = 0
nbins = 100 nbins = 100
size = len(self.fHeight) overburden = 0
self.fMassOverburden = np.zeros(size) for i in range(size-2)[::-1]:
self.fMassOverburden[size-1] = 0 delta_height = (self.height[i+1] - self.height[i]) / nbins
self.fLogMassOverburden = np.zeros(size)
self.fLogMassOverburden[size-1] = 0
for i in range(size - 2)[::-1]:
dheight = (self.fHeight[i+1] - fHeight[i]) / nbins
for j in range(nbins): for j in range(nbins):
overburden += exp(ss.Eval(self.fHeight[i] + dheight*j)) * (self.gkNs/self.gkLoschmidt) * dheight * 100 # calculate mass overburden in g/cm^2 h = self.height[i] + j*delta_height
self.fMassOverburden[i] = overburden overburden += self.get_density(h)*(self.gkNs/self.gkLoschmidt)*delta_height*100 # calculate mass overburden in g/cm^2
self.fLogMassOverburden[i] = log(overburden) self.mass_overburden[i] = overburden
self.log_mass_overburden[i] = log(overburden)
def get_density(self, height): def get_density(self, height):
if height < self.fHeight[0]: f = interp1d(self.height, self.log_density, kind='cubic')
return -1 return exp(f(height))
size = len(self.fHeight)
TSpline3 ss("spline", self.fHeight, self.fLogDensity, size)
return exp(ss.Eval(height))
def get_temperature(self, height): def get_temperature(self, height):
if height < self.fHeight[0]: f = interp1d(self.height, self.temperature, kind='cubic') # NOTE: Changed from linear
return -1 return f(height)
size = len(self.fHeight)
for i in range(size):
if height < self.fHeight[i]:
return self.fTemperature[i-1] + (self.fTemperature[i] - self.fTemperature[i-1])*(height - self.fHeight[i-1])/(self.fHeight[i] - fHeight[i-1])
return -1
def get_mass_overburden(self, height): def get_mass_overburden(self, height):
if height < self.fHeight[0]: """
return -1 Interpolate logarithm of mass overburden, which behaves more linearly
size = len(self.fHeight) than exponential decay.
"""
# interpolate logarithm of mass overburden, which behaves more linearly than exponential decay f = interp1d(self.height, self.log_mass_overburden, kind='cubic')
TSpline3 ss("spline", self.fHeight, self.fLogMassOverburden, size) return exp(f(height))
return exp(ss.Eval(height))
# TODO: This isn't working yet...
def height_at_mass_overburden(self, overburden): def height_at_mass_overburden(self, overburden):
size = len(self.fMassOverburden) size = len(self.mass_overburden)
if overburden > self.fMassOverburden[0]: if overburden > self.mass_overburden[0]:
return -1 return -1
if overburden < self.fMassOverburden[size-1]: if overburden < self.mass_overburden[size-1]:
return 99 return 99
# search for the next pivot point # search for the next pivot point
i = 0 # i = 0
for j in range(size): # for j in range(size):
if self.fMassOverburden[j] < overburden: # if self.mass_overburden[j] < overburden:
i = j # i = j
break # break
#
TSpline3 ss("spline", self.fHeight, self.fLogMassOverburden, size) # TSpline3 ss("spline", self.height, self.log_mass_overburden, size)
ss.GetCoeff(i-1, x, y, b, c, d) # ss.GetCoeff(i-1, x, y, b, c, d)
#
delta = self.fHeight[i] - self.fHeight[i-1] # delta = self.height[i] - self.height[i-1]
nbins = 1000 # nbins = 1000
for j in range(nbins): # for j in range(nbins):
dx = j * delta / nbins # dx = j * delta / nbins
res = exp(y + dx*(b + dx*(c + dx*d))) # res = exp(y + dx*(b + dx*(c + dx*d)))
if res < overburden: # if res < overburden:
return self.fHeight[i-1] + dx # return self.height[i-1] + dx
# return -1
return -1
# TODO: FIXME!
def plot(self, lowlim=2.2, uplim=19): # def plot(self, lowlim=2.2, uplim=19):
c = new TCanvas("c", "c", 1000, 500) # c = new TCanvas("c", "c", 1000, 500)
c.Divide(2, 1) # c.Divide(2, 1)
c.cd(1) # c.cd(1)
#
grho = new TGraph(len(self.fHeight) - 1, self.fHeight, self.fLogDensity) # grho = new TGraph(len(self.height) - 1, self.height, self.log_density)
grho.SetTitle("Density;altitude a.s.l. (km);ln (n / n_s)") # grho.SetTitle("Density;altitude a.s.l. (km);ln (n / n_s)")
grho.SetMarkerStyle(7) # grho.SetMarkerStyle(7)
grho.Draw("AP") # grho.Draw("AP")
#
flin = new TF1("flin", "[0]+[1]*x/10", 0, 19) # flin = new TF1("flin", "[0]+[1]*x/10", 0, 19)
flin.SetParameter(0, -6.69) # flin.SetParameter(0, -6.69)
flin.SetParameter(1, -1.06) # flin.SetParameter(1, -1.06)
grho.Fit(flin, "RW", "same", lowlim, uplim) # grho.Fit(flin, "RW", "same", lowlim, uplim)
#
self.fQuad.SetLineColor(kBlue) # self.fQuad.SetLineColor(kBlue)
self.fQuad.SetParameter(0, flin.GetParameter(0)) # self.fQuad.SetParameter(0, flin.GetParameter(0))
self.fQuad.SetParameter(1, -0.95) # self.fQuad.SetParameter(1, -0.95)
self.fQuad.SetParameter(2, -0.17) # self.fQuad.SetParameter(2, -0.17)
#
p0, p1, p2, chi2 = self.fit(grho, lowlim, uplim, "RWM+", "same") # p0, p1, p2, chi2 = self.fit(grho, lowlim, uplim, "RWM+", "same")
print(p0, " ", p1, " ", p2, " ", chi2) # print(p0, " ", p1, " ", p2, " ", chi2)
#
gPad.SetGridx() # gPad.SetGridx()
gPad.SetGridy() # gPad.SetGridy()
#
c.cd(2) # c.cd(2)
#
gtemp = new TGraph(len(self.fHeight)-1, self.fHeight, self.fTemperature) # gtemp = new TGraph(len(self.height)-1, self.height, self.temperature)
gtemp.SetTitle("Temperature;altitude a.s.l. (km);T (K)") # gtemp.SetTitle("Temperature;altitude a.s.l. (km);T (K)")
gtemp.SetMarkerStyle(7) # gtemp.SetMarkerStyle(7)
gtemp.Draw("AP") # gtemp.Draw("AP")
#
gtemp.self.fit(flin, "RW", "same", lowlim, uplim) # gtemp.self.fit(flin, "RW", "same", lowlim, uplim)
self.fQuad.SetParameter(0, flin.GetParameter(0)) # self.fQuad.SetParameter(0, flin.GetParameter(0))
self.fQuad.SetParameter(1, flin.GetParameter(1)) # self.fQuad.SetParameter(1, flin.GetParameter(1))
self.fQuad.SetParameter(2, -0.17) # self.fQuad.SetParameter(2, -0.17)
p0, p1, p2, chi2 = self.fit(gtemp, lowlim, uplim, "RWM+", "same") # p0, p1, p2, chi2 = self.fit(gtemp, lowlim, uplim, "RWM+", "same")
#
gPad.SetGridx() # gPad.SetGridx()
gPad.SetGridy() # gPad.SetGridy()
msis_reader.py 0 → 100644
View file @ b296cd78
import re
class MSISReader:
def __init__(self, filename):
fObservatory.SetLocation(MObservatory::kMagic1)
def PreProcess(self, pList):
"""Check for correct file names format."""
TRegexp reg("[a-z]+msis[a-z0]+_[0-9]+_[0-9]+_[0-9]+_[0-9]+.dat")
TObject o
TIter Next(fFileNames)
TString file
while o=Next():
filename = o.GetName()
if reg not in filename:
print(file, " is not a valid MSIS file, it does not contain the string msis_[0-9]+_[0-9]+_[0-9]+_[0-9]+.dat")
return False
fReport = (MSISEntry*)pList.FindCreateObj("MSISEntry")
if not fReport:
return False
return True
def Process(self):
"""Read the files line by line and fill the MSISEntry class."""
if not OpenNextFile():
return False
print("NEXT FILE:", GetFullFileName())
line = fIn.readline()
# get rid of the possible inherent java-scripts and HTML-code
is_script = False
while line.startswith("<") or is_script:
if line.startswith("<script"):
is_script = True
if line.startswith("</script"):
is_script = False
line = fIn.readline()
fReport.Clear()
fReport.SetReadyToSave(False)
if not fReport.FillDate(line):
return kCONTINUE
line = fIn.readline()
line = fIn.readline()
line = fIn.readline()
if not fReport.FillObservatory(line):
return kCONTINUE
line = fIn.readline()
line = fIn.readline()
while line[0] == "<" or line.isspace():
line = fIn.readline()
if not line.startswith("1 Height, km"):
print("File ", GetFullFileName(), " does not include height as a first parameter, cannot interpret")
print(dbg, line)
return kCONTINUE
line = fIn.readline()
if not line.startswith("2 Mass_density, g/cm-3"):
print("File ", GetFullFileName(), " does not include mass density as a second parameter, cannot interpret")
print(dbg, line)
return kCONTINUE
line = fIn.readline()
if not line.startswith("3 Temperature_neutral, K"):
print("File ", GetFullFileName(), " does not include temperature as a third parameter, cannot interpret")
print(dbg, line)
return kCONTINUE
line = fIn.readline()
while line.isspace():
line = fIn.readline()
line = fIn.readline()
while not line.isspace() and not fIn.eof() and not line.startswith("<"):
if not fReport.FillProfile(line,fObservatory):
return kCONTINUE
line = fIn.readline()
fReport.FitProfiles()
print("READY")
fReport.SetReadyToSave()
return True
Supports Markdown
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