Compare commits
No commits in common. "4ad238a4170fafc04f73d0998221d669fd103535" and "195344e774b79906c002f82f54c43ed7caf44d97" have entirely different histories.
4ad238a417
...
195344e774
8 changed files with 55 additions and 735 deletions
86
CHAOSfish.py
86
CHAOSfish.py
|
|
@ -1,86 +0,0 @@
|
|||
import numpy as np
|
||||
import fileinput
|
||||
import sys
|
||||
import pandas as pd
|
||||
from datetime import datetime
|
||||
from pathlib import Path
|
||||
import re
|
||||
|
||||
ch = np.arange(18)
|
||||
global u; u = np.arange(18)
|
||||
name=[None]*18
|
||||
thr=[12]*18
|
||||
u = [1.0]*18
|
||||
s = [1.0]*18
|
||||
|
||||
name[0] = "A1H"; ch[0] = 13; thr[0] = 10; s[0] = 15.0; u[0] = 4.55
|
||||
name[1] = "A1L"; ch[1] = 12; thr[1] = 10; s[1] = 1.0; u[1] = 1.0
|
||||
name[2] = "A2H"; ch[2] = 10; thr[2] = 15; s[2] = 15.0; u[2] = 4.55
|
||||
name[3] = "A2L"; ch[3] = 9; thr[3] = 10; s[3] = 1.0; u[3] = 1.0
|
||||
name[4] = "B"; ch[4] = 17; thr[4] = 50; s[4] = 15.0; u[4] = 1.0
|
||||
name[5] = "C1H"; ch[5] = 7; thr[5] = 10; s[5] = 15.0; u[5] = 4.55
|
||||
name[6] = "C1L"; ch[6] = 6; thr[6] = 10; s[6] = 1.0; u[6] = 1.0
|
||||
name[7] = "C2H"; ch[7] = 4; thr[7] = 15; s[7] = 15.0; u[7] = 4.55
|
||||
name[8] = "C2L"; ch[8] = 3; thr[8] = 10; s[8] = 1.0; u[8] = 1.0
|
||||
name[9] = "D1H"; ch[9] = 16; thr[9] = 20; s[9] = 15.0; u[9] = 1.577
|
||||
name[10] = "D1L"; ch[10] = 15; thr[10] = 20; s[10] = 1.0; u[10] = 1.0
|
||||
name[11] = "D2H"; ch[11] = 14; thr[11] = 20; s[11] = 15.0; u[11] = 1.577
|
||||
name[12] = "D2L"; ch[12] = 11; thr[12] = 20; s[12] = 1.0; u[12] = 1.0
|
||||
name[13] = "EH"; ch[13] = 1; thr[13] = 10; s[13] = 15.0; u[13] = 4.55
|
||||
name[14] = "EL"; ch[14] = 0; thr[14] = 10; s[14] = 1.0; u[14] = 1.0
|
||||
name[15] = "E1"; ch[15] = 8; thr[15] = 8; s[15] = 2.2; u[15] = 11.32
|
||||
name[16] = "E2"; ch[16] = 5; thr[16] = 8; s[16] = 2.2; u[16] = 11.32
|
||||
name[17] = "E3"; ch[17] = 2; thr[17] = 8; s[17] = 2.2; u[17] = 11.32
|
||||
|
||||
mV = 14000
|
||||
minV = -100
|
||||
maxV = 3500
|
||||
resV = 4*0.838214/2
|
||||
pol = [1.88909231e-04, -6.84421799e-03, -3.28911747e+00, 3.24630271e+02]
|
||||
|
||||
bin = int((maxV-minV)/resV)
|
||||
hist = np.zeros((bin+1,bin+2))
|
||||
|
||||
edges = [minV + resV * i for i in range(bin + 1)]
|
||||
for i, edge in enumerate(edges):
|
||||
hist[i][0] = edge
|
||||
|
||||
def main():
|
||||
with fileinput.input(files=("data_flight/2024-10-02-flight-2.EI"), encoding="utf-8", errors='ignore') as f:
|
||||
for line in f:
|
||||
if not line.strip():
|
||||
continue
|
||||
l = line.split()
|
||||
if l[0] == 'EI' and len(l)>3*18 and conditions(l):
|
||||
try:
|
||||
bx,by = get_values(l)
|
||||
except:
|
||||
bx = maxV; by = maxV
|
||||
|
||||
if minV <= bx <= maxV and minV <= by <= maxV:
|
||||
bxx = int((bx-minV)/resV)
|
||||
byy = int((by-minV)/resV)+1
|
||||
hist[bxx,byy] += 1
|
||||
|
||||
frame = pd.DataFrame(hist, columns = ['keV']+edges)
|
||||
frame.to_csv(f'histograms/fish.2dhist', sep = ' ', index = False)
|
||||
|
||||
def CmV(line, i):
|
||||
return int(line[3*i+5])/mV
|
||||
|
||||
def istrigger(l,channame):
|
||||
return CmV(l, ch[name.index(channame)])>thr[name.index(channame)]
|
||||
|
||||
def value(l,channame):
|
||||
return CmV(l, ch[name.index(channame)])* u[name.index(channame)]
|
||||
|
||||
def conditions(l):
|
||||
return any([istrigger(l,"C1H"), istrigger(l,"C2H")]) and all(istrigger(l,chans) for chans in ["D1H", "D2H", "EH"])
|
||||
|
||||
def get_values(l):
|
||||
y = max(value(l,"C1H"),value(l,"C2H"))/value(l,"EH")
|
||||
x = value(l,"D1H")+value(l,"D2H")
|
||||
return x,y
|
||||
|
||||
main()
|
||||
|
||||
193
CHAOSx.py
193
CHAOSx.py
|
|
@ -5,7 +5,6 @@ import sys
|
|||
import pandas as pd
|
||||
from datetime import datetime
|
||||
from pathlib import Path
|
||||
import re
|
||||
|
||||
parser = argparse.ArgumentParser(description="This script parses EI-data for AHBGO(C/D/S) and returns histogramms")
|
||||
parser.add_argument("-v", "--verbose", choices = ['0','1','2'], help = "when zero only emit serious error messages; when one also emit warnings, when two emit informational messages")
|
||||
|
|
@ -15,16 +14,12 @@ parser.add_argument('-Itime', action='store_true', help = "if chosen, Itime is c
|
|||
|
||||
parser.add_argument('-hist', action='store_true', help = "makes histogramm of all channels using EI-lines") # on/off flag
|
||||
parser.add_argument('-tdhist', type=str, nargs='+', help = "need two arguments: channely, channelx -> makes 2D histogramm for channely and channelx, operand + and d (divison) possible") # on/off flag
|
||||
parser.add_argument("-fish", action='store_true', help = "does 2d histogramm with y: max(C1H,C2H)/EH; x:D1H+D2H")
|
||||
|
||||
|
||||
parser.add_argument("-mV", action='store_true', help = "needs option -hist or -tdhist; if chosen, data will be kept in mV changed to keV otherwise") # on/off flag
|
||||
parser.add_argument("-trigger", type=str, nargs='+', help = "needs option -EI; letters of Triggerchannels to chose")
|
||||
parser.add_argument("-xray", action='store_true', help = "needs option -EI; only select xray data") # on/off flag
|
||||
parser.add_argument("-float", action='store_true', help = "needs option -EI; only select data from floating phase") # on/off flag
|
||||
|
||||
parser.add_argument("-deltaT", action='store_true', help = "saves deltatimes between events")
|
||||
|
||||
|
||||
parser.add_argument("-nameadd", type=str, nargs='+', help = "custom nameadd to add to filename")
|
||||
|
||||
|
|
@ -72,7 +67,7 @@ mV = 14000
|
|||
minV = -100
|
||||
maxV = 3500
|
||||
resV = 0.838214/2
|
||||
if args.tdhist or args.fish:
|
||||
if args.tdhist:
|
||||
resV=4*resV
|
||||
minV = 0
|
||||
maxV = 2500
|
||||
|
|
@ -95,43 +90,6 @@ def main():
|
|||
generate_2d_histogram(args.tdhist)
|
||||
print('.2dhist was created')
|
||||
|
||||
if args.deltaT:
|
||||
generate_deltaTimes()
|
||||
print('.deltaT was created')
|
||||
|
||||
if args.fish:
|
||||
generate_2d_fish()
|
||||
print('fish.2dhist was created')
|
||||
|
||||
def generate_deltaTimes():
|
||||
resT = 10
|
||||
minT=0
|
||||
maxT=100000
|
||||
bin = int((maxT-minT)/resT)
|
||||
dTh = np.zeros((bin+1,2))
|
||||
dTh[:, 0] = np.arange(0, bin+1)
|
||||
with fileinput.input(files=(file), encoding="utf-8", errors='ignore') as f:
|
||||
for line in f:
|
||||
if not line.strip() or len(line.split()) < 10:
|
||||
continue
|
||||
l = line.split()
|
||||
if not 'EOF' in l and not 'X' in l and not 'xF' in l:
|
||||
consider=True
|
||||
if l[0] == 'H':
|
||||
timestamp = int(l[1])
|
||||
if not args.float:
|
||||
consider=True
|
||||
elif 1727859780 < timestamp < 1727872980: # between 2024-10-02 11:23:00 and 2024-10-02 15:03:00 (Sweden time)
|
||||
consider=True
|
||||
else: consider = False
|
||||
if l[0] == 'EI' and len(l)>3*18 and consider:
|
||||
dT = int(l[4].split(":")[0])
|
||||
#x = int(dT/3)
|
||||
x=int((dT/3-minT)/resT)
|
||||
dTh[x, 1] += 1
|
||||
frame = pd.DataFrame(dTh, columns = ["deltaT","count"])
|
||||
frame.to_csv('histograms/'+filename+'.deltaT', sep = ' ', index = False)
|
||||
|
||||
|
||||
def generate_histogram(u):
|
||||
|
||||
|
|
@ -183,7 +141,9 @@ def generate_histogram(u):
|
|||
elif 1727859780 < timestamp < 1727872980: # between 2024-10-02 11:23:00 and 2024-10-02 15:03:00 (Sweden time)
|
||||
consider=True
|
||||
else: consider = False
|
||||
if l[0] == 'EI' and len(l)>3*18 and conditions(l,args) and consider:
|
||||
|
||||
T = adc_to_T(int(l[8]))
|
||||
elif l[0] == 'EI' and len(l)>18 and conditions(l,args) and consider:
|
||||
for i, chan in enumerate(chans):
|
||||
b = CmV(l,ch[chan])*u[chan]
|
||||
if minV <= b <= maxV:
|
||||
|
|
@ -193,142 +153,27 @@ def generate_histogram(u):
|
|||
frame = pd.DataFrame(hist, columns = [energy] + name)
|
||||
frame.to_csv(f'histograms/{filename}{nameadd}.hist', sep = ' ', index = False)
|
||||
###########################################################################################
|
||||
# def parse_channel_expression(expr):
|
||||
# if '+' in expr:
|
||||
# return expr.split('+'), 'add'
|
||||
# elif 'd' in expr:
|
||||
# return expr.split('d'), 'div'
|
||||
# else:
|
||||
# return [expr], 'single'
|
||||
|
||||
# def get_channel_value(l, channames,op):
|
||||
# value = CmV(l, ch[name.index(channames[0])]) * u[name.index(channames[0])]
|
||||
# if len(channames) > 1:
|
||||
# if op == 'add':
|
||||
# value += CmV(l, ch[name.index(channames[1])]) * u[name.index(channames[1])]
|
||||
# elif op == 'div':
|
||||
# value /= CmV(l, ch[name.index(channames[1])]) * u[name.index(channames[1])]
|
||||
# return value
|
||||
#
|
||||
# Hilfsfunktion für die max-Funktion
|
||||
|
||||
def apply_max(channames, l, name, chan, u):
|
||||
values = [CmV(l, ch[name.index(chan)]) * u[name.index(ch)] for ch in channames]
|
||||
return max(values)
|
||||
|
||||
def parse_channel_expression(expr):
|
||||
expr = expr.replace(' ', '') # Entferne Leerzeichen
|
||||
|
||||
# Fall: max(A1H, A2H) nach einem 'd'
|
||||
if 'd' in expr and 'max(' in expr:
|
||||
parts = expr.split('d')
|
||||
outside_expr = parts[0] # Z.B. 'EH'
|
||||
|
||||
# Suche nach 'max()' nur, wenn es existiert
|
||||
match = re.search(r'max\((.*?)\)', expr)
|
||||
if match:
|
||||
inside = match.group(1)
|
||||
channames = inside.split(',')
|
||||
return [parts[1]], 'div', 'max', channames
|
||||
else:
|
||||
raise ValueError(f"Invalid expression with 'd' but no 'max()': {expr}")
|
||||
|
||||
# Fall: max(A1H, A2H) ohne 'd'
|
||||
elif 'max(' in expr:
|
||||
match = re.search(r'max\((.*?)\)', expr)
|
||||
if match:
|
||||
inside = match.group(1)
|
||||
channames = inside.split(',')
|
||||
return channames, 'max', None, None
|
||||
|
||||
# Normaler Fall ohne max(), aber mit '+', 'd', oder einfachem Ausdruck
|
||||
elif '+' in expr:
|
||||
return expr.split('+'), 'add', None, None
|
||||
if '+' in expr:
|
||||
return expr.split('+'), 'add'
|
||||
elif 'd' in expr:
|
||||
return expr.split('d'), 'div', expr.split('d')[1], None
|
||||
return expr.split('d'), 'div'
|
||||
else:
|
||||
return [expr], 'single', None, None
|
||||
return [expr], 'single'
|
||||
|
||||
|
||||
def get_channel_value(l, channames, op, next_op=None, next_channel=None):
|
||||
value = None
|
||||
if op == 'max':
|
||||
value = apply_max(channames, l, name, ch, u)
|
||||
else:
|
||||
def get_channel_value(l, channames,op):
|
||||
value = CmV(l, ch[name.index(channames[0])]) * u[name.index(channames[0])]
|
||||
if len(channames) > 1 and op == 'add':
|
||||
if len(channames) > 1:
|
||||
if op == 'add':
|
||||
value += CmV(l, ch[name.index(channames[1])]) * u[name.index(channames[1])]
|
||||
|
||||
if next_op == 'div' and next_channel:
|
||||
value /= CmV(l, ch[name.index(next_channel)]) * u[name.index(next_channel)]
|
||||
elif next_op == 'div_max' and next_channel:
|
||||
max_value = apply_max(next_channel, l, name, ch, u)
|
||||
value /= max_value
|
||||
|
||||
elif op == 'div':
|
||||
value /= CmV(l, ch[name.index(channames[1])]) * u[name.index(channames[1])]
|
||||
return value
|
||||
|
||||
def generate_2d_fish():
|
||||
energy = 'keV'
|
||||
nameadd=''
|
||||
|
||||
if args.trigger:
|
||||
nameadd += '_trigger_' + '_'.join(args.trigger)
|
||||
|
||||
if args.nameadd:
|
||||
nameadd = args.nameadd[0]
|
||||
|
||||
global trigchan
|
||||
for t in trig:
|
||||
for i in range(len(name)):
|
||||
if t == name[i]:
|
||||
trigchan.append(i)
|
||||
#trigchan = [i for t in args.trigger for i, n in enumerate(name) if t == n]
|
||||
|
||||
bin = int((maxV-minV)/resV)
|
||||
hist = np.zeros((bin+1,bin+2))
|
||||
|
||||
edges = [minV + resV * i for i in range(bin + 1)]
|
||||
for i, edge in enumerate(edges):
|
||||
hist[i][0] = edge
|
||||
|
||||
|
||||
with fileinput.input(files=(file), encoding="utf-8", errors='ignore') as f:
|
||||
timestamp=0
|
||||
|
||||
for line in f:
|
||||
if not line.strip():
|
||||
continue
|
||||
l = line.split()
|
||||
cons=False
|
||||
if l[0] == 'H':
|
||||
timestamp = int(l[1])
|
||||
if not args.float:
|
||||
cons=True
|
||||
elif 1727859780 < timestamp < 1727872980: # between 2024-10-02 11:23:00 and 2024-10-02 15:03:00 (Sweden time)
|
||||
cons=True
|
||||
else: cons = False
|
||||
if l[0] == 'EI' and len(l)>3*18 and conditions(l,args) and cons:
|
||||
try:
|
||||
by = max(CmV(l, ch[name.index('C1H')]),CmV(l, ch[name.index('C2H')]))/CmV(l, ch[name.index('EH')])
|
||||
bx = CmV(l, ch[name.index('D1H')])+CmV(l, ch[name.index('D2H')])
|
||||
except:
|
||||
bx = maxV; by = maxV
|
||||
|
||||
if minV <= bx <= maxV and minV <= by <= maxV:
|
||||
bxx = int((bx-minV)/resV)
|
||||
byy = int((by-minV)/resV)+1
|
||||
hist[bxx,byy] += 1
|
||||
|
||||
frame = pd.DataFrame(hist, columns = ['keV']+edges)
|
||||
frame.to_csv(f'histograms/{filename}_fish{nameadd}~max(C1H,C2H)dEH-D1H+D2H.2dhist', sep = ' ', index = False)
|
||||
|
||||
def generate_2d_histogram(chans):
|
||||
"""Generate a 2D histogram for two specified channels."""
|
||||
#x_expr, x_op = parse_channel_expression(chans[1])
|
||||
#y_expr, y_op = parse_channel_expression(chans[0])
|
||||
x_expr, x_op, x_next_op, x_next_channel = parse_channel_expression(chans[1])
|
||||
y_expr, y_op, y_next_op, y_next_channel = parse_channel_expression(chans[0])
|
||||
|
||||
x_expr, x_op = parse_channel_expression(chans[1])
|
||||
y_expr, y_op = parse_channel_expression(chans[0])
|
||||
|
||||
nameadd='x'
|
||||
energy = 'keV'
|
||||
|
|
@ -369,10 +214,8 @@ def generate_2d_histogram(chans):
|
|||
l = line.split()
|
||||
if l[0] == 'EI' and len(l)>3*18 and conditions(l,args):
|
||||
try:
|
||||
#bx = get_channel_value(l,x_expr,x_op)
|
||||
#by = get_channel_value(l,y_expr,y_op)
|
||||
bx = get_channel_value(l, x_expr, x_op, x_next_op, x_next_channel)
|
||||
by = get_channel_value(l, y_expr, y_op, y_next_op, y_next_channel)
|
||||
bx = get_channel_value(l,x_expr,x_op)
|
||||
by = get_channel_value(l,y_expr,y_op)
|
||||
except:
|
||||
bx = maxV; by = maxV
|
||||
|
||||
|
|
@ -382,7 +225,7 @@ def generate_2d_histogram(chans):
|
|||
hist[bxx,byy] += 1
|
||||
|
||||
frame = pd.DataFrame(hist, columns = ['keV']+edges)
|
||||
frame.to_csv(f'histograms/{filename}_{nameadd}~{chans[0]}-{chans[1]}.2dhist', sep = ' ', index = False)
|
||||
frame.to_csv(f'histograms/{filename}{nameadd}~{chans[0]}-{chans[1]}.2dhist', sep = ' ', index = False)
|
||||
|
||||
def calculate_Itime():
|
||||
"""Calculate Itime and write to file."""
|
||||
|
|
|
|||
64
HK_flight.py
64
HK_flight.py
|
|
@ -11,7 +11,6 @@ import matplotlib.dates as md
|
|||
from datetime import datetime
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
from us_standard_atmosphere import *
|
||||
|
||||
file_hk = "hk_data/2024-10-02-flight-2"
|
||||
|
||||
|
|
@ -39,14 +38,13 @@ names_c64 = ["time", "A1H", "A1L", "A2H", "A2L", "B", "C1H", "C1L", "C2H", "C2L"
|
|||
data_c64 = pd.read_csv(file_hk + ".c64", sep="\s+", names=names_c64, skiprows=1, on_bad_lines='skip')
|
||||
|
||||
# Global font size
|
||||
fs = 30
|
||||
fs = 25
|
||||
|
||||
def main():
|
||||
# Plot all graphs
|
||||
#plot_and_layout(plot_allT,"","Esrange_HK_T.pdf")
|
||||
plot_and_layout(plot_allT, "","recap_HK_Ts")
|
||||
plot_and_layout(plot_pressures, "Pressures","HK_ps")
|
||||
plot_and_layout(plot_I_and_V, "","HK_VI")
|
||||
plot_and_layout(plot_allT, "All Temperatures","HK_Ts.pdf")
|
||||
plot_and_layout(plot_pressures, "Pressures","HK_ps.pdf")
|
||||
plot_and_layout(plot_I_and_V, "Current and Voltage","HK_VI.pdf")
|
||||
|
||||
# Channels grouped by type
|
||||
allchans = ["A1H", "A1L", "A2H", "A2L", "B", "C1H", "C1L", "C2H", "C2L", "D1H", "D1L", "D2H", "D2L", "EH", "EL", "E1", "E2", "E3"]
|
||||
|
|
@ -54,24 +52,21 @@ def main():
|
|||
AC = ["A1H", "A2H", "C1H", "C2H"]
|
||||
E = ["EH", "E1", "E2", "E3"]
|
||||
|
||||
plot_and_layout(lambda ax: plot_c64(ax, E), "","recap_HK_countsE")
|
||||
plot_and_layout(lambda ax: plot_c64(ax, AC), "","recap_HK_countsAC")
|
||||
plot_and_layout(lambda ax: plot_c64(ax, DB), "","recap_HK_countsDB")
|
||||
plot_and_layout(lambda ax: plot_c64(ax, E), "Counts channels E","plots/HK_countsE.pdf")
|
||||
plot_and_layout(lambda ax: plot_c64(ax, AC), "Counts channels A and C","plots/HK_countsAC.pdf")
|
||||
plot_and_layout(lambda ax: plot_c64(ax, DB), "Counts channels D and B","plots/HK_countsDB")
|
||||
|
||||
#plot_and_layout(plot_pressures_and_height, "","recap_HK_p")
|
||||
plot_six()
|
||||
|
||||
#plot_six()
|
||||
|
||||
plt.show()
|
||||
#plt.show()
|
||||
|
||||
|
||||
def plot_and_layout(plot_func, title, savefile):
|
||||
fig, ax = plt.subplots(figsize=(20, 13))
|
||||
fig, ax = plt.subplots(figsize=(20, 15))
|
||||
ax.set_title(title, fontsize=fs + 7, pad =fs-5)
|
||||
plot_func(ax)
|
||||
layout(ax)
|
||||
if savefile!= None:
|
||||
fig.tight_layout()
|
||||
fig.savefig("plots/"+savefile)
|
||||
else:plt.show()
|
||||
|
||||
|
|
@ -85,7 +80,7 @@ def plot_allT(ax):
|
|||
ax.plot(pd.to_datetime(data_h["time"]),data_h["T"],linewidth=3,label = "$T_{BGO}$")
|
||||
#ax.plot(pd.to_datetime(data_h["time"]),data_h["Tava"],label = "$T_{ava}$")
|
||||
ax.plot(pd.to_datetime(data_hc["time"]),data_hc["Tntc"],linewidth=3,label = "$T_{ntc}$")
|
||||
ax.set_ylabel("T in °C", fontsize = fs)
|
||||
ax.set_ylabel("T [°C]", fontsize = fs)
|
||||
|
||||
def plot_T_and_p(ax):
|
||||
ax2 = ax.twinx()
|
||||
|
|
@ -95,35 +90,13 @@ def plot_T_and_p(ax):
|
|||
|
||||
ax2.plot(pd.to_datetime(data_p["time1"]),data_p["p1"], label = names_p[2])
|
||||
ax2.plot(pd.to_datetime(data_p["time2"]),data_p["p2"], label = names_p[5])
|
||||
ax.set_ylabel("T in °C", fontsize = fs)
|
||||
ax2.set_ylabel("p in mbar", fontsize = fs)
|
||||
ax.set_ylabel("T [°C]", fontsize = fs)
|
||||
ax2.set_ylabel("p [mbar]", fontsize = fs)
|
||||
|
||||
def plot_pressures(ax):
|
||||
ax.plot(pd.to_datetime(data_p["time1"]),data_p["p1"], linewidth=3, label = names_p[2])
|
||||
ax.plot(pd.to_datetime(data_p["time2"]),data_p["p2"], linewidth=3, label = names_p[5])
|
||||
ax.set_ylabel("p in mbar", fontsize = fs)
|
||||
|
||||
def plot_pressures_and_height(ax):
|
||||
ax.plot(pd.to_datetime(data_p["time1"]),data_p["p1"], linewidth=3, label = names_p[2])
|
||||
ax.plot(pd.to_datetime(data_p["time2"]),data_p["p2"], linewidth=3, label = names_p[5])
|
||||
ax.set_ylabel("p in mbar", fontsize = fs)
|
||||
|
||||
print(data_p["p2"])
|
||||
|
||||
heights=[]
|
||||
for p in data_p["p2"]:
|
||||
heights.append(geometric_height(geopotential_height(float(p)*100))/1000)
|
||||
|
||||
ax2=ax.twinx()
|
||||
ax2.plot(pd.to_datetime(data_p["time2"][0]),heights[0], linewidth=0.1)
|
||||
ax2.plot(pd.to_datetime(data_p["time2"][0]),heights[0], linewidth=0.1)
|
||||
ax2.plot(pd.to_datetime(data_p["time2"]),heights, linewidth=3,label="estimated height")
|
||||
ax2.set_ylabel("height in km", fontsize = fs)
|
||||
ax2.legend(fontsize=fs-2, loc='upper right')
|
||||
ax2.tick_params(axis='y', labelsize=fs)
|
||||
ax.set_ylim(0,1200)
|
||||
ax2.set_ylim(0,33)
|
||||
#geopotential_height()
|
||||
ax.set_ylabel("p [mbar]", fontsize = fs)
|
||||
|
||||
def plot_I_and_V(ax):
|
||||
ax2 = ax.twinx()
|
||||
|
|
@ -135,16 +108,15 @@ def plot_I_and_V(ax):
|
|||
ax.plot(pd.to_datetime(data_hc["time"]),data_hc["Ibias"], linewidth=0.8, label = "$I_{biasHC}$", color='b')
|
||||
#ax.plot(pd.to_datetime(data_hc["time"]),data_hc["Idrv"], label = "$I_{drv}$")
|
||||
|
||||
ax2.plot(pd.to_datetime(data_h["time"][0]),data_h["Vbias"][0])
|
||||
#ax2.plot(pd.to_datetime(data_h["time"][0]),data_h["Vbias"][0])
|
||||
|
||||
ax2.plot(pd.to_datetime(data_hc["time"]),np.divide(data_hc["HV"],100), linewidth=3, label = "HV $\cdot 10^{-2}$")
|
||||
#ax2.plot(pd.to_datetime(data_h["time"]),data_h["Vbias"], label = "$V_{bias}$")
|
||||
ax2.plot(pd.to_datetime(data_hc["time"]),data_hc["Vprim"], linewidth=3, label = "$V_{prim}$")
|
||||
ax2.plot(pd.to_datetime(data_h["time"]),data_h["Vcc"], linewidth=3, label = "$V_{cc}$")
|
||||
ax2.plot(pd.to_datetime(data_h["time"]),data_h["Vss"], linewidth=3, label = "$V_{ss}$")
|
||||
ax2.plot(pd.to_datetime(data_hc["time"]),data_hc["Vprim"], linewidth=3, label = "$V_{prim}$")
|
||||
#ax2.plot(pd.to_datetime(data_hc["time"]),data_hc["Vana"], label = "$V_{ana}$")
|
||||
#ax2.plot(pd.to_datetime(data_hc["time"]),data_hc["Vdrv"], label = "$V_{drv}$")
|
||||
ax2.tick_params(axis='y', labelsize=fs)
|
||||
|
||||
ax.set_ylabel("Current I in A", fontsize = fs)
|
||||
ax2.set_ylabel("Voltage V in V", fontsize = fs)
|
||||
|
|
@ -160,10 +132,10 @@ def plot_c64(ax,chans):
|
|||
def add_highlighted_events(ax):
|
||||
for time, event in zip(highlighted_times, events):
|
||||
ax.axvline(x=time, color='k', linestyle='--')
|
||||
ax.text(time, ax.get_ylim()[1] * 0.98, event, rotation=90, verticalalignment='top', horizontalalignment='right', fontsize=fs - 2)
|
||||
ax.text(time, ax.get_ylim()[1] * 0.95, event, rotation=90, verticalalignment='top', horizontalalignment='right', fontsize=fs - 2)
|
||||
|
||||
def layout(ax):
|
||||
ax.legend(fontsize=fs-2, loc='upper left')
|
||||
ax.legend(fontsize=fs - 2, loc='upper left')
|
||||
ax.tick_params(axis='x', labelsize=fs - 3)
|
||||
ax.tick_params(axis='y', labelsize=fs)
|
||||
ax.xaxis.set_major_formatter(md.DateFormatter('%H:%M'))
|
||||
|
|
|
|||
|
|
@ -39,11 +39,7 @@ def main():
|
|||
|
||||
|
||||
file = args.filepath.split('.')[0]
|
||||
try:
|
||||
y,x = file.split('~')[-1].split('-')
|
||||
except:
|
||||
x='unknown'
|
||||
y='unknown'
|
||||
|
||||
end = '.2dhist'
|
||||
energy = 'keV'
|
||||
|
|
@ -84,7 +80,7 @@ def main():
|
|||
c = np.array(data)[:,1:]
|
||||
|
||||
#cb = plt.pcolormesh(data[energy],data[energy],np.multiply(c, 3600/Itime/resV), cmap = cmap, shading='nearest', vmin=0.0001)
|
||||
cb = plt.pcolormesh(data[energy],data[energy],c, cmap = cmap, shading='nearest', vmin=0.0001, norm='log')
|
||||
cb = plt.pcolormesh(data[energy],data[energy],c, cmap = cmap, shading='nearest', vmin=0.0001)
|
||||
|
||||
#scatter = plt.scatter(data[x], data[y], c=data[c], vmin = 0.1)
|
||||
cbar = plt.colorbar(cb, orientation = 'vertical', shrink = 0.95, pad = 0.1)
|
||||
|
|
@ -93,13 +89,13 @@ def main():
|
|||
fs = 15
|
||||
lim = [None, None]
|
||||
|
||||
# if not args.limits == None:
|
||||
# lim[0] = float(args.limits[0])
|
||||
# lim[1] = float(args.limits[1])
|
||||
# if lim[0] == None:
|
||||
# lim[0] = 0
|
||||
# if lim[1] == None:
|
||||
# lim[1] = 1000
|
||||
if not args.limits == None:
|
||||
lim[0] = float(args.limits[0])
|
||||
lim[1] = float(args.limits[1])
|
||||
if lim[0] == None:
|
||||
lim[0] = 0
|
||||
if lim[1] == None:
|
||||
lim[1] = 500
|
||||
|
||||
ylab = y + ' in ' + energy+'$_{Si}$'
|
||||
xlab = x + ' in ' + energy+'$_{Si}$'
|
||||
|
|
@ -115,11 +111,8 @@ def main():
|
|||
|
||||
plt.ylabel(ylab, fontsize = fs+5)
|
||||
plt.xlabel(xlab, fontsize=fs+5)
|
||||
plt.ylim(1,10000)
|
||||
plt.xlim(1,10000)
|
||||
|
||||
plt.yscale('log')
|
||||
plt.xscale('log')
|
||||
plt.ylim(lim[0],lim[1])
|
||||
plt.xlim(lim[0],lim[1])
|
||||
|
||||
cbar.set_label('counts', fontsize = fs+5)
|
||||
cbar.ax.tick_params(labelsize=fs)
|
||||
|
|
|
|||
|
|
@ -1,88 +0,0 @@
|
|||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Created on Wed Aug 21 13:33:03 2024
|
||||
|
||||
@author: ava
|
||||
"""
|
||||
|
||||
import argparse
|
||||
import sys
|
||||
import matplotlib.pyplot as plt
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
import os.path
|
||||
|
||||
def main():
|
||||
|
||||
parser = argparse.ArgumentParser(description="This script plots histogramms of given data ending on '.hist'")
|
||||
parser.add_argument("-v", "--verbose", choices = ['0','1','2'], help = "when zero only emit serious error messages; when one also emit warnings, when two emit informational messages")
|
||||
parser.add_argument("filepath", type = str, help="path/file were data is stored in")
|
||||
parser.add_argument("-t", "--title", type = str, help="title to use in plot; if not given a standard title will be chosen")
|
||||
parser.add_argument("-e", "--export", type = str, nargs='+', help = "first: filepath were figure shall be exported to (if f, filepath of histdata will be chosen with filepath_plot) \n second:format (png,pdf,...)")
|
||||
parser.add_argument("-l", "--limits", type=str, nargs='+',help = "limits of plot in order xmin,xmax, ymin,ymax,")
|
||||
|
||||
args = parser.parse_args()
|
||||
|
||||
if args.verbose == None:
|
||||
verbosity = 1
|
||||
else:
|
||||
verbosity = args.verbose
|
||||
|
||||
if verbosity == 2:
|
||||
sys.stderr.write("--verbostiy = '" + str(args.verbose)+"'\n")
|
||||
sys.stderr.write("--filepath = '" + str(args.filepath)+"'\n")
|
||||
sys.stderr.write("--title = '" + str(args.title)+"'\n")
|
||||
sys.stderr.write("--export = '" + str(args.export) +"'\n")
|
||||
sys.stderr.write("--limits = '" + str(args.limits) +"'\n")
|
||||
|
||||
file = args.filepath
|
||||
|
||||
if args.title == None:
|
||||
tit = 'Histplot, '+"'"+file.split('/')[-1]+"' "
|
||||
else:
|
||||
tit = args.title
|
||||
|
||||
data = pd.read_csv(file, sep='\s+', skiprows=1, names = ["dTime","counts"])
|
||||
|
||||
plt.figure(figsize=(15,10))
|
||||
|
||||
plt.step(data["dTime"], data["counts"], label='deltatimes')
|
||||
|
||||
fs = 15
|
||||
lim = [-100, 9000, 0.9, 30000]
|
||||
|
||||
if not args.limits == None:
|
||||
for i in range(0,4):
|
||||
if not args.limits[i] == None:
|
||||
lim[i] = float(args.limits[i])
|
||||
|
||||
xlab = 'deltatime in $\mu$s'
|
||||
ylab = 'counts'
|
||||
plt.legend(fontsize = fs, loc="upper right")
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
plt.title(tit, fontsize = fs+10, pad = 15)
|
||||
|
||||
plt.yscale('log')
|
||||
plt.ylabel(ylab, fontsize = fs+5)
|
||||
plt.xlabel(xlab, fontsize=fs+5)
|
||||
plt.ylim(lim[2],lim[3])
|
||||
plt.xlim(lim[0],lim[1])
|
||||
|
||||
|
||||
if args.export == None:
|
||||
plt.show()
|
||||
else:
|
||||
if args.export[0] == 'f':
|
||||
name = 'plots/dThist_' + file.split('/')[-1].split('.')[0]
|
||||
else:
|
||||
name = args.export[0]
|
||||
if len(args.export)>1:
|
||||
end = '.' + args.export[1]
|
||||
else:
|
||||
end=''
|
||||
plt.savefig(name+end)
|
||||
|
||||
main()
|
||||
|
||||
114
plotfish.py
114
plotfish.py
|
|
@ -1,114 +0,0 @@
|
|||
'''
|
||||
Autor: Ava Pohley
|
||||
'''
|
||||
|
||||
import argparse
|
||||
import sys
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib as mpl
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
import os.path
|
||||
|
||||
|
||||
|
||||
def main():
|
||||
|
||||
parser = argparse.ArgumentParser(description="This script plots histogramms of given data ending on '.2dhist'")
|
||||
parser.add_argument("-v", "--verbose", choices = ['0','1','2'], help = "when zero only emit serious error messages; when one also emit warnings, when two emit informational messages")
|
||||
parser.add_argument("filepath", type = str, help="path/file were data is stored in")
|
||||
parser.add_argument("-t", "--title", type = str, help="title to use in plot; if not given a standard title will be chosen")
|
||||
parser.add_argument("-e", "--export", type = str, nargs='+', help = "first: filepath were figure shall be exported to (if f, filepath of histdata will be chosen with filepath_plot) \n second:format (png,pdf,...)")
|
||||
parser.add_argument('-r', '--raw', action='store_true', help = "just prints data without considering itime") # on/off flag
|
||||
parser.add_argument("-l", "--limits", type=str, nargs='+',help = "limits of plot in order xmin,xmax, ymin,ymax,")
|
||||
|
||||
args = parser.parse_args()
|
||||
|
||||
if args.verbose == None:
|
||||
verbosity = 1
|
||||
else:
|
||||
verbosity = args.verbose
|
||||
|
||||
if verbosity == 2:
|
||||
sys.stderr.write("--verbostiy = '" + str(args.verbose)+"'\n")
|
||||
sys.stderr.write("--filepath = '" + str(args.filepath)+"'\n")
|
||||
sys.stderr.write("--title = '" + str(args.title)+"'\n")
|
||||
sys.stderr.write("--export = '" + str(args.export) +"'\n")
|
||||
sys.stderr.write("--raw = '" + str(args.raw) +"'\n")
|
||||
sys.stderr.write("--limits = '" + str(args.limits) +"'\n")
|
||||
|
||||
|
||||
file = args.filepath.split('.')[0]
|
||||
|
||||
end = '.2dhist'
|
||||
energy = 'keV'
|
||||
|
||||
if args.title == None:
|
||||
tit = 'Histplot2D, '+"'"+file.split('/')[-1].split('~')[0]+"' "
|
||||
else:
|
||||
tit = args.title
|
||||
|
||||
|
||||
plt.figure(figsize=(15,10))
|
||||
|
||||
|
||||
energy='keV'
|
||||
|
||||
data = pd.read_csv(file+end, sep='\s+',skiprows=0)
|
||||
|
||||
#cmap = cm.get_cmap()
|
||||
cmap = mpl.colormaps['viridis']
|
||||
cmap.set_under('white')
|
||||
|
||||
c = np.array(data)[:,1:]
|
||||
|
||||
#cb = plt.pcolormesh(data[energy],data[energy],np.multiply(c, 3600/Itime/resV), cmap = cmap, shading='nearest', vmin=0.0001)
|
||||
cb = plt.pcolormesh(data[energy],data[energy],c, cmap = cmap, shading='nearest', vmin=0.0001, norm='log')
|
||||
|
||||
#scatter = plt.scatter(data[x], data[y], c=data[c], vmin = 0.1)
|
||||
cbar = plt.colorbar(cb, orientation = 'vertical', shrink = 0.95, pad = 0.1)
|
||||
|
||||
|
||||
fs = 15
|
||||
lim = [None, None]
|
||||
|
||||
# if not args.limits == None:
|
||||
# lim[0] = float(args.limits[0])
|
||||
# lim[1] = float(args.limits[1])
|
||||
# if lim[0] == None:
|
||||
# lim[0] = 0
|
||||
# if lim[1] == None:
|
||||
# lim[1] = 1000
|
||||
|
||||
ylab = 'C/E'
|
||||
xlab = 'D' + ' in ' + energy+'$_{Si}$'
|
||||
|
||||
plt.xticks(fontsize=fs)
|
||||
plt.yticks(fontsize=fs)
|
||||
plt.title(tit, fontsize = fs+10, pad = 15)
|
||||
|
||||
plt.ylabel(ylab, fontsize = fs+5)
|
||||
plt.xlabel(xlab, fontsize=fs+5)
|
||||
#plt.ylim(1,10000)
|
||||
#plt.xlim(1,10000)
|
||||
|
||||
plt.yscale('log')
|
||||
plt.xscale('log')
|
||||
|
||||
cbar.set_label('counts', fontsize = fs+5)
|
||||
cbar.ax.tick_params(labelsize=fs)
|
||||
|
||||
if args.export == None:
|
||||
plt.show()
|
||||
else:
|
||||
if args.export[0] == 'f':
|
||||
name = 'plots/plot_' + file.split('/')[-1]
|
||||
else:
|
||||
name = args.export[0]
|
||||
if len(args.export)>1:
|
||||
end = '.' + args.export[1]
|
||||
else:
|
||||
end=''
|
||||
plt.savefig(name+end)
|
||||
|
||||
main()
|
||||
84
tryfish.py
84
tryfish.py
|
|
@ -1,84 +0,0 @@
|
|||
import numpy as np
|
||||
import fileinput
|
||||
import pandas as pd
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
ch = np.arange(18)
|
||||
global u; u = np.arange(18)
|
||||
name=[None]*18
|
||||
thr=[12]*18
|
||||
u = [1.0]*18
|
||||
s = [1.0]*18
|
||||
|
||||
name[0] = "A1H"; ch[0] = 13; thr[0] = 10; s[0] = 15.0; u[0] = 4.55
|
||||
name[1] = "A1L"; ch[1] = 12; thr[1] = 10; s[1] = 1.0; u[1] = 1.0
|
||||
name[2] = "A2H"; ch[2] = 10; thr[2] = 15; s[2] = 15.0; u[2] = 4.55
|
||||
name[3] = "A2L"; ch[3] = 9; thr[3] = 10; s[3] = 1.0; u[3] = 1.0
|
||||
name[4] = "B"; ch[4] = 17; thr[4] = 50; s[4] = 15.0; u[4] = 1.0
|
||||
name[5] = "C1H"; ch[5] = 7; thr[5] = 10; s[5] = 15.0; u[5] = 4.55
|
||||
name[6] = "C1L"; ch[6] = 6; thr[6] = 10; s[6] = 1.0; u[6] = 1.0
|
||||
name[7] = "C2H"; ch[7] = 4; thr[7] = 15; s[7] = 15.0; u[7] = 4.55
|
||||
name[8] = "C2L"; ch[8] = 3; thr[8] = 10; s[8] = 1.0; u[8] = 1.0
|
||||
name[9] = "D1H"; ch[9] = 16; thr[9] = 20; s[9] = 15.0; u[9] = 1.577
|
||||
name[10] = "D1L"; ch[10] = 15; thr[10] = 20; s[10] = 1.0; u[10] = 1.0
|
||||
name[11] = "D2H"; ch[11] = 14; thr[11] = 20; s[11] = 15.0; u[11] = 1.577
|
||||
name[12] = "D2L"; ch[12] = 11; thr[12] = 20; s[12] = 1.0; u[12] = 1.0
|
||||
name[13] = "EH"; ch[13] = 1; thr[13] = 10; s[13] = 15.0; u[13] = 4.55
|
||||
name[14] = "EL"; ch[14] = 0; thr[14] = 10; s[14] = 1.0; u[14] = 1.0
|
||||
name[15] = "E1"; ch[15] = 8; thr[15] = 8; s[15] = 2.2; u[15] = 11.32
|
||||
name[16] = "E2"; ch[16] = 5; thr[16] = 8; s[16] = 2.2; u[16] = 11.32
|
||||
name[17] = "E3"; ch[17] = 2; thr[17] = 8; s[17] = 2.2; u[17] = 11.32
|
||||
|
||||
# Voltage parameters
|
||||
mV = 14000
|
||||
minV = 1e-2 # Adjust the minV to be above 0 for log scale
|
||||
maxV = 3500
|
||||
|
||||
# Define logarithmic bin edges
|
||||
n_bins = 50 # Adjust number of bins if needed
|
||||
edges = np.logspace(np.log10(minV), np.log10(maxV), n_bins)
|
||||
|
||||
# Initialize histogram
|
||||
hist = np.zeros((n_bins, n_bins))
|
||||
|
||||
def main():
|
||||
with fileinput.input(files=("data_flight/2024-10-02-flight-2.EI"), encoding="utf-8", errors='ignore') as f:
|
||||
for line in f:
|
||||
if not line.strip():
|
||||
continue
|
||||
l = line.split()
|
||||
if l[0] == 'EI' and len(l)>3*18 and conditions(l):
|
||||
try:
|
||||
bx,by = get_values(l)
|
||||
except:
|
||||
bx = maxV; by = maxV
|
||||
|
||||
if minV <= bx <= maxV and minV <= by <= maxV:
|
||||
# Find appropriate bins for the logarithmic scale
|
||||
bxx = np.digitize(bx, edges) - 1
|
||||
byy = np.digitize(by, edges) - 1
|
||||
if 0 <= bxx < n_bins and 0 <= byy < n_bins:
|
||||
hist[bxx, byy] += 1
|
||||
|
||||
# Save the 2D histogram to a CSV file
|
||||
frame = pd.DataFrame(hist, columns = edges[:-1], index=edges[:-1])
|
||||
frame.to_csv(f'histograms/tryfish.2dhist', sep=' ', index=False)
|
||||
|
||||
def CmV(line, i):
|
||||
return int(line[3*i+5])/mV
|
||||
|
||||
def istrigger(l, channame):
|
||||
return CmV(l, ch[name.index(channame)]) > thr[name.index(channame)]
|
||||
|
||||
def value(l, channame):
|
||||
return CmV(l, ch[name.index(channame)]) * u[name.index(channame)]
|
||||
|
||||
def conditions(l):
|
||||
return any([istrigger(l,"C1H"), istrigger(l,"C2H")]) and all(istrigger(l, chans) for chans in ["D1H", "D2H", "EH"])
|
||||
|
||||
def get_values(l):
|
||||
y = max(value(l,"C1H"), value(l,"C2H")) / value(l,"EH")
|
||||
x = value(l,"D1H") + value(l,"D2H")
|
||||
return x, y
|
||||
|
||||
main()
|
||||
|
|
@ -1,116 +0,0 @@
|
|||
###########################################################
|
||||
### Implementation of the U.S. Standard Atmosphere 1976 ###
|
||||
###########################################################
|
||||
|
||||
import math
|
||||
import numpy as np
|
||||
|
||||
# Define constants
|
||||
R_star = 8.31432 # Gas constant for dry air, J/(mol·K)
|
||||
g0 = 9.80665 # Standard gravitational acceleration at sea-level, m/s^2
|
||||
g0_prime = 9.80665 # Geopotential constant, m^2/(s^2*m_prime), geopotential meter m_prime represents work done to lift unit mass 1 geometrical meter in standard gravitation
|
||||
P0 = 101325 # Standard atmospheric pressure at sea-level, Pa = N/m^2
|
||||
r0 = 6356766 # Effective earth's radius (corresponds to g0), m
|
||||
T0 = 288.15 # Standard sea-level temperature, K
|
||||
M0 = 28.9644e-3 # Mean molecular weight at sea-level, kg/mol
|
||||
|
||||
# Standard atmosphere constants for each layer
|
||||
# L corresponds to L_{M,b}
|
||||
# Heights are given in meters
|
||||
# The layers are defined as geopotential heights
|
||||
# The model is only valid up to the mesosphere
|
||||
layers = [
|
||||
{"h_base": 0, "T_base": T0, "P_base": P0, "L": -0.0065}, # Troposphere (0 - 11 km)
|
||||
{"h_base": 11000, "L": 0}, # Tropopause (11 - 20 km)
|
||||
{"h_base": 20000, "L": 0.001},# Lower Stratosphere (20 - 32 km)
|
||||
{"h_base": 32000, "L": 0.0028},# Upper Stratosphere (32 - 47 km)
|
||||
{"h_base": 47000, "L": 0}, # Stratopause (47 - 51 km)
|
||||
{"h_base": 51000, "L": -0.0028}, # Lower Mesosphere (51 - 71 km)
|
||||
{"h_base": 71000, "L": -0.002}, # Upper Mesosphere (71 - 84.8520 km)
|
||||
]
|
||||
|
||||
def molecular_scale_temp(H,H_base,T_base,L):
|
||||
return T_base + L * (H - H_base)
|
||||
|
||||
def pressure(H,H_base,T_base,P_base,L):
|
||||
if L == 0:
|
||||
return P_base * np.exp(-g0_prime*M0*(H-H_base)/(R_star*T_base))
|
||||
else:
|
||||
return P_base * (T_base/(T_base+L*(H-H_base)))**(g0_prime*M0/(R_star*L))
|
||||
|
||||
|
||||
for i in range(len(layers)):
|
||||
if i == 0:
|
||||
layers[i]["T_base"] = T0
|
||||
layers[i]["P_base"] = P0
|
||||
else:
|
||||
layers[i]["T_base"] = molecular_scale_temp(layers[i]["h_base"],layers[i-1]["h_base"],layers[i-1]["T_base"],layers[i-1]["L"])
|
||||
layers[i]["P_base"] = pressure(layers[i]["h_base"],layers[i-1]["h_base"],layers[i-1]["T_base"],layers[i-1]["P_base"],layers[i-1]["L"])
|
||||
|
||||
def temp_pressure(height):
|
||||
'''Calculates the temperature and pressure profiles in the atmosphere for given geometric heights.'''
|
||||
t, p = np.zeros(len(height)), np.zeros(len(height))
|
||||
for i,h_geom in enumerate(height):
|
||||
h = geometric_to_geopotential(h_geom)
|
||||
for j,layer in enumerate(layers):
|
||||
if j < len(layers) - 1 and h >= layer["h_base"] and h < layers[j+1]["h_base"]:
|
||||
t[i] = molecular_scale_temp(h,layer["h_base"],layer["T_base"],layer["L"])
|
||||
p[i] = pressure(h,layer["h_base"],layer["T_base"],layer["P_base"],layer["L"])
|
||||
break
|
||||
elif j == len(layers) - 1:
|
||||
t[i] = molecular_scale_temp(h,layer["h_base"],layer["T_base"],layer["L"])
|
||||
p[i] = pressure(h,layer["h_base"],layer["T_base"],layer["P_base"],layer["L"])
|
||||
return np.array(t), np.array(p)
|
||||
|
||||
def geopotential_height(P):
|
||||
"""
|
||||
Calculate the geopotential height from the given pressure using the US Standard Atmosphere 1976 model.
|
||||
|
||||
Parameters:
|
||||
P (float): Pressure at the desired height, in Pascals (Pa), takes no numpy-arrays
|
||||
|
||||
Returns:
|
||||
float: Geopotential height in meters (m)
|
||||
"""
|
||||
for i,layer in enumerate(layers):
|
||||
h_base = layer["h_base"]
|
||||
T_base = layer["T_base"]
|
||||
P_base = layer["P_base"]
|
||||
L = layer["L"]
|
||||
|
||||
# Checks for correct layer (first layer where the given pressure is smaller than the base pressure)
|
||||
# print(P,P_base,layers[i+1]["P_base"])
|
||||
if i < len(layers) - 1 and P <= P_base and P > layers[i+1]["P_base"]:
|
||||
if L == 0: # Isothermal layer
|
||||
H = R_star * T_base / g0_prime / M0
|
||||
H_diff = H * math.log(P_base/P)
|
||||
else: # Non-isothermal layer
|
||||
exponential = (P / P_base) ** (- R_star * L / (g0_prime * M0))
|
||||
H_diff = (T_base * (exponential - 1)) / L
|
||||
break
|
||||
elif i == len(layers) - 1:
|
||||
if L == 0: # Isothermal layer
|
||||
H = R_star * T_base / g0_prime / M0
|
||||
H_diff = H * math.log(P_base/P)
|
||||
else: # Non-isothermal layer
|
||||
exponential = (P / P_base) ** (- R_star * L / (g0_prime * M0))
|
||||
H_diff = (T_base * (exponential - 1)) / L
|
||||
|
||||
# print(h_base,H_diff)
|
||||
geopotential_height = h_base + H_diff
|
||||
return geopotential_height
|
||||
|
||||
def geometric_height(geopotential_height):
|
||||
"Calculates geometric height from geopotential height"
|
||||
return r0 * geopotential_height / (g0/g0_prime * r0 - geopotential_height)
|
||||
|
||||
def geometric_to_geopotential(geometric_height):
|
||||
"Calculates geopotential height from geometric_height"
|
||||
return g0/g0_prime * (r0 * geometric_height / (r0 + geometric_height))
|
||||
|
||||
# Example usage:
|
||||
# pressure = 1500 # Pressure in Pa (for example, at some altitude in the upper mesosphere)
|
||||
# geo_height = geopotential_height(pressure)
|
||||
# height = geometric_height(geo_height)
|
||||
# print(f"Geopotential height: {geo_height} meters")
|
||||
# print(f"Geometric height: {height} meters")
|
||||
Loading…
Add table
Add a link
Reference in a new issue