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()
-

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):
 
@@ -178,12 +136,14 @@ def generate_histogram(u):
                 consider=False
                 if l[0] == 'H':
                     timestamp = int(l[1])
-                if not args.float:
+                    if not args.float:
-                    consider=True
+                        consider=True
-                elif 1727859780 < timestamp < 1727872980:  # between 2024-10-02 11:23:00 and 2024-10-02 15:03:00 (Sweden time)
+                    elif 1727859780 < timestamp < 1727872980:  # between 2024-10-02 11:23:00 and 2024-10-02 15:03:00 (Sweden time)
-                    consider=True
+                        consider=True
-                else: consider = False
+                    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
+    if '+' in expr:
-
+        return expr.split('+'), 'add'
-    # 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
+        return expr.split('d'), 'div'
     else:
-        return [expr], 'single', None, None
+        return [expr], 'single'  
 
-
+def get_channel_value(l, channames,op):
-def get_channel_value(l, channames, op, next_op=None, next_channel=None):
+    value = CmV(l, ch[name.index(channames[0])]) * u[name.index(channames[0])]
-    value = None 
+    if len(channames) > 1:
-    if op == 'max':
+        if op == 'add':
-        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 and op == 'add':
             value += CmV(l, ch[name.index(channames[1])]) * u[name.index(channames[1])]
-
+        elif op == 'div':
-    if next_op == 'div' and next_channel:
+            value /= CmV(l, ch[name.index(channames[1])]) * u[name.index(channames[1])]
-        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])
+    x_expr, x_op = parse_channel_expression(chans[1])
-    #y_expr, y_op = parse_channel_expression(chans[0])
+    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'
@@ -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)
+                    bx = get_channel_value(l,x_expr,x_op)
-                    #by = get_channel_value(l,y_expr,y_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
 
@@ -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."""

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, "All Temperatures","HK_Ts.pdf")
-    plot_and_layout(plot_allT, "","recap_HK_Ts")
+    plot_and_layout(plot_pressures, "Pressures","HK_ps.pdf")
-    plot_and_layout(plot_pressures, "Pressures","HK_ps")
+    plot_and_layout(plot_I_and_V, "Current and Voltage","HK_VI.pdf")
-    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"]
@@ -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, E), "Counts channels E","plots/HK_countsE.pdf")
-    plot_and_layout(lambda ax: plot_c64(ax, AC), "","recap_HK_countsAC")
+    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), "","recap_HK_countsDB")
+    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)
+    ax.set_ylabel("T [°C]", fontsize = fs)
-    ax2.set_ylabel("p in mbar", 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)
+    ax.set_ylabel("p [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()
@@ -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'))

histplot2d.py

@@ -39,11 +39,7 @@ def main():
 
    
     file = args.filepath.split('.')[0]
-    try:
+    y,x = file.split('~')[-1].split('-')
-        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:
+    if not args.limits == None:
-    #     lim[0] = float(args.limits[0])
+        lim[0] = float(args.limits[0])
-    #     lim[1] = float(args.limits[1])    
+        lim[1] = float(args.limits[1])    
-    # if lim[0] == None:
+    if lim[0] == None:
-    #     lim[0] = 0
+        lim[0] = 0
-    # if lim[1] == None:
+    if lim[1] == None:
-    #     lim[1] = 1000
+        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.ylim(lim[0],lim[1])
-    plt.xlim(1,10000)
+    plt.xlim(lim[0],lim[1])
-
-    plt.yscale('log')
-    plt.xscale('log')
     
     cbar.set_label('counts', fontsize = fs+5)
     cbar.ax.tick_params(labelsize=fs)

plot_deltaT.py

@@ -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()    
-

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()    

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()

us_standard_atmosphere.py

@@ -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")