ADD new files from Kiruna

CHAOSfish.py

@@ -0,0 +1,86 @@
+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()
+

CHAOSx.py

@@ -5,6 +5,7 @@ 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")
@@ -14,12 +15,16 @@ 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")
 
@@ -67,7 +72,7 @@ mV = 14000
 minV = -100
 maxV = 3500
 resV = 0.838214/2
-if args.tdhist:
+if args.tdhist or args.fish:
     resV=4*resV
     minV = 0
     maxV = 2500
@@ -90,6 +95,43 @@ 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):
 
@@ -141,9 +183,7 @@ 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
-                    
-                    T = adc_to_T(int(l[8]))
-                elif l[0] == 'EI' and len(l)>18 and conditions(l,args) and consider:  
+                if l[0] == 'EI' and len(l)>3*18 and conditions(l,args) and consider:  
                     for i, chan in enumerate(chans):
                         b = CmV(l,ch[chan])*u[chan] 
                         if minV <= b <= maxV:
@@ -153,27 +193,142 @@ 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 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):
+# 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
+    elif 'd' in expr:
+        return expr.split('d'), 'div', expr.split('d')[1], None
+    else:
+        return [expr], 'single', None, None
+
+
+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:
         value = CmV(l, ch[name.index(channames[0])]) * u[name.index(channames[0])]
-    if len(channames) > 1:
-        if op == 'add':
+        if len(channames) > 1 and 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])]
+
+    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
+
     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 = 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])
+
 
     nameadd='x'
     energy = 'keV'
@@ -214,8 +369,10 @@ 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)
+                    #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)
                 except:
                     bx = maxV; by = maxV
 
@@ -225,7 +382,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."""

HK_flight.py

@@ -11,6 +11,7 @@ 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"
 
@@ -38,13 +39,14 @@ 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 = 25
+fs = 30
 
 def main():
     # Plot all graphs
-    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")
+    #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")
     
     # Channels grouped by type
     allchans = ["A1H", "A1L", "A2H", "A2L", "B", "C1H", "C1L", "C2H", "C2L", "D1H", "D1L", "D2H", "D2L", "EH", "EL", "E1", "E2", "E3"]
@@ -52,21 +54,24 @@ def main():
     AC = ["A1H", "A2H", "C1H", "C2H"]
     E = ["EH", "E1", "E2", "E3"]
 
-    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(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_six()
+    #plot_and_layout(plot_pressures_and_height, "","recap_HK_p")
 
-    #plt.show()
+    #plot_six()
+
+    plt.show()
 
 
 def plot_and_layout(plot_func, title, savefile):
-    fig, ax = plt.subplots(figsize=(20, 15))
+    fig, ax = plt.subplots(figsize=(20, 13))
     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()
 
@@ -80,7 +85,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 [°C]", fontsize = fs)
+    ax.set_ylabel("T in °C", fontsize = fs)
 
 def plot_T_and_p(ax):
     ax2 = ax.twinx()
@@ -90,13 +95,35 @@ 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 [°C]", fontsize = fs)
-    ax2.set_ylabel("p [mbar]", fontsize = fs)
+    ax.set_ylabel("T in °C", fontsize = fs)
+    ax2.set_ylabel("p in 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 [mbar]", fontsize = fs)
+    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()
 
 def plot_I_and_V(ax):
     ax2 = ax.twinx()
@@ -108,15 +135,16 @@ 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)
@@ -132,10 +160,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.95, event, rotation=90, verticalalignment='top', horizontalalignment='right', fontsize=fs - 2)
+        ax.text(time, ax.get_ylim()[1] * 0.98, 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'))

histplot2d.py

@@ -39,7 +39,11 @@ def main():
 
    
     file = args.filepath.split('.')[0]
+    try:
         y,x = file.split('~')[-1].split('-')
+    except:
+        x='unknown'
+        y='unknown'
     
     end = '.2dhist'
     energy = 'keV'
@@ -80,7 +84,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)
+    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)
@@ -89,13 +93,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] = 500
+    # 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 = y + ' in ' + energy+'$_{Si}$'
     xlab = x + ' in ' + energy+'$_{Si}$'
@@ -111,8 +115,11 @@ def main():
 
     plt.ylabel(ylab, fontsize = fs+5)
     plt.xlabel(xlab, fontsize=fs+5)
-    plt.ylim(lim[0],lim[1])
-    plt.xlim(lim[0],lim[1])
+    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)

plot_deltaT.py

@@ -0,0 +1,88 @@
+#!/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()    
+

plotfish.py

@@ -0,0 +1,114 @@
+'''
+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()    

tryfish.py

@@ -0,0 +1,84 @@
+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()

us_standard_atmosphere.py

@@ -0,0 +1,116 @@
+###########################################################
+### 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")