modified to suit current project

This commit is contained in:
Abyss Halley 2025-08-22 12:35:52 -04:00
parent 002fe62744
commit fd91791970

View file

@ -3,6 +3,8 @@
""" """
@author: Shawaiz Tabassum @author: Shawaiz Tabassum
This script uses code from YOUR. This script uses code from YOUR.
Modified by Abyss Halley.
""" """
import h5py import h5py
import numpy as np import numpy as np
@ -16,15 +18,20 @@ import pandas as pd
from scipy.signal import find_peaks from scipy.signal import find_peaks
import smplotlib #plot style, not needed import smplotlib #plot style, not needed
# from pygdsm import GlobalSkyModel from pygdsm import GlobalSkyModel
# from astropy.coordinates import SkyCoord from astropy.coordinates import SkyCoord
# from astropy.coordinates import ICRS, Galactic, FK4, FK5 from astropy.coordinates import ICRS, Galactic, FK4, FK5
# from astropy.coordinates import Angle, Latitude, Longitude from astropy.coordinates import Angle, Latitude, Longitude
# import astropy.units as u import astropy.units as u
pulsar_dir = "/home/st3105/greenburst/J1944+0907/data_2025-02-27_17-41-33_composite/" #AH: Modifying this to support multiple directories
pulsar_period = 5.185201908798642 #period of the pulsar. main_dir = os.path.join("ldata","trunk")
pulsar_dir_list = [os.path.join(main_dir, name, "") for name in
["data_2025-04-24_07-36-04", "data_2025-04-29_07-50-16", "data_2025-04-30_07-53-07", "data_2025-04-30_08-18-17", "data_2025-05-01_07-47-34"]]
#AH: old variable: pulsar_dir
pulsar_period = 1.34789947290 * 1000 #period of the pulsar. AH: in ms
#AH: Leaving these as default for now.
flux_threshold = 35 #minimum flux in mJy flux_threshold = 35 #minimum flux in mJy
distance_peaks = 6 #Number of samples/bins each pulse needs to be away from another pulse to be considered a detection. distance_peaks = 6 #Number of samples/bins each pulse needs to be away from another pulse to be considered a detection.
time_threshold = 0.1 #error interval for matching pulses only at pulsar spin period apart. Set to 1 if don't want to apply this filter time_threshold = 0.1 #error interval for matching pulses only at pulsar spin period apart. Set to 1 if don't want to apply this filter
@ -32,21 +39,22 @@ time_threshold = 0.1 #error interval for matching pulses only at pulsar spin per
#Set options #Set options
detrend_ft = True detrend_ft = True
mad_filter = True mad_filter = True
save_path = "/home/st3105/greenburst/J1944+0907/data_2025-02-27_17-41-33_composite/calib-pngs/take3/" save_path_list = [os.path.join(name, "calib-pngs", "") for name in pulsar_dir_list]
#AH: old variable: save_path
#%% #%%
#Need to get radio sky temparature at the location of the pulsar. Get it using pygdsm. #Need to get radio sky temparature at the location of the pulsar. Get it using pygdsm.
t_sys = 18 #Fig 3 of the GBT document: https://greenbankobservatory.org/wp-content/uploads/2017/07/GBTpg.pdf) t_sys = 18 #Fig 3 of the GBT document: https://greenbankobservatory.org/wp-content/uploads/2017/07/GBTpg.pdf)
# freq_mhz = 1400 freq_mhz = 1400
# sky_coord = SkyCoord(47.160175 * u.deg, -7.357055 * u.deg, frame="galactic") #For J1944+0907 sky_coord = SkyCoord(339.492 * u.deg, 0.406 * u.deg, frame="galactic") #AH: For J1643-4522
# gsm = GlobalSkyModel() gsm = GlobalSkyModel()
# t_sky = gsm.get_sky_temperature(sky_coord, freq_mhz) t_sky = gsm.get_sky_temperature(sky_coord, freq_mhz)
t_sky = 4.9805965 #t_sky = 4.9805965
t = t_sys + t_sky t = t_sys + t_sky
gain = 2.0 #same doc as above gain = 2.0 #same doc as above
n_p = 2.0 #number of polarizations n_p = 2.0 #number of polarizations
bw = 800*10**6 #Hz bw = 960*10**6 #Hz
def radiometer(T, G, np, bandwidth, W, sigma): def radiometer(T, G, np, bandwidth, W, sigma):
""" """
@ -317,10 +325,11 @@ def compute_snr_from_h5(h5_file, detrend_ft=True, mad_filter=True, save=False, s
#%% #%%
h5_dir = pulsar_dir + "h5/" for i, pulsar_dir in enumerate(pulsar_dir_list):
save_path = save_path_list[i] #AH: link directory and save path in a bit of a messy way oops
h5_dir = os.join(pulsar_dir,"cands","")
h5_files = sorted(glob.glob(f"{pulsar_dir}*.h5")) h5_files = sorted(glob.glob(f"{pulsar_dir}*.h5"))
#%%
results = [] results = []
for filepath in h5_files: for filepath in h5_files: