Compare commits
7 commits
fff11e4648
...
58d5188c11
| Author | SHA1 | Date | |
|---|---|---|---|
| 58d5188c11 | |||
| ae4e8ec66a | |||
| 129f3cb7e7 | |||
| 1626446c72 | |||
| f334161b70 | |||
| cd9ced7b95 | |||
|
|
0da36eb56e |
2 changed files with 76 additions and 27 deletions
|
|
@ -27,7 +27,7 @@ fp_out = args.outfile
|
|||
current_time = timeit.default_timer()
|
||||
filename, rest = fp.split(".")
|
||||
#filename = input("Filename (without extension): ")
|
||||
file = open(filename+'.AHA')
|
||||
file = open(filename+'.AHA', encoding='utf-8', errors='ignore')
|
||||
content = file.readlines()
|
||||
file.close()
|
||||
|
||||
|
|
@ -75,8 +75,11 @@ for i in range(len(EDB_Array)):
|
|||
# EDBnpf_array[i][j-7] = float(parts[j+1])
|
||||
#else:
|
||||
# EDBnp_array[i][j] = int(parts[j+1])
|
||||
try:
|
||||
if not(int(parts[3])>= 24) and float(parts[5])/10000 > min_value and float(parts[5])/10000 < max_value:
|
||||
len_array[int(parts[2])][int(parts[3])] += 1
|
||||
except:
|
||||
pass
|
||||
#EDB_File.close()
|
||||
print("Element counting completed", timeit.default_timer()-current_time)
|
||||
current_time = timeit.default_timer()
|
||||
|
|
|
|||
86
seth_hk.py
86
seth_hk.py
|
|
@ -23,6 +23,8 @@ fp = args.file
|
|||
plot_flag = args.plot
|
||||
live = args.live
|
||||
|
||||
|
||||
|
||||
global time
|
||||
global data
|
||||
time = []
|
||||
|
|
@ -64,8 +66,12 @@ def counter(l):
|
|||
lower_count = (math.fsum(sl0)-sl0[3]-sl0[11]-sl0[20]-sl0[11])/20
|
||||
if lower_count <= 0:
|
||||
lower_count = 'nan'
|
||||
if upper_count == 'nan' or lower_count == 'nan':
|
||||
if upper_count == 'nan' and lower_count == 'nan':
|
||||
return [upper_count, BGO_upper/12, HETA/12, HETB/12, BGO_lower/12, lower_count] # divide by 12 to get counts per second
|
||||
elif upper_count == 'nan':
|
||||
return [upper_count, BGO_upper/12, HETA/12, HETB/12, BGO_lower/12, lower_count/12] # divide by 12 to get counts per second
|
||||
elif lower_count == 'nan':
|
||||
return [upper_count/12, BGO_upper/12, HETA/12, HETB/12, BGO_lower/12, lower_count] # divide by 12 to get counts per second
|
||||
else:
|
||||
return [upper_count/12, BGO_upper/12, HETA/12, HETB/12, BGO_lower/12, lower_count/12] # divide by 12 to get counts per second
|
||||
|
||||
|
|
@ -136,14 +142,6 @@ def calibrate_vectors(acc_raw, mag_raw):
|
|||
"""
|
||||
mag_sens = 6842 # LSB/gauss
|
||||
acc_sens = 16000 # 1000 LSB/g and it is only a 12 bit number so divide by 16 (1hex bit) equivalent to shifting 4 bits
|
||||
mag_raw[0] = mag_raw[0] + 700 #Values for 360-29e (guessed)
|
||||
mag_raw[1] = mag_raw[1] +4000 # For tvw-36e use -500 and +2500
|
||||
if mag_raw[0] < 0:
|
||||
mag_raw[0] = mag_raw[0]+65535
|
||||
if mag_raw[0] > 65535:
|
||||
mag_raw[0] = mag_raw[0]-65535
|
||||
if mag_raw[1] > 65535:
|
||||
mag_raw[1] = mag_raw[1]-65535
|
||||
mag_calib = [1.02305, 1.13811, 1.16184, 0.0265714, -0.188856, -0.177294] # fit coefficients from BA thesis
|
||||
acc_calib = [0.972148, 0.968608, 0.974061, -0.00109779, -0.0154984, 0.0135722]
|
||||
acc_raw = [struct.unpack('<h', struct.pack('<H', int(num)))[0] for num in acc_raw] # turns int into hex, turns hex into signed int
|
||||
|
|
@ -178,16 +176,14 @@ def calculate_heading(mag, phi, theta):
|
|||
thus, the heading is simply the angle between x and north. Translate this to compass.
|
||||
"""
|
||||
declination = 0 # Kiruna
|
||||
#phi = 0
|
||||
theta = 0 # Set zero because it works...
|
||||
mag_world = np.matmul(Rz(-phi), np.matmul(Ry(theta-np.pi), np.matmul(Rz(phi), mag)))
|
||||
mag_world_psi = np.arctan2(mag_world[1], mag_world[0]) + declination # angle mag vector in wf, x axis in wf + delta because TN is delta left of the MN (=mag vector)
|
||||
if mag_world_psi <= 0:
|
||||
mag_world_psi = 2*np.pi + mag_world_psi
|
||||
#mag_world_psi = abs(mag_world_psi) # because compass counts clock-wise and atan2 counter CW
|
||||
#mag_world_psi = 2*np.pi + mag_world_psi
|
||||
mag_world_psi = abs(mag_world_psi) # because compass counts clock-wise and atan2 counter CW
|
||||
else:
|
||||
pass
|
||||
#mag_world_psi = 2*np.pi - mag_world_psi
|
||||
#pass
|
||||
mag_world_psi = 2*np.pi - mag_world_psi
|
||||
heading = mag_world_psi * 180/np.pi # heading in degrees (0deg = True North)
|
||||
return heading
|
||||
|
||||
|
|
@ -206,14 +202,20 @@ def altitude(pressure):
|
|||
'''
|
||||
Uses the barometric formula to solve for the height
|
||||
'''
|
||||
if type(pressure) is int or type(pressure) is float:
|
||||
H = 8400 # scale height for T=288K
|
||||
return -math.log(pressure/1013)*H
|
||||
else:
|
||||
return 'nan'
|
||||
|
||||
def main(input):
|
||||
i = -1 # HK comes in 12s intervals i ended by C64 lines: ... C64 P H P HD HD P C64 P H P HD HD P C64 ...
|
||||
global time
|
||||
global data
|
||||
global live
|
||||
time_wrong_begin = False
|
||||
time_wrong_middle = False
|
||||
time_changed = False
|
||||
dorn_temps = [] # [Tadc0, Tpa0, Tbgo01, Tbgo02, Tpwr, Tadc1, Tpa1, Tbgo11, Tbgo12, Text] slice 0 and slice1
|
||||
temp = []
|
||||
pres = []
|
||||
|
|
@ -277,18 +279,49 @@ def main(input):
|
|||
|
||||
if l[:2] == 'H ': # One every 12s
|
||||
ts = int(l.split()[1])
|
||||
if fp and fp.split('/')[-1] == "2025-08-14-longtime-23e.AHA" and ts > 1800000000: # timestamps wrong in this file
|
||||
ts -= 1963615744
|
||||
if i < 2000 and ts >= 3718824668-1963615744:
|
||||
ts -= 12*3600
|
||||
#if fp and fp.split('/')[-1] == "2025-08-14-longtime-23e.AHA" and ts > 1800000000: # timestamps wrong in this file
|
||||
# ts -= 1963615744
|
||||
# if i < 2000 and ts >= 3718824668-1963615744:
|
||||
# ts -= 12*3600
|
||||
if fp and fp.split('/')[-1] == '2025-08-08-seth-17e.AHA' and ts > 1800000000:
|
||||
ts -= 2138848512
|
||||
if i < 500 and ts > 3893525631-2138848512:
|
||||
ts -= 16404
|
||||
|
||||
#if live and (ts > 1800000000 or ts < 1700000000) and not time_wrong_begin and time == []:
|
||||
# time_wrong_begin = True
|
||||
#if live and (ts < 1800000000 and ts > 1700000000) and time_wrong_begin:
|
||||
# for i in range(len(time)):
|
||||
# time[i] += ts - time[-1] -12
|
||||
# time_wrong_begin = False
|
||||
# time_changed = True
|
||||
# if live and (ts > 1800000000 or ts < 1700000000) and not time_wrong_middle and not time == []:
|
||||
# time_wrong_middle == len(time)
|
||||
# if live and (ts < 1800000000 and ts > 1700000000) and time_wrong_middle:
|
||||
# for i in range(i, len(time)):
|
||||
# time[i] += ts - time[-1] -12
|
||||
# time_wrong_middle = False
|
||||
# time_changed = True
|
||||
#if live and len(time) and abs(ts - time[-1]) > 120:
|
||||
#print(ts-time[-1])
|
||||
#print(time, "time (pre)")
|
||||
#print(len(time))
|
||||
#print(list(data), "data")
|
||||
#print(len(time)==len(list(data)), len(time), len(data))
|
||||
# for i in range(len(time)-1):
|
||||
# data[time[i] + ts - time[-1] -12] = data.pop(time[i])
|
||||
# time[i] += ts - time[-1] -12
|
||||
#print(i)
|
||||
#time[i] += ts - time[-1] -12
|
||||
# time_changed = True
|
||||
#print(time, "time (post)")
|
||||
#print(len(time))
|
||||
if len(time) > i: # i is zero-indexed sometimes there are >1 H lines in one interval, use the last
|
||||
time[-1] = ts
|
||||
else:
|
||||
time.append(ts) # [time0, time1, ...] times of the intervals
|
||||
if time_changed:
|
||||
print(time) #debugging code
|
||||
|
||||
if l[:6] == 'HDORN ': # two lines, sl0 and sl1 per interval
|
||||
sl = l.split()[1]
|
||||
|
|
@ -302,8 +335,14 @@ def main(input):
|
|||
|
||||
if l[:4] == 'C64 ':
|
||||
heading, roll, pitch = calculate_ads(acc_raw, mag_raw)
|
||||
if pres:
|
||||
mean_pres = math.fsum(pres)/len(pres)
|
||||
else:
|
||||
mean_pres = 'nan'
|
||||
if temp:
|
||||
mean_temp = math.fsum(temp)/len(temp)
|
||||
else:
|
||||
mean_temp = 'nan'
|
||||
if len(dorn_temps) < 10:
|
||||
dorn_temps = ['nan' for i in range(10)]
|
||||
data[time[-1]] = [mean_pres, mean_temp, counter(l), dorn_temps, [heading, roll, pitch]]
|
||||
|
|
@ -314,6 +353,13 @@ def main(input):
|
|||
# Plotting
|
||||
if plot_flag and live: # only plot if -pl option is True
|
||||
### Date formatting
|
||||
if time_changed:
|
||||
#print(x, "x (pre)")
|
||||
#print(data)
|
||||
for i in range(len(time)-1):
|
||||
x[i] = datetime.datetime.fromtimestamp(list(data)[i], datetime.timezone.utc)
|
||||
time_changed = False
|
||||
#print(x, "x (post)")
|
||||
x.append(datetime.datetime.fromtimestamp(list(data)[-1], datetime.timezone.utc))
|
||||
fmtdate = md.DateFormatter('%d-%H:%M:%S')
|
||||
ax1.xaxis.set_major_formatter(fmtdate)
|
||||
|
|
@ -383,7 +429,7 @@ def main(input):
|
|||
|
||||
|
||||
if fp:
|
||||
with open(fp) as file:
|
||||
with open(fp, encoding='utf-8', errors='ignore') as file:
|
||||
input = file
|
||||
main(input)
|
||||
if plot_flag and not live:
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue