diff --git a/singlepulse.py b/singlepulse.py index 68cc155..3df4840 100644 --- a/singlepulse.py +++ b/singlepulse.py @@ -3,6 +3,8 @@ """ @author: Shawaiz Tabassum This script uses code from YOUR. + +Modified by Abyss Halley. """ import h5py import numpy as np @@ -16,15 +18,20 @@ import pandas as pd from scipy.signal import find_peaks import smplotlib #plot style, not needed -# from pygdsm import GlobalSkyModel -# from astropy.coordinates import SkyCoord -# from astropy.coordinates import ICRS, Galactic, FK4, FK5 -# from astropy.coordinates import Angle, Latitude, Longitude -# import astropy.units as u +from pygdsm import GlobalSkyModel +from astropy.coordinates import SkyCoord +from astropy.coordinates import ICRS, Galactic, FK4, FK5 +from astropy.coordinates import Angle, Latitude, Longitude +import astropy.units as u -pulsar_dir = "/home/st3105/greenburst/J1944+0907/data_2025-02-27_17-41-33_composite/" -pulsar_period = 5.185201908798642 #period of the pulsar. +#AH: Modifying this to support multiple directories +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 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 @@ -32,21 +39,22 @@ time_threshold = 0.1 #error interval for matching pulses only at pulsar spin per #Set options detrend_ft = 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. t_sys = 18 #Fig 3 of the GBT document: https://greenbankobservatory.org/wp-content/uploads/2017/07/GBTpg.pdf) -# freq_mhz = 1400 -# sky_coord = SkyCoord(47.160175 * u.deg, -7.357055 * u.deg, frame="galactic") #For J1944+0907 -# gsm = GlobalSkyModel() -# t_sky = gsm.get_sky_temperature(sky_coord, freq_mhz) -t_sky = 4.9805965 +freq_mhz = 1400 +sky_coord = SkyCoord(339.492 * u.deg, 0.406 * u.deg, frame="galactic") #AH: For J1643-4522 +gsm = GlobalSkyModel() +t_sky = gsm.get_sky_temperature(sky_coord, freq_mhz) +#t_sky = 4.9805965 t = t_sys + t_sky gain = 2.0 #same doc as above 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): """ @@ -317,36 +325,37 @@ def compute_snr_from_h5(h5_file, detrend_ft=True, mad_filter=True, save=False, s #%% -h5_dir = pulsar_dir + "h5/" -h5_files = sorted(glob.glob(f"{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")) -#%% -results = [] + results = [] -for filepath in h5_files: - try: - fluxes_mjy, flux_stds, widths_t, peak_mjds, time_peaks, time_pulse, dm = compute_snr_from_h5(filepath, - detrend_ft=True, - mad_filter=True, - save=True, - show=False, - save_path=save_path) - - for i in range(0,(len(fluxes_mjy))): - results.append({ - "mjd": peak_mjds[i], - "flux_std": flux_stds[i], - "flux_mJy": fluxes_mjy[i], - "width_s": widths_t[i], - "Pulse time (ms)": time_peaks[i], - "tcand": time_pulse, - "dm": dm, - "filename": os.path.basename(filepath) - }) - except Exception as e: - print(f"Error processing {filepath}: {e}") + for filepath in h5_files: + try: + fluxes_mjy, flux_stds, widths_t, peak_mjds, time_peaks, time_pulse, dm = compute_snr_from_h5(filepath, + detrend_ft=True, + mad_filter=True, + save=True, + show=False, + save_path=save_path) + + for i in range(0,(len(fluxes_mjy))): + results.append({ + "mjd": peak_mjds[i], + "flux_std": flux_stds[i], + "flux_mJy": fluxes_mjy[i], + "width_s": widths_t[i], + "Pulse time (ms)": time_peaks[i], + "tcand": time_pulse, + "dm": dm, + "filename": os.path.basename(filepath) + }) + except Exception as e: + print(f"Error processing {filepath}: {e}") -df = pd.DataFrame(results) -#%% -df['pulse_energy_jyus'] = (df['flux_mJy']/1000) * (df['width_s']*10**6) -df.to_csv(save_path+"flux.csv", index=False) + df = pd.DataFrame(results) + #%% + df['pulse_energy_jyus'] = (df['flux_mJy']/1000) * (df['width_s']*10**6) + df.to_csv(save_path+"flux.csv", index=False)