diff --git a/.ruff.toml b/.ruff.toml index 0ee0427..00ac998 100644 --- a/.ruff.toml +++ b/.ruff.toml @@ -5,6 +5,7 @@ exclude = [ "__pycache__", "build", "bella/version.py", + "bella/type_III_fitter/read_wi_wa_l2_60s_data.py", ] [lint] @@ -32,6 +33,8 @@ extend-ignore = [ ] "__init__.py" = ["E402", "F401", "F403"] "test_*.py" = ["B011", "D", "E402", "PGH001", "S101"] +"bella/multilaterate/bayes_positioner.py" = ["F841"] + [lint.pydocstyle] convention = "numpy" diff --git a/bella/multilaterate/bayes_positioner.py b/bella/multilaterate/bayes_positioner.py index dafb80d..6a22c81 100755 --- a/bella/multilaterate/bayes_positioner.py +++ b/bella/multilaterate/bayes_positioner.py @@ -61,7 +61,7 @@ def sample(self, toa, draws=2000, tune=2000, chains=4,cores=4, init='jitter+adap # assert correct number of observations if len(toa) != len(stations): - raise Exception("ERROR: number of observations must match number of stations! (%i, %i)"%(len(toa), len(stations))) + raise Exception(f"ERROR: number of observations must match number of stations! ({len(toa)}, {len(stations)})") # assert max toa is not larger than t_lim if np.max(toa) > t_lim: @@ -167,7 +167,7 @@ def triangulate(coords, times,t_cadence=60,v_sd=3E5, chains=4, cores=4, N_SAMPLE print(summary) - if traceplot==True: + if traceplot is True: # trace plot left = 0.04 # the left side of the subplots of the figure right = 0.972 # the right side of the subplots of the figure @@ -331,7 +331,7 @@ def triangulate(coords, times,t_cadence=60,v_sd=3E5, chains=4, cores=4, N_SAMPLE ax[3, 1].set_ylim([t_sampled-2.9*t_sd_sampled, t_sampled+2.9*t_sd_sampled]) # plt.show(block=False) - if savetraceplot == True: + if savetraceplot is True: t0 = times[0] if traceplotdir=="": direct = mkdirectory("./Traceplots/") @@ -345,7 +345,7 @@ def triangulate(coords, times,t_cadence=60,v_sd=3E5, chains=4, cores=4, N_SAMPLE print(f"saved {direct}/traceplot_{traceplotfn}") - if showplot==True: + if showplot is True: plt.show(block=False) else: plt.close() @@ -370,7 +370,7 @@ def triangulate(coords, times,t_cadence=60,v_sd=3E5, chains=4, cores=4, N_SAMPLE plt.xlabel(r"'HEE - X / $R_{\odot}$'") plt.ylabel(r"'HEE - Y / $R_{\odot}$'") # plt.savefig("bayes_positioner_result2.jpg", bbox_inches='tight', pad_inches=0.01, dpi=300) - if showplot==True: + if showplot is True: plt.show(block=False) else: plt.close() diff --git a/bella/multilaterate/bayesian_tracker.py b/bella/multilaterate/bayesian_tracker.py index 5b4c62b..9abab77 100755 --- a/bella/multilaterate/bayesian_tracker.py +++ b/bella/multilaterate/bayesian_tracker.py @@ -1,7 +1,6 @@ # Standard Library imports import os -import sys # General import logging import datetime as dt @@ -17,12 +16,15 @@ import numpy as np import pymc3 as pm import solarmap -# Local imports -from bayes_positioner import * from scipy.ndimage import median_filter from astropy.constants import R_sun, au, c +# Local imports +import sys +sys.path.insert(1, '/Users/canizares/Library/CloudStorage/OneDrive-Personal/Work/0_PhD/Projects/BELLA_Projects/TCDSolarBELLA') +from bella.multilaterate.bayes_positioner import triangulate + @contextmanager def suppress_stdout(): @@ -52,7 +54,6 @@ def parallel_pos_map(x1,y1,stations,xrange,yrange,xres,yres, cores=1, traceplots x_true = np.array([x1 * R_sun.value, y1 * R_sun.value]) # true source position (m) v_true = c.value # speed of light (m/s) - t0_true = 100 # source time. can be any constant, as long as it is within the uniform distribution prior on t0 d_true = np.linalg.norm(stations - x_true, axis=1) t1_true = d_true / v_true # true time of flight values # t_obs = t1_true-t0_true# true time difference of arrival values @@ -72,7 +73,7 @@ def parallel_pos_map(x1,y1,stations,xrange,yrange,xres,yres, cores=1, traceplots # print(f"t1_pred: {t1_pred}") # print(f"stations: {stations}") - if traceplotsave == True: + if traceplotsave is True: traceplotpath = f"./traceplots/{date_str}/traceplot_{xrange[0]}_{xrange[1]}_{yrange[0]}_{yrange[1]}_{xres}_{yres}" mkdirectory(traceplotpath) @@ -96,12 +97,11 @@ def parallel_pos_map(x1,y1,stations,xrange,yrange,xres,yres, cores=1, traceplots res = np.array([delta_obs, np.nan, np.nan]) tloop1 = dt.datetime.now() - except: - delta_obs = 200 - print(f"SIM FAILED at P({x1},{y1}), SKIPPED") - res = np.array([delta_obs, x1, y1]) + + except Exception as e: + print(f"Error at P({x1},{y1}): {str(e)}") + res = np.array([600, x1, y1]) # Fallback error handling tloop1 = dt.datetime.now() - pass tloopinseconds = (tloop1 - tloop0).total_seconds() print(f"Time Loop : {tloop1 - tloop0} : {tloopinseconds}s ") @@ -128,8 +128,9 @@ def parallel_tracker(freq,t_obs, stations ): res = np.array([xy, np.nan, np.nan]) tloop1 = dt.datetime.now() - except: + except Exception as e: xy = 0 + print(f"ERROR: {e}") print(f"SIM FAILED at Freq = {freq} MHz, SKIPPED") res = np.array([xy, freq]) tloop1 = dt.datetime.now() @@ -178,12 +179,12 @@ def plot_map_simple(delta_obs, xmapaxis, ymapaxis, stations,vmin=0,vmax=30, save ax.set_ylabel(r"'HEE - Y / $R_{\odot}$'", fontsize=22) ax.set_title(title, fontsize=22) - if savefigure == True: + if savefigure is True: figdir = f"{figdir}/{date_str}" mkdirectory(figdir) plt.savefig(figdir+f"/bayes_positioner_map_{xmapaxis[0]}_{xmapaxis[-1]}_{ymapaxis[0]}_{ymapaxis[-1]}_{xres}_{yres}_{N_STATIONS}.jpg", bbox_inches='tight', pad_inches=0.01, dpi=300) - if showfigure == True: + if showfigure is True: plt.show(block=False) else: plt.close(fig) @@ -208,7 +209,7 @@ def savetrackedtypeiii(results, dir="./Data/", date_str="date", profile="PROFILE title = f"TRACKING_{date_str}_results_{N_STATIONS}stations_{profile}.pkl" - direct = mkdirectory(dir) + mkdirectory(dir) stringpath = dir+f'{title}' with open(stringpath, 'wb') as outp: diff --git a/bella/multilaterate/bella_plotter.py b/bella/multilaterate/bella_plotter.py index eebcdb5..3e17473 100755 --- a/bella/multilaterate/bella_plotter.py +++ b/bella/multilaterate/bella_plotter.py @@ -1,5 +1,4 @@ import os -import sys import math import logging import datetime as dt @@ -12,7 +11,6 @@ import numpy as np import scipy.io import solarmap -from bayes_positioner import * from matplotlib import cm from matplotlib.ticker import FormatStrFormatter, LogFormatter # import pymc3 as pm @@ -22,6 +20,10 @@ # import arviz as az from astropy.constants import R_sun, au +import sys +# sys.path.insert(1, '/Users/canizares/Library/CloudStorage/OneDrive-Personal/Work/0_PhD/Projects/BELLA_Projects/TCDSolarBELLA') +# from bella.multilaterate.bayes_positioner import * + plt.rcParams.update({'font.size': 18}) plt.rcParams["font.family"] = "Times New Roman" @@ -81,9 +83,9 @@ def plot_map_simple_withTracked(delta_obs, trackedtypeIII, xmapaxis, ymapaxis, s sd = sd_vals/R_sun.value - xres = xmapaxis[1]-xmapaxis[0] - yres = ymapaxis[1]-ymapaxis[0] - N_STATIONS = len(stations) + xmapaxis[1]-xmapaxis[0] + ymapaxis[1]-ymapaxis[0] + len(stations) # PARKER SPIRAL @@ -128,7 +130,7 @@ def plot_map_simple_withTracked(delta_obs, trackedtypeIII, xmapaxis, ymapaxis, s i+=1 - if parker_spiral == True: + if parker_spiral is True: ax.plot(x_parker,y_parker,"k--") ax.plot(x_parker_plus,y_parker_plus,"k--", markersize=0.5) ax.plot(x_parker_minus,y_parker_minus,"k--") @@ -144,7 +146,7 @@ def plot_map_simple_withTracked(delta_obs, trackedtypeIII, xmapaxis, ymapaxis, s # # ELLIPSES colors = cm.turbo(list(np.linspace(0,1.0,len(xy)))) - if confidence ==True: + if confidence is True: i = 0 ell_track_uncertainty = matplotlib.patches.Ellipse(xy=(xy[i, 0], xy[i, 1]), width=2 * sd[i, 0], height=2 * sd[i, 1], @@ -204,12 +206,12 @@ def plot_map_simple_withTracked(delta_obs, trackedtypeIII, xmapaxis, ymapaxis, s ax.set_ylabel(r"HEE - Y (R$_\odot$)", fontsize=22) ax.set_title(title, fontsize=22) - if savefigure == True: + if savefigure is True: figdir = f"{figdir}/{date_str}" mkdirectory(figdir) plt.savefig(figdir+filename, bbox_inches='tight', pad_inches=0.01, dpi=300) - if showfigure == True: + if showfigure is True: plt.show(block=False) else: plt.close(fig) @@ -225,9 +227,9 @@ def plot_bella_map(fig,ax,delta_obs, xmapaxis, ymapaxis, stations, linecolor = "k"): - xres = xmapaxis[1]-xmapaxis[0] - yres = ymapaxis[1]-ymapaxis[0] - N_STATIONS = len(stations) + xmapaxis[1]-xmapaxis[0] + ymapaxis[1]-ymapaxis[0] + len(stations) # fig, ax = plt.subplots(1,1,figsize=(11,11)) # plt.subplots_adjust(top=1, bottom=0) @@ -258,7 +260,7 @@ def plot_bella_map(fig,ax,delta_obs, xmapaxis, ymapaxis, stations, ax.set_ylim(ymapaxis[0], ymapaxis[-1]) # fig.subplots_adjust(right=0.9) - if showcolorbar == True: + if showcolorbar is True: cbar_ax = fig.add_axes(cbar_axes) # [left, bottom, width, height] fig.colorbar(im_0, cax=cbar_ax) cbar_ax.set_ylabel(r'BELLA uncertainty (R$_\odot$)', fontsize=18) @@ -269,7 +271,7 @@ def plot_bella_map(fig,ax,delta_obs, xmapaxis, ymapaxis, stations, if "sun" in objects: ax.plot(0, 0, 'yo', label="Sun", markersize=10, markeredgecolor ='k') - if showlegend == True: + if showlegend is True: ax.legend(loc=1) ax.set_xlabel(r"HEE - X (R$_\odot$)", fontsize=20) @@ -299,7 +301,7 @@ def plot_tracked_typeIII(fig, ax, trackedtypeIII, confidence=False, showcolorbar # colors = cm.turbo(list(np.linspace(0,1.0,len(xy)))) - if confidence ==True: + if confidence is True: i = 0 ell_track_uncertainty = matplotlib.patches.Ellipse(xy=(xy[i, 0], xy[i, 1]), width=2*sd[i, 0], height=2*sd[i, 1], @@ -316,7 +318,7 @@ def plot_tracked_typeIII(fig, ax, trackedtypeIII, confidence=False, showcolorbar im_track = ax.scatter(xy[:,0], xy[:,1],c = tracked_freqs, cmap=cmap, marker=marker, edgecolors=edgecolors, s=s, label=label, norm=norm,zorder=zorder) - if showcolorbar == True: + if showcolorbar is True: if cbar_sep_ax: cbar_ax2 = fig.add_axes(cbar_axes) # [left, bottom, width, height] formatter = LogFormatter(10, labelOnlyBase=False) @@ -354,7 +356,7 @@ def plot_typeIII_sim(fig, ax, trackedtypeIII, confidence=False, showcolorbar=Tru colors = cmap_func(np.linspace(0, 1.0, len(xy))) # colors = cm.turbo(list(np.linspace(0,1.0,len(xy)))) - if confidence ==True: + if confidence is True: i = 0 ell_track_uncertainty = matplotlib.patches.Ellipse(xy=(xy[i, 0], xy[i, 1]), width=2*sd[i, 0], height=2*sd[i, 1], @@ -373,10 +375,10 @@ def plot_typeIII_sim(fig, ax, trackedtypeIII, confidence=False, showcolorbar=Tru ax.add_patch(ell_track_uncertainty) if showtruesources: - im_true = ax.scatter(xy_true[:,0], xy_true[:,1], cmap=cmap, marker="o", edgecolors=edgecolors, c="yellow", s=s, label="True Sources", norm=norm,zorder=zorder-1) + ax.scatter(xy_true[:,0], xy_true[:,1], cmap=cmap, marker="o", edgecolors=edgecolors, c="yellow", s=s, label="True Sources", norm=norm,zorder=zorder-1) im_track = ax.scatter(xy[:,0], xy[:,1], cmap=cmap, marker=marker, edgecolors=edgecolors, s=s, label=label, norm=norm,zorder=zorder) - if showcolorbar == True: + if showcolorbar is True: cbar_ax2 = fig.add_axes(cbar_axes) # [left, bottom, width, height] formatter = LogFormatter(10, labelOnlyBase=False) fig.colorbar(im_track, cax=cbar_ax2, format=formatter) @@ -720,7 +722,7 @@ def __init__(self, xy_true, xy_detected, sd_detected): ax.tick_params(axis='both', which='major', labelsize=18) plt.show(block = False) - if savefigs == True: + if savefigs is True: dir = mkdirectory("./Figures/") if trackedfile[-11:-4] != 'SCATTER': plt.savefig(dir+'BELLA_map0.png', dpi=300) @@ -817,7 +819,7 @@ def __init__(self, xy_true, xy_detected, sd_detected): # plt.show(block = False) # # savefigs = True - # if savefigs == True: + # if savefigs is True: # plt.savefig('/Users/canizares/OneDrive/Work/0_PhD/Projects/2012_06_07_WSTASTB/BELLA_map_mosaic.png', dpi=300) # # diff --git a/bella/multilaterate/density_models.py b/bella/multilaterate/density_models.py index cde4584..a646c9e 100755 --- a/bella/multilaterate/density_models.py +++ b/bella/multilaterate/density_models.py @@ -4,6 +4,7 @@ import astropy.constants as const import astropy.units as u +from scipy.optimize import root_scalar matplotlib.rcParams.update({'font.size': 26}) ### @@ -14,7 +15,7 @@ def Ne_leblanc(rho) : In: rho - heliocentric distance in R_sun Out : Density in per cc ''' - alpha_1,alpha_2,alpha_3 = 3.3e5,4.1e6,8.0e7 + alpha_1, alpha_2, alpha_3 = 3.3e5, 4.1e6, 8.0e7 return alpha_1*rho**-2 + alpha_2*rho**-4 + alpha_3*rho**-6 def Ne_leblanc_deriv(rho) : @@ -22,7 +23,7 @@ def Ne_leblanc_deriv(rho) : In: rho - heliocentric distance in R_sun Out : delta Density in per cc/ delta rho in R_sun ''' - alpha_1,alpha_2,alpha_3 = 3.3e5,4.1e6,8.0e7 + alpha_1, alpha_2, alpha_3 = 3.3e5, 4.1e6, 8.0e7 return -2*alpha_1*rho**-3 + -4*alpha_2*rho**-5 + -6*alpha_3*rho**-7 def Ne_saito(r) : @@ -54,17 +55,17 @@ def Ne_parker(r): In : r # Heliocentric radius in R_S Out : n_e : heliospheric electorn density in cm^-3 ''' - m1,m2,m3 = 1.39e6,3e8,4.8e9 - alpha,beta,gamma = 2.3,6,14 + m1, m2, m3 = 1.39e6, 3e8, 4.8e9 + alpha, beta, gamma = 2.3, 6, 14 return m1/r**alpha+m2/r**beta+m3/r**gamma -def Ne_mann(r,N0 = 8.775e8, D = 9.88): +def Ne_mann(r, N0=8.775e8, D=9.88): # r in rsun n_e_model_MANN = N0 * np.exp(D * (1 / r - 1)) # Mann & Klassen (2005) return n_e_model_MANN -def Ne_allen(r, a0=2.99,a1=1.55,a2=0.036): +def Ne_allen(r, a0=2.99, a1=1.55, a2=0.036): # r in rsun n_e_model_ALLEN = 1e8 * ( a0*np.power(r,-16.) + a1*np.power(r,-6.) + a2*np.power(r,-1.5)) # Allen return n_e_model_ALLEN @@ -82,7 +83,6 @@ def parker(hr, f, a): """ [1] for fundamental backbone, [2] for Harmonic Backbone. f: frequency (MHz). a: fold (1 - 4, [1] for quiet Sun, [4] for active regions). """ - hr=0 n_e = Ne(f)*a r = R_parker(n_e, calibration_factor=1) @@ -102,38 +102,40 @@ def Ne_nk(r,fold=1): def Ne_parker_deriv(r) : - m1,m2,m3 = 1.39e6,3e8,4.8e9 - alpha,beta,gamma = 2.3,6,14 + m1, m2, m3 = 1.39e6, 3e8, 4.8e9 + alpha, beta, gamma = 2.3, 6, 14 return -alpha*m1/r**(alpha+1) - beta*m2/r**(beta+1) - gamma*m3/r**(gamma+1) -def Ne_helios(r) : +def Ne_helios(r): # https://link.springer.com/content/pdf/10.1007%2FBF00173965.pdf n0 = 6.1*u.cm**-3 - if not isinstance(r,u.quantity.Quantity) : r *= u.R_sun + if not isinstance(r, u.quantity.Quantity): + r *= u.R_sun return (n0 * (const.au.to('m')/r.to('m'))**2.1).to("cm^-3") -def Ne_moncuquet(r) : +def Ne_moncuquet(r): # https://iopscience.iop.org/article/10.3847/1538-4365/ab5a84/pdf n0 = 10*u.cm**-3 - return (n0 * (const.au.to('m')/r.to('m'))**2).to("cm^-3") # Plasma Frequency in MHz def f_pe(Ne) : - if not isinstance(Ne,u.quantity.Quantity) : Ne *= u.cm**-3 - return ((80.6*Ne.to("cm^-3").value)**0.5/1e3 )*u.MHz + if not isinstance(Ne, u.quantity.Quantity): + Ne *= u.cm**-3 + return ((80.6*Ne.to("cm^-3").value)**0.5/1e3)*u.MHz def Ne(f_pe) : - if not isinstance(f_pe,u.quantity.Quantity) : f_pe *= u.MHz + if not isinstance(f_pe,u.quantity.Quantity): + f_pe *= u.MHz return f_pe.to("kHz").value**2/80.6*u.cm**-3 def R_moncuquet(Ne) : - if not isinstance(Ne,u.quantity.Quantity) : Ne *= u.cm**-3 + if not isinstance(Ne,u.quantity.Quantity): + Ne *= u.cm**-3 n0 = 10*u.cm**-3 return (n0/Ne)**0.5 * u.au -from scipy.optimize import root_scalar def find_radius_parker(target_ne, calibration_factor=1): @@ -152,7 +154,7 @@ def difference(r): def R_parker(target_ne, calibration_factor=1): target_ne = target_ne.value - sol = root_scalar(find_radius_parker(target_ne,calibration_factor), bracket=[1, 215], method='brentq') + sol = root_scalar(find_radius_parker(target_ne, calibration_factor), bracket=[1, 215], method='brentq') if sol.converged: return sol.root else: @@ -168,7 +170,7 @@ def R_saito(target_ne, calibration_factor=1): def R_leblanc(target_ne, calibration_factor=1): target_ne = target_ne.value - sol = root_scalar(find_radius_leblanc(target_ne,calibration_factor), bracket=[1, 215], method='brentq') + sol = root_scalar(find_radius_leblanc(target_ne, calibration_factor), bracket=[1, 215], method='brentq') if sol.converged: return sol.root else: @@ -240,4 +242,4 @@ def R_leblanc(target_ne, calibration_factor=1): print(f"Parker : {Ne_parker(r)} cm^{-3}") print(f"Allen : {Ne_allen(r)} cm^{-3}") - print(f"BELLA : {Ne_allen(r, *popt)} cm^{-3}") + # print(f"BELLA : {Ne_allen(r, *popt)} cm^{-3}") diff --git a/bella/multilaterate/positioner_mapping_parallel.py b/bella/multilaterate/positioner_mapping_parallel.py index 076ab4a..fe4fd58 100755 --- a/bella/multilaterate/positioner_mapping_parallel.py +++ b/bella/multilaterate/positioner_mapping_parallel.py @@ -1,5 +1,4 @@ import os -import sys import math import logging import argparse @@ -9,10 +8,6 @@ from contextlib import contextmanager import arviz as az -# caution: path[0] is reserved for script path (or '' in REPL) -# IF running from a different directory, point to the multilaterate directory with bayes_positioner.py file -#sys.path.insert(1, 'PATH/TO/Multilateratefolder') -import bayes_positioner as bp import matplotlib import matplotlib.patheffects as PathEffects import matplotlib.pyplot as plt @@ -28,6 +23,15 @@ from astropy.constants import R_sun, au, c +# caution: path[0] is reserved for script path (or '' in REPL) +# IF running from a different directory, point to the multilaterate directory with bayes_positioner.py file +#sys.path.insert(1, 'PATH/TO/Multilateratefolder') + +import sys +sys.path.insert(1, '/Users/canizares/Library/CloudStorage/OneDrive-Personal/Work/0_PhD/Projects/BELLA_Projects/TCDSolarBELLA') + +import bella.multilaterate.bayes_positioner as bp + @contextmanager def suppress_stdout(): @@ -68,7 +72,6 @@ def parallel_pos_map(x1,y1,stations,xrange,yrange,xres,yres,t_cadence=60, cores= x_true = np.array([x1 * R_sun.value, y1 * R_sun.value]) # true source position (m) v_true = c.value # speed of light (m/s) - t0_true = 100 # source time. can be any constant, as long as it is within the uniform distribution prior on t0 d_true = np.linalg.norm(stations - x_true, axis=1) t1_true = d_true / v_true # true time of flight values # t_obs = t1_true-t0_true# true time difference of arrival values @@ -82,11 +85,11 @@ def parallel_pos_map(x1,y1,stations,xrange,yrange,xres,yres,t_cadence=60, cores= mu = results_out[0] sd = results_out[1] - t1_pred = results_out[2] + results_out[2] trace = results_out[3] - summary = results_out[4] - t_emission = results_out[5] - v_analysis = results_out[6] + results_out[4] + results_out[5] + results_out[6] @@ -100,7 +103,7 @@ def parallel_pos_map(x1,y1,stations,xrange,yrange,xres,yres,t_cadence=60, cores= # print(f"t1_pred: {t1_pred}") # print(f"stations: {stations}") - if traceplotsave == True: + if traceplotsave is True: traceplotpath = f"./traceplots/{date_str}/traceplot_{xrange[0]}_{xrange[1]}_{yrange[0]}_{yrange[1]}_{xres}_{yres}" mkdirectory(traceplotpath) @@ -117,7 +120,7 @@ def parallel_pos_map(x1,y1,stations,xrange,yrange,xres,yres,t_cadence=60, cores= # difference detected observed x_true = x_true / R_sun.value - xy = mu[0] / R_sun.value, mu[1] / R_sun.value + mu[0] / R_sun.value, mu[1] / R_sun.value #delta_obs = sqrt((xy[0] - x_true[0]) ** 2 + (xy[1] - x_true[1]) ** 2) # Deltaobs = xy - x_true Deprecated delta_obs = np.amax([sd[0]/R_sun.value , sd[1]/R_sun.value]) # Bayesian uncertainty. Pick the largest one. @@ -129,8 +132,9 @@ def parallel_pos_map(x1,y1,stations,xrange,yrange,xres,yres,t_cadence=60, cores= print(f" P({x1},{y1}), {colored('SUCCESS', 'green')}: {res}") - except: + except Exception as e: delta_obs = 600 + print(f"SIM ERROR: {e}") print(f"SIM {colored('FAILED', 'red')} at P({x1},{y1}), SKIPPED") res = np.array([delta_obs, x1, y1]) tloop1 = dt.datetime.now() @@ -146,9 +150,9 @@ def parallel_pos_map(x1,y1,stations,xrange,yrange,xres,yres,t_cadence=60, cores= return res def plot_map_simple(delta_obs, xmapaxis, ymapaxis, stations,vmin=0,vmax=30, savefigure=False, showfigure=True, title="",figdir="./MapFigures", date_str="date", filename="fig.jpg"): - xres = xmapaxis[1]-xmapaxis[0] - yres = ymapaxis[1]-ymapaxis[0] - N_STATIONS = len(stations) + xmapaxis[1]-xmapaxis[0] + ymapaxis[1]-ymapaxis[0] + len(stations) fig, ax = plt.subplots(1,1,figsize=(8,8)) plt.subplots_adjust(top=1, bottom=0) @@ -179,12 +183,12 @@ def plot_map_simple(delta_obs, xmapaxis, ymapaxis, stations,vmin=0,vmax=30, save ax.set_ylabel(r"'HEE - Y / $R_{\odot}$'", fontsize=22) ax.set_title(title, fontsize=22) - if savefigure == True: + if savefigure is True: figdir = f"{figdir}/{date_str}" mkdirectory(figdir) plt.savefig(figdir+filename, bbox_inches='tight', pad_inches=0.01, dpi=300) - if showfigure == True: + if showfigure is True: plt.show(block=False) else: plt.close(fig) @@ -217,9 +221,9 @@ def plot_map_simple_withTracked(delta_obs, trackedtypeIII, xmapaxis, ymapaxis, s sd = sd_vals/R_sun.value - xres = xmapaxis[1]-xmapaxis[0] - yres = ymapaxis[1]-ymapaxis[0] - N_STATIONS = len(stations) + xmapaxis[1]-xmapaxis[0] + ymapaxis[1]-ymapaxis[0] + len(stations) # PARKER SPIRAL @@ -264,7 +268,7 @@ def plot_map_simple_withTracked(delta_obs, trackedtypeIII, xmapaxis, ymapaxis, s i+=1 - if parker_spiral == True: + if parker_spiral is True: ax.plot(x_parker,y_parker,"k--") ax.plot(x_parker_plus,y_parker_plus,"k--", markersize=0.5) ax.plot(x_parker_minus,y_parker_minus,"k--") @@ -280,7 +284,7 @@ def plot_map_simple_withTracked(delta_obs, trackedtypeIII, xmapaxis, ymapaxis, s # # ELLIPSES colors = cm.turbo(list(np.linspace(0,1.0,len(xy)))) - if confidence ==True: + if confidence is True: i = 0 ell_track_uncertainty = matplotlib.patches.Ellipse(xy=(xy[i, 0], xy[i, 1]), width=2 * sd[i, 0], height=2 * sd[i, 1], @@ -307,12 +311,6 @@ def plot_map_simple_withTracked(delta_obs, trackedtypeIII, xmapaxis, ymapaxis, s # windlab.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='k')]) - - - - - - ax.set_aspect('equal') ax.set_xlim(xmapaxis[0], xmapaxis[-1]) @@ -340,12 +338,12 @@ def plot_map_simple_withTracked(delta_obs, trackedtypeIII, xmapaxis, ymapaxis, s ax.set_ylabel(r"'HEE - Y / $R_{\odot}$'", fontsize=22) ax.set_title(title, fontsize=22) - if savefigure == True: + if savefigure is True: figdir = f"{figdir}/{date_str}" mkdirectory(figdir) plt.savefig(figdir+filename, bbox_inches='tight', pad_inches=0.01, dpi=300) - if showfigure == True: + if showfigure is True: plt.show(block=False) else: plt.close(fig) @@ -361,12 +359,12 @@ def plot_map_simple_withTracked_zoom(delta_obs, trackedtypeIII, xmapaxis, ymapax xy = xy_vals/R_sun.value tracked_freqs = trackedtypeIII.freq sd_vals = np.array(list(trackedtypeIII.sd)) - sd = sd_vals/R_sun.value + sd_vals/R_sun.value - xres = xmapaxis[1]-xmapaxis[0] - yres = ymapaxis[1]-ymapaxis[0] - N_STATIONS = len(stations) + xmapaxis[1]-xmapaxis[0] + ymapaxis[1]-ymapaxis[0] + len(stations) # PARKER SPIRAL @@ -418,7 +416,6 @@ def plot_map_simple_withTracked_zoom(delta_obs, trackedtypeIII, xmapaxis, ymapax im_track = ax.scatter(xy[:,0], xy[:,1],c = tracked_freqs, marker="s",edgecolors="w",cmap='plasma', s=250, label="TrackedBeam") - i = 0 # ell_track_uncertainty = matplotlib.patches.Ellipse(xy=(xy[i, 1], xy[i, 0]), # width=4 * sd[i, 1], height=4 * sd[i, 0], # angle=0., edgecolor=tracked_freqs, cmap='plasma', label="Posterior ($2\sigma$)", lw=1.5) @@ -482,12 +479,12 @@ def plot_map_simple_withTracked_zoom(delta_obs, trackedtypeIII, xmapaxis, ymapax ax.indicate_inset_zoom(axins, edgecolor="black") - if savefigure == True: + if savefigure is True: figdir = f"{figdir}/{date_str}" mkdirectory(figdir) plt.savefig(figdir+filename, bbox_inches='tight', pad_inches=0.01, dpi=300) - if showfigure == True: + if showfigure is True: plt.show(block=False) else: plt.close(fig) @@ -502,9 +499,9 @@ def savedata(delta_obs, xmapaxis, ymapaxis, stations, filename="output.pkl"): mkdirectory(directory) - xres = xmapaxis[1]-xmapaxis[0] - yres = ymapaxis[1]-ymapaxis[0] - N_STATIONS = len(stations) + xmapaxis[1]-xmapaxis[0] + ymapaxis[1]-ymapaxis[0] + len(stations) results = [xmapaxis, ymapaxis, delta_obs, stations] with open(filename, 'wb') as outp: @@ -640,11 +637,12 @@ def __init__(self, freq, xy, sd): year = observation_date.year except ValueError: print("Invalid date format. Please use YYYY-MM-DD format.") - except: - print("No date provided.") + except Exception as e: + print(f"Warning. {e}, Date assigned to manual") day = 4 month = 12 year = 2021 + print( f"Date assigned: {year}-{month:02d}-{day:02d}") date_str = f"{year}_{month:02d}_{day:02d}" @@ -680,32 +678,43 @@ def __init__(self, freq, xy, sd): if "wind" in spacecraft: windsc=True sc_buff.append("wind") - else:windsc=False + else: + windsc=False + if "stereo_a" in spacecraft: steasc=True sc_buff.append("stereo_a") - else:steasc=False + else: + steasc=False + if "stereo_b" in spacecraft: stebsc=True sc_buff.append("stereo_b") - else:stebsc=False + else: + stebsc=False + if "solo" in spacecraft: solosc=True sc_buff.append("solo") - else:solosc=False + else: + solosc=False + if "psp" in spacecraft: pspsc=True sc_buff.append("psp") - else:pspsc=False + else: + pspsc=False + if "mex" in spacecraft: mexsc=True sc_buff.append("mex") - else:mexsc=False + else: + mexsc=False spacecraft = sc_buff print(f"Spacecraft selected: {spacecraft}") - sc_str="" + sc_str = "" if windsc: sc_str += '_WIND' if steasc: @@ -725,9 +734,9 @@ def __init__(self, freq, xy, sd): theta_sc = int(sys.argv[1]) print(f"theta_sc: {theta_sc}") - L1 = [0.99*(au/R_sun),0] - L4 = [(au/R_sun)*np.cos(radians(60)),(au/R_sun)*np.sin(radians(60))] - L5 = [(au/R_sun)*np.cos(radians(60)),-(au/R_sun)*np.sin(radians(60))] + L1 = [0.99*(au/R_sun), 0] + L4 = [(au/R_sun)*np.cos(radians(60)), (au/R_sun)*np.sin(radians(60))] + L5 = [(au/R_sun)*np.cos(radians(60)), -(au/R_sun)*np.sin(radians(60))] # ahead = [(au/R_sun)*np.cos(radians(theta_sc)),(au/R_sun)*np.sin(radians(theta_sc))] # behind = [(au/R_sun)*np.cos(radians(theta_sc)),-(au/R_sun)*np.sin(radians(theta_sc))] @@ -744,7 +753,7 @@ def __init__(self, freq, xy, sd): elif date_str == "test": stations_rsun = np.array([[200, 200], [-200, -200], [-200, 200], [200, -200]]) elif date_str == "manual": - stations_rsun = np.array([[45.27337378, 9.90422281],[-24.42715218,-206.46280171],[ 212.88183411,0.]]) + stations_rsun = np.array([[45.27337378, 9.90422281], [-24.42715218, -206.46280171], [212.88183411, 0.]]) date_str = f"{year}_{month:02d}_{day:02d}" else: solarsystem = solarmap.get_HEE_coord(date=[year, month, day], objects=spacecraft) @@ -801,7 +810,7 @@ def __init__(self, freq, xy, sd): else: runserver = False ############################################################### - if runserver == True: + if runserver is True: # Settings for server (Make sure runserver=True) runsimulation = True run_failed_again = False # Keep False. NEEDS DEBUGGING> FAILS @@ -834,10 +843,10 @@ def __init__(self, freq, xy, sd): filename = f"./Data/BG_Maps/{date_str}/results_{xrange[0]}_{xrange[-1]}_{yrange[0]}_{yrange[-1]}_{xres}_{yres}_{N_STATIONS}stations{sc_str}_{cadence}s.pkl" # Does not make sense to save data or run flagged points if simulation is off. - if runsimulation == False: + if runsimulation is False: run_savedata = False run_failed_again = False - if run_loaddata == True: + if run_loaddata is True: #Check if data exists isExist = os.path.exists(filename) if not isExist: @@ -853,7 +862,7 @@ def __init__(self, freq, xy, sd): run_savedata = True run_loaddata = False - if runsimulation == True: + if runsimulation is True: # Check if data exists isExist = os.path.exists(filename) if isExist: @@ -875,7 +884,7 @@ def __init__(self, freq, xy, sd): """ --------------------------------------------------------------------- """ """ ------------------------- Simulation -------------------------------- """ """ --------------------------------------------------------------------- """ - if runsimulation == True: + if runsimulation is True: compcost0 = dt.datetime.now() coresforloop = multiprocessing.cpu_count() @@ -908,26 +917,26 @@ def __init__(self, freq, xy, sd): """ ------------------------ Save Data ---------------------------------- """ - if run_savedata == True: + if run_savedata is True: savedata(delta_obs, xmapaxis, ymapaxis, stations_rsun, filename=filename) """ --------------------------------------------------------------------- """ """ ------------------------ Load Data ---------------------------------- """ - if run_loaddata == True: + if run_loaddata is True: xmapaxis, ymapaxis, delta_obs, stations_rsun = loaddata(filename) """ --------------------------------------------------------------------- """ """ ------------------------ PLOT Data ---------------------------------- """ - if run_plotdata == True: + if run_plotdata is True: fname = f"/bayes_positioner_map_{xmapaxis[0]}_{xmapaxis[-1]}_{ymapaxis[0]}_{ymapaxis[-1]}_{xres}_{yres}_{N_STATIONS}_30s.jpg" plot_map_simple(delta_obs*2, xmapaxis, ymapaxis, stations_rsun, vmin=np.min(delta_obs*2), vmax=np.max(delta_obs*2), showfigure=True, savefigure=True, date_str=date_str, filename=fname) """ --------------------------------------------------------------------- """ """ ------------------------ Median Filter --------------------------- """ - if run_median_filter == True: + if run_median_filter is True: # from scipy.ndimage import median_filter as medfil delta_obs2 = delta_obs*2 # 2std for 95% confidence intervalm median_filter_image = medfil(delta_obs2, size=(6,6)) diff --git a/bella/type_III_fitter/dynspec.py b/bella/type_III_fitter/dynspec.py index d180c5a..e33ddcc 100644 --- a/bella/type_III_fitter/dynspec.py +++ b/bella/type_III_fitter/dynspec.py @@ -16,12 +16,12 @@ from astropy.constants import R_sun, au from astropy.time import Time -r_sun = R_sun.value -AU=au.value from radiospectra.spectrogram import Spectrogram # in the process of updating old spectrogram from sunpy.net import attrs as a +r_sun = R_sun.value +AU=au.value ########################################################################## # Functions to LOAD SC dyn spectra @@ -91,7 +91,8 @@ def waves_spec(start, endt,datatype="RAD1", bg_subtraction=False, lighttravelshi return waves_spec def local_waves_spec_l2_60s(file, datatype='RAD1', kind='SMEAN', bg_subtraction=False, lighttravelshift=0): - from read_wi_wa_l2_60s_data import read_l2_60s # This file can be found at https://cdpp-archive.cnes.fr/ + from bella.type_III_fitter.read_wi_wa_l2_60s_data import \ + read_l2_60s # This file can be found at https://cdpp-archive.cnes.fr/ header, data = read_l2_60s(file) frequencies = data['FKHZ'] @@ -320,13 +321,13 @@ def solo_rpw_hfr(filepath): # V3-V1=6, B_MF=7, HF_V1-V2=9, HF_V2-V3=10, HF_V3-V1=11) sensor = rpw_l2_hfr.varget('SENSOR_CONFIG') freq_uniq = np.unique(rpw_l2_hfr.varget('FREQUENCY')) # frequency channels list - sample_time = rpw_l2_hfr.varget('SAMPLE_TIME') + rpw_l2_hfr.varget('SAMPLE_TIME') agc1 = rpw_l2_hfr.varget('AGC1') agc2 = rpw_l2_hfr.varget('AGC2') - flux_density1 = rpw_l2_hfr.varget('FLUX_DENSITY1') - flux_density2 = rpw_l2_hfr.varget('FLUX_DENSITY2') + rpw_l2_hfr.varget('FLUX_DENSITY1') + rpw_l2_hfr.varget('FLUX_DENSITY2') rpw_l2_hfr.close() # l2_cdf_file.close() @@ -420,7 +421,7 @@ def check_cadence(times, plot=True, method="mode", title=""): time_idx.append(times[i].datetime) - if plot==True: + if plot is True: plt.figure() ax = plt.gca() plt.plot_date(time_idx, dtime, ".") diff --git a/bella/type_III_fitter/openmarsis.py b/bella/type_III_fitter/openmarsis.py index 3adc92e..cfba3c4 100644 --- a/bella/type_III_fitter/openmarsis.py +++ b/bella/type_III_fitter/openmarsis.py @@ -188,7 +188,7 @@ def marsis_spectra(file_path, quickplot=False, histogram=[], lighttravelshift=0) } mars_spec = Spectrogram(mars_spec_bg.T, meta) # Generates radiospectra spectrogram - if quickplot == True: + if quickplot is True: mars_mm_h = np.percentile(mars_spec.data, [histogram[0], histogram[1]]) # histogram levels fig, axes = plt.subplots(1, 1, sharex=True, figsize=(10, 10)) diff --git a/bella/type_III_fitter/psp_quicklook.py b/bella/type_III_fitter/psp_quicklook.py index 177e7b3..dd513f7 100644 --- a/bella/type_III_fitter/psp_quicklook.py +++ b/bella/type_III_fitter/psp_quicklook.py @@ -8,8 +8,6 @@ import astropy.units as u from astropy.time import Time -plt.rcParams.update({'font.size': 22}) -plt.rcParams.update({'font.family': "Times New Roman"}) # import pickle import argparse @@ -29,6 +27,13 @@ from sunpy.net import attrs as a + + +plt.rcParams.update({'font.size': 22}) +plt.rcParams.update({'font.family': "Times New Roman"}) + + + def backSub(data, percentile=1): """ Background subtraction: This function has been modified from Eoin Carley's backsub function diff --git a/bella/type_III_fitter/solo_quicklook.py b/bella/type_III_fitter/solo_quicklook.py index 22f7163..3d6a3fb 100644 --- a/bella/type_III_fitter/solo_quicklook.py +++ b/bella/type_III_fitter/solo_quicklook.py @@ -9,9 +9,6 @@ import astropy.units as u from astropy.time import Time -plt.rcParams.update({'font.size': 22}) -plt.rcParams.update({'font.family': "Times New Roman"}) -# import pickle import argparse import pyspedas # pip install git+https://github.com/STBadman/pyspedas @@ -29,6 +26,8 @@ # from sunpy.net import attrs as a +plt.rcParams.update({'font.size': 22}) +plt.rcParams.update({'font.family': "Times New Roman"}) def backSub(data, percentile=1): """ Background subtraction: diff --git a/bella/type_III_fitter/solo_quicklook_L3_data.py b/bella/type_III_fitter/solo_quicklook_L3_data.py index 27bf27b..b102c0b 100644 --- a/bella/type_III_fitter/solo_quicklook_L3_data.py +++ b/bella/type_III_fitter/solo_quicklook_L3_data.py @@ -1,6 +1,5 @@ import datetime as dt -import dynspec import matplotlib as mpl from matplotlib import pyplot as plt from radiospectra.spectrogram import Spectrogram @@ -11,11 +10,11 @@ from sunpy.net import attrs as a +from bella.type_III_fitter import dynspec + plt.rcParams.update({'font.size': 22}) plt.rcParams.update({'font.family': "Times New Roman"}) - - # def open_rpw_l3(cdf_file_path, bg_subtraction=False, lighttravelshift=0): # #CDFLIB reads Epoch inccorrectly. DO NOT USE CDFLIB with solo data # # Open the CDF file @@ -70,18 +69,18 @@ def open_rpw_l3(cdf_file_path, bg_subtraction=False, lighttravelshift=0): cdf = pycdf.CDF(cdf_file_path) epoch = cdf["Epoch"] frequency = cdf["FREQUENCY"] - background = cdf["BACKGROUND"] - sensor_config = cdf["SENSOR_CONFIG"] - channel = cdf["CHANNEL"] - timing = cdf["TIMING"] - quality_flag = cdf["QUALITY_FLAG"] - interpol_flags = cdf["INTERPOL_FLAG"] - psd_v2 = cdf["PSD_V2"] - psd_flux = cdf["PSD_FLUX"] + cdf["BACKGROUND"] + cdf["SENSOR_CONFIG"] + cdf["CHANNEL"] + cdf["TIMING"] + cdf["QUALITY_FLAG"] + cdf["INTERPOL_FLAG"] + cdf["PSD_V2"] + cdf["PSD_FLUX"] psd_sfu = cdf["PSD_SFU"] - lbl1_sc_pos_hci = cdf["LBL1_SC_POS_HCI"] - sc_pos_hci = cdf["SC_POS_HCI"] - rep1_sc_pos_hci = cdf["REP1_SC_POS_HCI"] + cdf["LBL1_SC_POS_HCI"] + cdf["SC_POS_HCI"] + cdf["REP1_SC_POS_HCI"] epoch_dt = epoch[:] diff --git a/bella/type_III_fitter/solo_quicklook_solo_env.py b/bella/type_III_fitter/solo_quicklook_solo_env.py index d5a5ef9..5264e30 100644 --- a/bella/type_III_fitter/solo_quicklook_solo_env.py +++ b/bella/type_III_fitter/solo_quicklook_solo_env.py @@ -204,7 +204,7 @@ def mkdirectory(directory): plt.show(block=False) data = [rpw_AGC1, rpw_AGC2] - if savedata == True: + if savedata is True: directory = mkdirectory("solo_data/RPW/") with open(f'{directory}rpw_extracted_radiospectra_{YYYY}{MM:02}{dd:02}.pkl', 'wb') as output_file: pickle.dump(data, output_file) diff --git a/bella/type_III_fitter/typeIIIfitter.py b/bella/type_III_fitter/typeIIIfitter.py index 2afe334..1bce0f5 100755 --- a/bella/type_III_fitter/typeIIIfitter.py +++ b/bella/type_III_fitter/typeIIIfitter.py @@ -1,626 +1,18 @@ #!/usr/bin/env python -# ------------------------------------------------------------ -# Script which demonstrates how to find the best-fit -# parameters of a Gauss-Hermite line shape model -# -# Vog, 26 Mar 2012 -# ------------------------------------------------------------ + import os -import matplotlib.dates as mdates import matplotlib.pyplot as plt import numpy as np -from kapteyn import kmpfit # https://www.astro.rug.nl/software/kapteyn/intro.html#installinstructions -from scipy.optimize import fsolve -from scipy.special import wofz -ln2 = np.log(2) -PI = np.pi -from math import sqrt from datetime import datetime from scipy.optimize import curve_fit -########################################################################## -# GH fitting functions -########################################################################## -def gausshermiteh3h4(x, A, x0, s, h3, h4): - # ------------------------------------------------------------ - # The Gauss-Hermite function is a superposition of functions of the form - # F = (x-xc)/s - # E = A.Exp[-1/2.F^2] * {1 + h3[c1.F+c3.F^3] + h4[c5+c2.F^2+c4.F^4]} - # ------------------------------------------------------------ - c0 = sqrt(6.0) / 4.0 - c1 = -sqrt(3.0) - c2 = -sqrt(6.0) - c3 = 2.0 * sqrt(3.0) / 3.0 - c4 = sqrt(6.0) / 3.0 - - F = (x - x0) / s - E = A * np.exp(-0.5 * F * F) * (1.0 + h3 * F * (c3 * F * F + c1) + h4 * (c0 + F * F * (c2 + c4 * F * F))) - return E - -def hermite2gauss(par, dpar): - # ------------------------------------------------------------ - # Convert Gauss-Hermite parameters to Gauss(like)parameters. - # - # We use the first derivative of the Gauss-Hermite function - # to find the maximum, usually around 'x0' which is the center - # of the (pure) Gaussian part of the function. - # If F = (x-x0)/s then the function for which we want the - # the zero's is A0+A1*F+A2*F^2+A3*F^3+A4*F^4+A5*F^5 = 0 - # c0 = 1/4sqrt(6) c1 = -sqrt(3) c2 = -sqrt(6) - # c3 = 2/3sqrt(3) c4 = 1/3sqrt(6) - # ------------------------------------------------------------ - sqrt2pi = sqrt(2.0 * PI) - amp, x0, s, h3, h4 = par - damp, dx0, ds, dh3, dh4 = dpar # The errors in those parameters - c0 = sqrt(6.0) / 4.0 - c1 = -sqrt(3.0) - c2 = -sqrt(6.0) - c3 = 2.0 * sqrt(3.0) / 3.0 - c4 = sqrt(6.0) / 3.0 - - A = np.zeros(6) - A[0] = -c1 * h3 - A[1] = h4 * (c0 - 2.0 * c2) + 1.0 - A[2] = h3 * (c1 - 3.0 * c3) - A[3] = h4 * (c2 - 4.0 * c4) - A[4] = c3 * h3 - A[5] = c4 * h4 - - # Define the function that represents the derivative of - # the GH function. You need it to find the position of the maximum. - fx = lambda x: A[0] + x * (A[1] + x * (A[2] + x * (A[3] + x * (A[4] + x * A[5])))) - xr = fsolve(fx, 0, full_output=True) - xm = s * xr[0] + x0 - ampmax = gausshermiteh3h4(xm, amp, x0, s, h3, h4) - - # Get line strength - f = 1.0 + h4 * sqrt(6.0) / 4.0 - area = amp * s * f * sqrt2pi - d_area = sqrt2pi * sqrt(s * s * f * f * damp * damp + \ - amp * amp * f * f * ds * ds + \ - 3.0 * amp * amp * s * s * dh4 * dh4 / 8.0) - - # Get mean - mean = x0 + sqrt(3.0) * h3 * s - d_mean = sqrt(dx0 * dx0 + 3.0 * h3 * h3 * ds * ds + 3.0 * s * s * dh3 * dh3) - - # Get dispersion - f = 1.0 + h4 * sqrt(6.0) - dispersion = abs(s * f) - d_dispersion = sqrt(f * f * ds * ds + 6.0 * s * s * dh4 * dh4) - - # Skewness - f = 4.0 * sqrt(3.0) - skewness = f * h3 - d_skewness = f * dh3 - - # Kurtosis - f = 8.0 * sqrt(6.0) - kurtosis = f * h4 - d_kurtosis = f * dh4 - - res = dict(xmax=xm, amplitude=ampmax, area=area, mean=mean, dispersion=dispersion, \ - skewness=skewness, kurtosis=kurtosis, d_area=d_area, d_mean=d_mean, \ - d_dispersion=d_dispersion, d_skewness=d_skewness, d_kurtosis=d_kurtosis) - return res - -def voigt(x, y): - # The Voigt function is also the real part of - # w(z) = exp(-z^2) erfc(iz), the complex probability function, - # which is also known as the Faddeeva function. Scipy has - # implemented this function under the name wofz() - z = x + 1j * y - I = wofz(z).real - return I - -def Voigt(nu, alphaD, alphaL, nu_0, A): - # The Voigt line shape in terms of its physical parameters - f = np.sqrt(ln2) - x = (nu - nu_0) / alphaD * f - y = alphaL / alphaD * f - V = A * f / (alphaD * np.sqrt(np.pi)) * voigt(x, y) - return V - -def funcV(p, x): - # Compose the Voigt line-shape - alphaD, alphaL, nu_0, I, z0 = p - return Voigt(x, alphaD, alphaL, nu_0, I) + z0 - -def funcG(p, x): - # Model function is a gaussian - A, mu, sigma, zerolev = p - return (A * np.exp(-(x - mu) * (x - mu) / (2 * sigma * sigma)) + zerolev) - -def funcGH(p, x): - # Model is a Gauss-Hermite function - A, xo, s, h3, h4, zerolev = p - return gausshermiteh3h4(x, A, xo, s, h3, h4) + zerolev - -def residualsV(p, data): - # Return weighted residuals of Voigt - x, y, err = data - return (y - funcV(p, x)) / err - -def residualsG(p, data): - # Return weighted residuals of Gauss - x, y, err = data - return (y - funcG(p, x)) / err - -def residualsGH(p, data): - # Return weighted residuals of Gauss-Hermite - x, y, err = data - return (y - funcGH(p, x)) / err - -# def fit_lc(x,y,freq,sigma=1, sdauto=2, h3guess=0.1, h4guess=0, z0guess = 1,method="sigma", plotresults=False, saveplots=False, dir="/lightcurves", fname="freq"): -# """sdauto picks n number of points from max point as initial guess for standard deviation""" -# """ Recommended sigma 1 or 2""" -# """ h3 - skewness , h4 - kurtosis""" -# N = len(y) -# err = np.ones(N) -# -# if not type(sdauto) is int: -# raise TypeError("sdauto: Only integers are allowed") -# # Fit the Gauss-Hermite model -# # Initial estimates for A, xo, s, h3, h4, z0 -# -# Aguess = np.max(y) -# xguess = x[np.where(y==Aguess)[0][0]] -# sguess = xguess - x[np.where(y==Aguess)[0][0]-sdauto] -# # h3guess = h3guess -# # h4guess = h4guess -# # z0guess = z0guess -# p0 = [Aguess, xguess, sguess, h3guess, h4guess, z0guess] -# fitterGH = kmpfit.Fitter(residuals=residualsGH, data=(x, y, err)) -# # fitterGH.parinfo = [{}, {}, {}, {}, {}] # Take zero level fixed in fit -# fitterGH.fit(params0=p0) -# print("\n========= Fit results Gaussian profile ==========") -# print("Initial params:", fitterGH.params0) -# print("Params: ", fitterGH.params) -# print("Iterations: ", fitterGH.niter) -# print("Function ev: ", fitterGH.nfev) -# print("Uncertainties: ", fitterGH.xerror) -# print("dof: ", fitterGH.dof) -# print("chi^2, rchi2: ", fitterGH.chi2_min, fitterGH.rchi2_min) -# print("stderr: ", fitterGH.stderr) -# print("Status: ", fitterGH.status) -# -# A, x0, s, h3, h4, z0GH = fitterGH.params -# -# -# # xm, ampmax, area, mean, dispersion, skewness, kurtosis -# res = hermite2gauss(fitterGH.params[:-1], fitterGH.stderr[:-1]) -# print("Gauss-Hermite max=%g at x=%g" % (res['amplitude'], res['xmax'])) -# print("Area :", res['area'], '+-', res['d_area']) -# print("Mean (X0) :", res['mean'], '+-', res['d_mean']) -# print("Dispersion:", res['dispersion'], '+-', res['d_dispersion']) -# print("Skewness :", res['skewness'], '+-', res['d_skewness']) -# print("Kurtosis :", res['kurtosis'], '+-', res['d_kurtosis']) -# -# xd = np.linspace(x.min(), x.max(), 500) -# -# xd_dt = timestamp2datetime(xd) -# x_date = timestamp2datetime(x) -# -# -# yfit = funcGH(fitterGH.params, xd) # signal model fit -# ybg = [z0GH] * len(xd) # background -# -# maxfit_idx = np.where(yfit==np.max(yfit))[0][0] -# -# supported_methods = ["sigma", "hwhm", "thirdmax", "peak"] -# if method == "sigma": -# # Get an n sigma (standard distributions) off the peak -# xrise = xd[maxfit_idx] - sigma*s -# elif method == "hwhm": -# # half width half maximum LHS -# # need to find nearest point in model because y is output not input, thus can't import into funcGH. -# # this is ok because we are under the resolution of the instrument -# ycheck = yfit[0:maxfit_idx] -# halfmax = (yfit[maxfit_idx] - z0GH)/2 + z0GH -# y_hwhm,y_hwhm_idx = find_nearest(ycheck,halfmax) -# xrise = xd[y_hwhm_idx] -# elif method == "thirdmax": -# # half width third maximum LHS -# # need to find nearest point in model because y is output not input, thus can't import into funcGH. -# # this is ok because we are under the resolution of the instrument -# ycheck = yfit[0:maxfit_idx] -# halfmax = (yfit[maxfit_idx] - z0GH)/3 + z0GH -# y_hwhm,y_hwhm_idx = find_nearest(ycheck,halfmax) -# xrise = xd[y_hwhm_idx] -# -# elif method == "peak": -# xrise = xd[maxfit_idx] -# -# else: -# raise AttributeError(f"method not supported. supported methods: {supported_methods}") -# -# -# -# yrise = funcGH(fitterGH.params, xrise) -# -# xrise = datetime.fromtimestamp(xrise) -# -# -# -# # Plot the result -# showparamstitle = False -# -# plt.rcParams.update({'font.size': 15}) -# plt.rc('legend', fontsize=15) -# fig1 = plt.figure(figsize=(9, 6), dpi=150) -# frame1 = fig1.add_subplot(1, 1, 1) -# frame1.plot(x_date, y, 'bo', label="data") -# label = "G-H Model" -# frame1.plot(xd_dt, yfit, 'c', ls='--', label=label) -# frame1.plot(xd_dt,ybg , "y", ls="--", label='G-H Background') -# frame1.set_xlabel(f"{x_date[0].year}/{x_date[0].month:02}/{x_date[0].day:02} [UTC]",fontsize=15) -# frame1.set_ylabel("Relative power",fontsize=15) -# title = ""#"Profile with Gauss-Hermite model\n" -# t = (res['area'], res['mean'], res['dispersion'], res['skewness'], res['kurtosis']) -# title += "GH: $\gamma_{gh}$=%.1f $x_{0_{gh}}$=%.1f $\sigma_{gh}$ = %.2f $\\xi_1$=%.2f $\\xi_f$=%.2f\n" % t -# title += f"freq: {freq}MHz" -# frame1.plot(xrise, yrise, 'r*', label="detection") -# -# frame1.grid(False) -# leg = plt.legend(loc='upper right') -# frame1.add_artist(leg) -# -# if showparamstitle == True: -# frame1.set_title(title, fontsize=15) -# else: -# plt.text(.01, .99, f"freq: {freq:.2f} MHz", ha='left', va='top', transform=frame1.transAxes) -# -# frame1.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M')) -# plt.gcf().autofmt_xdate() -# -# -# -# -# if saveplots==True: -# mkdirectory(dir) -# fid = f"{dir}/{fname}.jpg" -# plt.gcf().savefig(fid, dpi=150) -# print(f"{fid} saved") -# -# if plotresults == True: -# plt.show(block=False) -# else: -# plt.close() -# -# return [xrise, yrise] - -def fit_lc(x,y,freq,sigma=1, sdauto=2, h3guess=0.1, h4guess=0, z0guess = 1,method="sigma",profile="LE", plotresults=False, saveplots=False, dir="/lightcurves", fname="freq"): - """sdauto picks n number of points from max point as initial guess for standard deviation""" - """ Recommended sigma 1 or 2""" - """ h3 - skewness , h4 - kurtosis""" - N = len(y) - err = np.ones(N) - - if type(sdauto) is not int: - raise TypeError("sdauto: Only integers are allowed") - # Fit the Gauss-Hermite model - # Initial estimates for A, xo, s, h3, h4, z0 - - Aguess = np.max(y) - xguess = x[np.where(y==Aguess)[0][0]] - sguess = xguess - x[np.where(y==Aguess)[0][0]-sdauto] - # h3guess = h3guess - # h4guess = h4guess - # z0guess = z0guess - p0 = [Aguess, xguess, sguess, h3guess, h4guess, z0guess] - fitterGH = kmpfit.Fitter(residuals=residualsGH, data=(x, y, err)) - # fitterGH.parinfo = [{}, {}, {}, {}, {}] # Take zero level fixed in fit - fitterGH.fit(params0=p0) - print("\n========= Fit results Gaussian profile ==========") - print("Initial params:", fitterGH.params0) - print("Params: ", fitterGH.params) - print("Iterations: ", fitterGH.niter) - print("Function ev: ", fitterGH.nfev) - print("Uncertainties: ", fitterGH.xerror) - print("dof: ", fitterGH.dof) - print("chi^2, rchi2: ", fitterGH.chi2_min, fitterGH.rchi2_min) - print("stderr: ", fitterGH.stderr) - print("Status: ", fitterGH.status) - - A, x0, s, h3, h4, z0GH = fitterGH.params - - - # xm, ampmax, area, mean, dispersion, skewness, kurtosis - res = hermite2gauss(fitterGH.params[:-1], fitterGH.stderr[:-1]) - print("Gauss-Hermite max=%g at x=%g" % (res['amplitude'], res['xmax'])) - print("Area :", res['area'], '+-', res['d_area']) - print("Mean (X0) :", res['mean'], '+-', res['d_mean']) - print("Dispersion:", res['dispersion'], '+-', res['d_dispersion']) - print("Skewness :", res['skewness'], '+-', res['d_skewness']) - print("Kurtosis :", res['kurtosis'], '+-', res['d_kurtosis']) - - xd = np.linspace(x.min(), x.max(), 500) - - xd_dt = timestamp2datetime(xd) - x_date = timestamp2datetime(x) - - - yfit = funcGH(fitterGH.params, xd) # signal model fit - ybg = [z0GH] * len(xd) # background - - maxfit_idx = np.where(yfit==np.max(yfit))[0][0] - - supported_methods = ["sigma", "hwhm", "thirdmax","quartermax", "peak"] - if method == "sigma": - # Get an n sigma (standard distributions) off the peak - xrise = xd[maxfit_idx] - sigma*s - - elif method == "hwhm": - # half width half maximum LHS (LE) and RHS(TE) - # need to find nearest point in model because y is output not input, thus can't import into funcGH. - if profile == "LE": - ycheck = yfit[0:maxfit_idx] - ydetection = (yfit[maxfit_idx] - z0GH)/2 + z0GH - y_quarter, y_quarter_idx = find_nearest(ycheck, ydetection) - xrise = xd[y_quarter_idx] - elif profile == "TE": - ycheck = yfit[maxfit_idx:] - ydetection = (yfit[maxfit_idx] - z0GH)/2 + z0GH - y_quarter, y_quarter_idx = find_nearest(ycheck, ydetection) - xrise = xd[maxfit_idx+y_quarter_idx] - elif method == "thirdmax": - # half width half maximum LHS (LE) and RHS(TE) - # need to find nearest point in model because y is output not input, thus can't import into funcGH. - if profile == "LE": - ycheck = yfit[0:maxfit_idx] - ydetection = (yfit[maxfit_idx] - z0GH)/3 + z0GH - y_quarter, y_quarter_idx = find_nearest(ycheck, ydetection) - xrise = xd[y_quarter_idx] - elif profile == "TE": - ycheck = yfit[maxfit_idx:] - ydetection = (yfit[maxfit_idx] - z0GH)/3 + z0GH - y_quarter, y_quarter_idx = find_nearest(ycheck, ydetection) - xrise = xd[maxfit_idx+y_quarter_idx] - elif method == "quartermax": - # half width half maximum LHS (LE) and RHS(TE) - # need to find nearest point in model because y is output not input, thus can't import into funcGH. - if profile == "LE": - ycheck = yfit[0:maxfit_idx] - ydetection = (yfit[maxfit_idx] - z0GH)/4 + z0GH - y_quarter, y_quarter_idx = find_nearest(ycheck, ydetection) - xrise = xd[y_quarter_idx] - elif profile == "TE": - ycheck = yfit[maxfit_idx:] - ydetection = (yfit[maxfit_idx] - z0GH)/4 + z0GH - y_quarter, y_quarter_idx = find_nearest(ycheck, ydetection) - xrise = xd[maxfit_idx+y_quarter_idx] - - elif method == "peak": - xrise = xd[maxfit_idx] - - else: - raise AttributeError(f"method not supported. supported methods: {supported_methods}") - - - - yrise = funcGH(fitterGH.params, xrise) - - xrise = datetime.fromtimestamp(xrise) - - - - # Plot the result - showparamstitle = False - - plt.rcParams.update({'font.size': 15}) - plt.rc('legend', fontsize=15) - fig1 = plt.figure(figsize=(7, 8), dpi=150) - frame1 = fig1.add_subplot(1, 1, 1) - frame1.plot(x_date, y, 'bo')#, label="Data") - label = "G-H Model" - frame1.plot(xd_dt, yfit, 'c', ls='--', label=label) - frame1.plot(xd_dt,ybg , "y", ls="--", label='G-H Background') - frame1.set_xlabel(f"{x_date[0].year}/{x_date[0].month:02}/{x_date[0].day:02} [UTC]",fontsize=15) - frame1.set_ylabel("Relative power",fontsize=15) - title = ""#"Profile with Gauss-Hermite model\n" - t = (res['area'], res['mean'], res['dispersion'], res['skewness'], res['kurtosis']) - title += "GH: $\\gamma_{gh}$=%.1f $x_{0_{gh}}$=%.1f $\\sigma_{gh}$ = %.2f $\\xi_1$=%.2f $\\xi_f$=%.2f\n" % t - - title += f"{freq:.2f} MHz" - - frame1.plot(xrise, yrise, 'r*')#, label="Detection") - frame1.axvline(x=xrise, color='red', linestyle='--') - # plt.text(.11, .45, f"Detection", ha='left', va='top', transform=frame1.transAxes) - - # Add an arrow annotation - arrow_props = { - 'arrowstyle': '->', # Arrow style - 'color': 'red', # Arrow color - 'linewidth': 2, # Arrow linewidth - } - - ymidpoint = (yfit.max() + yfit.min())/2 - frame1.annotate('Detection', xy=(xrise, ymidpoint), xytext=(datetime.fromtimestamp(xd[3]), ymidpoint), arrowprops=arrow_props) - - frame1.grid(False) - leg = plt.legend(loc='upper right') - frame1.add_artist(leg) - - if showparamstitle == True: - frame1.set_title(title, fontsize=15) - else: - freq_lab = f"{freq:.2f}" - # plt.text(.01, .95, f"{freq:.2f} MHz", ha='left', va='top', transform=frame1.transAxes) - - if freq_lab == "0.28": - plt.text(.01, .95, f"(a) {freq:.2f} MHz", ha='left', va='top', transform=frame1.transAxes) - elif freq_lab == "1.42": - plt.text(.01, .95, f"(b) {freq:.2f} MHz", ha='left', va='top', transform=frame1.transAxes) - elif freq_lab == "1.98": - plt.text(.01, .95, f"(c) {freq:.2f} MHz", ha='left', va='top', transform=frame1.transAxes) - elif freq_lab == "0.12": - plt.text(.01, .95, f"(d) {freq:.2f} MHz", ha='left', va='top', transform=frame1.transAxes) - else: - plt.text(.01, .95, f" {freq:.2f} MHz", ha='left', va='top', transform=frame1.transAxes) - - frame1.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M')) - plt.gcf().autofmt_xdate() - - if saveplots==True: - mkdirectory(dir) - fid = f"{dir}{fname}.jpg" - plt.gcf().savefig(fid, dpi=150) - print(f"{fid} saved") - - if plotresults == True: - plt.show(block=False) - else: - plt.close() - - return [xrise, yrise] - -def auto_rise_times(spectra, freqlims=[], timelims=[],sigma=1, sdauto=2, h3guess=0.1, h4guess=0, z0guess=1, method="sigma",profile="LE", plotresults=False, saveplots=False): - - freqs = spectra.frequencies - times = spectra.times - data = spectra.data - - if freqlims==[]: - freqlimmin = freqs[0] - freqlimmax = freqs[-1] - else: - freqlimmin, freqlimmax = freqlims - - minfreq, minfreqidx = find_nearest(freqs, freqlimmin) - maxfreq, maxfreqidx = find_nearest(freqs, freqlimmax) - - if timelims==[]: - mintime = times[0] - maxtime = times[-1] - else: - mintime, maxtime = timelims - - time_idx = np.where(np.logical_and(times >= mintime, times <= maxtime)) - - times_dt = [] - for each in times[time_idx]: - times_dt.append(each.datetime) - - - dir = f"lightcurves/lc_{times_dt[0].year}_{times_dt[0].month:02}_{times_dt[0].day:02}/{spectra.observatory}/{method}_{profile}/" - - risetimes = [] - riseval = [] - testfreq = [] - x = np.array(datetime2timestamp(times_dt)) - for f in range(minfreqidx,maxfreqidx): #f = 1 - testdata = data[f, time_idx][0] - normfact = np.abs(testdata[0]) - if normfact != 0: - y = np.array(testdata / normfact) # normalised - fname = f"{spectra.observatory.replace(' ','')}_{method}_{freqs[f].value:.2f}" - risex, risey = fit_lc(x, y,freq=freqs[f].value, sigma=sigma,profile=profile, sdauto=sdauto, h3guess=h3guess, h4guess=h4guess, z0guess=z0guess,method=method, plotresults=plotresults, saveplots=saveplots, dir=dir, fname=fname ) - risetimes.append(risex) - riseval.append(risey) - testfreq.append(freqs[f].value) - - return risetimes, riseval, testfreq - - -def rise_times_atmax(spectra, freqlims=[], timelims=[]): - - freqs = spectra.frequencies - times = spectra.times - data = spectra.data - - if freqlims==[]: - freqlimmin = freqs[0] - freqlimmax = freqs[-1] - else: - freqlimmin, freqlimmax = freqlims - - minfreq, minfreqidx = find_nearest(freqs, freqlimmin) - maxfreq, maxfreqidx = find_nearest(freqs, freqlimmax) - - if timelims==[]: - mintime = times[0] - maxtime = times[-1] - else: - mintime, maxtime = timelims - - time_idx = np.where(np.logical_and(times >= mintime, times <= maxtime)) - - times_dt = [] - for each in times[time_idx]: - times_dt.append(each.datetime) - - risetimes = [] - riseval = [] - testfreq = [] - x = np.array(datetime2timestamp(times_dt)) - for f in range(minfreqidx,maxfreqidx): #f = 1 - testdata = data[f, time_idx][0] - y_max = testdata.max() - t = np.where(testdata==y_max) - risex = timestamp2datetime(x[t])[0] - risetimes.append(risex) - riseval.append(y_max) - testfreq.append(freqs[f].value) - - return risetimes, riseval, testfreq - - -def testmethod(spectra, freqlims=[], timelims=[],sigma=1, sdauto=2, h3guess=0.1, h4guess=0, z0guess=1, plotresults=False, saveplots=False): - - freqs = spectra.frequencies - times = spectra.times - data = spectra.data - - if freqlims==[]: - freqlimmin = freqs[0] - freqlimmax = freqs[-1] - else: - freqlimmin, freqlimmax = freqlims - - minfreq, minfreqidx = find_nearest(freqs, freqlimmin) - maxfreq, maxfreqidx = find_nearest(freqs, freqlimmax) - - if timelims==[]: - mintime = times[0] - maxtime = times[-1] - else: - mintime, maxtime = timelims - - time_idx = np.where(np.logical_and(times >= mintime, times <= maxtime)) - - times_dt = [] - for each in times[time_idx]: - times_dt.append(each.datetime) - - - - dir=f"lightcurves/{spectra.observatory}" - - - - risetimes = [] - riseval = [] - testfreq = [] - x = np.array(datetime2timestamp(times_dt)) - for f in range(minfreqidx,maxfreqidx): #f = 1 - testdata = data[f, time_idx][0] - normfact = np.abs(testdata[0]) - if normfact != 0: - y = np.array(testdata / normfact) # normalised - risex, risey = fit_lc(x, y,freq=freqs[f].value, sigma=sigma, sdauto=sdauto, h3guess=h3guess, h4guess=h4guess, z0guess=z0guess,plotresults=plotresults, saveplots=saveplots, dir=dir, fname=f"{f}" ) - risetimes.append(risex) - riseval.append(risey) - testfreq.append(freqs[f].value) - - return risetimes, riseval, testfreq - - - +ln2 = np.log(2) +PI = np.pi ########################################################################## # Type III fitting functions ########################################################################## @@ -647,15 +39,12 @@ def typeIII_func(times, popt, pcov, xref, num=100): def log_func(x,a,b,c): return (-1/b) * np.log((x-c)/a) -def exponential_func(x, a, b, c): - return a * np.exp(-b * x) + c - -def exponential_func2(x, a, b, c,d): - return a * np.exp((-b * x) + d) + c - def log_func2(x,a,b,c,d): return (1/b) * (d - np.log((x-c)/a)) +def exponential_func(x, a, b, c): + return a * np.exp(-b * x) + c + def exponential_func2(x, a, b, c,d): return a * np.exp((-b * x) + d) + c @@ -716,7 +105,7 @@ def typeIIIfitting(risetimes,testfreq, fitfreqs,freqs4tri, plot_residuals=False) fittimes_for_residuals = reciprocal_2ndorder(np.array(testfreq), *popt) residuals = np.subtract(xdata, fittimes_for_residuals) - if plot_residuals == True: + if plot_residuals is True: plt.figure() plt.plot(residuals, "r.") plt.title("residuals WAVES LEADING EDGE") @@ -774,57 +163,3 @@ def find_nearest(array, value): def f_to_angs(f_mhz,c=299792458): angstrom = (c / (f_mhz * 10 ** 6)) * 10 ** 10 return angstrom - - - - - - - -if __name__=="__main__": - - - spectra = solo_spec # spectra from radiospectra - freqs = spectra.frequencies - times = spectra.times - data = spectra.data - - minfreq, minfreqidx = find_nearest(freqs, freqlimmin) - maxfreq, maxfreqidx = find_nearest(freqs, freqlimmax) - - mintime = datetime(YYYY, MM, dd, HH_0,mm_0) - maxtime = datetime(YYYY, MM, dd, 3,30) - time_idx = np.where(np.logical_and(times >= mintime, times <= maxtime)) - - times_dt = [] - for each in times[time_idx]: - times_dt.append(each.datetime) - - - risetimes = [] - riseval = [] - testfreq = [] - x = np.array(datetime2timestamp(times_dt)) - for f in range(minfreqidx,maxfreqidx): #f = 1 - testdata = data[f, time_idx][0] - normfact = np.abs(testdata[0]) - if normfact != 0: - y = np.array(testdata / normfact) # normalised - risex, risey = fit_lc(x, y, sdauto=1, h3guess=0.1, h4guess=0, z0guess=0, plotresults=False,saveplots=True,dir=f"lightcurves/{spectra.observatory}", fname=f"{f}" ) - risetimes.append(risex) - riseval.append(risey) - testfreq.append(freqs[f].value) - - - - mm = np.percentile(spectra.data, [10, 99]) - - fig, axes = plt.subplots(1, 1, sharex=True, figsize=(10, 10)) - spectra.plot(axes=axes, vmin=mm[0], vmax=mm[1]) - axes.plot(risetimes, testfreq, 'r*') - axes.set_ylim(reversed(axes.get_ylim())) - axes.set_yscale('log') - axes.set_ylim([freqlimmax, freqlimmin]) - axes.set_xlim(datetime(YYYY, MM, dd, HH_0,mm_0), datetime(YYYY, MM, dd, 3,0)) - plt.subplots_adjust(hspace=0.31) - plt.show(block=False) diff --git a/examples/bella_triangulation_2012_06_07.py b/examples/bella_triangulation_2012_06_07.py index 8a39520..e7c6f36 100755 --- a/examples/bella_triangulation_2012_06_07.py +++ b/examples/bella_triangulation_2012_06_07.py @@ -11,8 +11,11 @@ import pymc3 as pm # General import solarmap -from bayes_positioner import * -from bayesian_tracker import * +import sys +sys.path.insert(1, '/Users/canizares/Library/CloudStorage/OneDrive-Personal/Work/0_PhD/Projects/BELLA_Projects/TCDSolarBELLA') + +from bella.multilaterate.bayes_positioner import triangulate +from bella.multilaterate.bayesian_tracker import loadtypeiii, savetrackedtypeiii # Standard Library imports from astropy.constants import R_sun @@ -131,14 +134,14 @@ tloop0 = dt.datetime.now() check = False - while check == False: + while check is False: try: # connect mu, sd, t1_pred, trace, summary, t_emission_fromtrace, v_analysis = triangulate(stations, typeIII_times[i_freq], t_cadence=60, cores=4, progressbar=True, report=0, plot=0,traceplot=True, savetraceplot=True, traceplotdir=f'{date_str}', traceplotfn=f'{i_freq}.jpg') check = True - except: - print("MULTILATERATION FAILED, most likely by divergance, try again") - pass + except Exception as e: + print(f"ERROR: {e}, try again") + tloop1 = dt.datetime.now() dtloop = tloop1-tloop0 diff --git a/examples/bella_triangulation_2012_06_07_SCATTER.py b/examples/bella_triangulation_2012_06_07_SCATTER.py index 33e626b..0e2bc32 100755 --- a/examples/bella_triangulation_2012_06_07_SCATTER.py +++ b/examples/bella_triangulation_2012_06_07_SCATTER.py @@ -11,8 +11,11 @@ import pymc3 as pm # General import solarmap -from bayes_positioner import * -from bayesian_tracker import * +import sys +sys.path.insert(1, '/Users/canizares/Library/CloudStorage/OneDrive-Personal/Work/0_PhD/Projects/BELLA_Projects/TCDSolarBELLA') + +from bella.multilaterate.bayes_positioner import triangulate +from bella.multilaterate.bayesian_tracker import loadtypeiii, savetrackedtypeiii # Standard Library imports from astropy.constants import R_sun @@ -131,14 +134,13 @@ tloop0 = dt.datetime.now() check = False - while check == False: + while check is False: try: # connect mu, sd, t1_pred, trace, summary, t_emission_fromtrace, v_analysis = triangulate(stations, typeIII_times[i_freq], t_cadence=60, cores=4, progressbar=True, report=0, plot=0,traceplot=True, savetraceplot=True, traceplotdir=f'{date_str}', traceplotfn=f'{i_freq}.jpg') check = True - except: - print("MULTILATERATION FAILED, most likely by divergance, try again") - pass + except Exception as e: + print(f"ERROR: {e}, try again") tloop1 = dt.datetime.now() dtloop = tloop1-tloop0 diff --git a/examples/bella_triangulation_2021_12_4_LE_BURST_A.py b/examples/bella_triangulation_2021_12_4_LE_BURST_A.py index f2edaa3..ed058f3 100755 --- a/examples/bella_triangulation_2021_12_4_LE_BURST_A.py +++ b/examples/bella_triangulation_2021_12_4_LE_BURST_A.py @@ -13,8 +13,12 @@ # Parallel processing imports import multiprocessing -import bayes_positioner as bp -import bayesian_tracker as btrack +import sys +sys.path.insert(1, '/Users/canizares/Library/CloudStorage/OneDrive-Personal/Work/0_PhD/Projects/BELLA_Projects/TCDSolarBELLA') + +import bella.multilaterate.bayes_positioner as bp +import bella.multilaterate.bayesian_tracker as btrack + import numpy as np # Third Party imports # BAYESIAN IMPORTS @@ -86,18 +90,36 @@ for sc in args.spacecraft: spacecraft.append(sc.lower()) - if "wind" in spacecraft:windsc=True - else:windsc=False - if "stereo_a" in spacecraft:steasc=True - else:steasc=False - if "stereo_b" in spacecraft:stebsc=True - else:stebsc=False - if "solo" in spacecraft:solosc=True - else:solosc=False - if "psp" in spacecraft:pspsc=True - else:pspsc=False - if "mex" in spacecraft:mexsc=True - else:mexsc=False + if "wind" in spacecraft: + windsc=True + else: + windsc=False + + if "stereo_a" in spacecraft: + steasc=True + else: + steasc=False + + if "stereo_b" in spacecraft: + stebsc=True + else: + stebsc=False + + if "solo" in spacecraft: + solosc=True + else: + solosc=False + + if "psp" in spacecraft: + pspsc=True + else: + pspsc=False + + if "mex" in spacecraft: + mexsc=True + else: + mexsc=False + print(f"Spacecraft selected: {spacecraft}") @@ -285,14 +307,14 @@ tloop0 = dt.datetime.now() check = False - while check == False: + while check is False: try: # connect mu, sd, t1_pred, trace, summary, t_emission_fromtrace, v_analysis = bp.triangulate(stations, typeIII_times[i_freq], t_cadence=60, v_sd=0.0001*c.value, cores=4, progressbar=True, report=0, plot=0,traceplot=True, savetraceplot=True, traceplotdir=f'{sc_str}_{profile}_{note}', traceplotfn=f'{i_freq}.jpg') check = True - except: - print("MULTILATERATION FAILED, most likely by divergance, try again") - pass + except Exception as e: + print(f"ERROR: {e}, try again") + tloop1 = dt.datetime.now() dtloop = tloop1-tloop0 diff --git a/examples/plot_2021_12_04.py b/examples/plot_2021_12_04.py index a8035dc..5035395 100644 --- a/examples/plot_2021_12_04.py +++ b/examples/plot_2021_12_04.py @@ -7,7 +7,11 @@ # from joblib import Parallel, delayed import multiprocessing -import bella_plotter as bplot + +import sys +sys.path.insert(1, '/Users/canizares/Library/CloudStorage/OneDrive-Personal/Work/0_PhD/Projects/BELLA_Projects/TCDSolarBELLA') + +import bella.multilaterate.bella_plotter as bplot import matplotlib.pyplot as plt import numpy as np import pymc3 as pm @@ -236,6 +240,6 @@ ax.tick_params(axis='both', which='major', labelsize=18) plt.show(block = False) - if savefigs == True: + if savefigs is True: dir = bplot.mkdirectory("./Figures/") plt.savefig(dir+'BELLA_map_4spacecraft.png', dpi=300) diff --git a/examples/stacked_dyn_spectra_2012_06_07_HIGHRES.py b/examples/stacked_dyn_spectra_2012_06_07_HIGHRES.py deleted file mode 100755 index 7949927..0000000 --- a/examples/stacked_dyn_spectra_2012_06_07_HIGHRES.py +++ /dev/null @@ -1,1031 +0,0 @@ -# Author: L Alberto Canizares canizares (at) cp.dias.ie -from datetime import datetime, timedelta - -import numpy as np -import pyspedas -from matplotlib import pyplot as plt -from matplotlib.colors import LogNorm -from numpy import interp -from pytplot import get_data -from scipy import stats -from scipy.interpolate import CubicSpline -from scipy.optimize import curve_fit - -import astropy.units as u -from astropy.time import Time - -plt.rcParams.update({'font.size': 22}) -from matplotlib.ticker import FormatStrFormatter - -from astropy.constants import R_sun, au - -r_sun = R_sun.value -AU=au.value - -import pickle - -import cdflib -from radiospectra.spectrogram import Spectrogram # in the process of updating old spectrogram -from typeIIIfitter import * - -from sunpy.net import attrs as a - -# from spacepy import pycdf - -# from rpw_mono.thr.hfr.reader import read_hfr_autoch - - -def f_to_angs(f_mhz,c=299792458): - angstrom = (c / (f_mhz * 10 ** 6)) * 10 ** 10 - return angstrom - -def backSub(data, percentile=1): - """ Background subtraction: - This function has been modified from Eoin Carley's backsub function - https://github.com/eoincarley/ilofar_scripts/blob/master/Python/bst/plot_spectro.py - - data: numpy 2D matrix of floating values for dynamic spectra - percentile: integer value def = 1. bottom X percentile of time slices - - - - METHOD ---- - * This function takes the bottom x percentile of time slices - * Averages those time slices. - * Subtracts the average value from the whole dataset - - - """ - # Get time slices with standard devs in the bottom nth percentile. - # Get average spectra from these time slices. - # Divide through by this average spec. - # Expects (row, column) - - print("Start of Background Subtraction of data") - dat = data.T - # dat = np.log10(dat) - # dat[np.where(np.isinf(dat) == True)] = 0.0 - dat_std = np.std(dat, axis=0) - dat_std = dat_std[np.nonzero(dat_std)] - min_std_indices = np.where(dat_std < np.percentile(dat_std, percentile))[0] - min_std_spec = dat[:, min_std_indices] - min_std_spec = np.mean(min_std_spec, axis=1) - dat = np.transpose(np.divide(np.transpose(dat), min_std_spec)) - - data = dat.T - - # Alternative: Normalizing frequency channel responses using median of values. - # for sb in np.arange(data.shape[0]): - # data[sb, :] = data[sb, :]/np.mean(data[sb, :]) - - print("Background Subtraction of data done") - return data - -def solo_rpw_hfr(filepath): - rpw_l2_hfr = cdflib.CDF(filepath) - # l2_cdf_file = pycdf.CDF(filepath) - - # times = l2_cdf_file['Epoch'] - # times = times[:] - - times = rpw_l2_hfr.varget('EPOCH') - - - freqs = rpw_l2_hfr.varget('FREQUENCY') - - # Indicates the THR sensor configuration (V1=1, V2=2, V3=3, V1-V2=4, V2-V3=5, - # V3-V1=6, B_MF=7, HF_V1-V2=9, HF_V2-V3=10, HF_V3-V1=11) - sensor = rpw_l2_hfr.varget('SENSOR_CONFIG') - freq_uniq = np.unique(rpw_l2_hfr.varget('FREQUENCY')) # frequency channels list - sample_time = rpw_l2_hfr.varget('SAMPLE_TIME') - - agc1 = rpw_l2_hfr.varget('AGC1') - agc2 = rpw_l2_hfr.varget('AGC2') - - flux_density1 = rpw_l2_hfr.varget('FLUX_DENSITY1') - flux_density2 = rpw_l2_hfr.varget('FLUX_DENSITY2') - - rpw_l2_hfr.close() - # l2_cdf_file.close() - - # For CH1 extract times, freqs and data points - slices1 = [] - times1 = [] - freq1 = [] - for cfreq in freq_uniq: - search = np.argwhere((freqs == cfreq) & (sensor[:, 0] == 9) & (agc1 != 0)) - if search.size > 0: - slices1.append(agc1[search]) - times1.append(times[search]) - freq1.append(cfreq) - - # For CH1 extract times, freqs and data points - slices2 = [] - times2 = [] - freq2 = [] - for cfreq in freq_uniq: - search = np.argwhere((freqs == cfreq) & (sensor[:, 1] == 9) & (agc2 != 0)) - if search.size > 0: - slices2.append(agc2[search]) - times2.append(times[search]) - freq2.append(cfreq) - - # Kinda arb but pick a time near middle of freq sweep - tt1 = np.hstack(times1)[:, -1]#160] - tt2 = np.hstack(times2)[:, -1]#50] - - spec1 = np.hstack(slices1) - spec2 = np.hstack(slices2) - - return tt1, freq1, spec1, tt2, freq2, spec2 - -def check_cadence(times, plot=True, method="mode", title=""): - dtime = [] - times = np.array(times) - time_idx = [] - for i in range(1, len(times)): - ti = times[i].datetime - ti_1 = times[i-1].datetime - diff = ti - ti_1 - dtime.append(diff.total_seconds()) - time_idx.append(times[i].datetime) - - - if plot==True: - from matplotlib import dates as mdates - from matplotlib import pyplot as plt - - plt.figure() - ax = plt.gca() - plt.plot_date(time_idx, dtime, ".") - plt.xlabel(f"Time {time_idx[0].year}/{time_idx[0].month:02}/{time_idx[0].day:02} ") - plt.ylabel("cadence") - plt.title(title) - formatter = mdates.DateFormatter("%H:%M") - ax.xaxis.set_major_formatter(formatter) - - plt.show(block=False) - - if method=="average": - return np.mean(dtime) - elif method=="max": - return np.amax(dtime) - elif method=="min": - return np.amin(dtime) - elif method=="mode": - return stats.mode(dtime)[0][0] - else: - print("method can only be mode, average, max or min") - -def find_nearest(array, value): - array = np.asarray(array) - idx = (np.abs(array - value)).argmin() - return array[idx], idx -def find_near_dt_idx(datetime_array, target_datetime): - return min(range(len(datetime_array)), key=lambda i: abs(datetime_array[i] - target_datetime)) - - - -def swaves_highres_spec(start, endt, probe='a', datatype='hfr', bg_subtraction=False): - """ - Downloads and processes high-resolution SWAVES spectrogram data from the STEREO mission. - - Parameters - ---------- - start : datetime - The start time of the data range to be downloaded. - endt : datetime - The end time of the data range to be downloaded. - probe : str, optional - The STEREO spacecraft probe to use. Default is 'a'. Alternative is 'b'. - datatype : str, optional - The type of SWAVES data to use. Default is 'hfr'. Alternative is 'lfr' - bg_subtraction : bool, optional - Whether or not to perform background subtraction on the spectrogram. Default is False. - - Returns - ------- - swaves_spec : `~sunpy.timeseries.Spectrogram` - A `sunpy.timeseries.Spectrogram` object containing the SWAVES spectrogram data. please see radiospectra https://github.com/sunpy/radiospectra - - Notes - ----- - This function uses the `pyspedas` and sunpy's 'radiospectra' libraries to download and process the SWAVES data. The resulting - spectrogram object includes metadata such as the observatory, instrument, detector, frequencies, and time range. - """ - stereo_variables = pyspedas.stereo.waves([start.strftime("%m/%d/%Y %H:%M:%S"), endt.strftime("%m/%d/%Y %H:%M:%S")], probe = probe, datatype = datatype) - print(f"Stereo {probe.upper()} {datatype.upper()} downloaded: stereo_variables:") - print(stereo_variables) - swaves_data = get_data('PSD_FLUX') - - - swaves_freqs_MHz = swaves_data.v/ 1e6 * u.MHz - - swaves_timestamps = swaves_data.times - swaves_times = Time(Time(np.array([datetime.utcfromtimestamp(t_) for t_ in swaves_timestamps])).isot) - - - swaves_psdarray = swaves_data.y - - meta = { - 'observatory': f"STEREO {probe.upper()}", - 'instrument': "SWAVES", - 'detector': datatype.upper(), - 'freqs': swaves_freqs_MHz, - 'times': swaves_times, - 'wavelength': a.Wavelength(swaves_freqs_MHz[0], swaves_freqs_MHz[-1]), - 'start_time': swaves_times[0], - 'end_time': swaves_times[-1] - } - swaves_spec = Spectrogram(swaves_psdarray.T, meta) - - if bg_subtraction: - swaves_spec.data = backSub(swaves_spec.data.T).T - - return swaves_spec - -def waves_spec(start, endt,datatype="RAD1", bg_subtraction=False): - """ - Downloads and processes high-resolution WAVES spectrogram data from the WIND mission. - - Parameters - ---------- - start : datetime - The start time of the data range to be downloaded. - endt : datetime - The end time of the data range to be downloaded. - - datatype : str, optional - The type of WAVES data to use. Default is 'RAD1'. Alternative is 'RAD2' - - bg_subtraction : bool, optional - Whether or not to perform background subtraction on the spectrogram. Default is False. - - Returns - ------- - wwaves_spec : `~sunpy.timeseries.Spectrogram` - A `sunpy.timeseries.Spectrogram` object containing the WAVES spectrogram data. please see radiospectra https://github.com/sunpy/radiospectra - - Notes - ----- - This function uses the `pyspedas` and sunpy's 'radiospectra' libraries to download and process the WAVES data. The resulting - spectrogram object includes metadata such as the observatory, instrument, detector, frequencies, and time range. - """ - wind_variables = pyspedas.wind.waves([start.strftime("%m/%d/%Y %H:%M:%S"), endt.strftime("%m/%d/%Y %H:%M:%S")]) - print("Wind WAVES downloaded: wind_variables:") - print(wind_variables) - waves_data = get_data(f'E_VOLTAGE_{datatype}') - - - - waves_freqs_MHz = waves_data.v/ 1e3 * u.MHz - - waves_timestamps = waves_data.times - waves_times = Time(Time(np.array([datetime.utcfromtimestamp(t_) for t_ in waves_timestamps])).isot) - - - waves_psdarray = waves_data.y - - meta = { - 'observatory': "WIND", - 'instrument': "WAVES", - 'detector': datatype, - 'freqs': waves_freqs_MHz, - 'times': waves_times, - 'wavelength': a.Wavelength(waves_freqs_MHz[0], waves_freqs_MHz[-1]), - 'start_time': waves_times[0], - 'end_time': waves_times[-1] - } - waves_spec = Spectrogram(waves_psdarray.T, meta) - - if bg_subtraction: - waves_spec.data = backSub(waves_spec.data.T).T - - return waves_spec - -def reciprocal_3rdorder(x, a0, a1, a2, a3): - return a0 + a1 / x + a2 / x ** 2+ a3 / x ** 3 - -def reciprocal_2ndorder(x, a0, a1, a2): - return a0 + a1 / x + a2 / x ** 2 - -def loadpickle(filenamef): - import pickle - with open(filenamef, 'rb') as inp: - results = pickle.load(inp) - - inp.close() - return results - - - - -if __name__=="__main__": - # Fido - YYYY = 2012 - MM = 6 - dd = 7 - HH_0 = 19 - mm_0 = 20 - HH_1 = 20 - mm_1 = 00 - # - background_subtraction = True - leadingedge = True - backbone = False - plot_residuals = False - figdir = mkdirectory(f"Figures/{YYYY}_{MM:02}_{dd:02}") - - mintime = datetime(YYYY, MM, dd, HH_0,mm_0) - maxtime = datetime(YYYY, MM, dd, HH_1,mm_1) - timelims = [mintime,maxtime] - - - - - - # Waves RAD 1 ( low freqs) - waves_spec_lfr = waves_spec(mintime, maxtime, datatype='RAD1', bg_subtraction=True) - # Waves RAD 2 (High freqs) - waves_spec_hfr = waves_spec(mintime, maxtime, datatype='RAD2', bg_subtraction=True) - - - # SWAVES A HFR - swaves_a_spec_hfr = swaves_highres_spec(mintime, maxtime, probe='a', datatype='hfr', bg_subtraction=True) - - # SWAVES B HFR - swaves_b_spec_hfr = swaves_highres_spec(mintime, maxtime, probe='b', datatype='hfr', bg_subtraction=True) - - - - - # Histogram levels - waves_mm_l = np.percentile(waves_spec_lfr.data, [20,100]) - waves_mm_h = np.percentile(waves_spec_hfr.data, [20,100]) - - # swaves_a_mm_l = np.percentile(swaves_a_spec_lfr.data, [70,99.99]) - swaves_a_mm_h = np.percentile(swaves_a_spec_hfr.data, [1,99.999]) - - # swaves_b_mm_l = np.percentile(swaves_b_spec_lfr.data, [70,99.99]) - swaves_b_mm_h = np.percentile(swaves_b_spec_hfr.data, [1,99.999]) - - - histograms_check = np.concatenate((waves_mm_l,waves_mm_h ,swaves_a_mm_h,swaves_b_mm_h)) - if any(check <= 0 for check in histograms_check): - print('WARNING: One or more histogram levels is 0. Increase this value to avoid error in plots. ') - - - # Dynamic spectra frequency range. - # LEADING EDGE - freqlimmin_LE = 0.02 #0.18 - freqlimmax_LE = 3 - freqlims_LE = [freqlimmin_LE, freqlimmax_LE] - - freqlimmin_BB = 0.1#0.18 - freqlimmax_BB = 3 - freqlims_BB = [freqlimmin_BB, freqlimmax_BB] - - - # Input for fit - freqfitmin = 0.02# 0.15 - freqfitmax = 13 - - # Frequencies to be extracted for Multilateration - freq4trimin = 0.15 - freq4trimax = 2 - - - - if leadingedge ==True: - # ------------------------- LEADING EDGE ----------------------------------------- # - # WAVES low - waves_risetimes_l_LE, waves_riseval_l_LE, waves_testfreq_l_LE = auto_rise_times(waves_spec_lfr, freqlims=[0.1, 3], timelims=timelims, sigma=1, sdauto=5, h3guess=0, h4guess=0,method="sigma", saveplots=True) - # WAVES high - waves_risetimes_h_LE, waves_riseval_h_LE, waves_testfreq_h_LE = auto_rise_times(waves_spec_hfr, freqlims=[0.15, 3], timelims=timelims, sigma=1, sdauto=2, h3guess=0.1, h4guess=0,method="sigma", saveplots=True) - - swaves_a_risetimes_h_LE, swaves_a_riseval_h_LE, swaves_a_testfreq_h_LE = auto_rise_times(swaves_a_spec_hfr, freqlims=freqlims_LE, timelims=timelims, sigma=1.5, sdauto=2, h3guess=0.1, h4guess=0,method="sigma", saveplots=True) - swaves_b_risetimes_h_LE, swaves_b_riseval_h_LE, swaves_b_testfreq_h_LE = auto_rise_times(swaves_b_spec_hfr, freqlims=freqlims_LE, timelims=timelims, sigma=1.5, sdauto=2, h3guess=0.1, h4guess=0,method="sigma", saveplots=True) - - - # Fine tuning. - # # MANUAL EDITS BASED ON VISUAL INSPECTION OF THE DYNSPECTRA. - # Fix Outliers and to prevent unphysical morphologies. - - waves_risetimes_l_LE[0] = datetime(2012, 6, 7, 19, 38) - waves_risetimes_l_LE[1] = datetime(2012, 6, 7, 19, 37) - waves_risetimes_l_LE[2] = datetime(2012, 6, 7, 19, 36, 50) - waves_risetimes_l_LE[3] = datetime(2012, 6, 7, 19, 36, 40) - waves_risetimes_l_LE[4] = datetime(2012, 6, 7, 19, 36, 30) - waves_risetimes_l_LE[5] = datetime(2012, 6, 7, 19, 36, 20) - waves_risetimes_l_LE[6] = datetime(2012, 6, 7, 19, 35, 50) - waves_risetimes_l_LE[7] = datetime(2012, 6, 7, 19, 35, 30) - waves_risetimes_l_LE[8] = datetime(2012, 6, 7, 19, 35, 50) - waves_risetimes_l_LE[9] = datetime(2012, 6, 7, 19, 35, 50) - waves_risetimes_l_LE[10] = datetime(2012, 6, 7, 19, 35,0) - waves_risetimes_l_LE[11] = datetime(2012, 6, 7, 19, 35,0) - waves_risetimes_l_LE[12] = datetime(2012, 6, 7, 19, 35,0) - - swaves_a_risetimes_h_LE[0] = datetime(2012, 6, 7, 19, 50, 0) - swaves_a_risetimes_h_LE[1] = datetime(2012, 6, 7, 19, 39, 42) - - swaves_b_risetimes_h_LE[0] = datetime(2012, 6, 7, 19, 46, 50) - swaves_b_risetimes_h_LE[1] = datetime(2012, 6, 7, 19, 40, 38) - - - - if backbone ==True: - # ------------------------- BACKBONE ----------------------------------------- # - # WAVES low - waves_risetimes_l_BB, waves_riseval_l_BB, waves_testfreq_l_BB = auto_rise_times(waves_spec_lfr, freqlims=freqlims_BB, timelims=timelims, sdauto=5, h3guess=0, h4guess=0,method="peak", saveplots=True) - # WAVES high - waves_risetimes_h_BB, waves_riseval_h_BB, waves_testfreq_h_BB = auto_rise_times(waves_spec_hfr, freqlims=freqlims_BB, timelims=timelims, sdauto=2, h3guess=0.1, h4guess=0,method="peak", saveplots=True) - - swaves_a_risetimes_h_BB, swaves_a_riseval_h_BB, swaves_a_testfreq_h_BB = auto_rise_times(swaves_a_spec_hfr, freqlims=freqlims_BB, timelims=timelims, sdauto=2, h3guess=0.1, h4guess=0,method="peak", saveplots=True) - swaves_b_risetimes_h_BB, swaves_b_riseval_h_BB, swaves_b_testfreq_h_BB = auto_rise_times(swaves_b_spec_hfr, freqlims=freqlims_BB, timelims=timelims, sdauto=2, h3guess=0.1, h4guess=0,method="peak", saveplots=True) - - - # -------------------------------------------------------------------------------------------------------------------- # - # FITTING TYPE IIIs - # -------------------------------------------------------------------------------------------------------------------- # - - - - - fitfreqs = np.logspace(np.log10(freqfitmin), np.log10(freqfitmax), num=200) - freqs4tri = np.logspace(np.log10(freq4trimin), np.log10(freq4trimax), num=200) - # fitfreqs2 = np.linspace(freqlimmin,freqlimmax, num=300 ) - - - # running this shows why fitfreqs needs to be in logspace. - # if fitfreqs is in logspace then there's an even distribution of points between low freqs and highfreqs - # xx = np.logspace(np.log10(freqfitmin), np.log10(freqfitmax), num=300) - # xx2 = np.linspace(freqfitmin,freqfitmax, num=300 ) - # - # plt.figure() - # plt.plot(xx, "r.") - # plt.plot(xx2, "b.") - # plt.yscale('log') - # plt.show(block=False) - - - - # ---------------------------------------------------------------- # - # LEADING EDGE - # ---------------------------------------------------------------- # - if leadingedge == True: - # waves_risetimes_l_LE_manual = datetime(2012, 6, 7, 19, 40, 0) - # waves_riseval_l_LE_manual = np.average(waves_riseval_l_LE) # just add arbitrary value to keep all matrices the same size - # waves_testfreq_l_LE_manual = - - # Wind - wavesrisetimes_LE = list(waves_risetimes_l_LE) + list(waves_risetimes_h_LE) - wavesriseval_LE = list(waves_riseval_l_LE) + list(waves_riseval_h_LE) - wavestestfreq_LE = list(waves_testfreq_l_LE) + list(waves_testfreq_h_LE) - - xdata_waves_LE, xref_waves_LE = epoch2time(wavesrisetimes_LE) - ydata_waves_LE = wavestestfreq_LE - - - - popt_waves_LE, pcov_waves_LE = curve_fit(reciprocal_2ndorder, ydata_waves_LE, xdata_waves_LE) - - fittimes_waves_LE = reciprocal_2ndorder(fitfreqs, *popt_waves_LE) - times4tri_waves_LE = reciprocal_2ndorder(freqs4tri, *popt_waves_LE) - - - - notnan = np.where(~np.isnan(fittimes_waves_LE)) - fitfreqs_waves_LE = fitfreqs[notnan] - fittimes_waves_notnan_LE = fittimes_waves_LE[notnan] - - times4tri_waves_LE_dt = time2epoch(times4tri_waves_LE, xref_waves_LE) - fittimes_corrected_waves_LE = time2epoch(fittimes_waves_notnan_LE, xref_waves_LE) - - # residuals waves Leading edge - fittimes_waves_for_residuals_LE = reciprocal_2ndorder(np.array(wavestestfreq_LE), *popt_waves_LE) - residuals_waves_LE = np.subtract(xdata_waves_LE, fittimes_waves_for_residuals_LE) - - if plot_residuals == True: - plt.figure() - plt.plot(residuals_waves_LE, "r.") - plt.title("residuals WAVES LEADING EDGE") - plt.xlabel("index") - plt.ylabel("difference") - plt.show(block=False) - - - - # Stereo A - swaves_a_risetimes_LE = list(swaves_a_risetimes_h_LE)# list([datetime(2012, 6, 7, 19, 58, 8)]) + - swaves_a_riseval_LE = list(swaves_a_riseval_h_LE)#list([np.average(swaves_a_riseval_h_LE)]) + - swaves_a_testfreq_LE = list(swaves_a_testfreq_h_LE)#list([0.15]) + - - xdata_swaves_a_LE, xref_swaves_a_LE = epoch2time(swaves_a_risetimes_LE) - ydata_swaves_a_LE = swaves_a_testfreq_LE - popt_swaves_a_LE, pcov_swaves_a_LE = curve_fit(reciprocal_2ndorder, ydata_swaves_a_LE, xdata_swaves_a_LE) - - - fittimes_swaves_a_LE = reciprocal_2ndorder(fitfreqs, *popt_swaves_a_LE) - times4tri_swaves_a_LE_dt = reciprocal_2ndorder(freqs4tri, *popt_swaves_a_LE) - - notnan = np.where(~np.isnan(fittimes_swaves_a_LE)) - fitfreqs_swaves_a_LE = fitfreqs[notnan] - fittimes_swaves_a_notnan_LE = fittimes_swaves_a_LE[notnan] - - - - fittimes_corrected_swaves_a_LE = time2epoch(fittimes_swaves_a_notnan_LE, xref_swaves_a_LE) - times4tri_swaves_a_LE_dt = time2epoch(times4tri_swaves_a_LE_dt, xref_swaves_a_LE) - - # residuals swaves a leading edge - fittimes_swavesa_for_residuals_LE = reciprocal_2ndorder(np.array(swaves_a_testfreq_LE), *popt_swaves_a_LE) - residuals_swaves_a_LE = np.subtract(xdata_swaves_a_LE, fittimes_swavesa_for_residuals_LE) - - if plot_residuals == True: - plt.figure() - plt.plot(residuals_swaves_a_LE, "r.") - plt.title("residuals SWAVES A LEADING EDGE") - plt.xlabel("index") - plt.ylabel("difference") - plt.show(block=False) - - # Stereo B - swaves_b_risetimes_LE = list(swaves_b_risetimes_h_LE) # list( [datetime(2012, 6, 7, 19, 59, 29)]) + - swaves_b_riseval_LE = list(swaves_b_riseval_h_LE) #list([np.average(swaves_a_riseval_h_LE)]) + - swaves_b_testfreq_LE = list(swaves_b_testfreq_h_LE) #list([0.15]) + - - xdata_swaves_b_LE, xref_swaves_b_LE = epoch2time(swaves_b_risetimes_LE) - ydata_swaves_b_LE = swaves_b_testfreq_LE - popt_swaves_b_LE, pcov_swaves_b_LE = curve_fit(reciprocal_2ndorder, ydata_swaves_b_LE, xdata_swaves_b_LE) - - - fittimes_swaves_b_LE = reciprocal_2ndorder(fitfreqs, *popt_swaves_b_LE) - times4tri_swaves_b_LE_dt = reciprocal_2ndorder(freqs4tri, *popt_swaves_b_LE) - - notnan = np.where(~np.isnan(fittimes_swaves_b_LE)) - fitfreqs_swaves_b_LE = fitfreqs[notnan] - fittimes_swaves_b_notnan_LE = fittimes_swaves_b_LE[notnan] - fittimes_corrected_swaves_b_LE = time2epoch(fittimes_swaves_b_notnan_LE, xref_swaves_b_LE) - times4tri_swaves_b_LE_dt = time2epoch(times4tri_swaves_b_LE_dt, xref_swaves_b_LE) - - # residuals swaves b leading edge - fittimes_swavesb_for_residuals_LE = reciprocal_2ndorder(np.array(swaves_b_testfreq_LE), *popt_swaves_b_LE) - residuals_swaves_b_LE = np.subtract(xdata_swaves_b_LE, fittimes_swavesb_for_residuals_LE) - if plot_residuals == True: - plt.figure() - plt.plot(residuals_swaves_b_LE, "r.") - plt.title("residuals SWAVES B LEADING EDGE") - plt.xlabel("index") - plt.ylabel("difference") - plt.show(block=False) - - - - - - # ---------------------------------------------------------------- # - # BACKBONE - # ---------------------------------------------------------------- # - if backbone == True: - - # Wind - wavesrisetimes_BB = list(waves_risetimes_l_BB) + list(waves_risetimes_h_BB) - wavesriseval_BB = list(waves_riseval_l_BB) + list(waves_riseval_h_BB) - wavestestfreq_BB = list(waves_testfreq_l_BB) + list(waves_testfreq_h_BB) - - xdata_waves_BB, xref_waves_BB = epoch2time(wavesrisetimes_BB) - ydata_waves_BB = wavestestfreq_BB - - - - popt_waves_BB, pcov_waves_BB = curve_fit(reciprocal_2ndorder, ydata_waves_BB, xdata_waves_BB) - - fittimes_waves_BB = reciprocal_2ndorder(fitfreqs, *popt_waves_BB) - times4tri_waves_BB = reciprocal_2ndorder(freqs4tri, *popt_waves_BB) - - - - notnan = np.where(~np.isnan(fittimes_waves_BB)) - fitfreqs_waves_BB = fitfreqs[notnan] - fittimes_waves_notnan_BB = fittimes_waves_BB[notnan] - - times4tri_waves_BB_dt = time2epoch(times4tri_waves_BB, xref_waves_BB) - fittimes_corrected_waves_BB = time2epoch(fittimes_waves_notnan_BB, xref_waves_BB) - - # residuals waves Leading edge - fittimes_waves_for_residuals_BB = reciprocal_2ndorder(np.array(wavestestfreq_BB), *popt_waves_BB) - residuals_waves_BB = np.subtract(xdata_waves_BB, fittimes_waves_for_residuals_BB) - - if plot_residuals == True: - plt.figure() - plt.plot(residuals_waves_BB, "r.") - plt.title("residuals WAVES BACKBONE") - plt.xlabel("index") - plt.ylabel("difference") - plt.show(block=False) - - - - # Stereo A - swaves_a_risetimes_BB = list(swaves_a_risetimes_h_BB)# list([datetime(2012, 6, 7, 19, 58, 8)]) + - swaves_a_riseval_BB = list(swaves_a_riseval_h_BB)#list([np.average(swaves_a_riseval_h_BB)]) + - swaves_a_testfreq_BB = list(swaves_a_testfreq_h_BB)#list([0.15]) + - - xdata_swaves_a_BB, xref_swaves_a_BB = epoch2time(swaves_a_risetimes_BB) - ydata_swaves_a_BB = swaves_a_testfreq_BB - popt_swaves_a_BB, pcov_swaves_a_BB = curve_fit(reciprocal_2ndorder, ydata_swaves_a_BB, xdata_swaves_a_BB) - - - fittimes_swaves_a_BB = reciprocal_2ndorder(fitfreqs, *popt_swaves_a_BB) - times4tri_swaves_a_BB_dt = reciprocal_2ndorder(freqs4tri, *popt_swaves_a_BB) - - notnan = np.where(~np.isnan(fittimes_swaves_a_BB)) - fitfreqs_swaves_a_BB = fitfreqs[notnan] - fittimes_swaves_a_notnan_BB = fittimes_swaves_a_BB[notnan] - - - - fittimes_corrected_swaves_a_BB = time2epoch(fittimes_swaves_a_notnan_BB, xref_swaves_a_BB) - times4tri_swaves_a_BB_dt = time2epoch(times4tri_swaves_a_BB_dt, xref_swaves_a_BB) - - # residuals swaves a leading edge - fittimes_swavesa_for_residuals_BB = reciprocal_2ndorder(np.array(swaves_a_testfreq_BB), *popt_swaves_a_BB) - residuals_swaves_a_BB = np.subtract(xdata_swaves_a_BB, fittimes_swavesa_for_residuals_BB) - - if plot_residuals == True: - plt.figure() - plt.plot(residuals_swaves_a_BB, "r.") - plt.title("residuals SWAVES A BACKBONE") - plt.xlabel("index") - plt.ylabel("difference") - plt.show(block=False) - - # Stereo B - swaves_b_risetimes_BB = list(swaves_b_risetimes_h_BB) # list( [datetime(2012, 6, 7, 19, 59, 29)]) + - swaves_b_riseval_BB = list(swaves_b_riseval_h_BB) #list([np.average(swaves_a_riseval_h_BB)]) + - swaves_b_testfreq_BB = list(swaves_b_testfreq_h_BB) #list([0.15]) + - - xdata_swaves_b_BB, xref_swaves_b_BB = epoch2time(swaves_b_risetimes_BB) - ydata_swaves_b_BB = swaves_b_testfreq_BB - popt_swaves_b_BB, pcov_swaves_b_BB = curve_fit(reciprocal_2ndorder, ydata_swaves_b_BB, xdata_swaves_b_BB) - - - fittimes_swaves_b_BB = reciprocal_2ndorder(fitfreqs, *popt_swaves_b_BB) - times4tri_swaves_b_BB_dt = reciprocal_2ndorder(freqs4tri, *popt_swaves_b_BB) - - notnan = np.where(~np.isnan(fittimes_swaves_b_BB)) - fitfreqs_swaves_b_BB = fitfreqs[notnan] - fittimes_swaves_b_notnan_BB = fittimes_swaves_b_BB[notnan] - fittimes_corrected_swaves_b_BB = time2epoch(fittimes_swaves_b_notnan_BB, xref_swaves_b_BB) - times4tri_swaves_b_BB_dt = time2epoch(times4tri_swaves_b_BB_dt, xref_swaves_b_BB) - - # residuals swaves b leading edge - fittimes_swavesb_for_residuals_BB = reciprocal_2ndorder(np.array(swaves_b_testfreq_BB), *popt_swaves_b_BB) - residuals_swaves_b_BB = np.subtract(xdata_swaves_b_BB, fittimes_swavesb_for_residuals_BB) - if plot_residuals == True: - plt.figure() - plt.plot(residuals_swaves_b_BB, "r.") - plt.title("residuals SWAVES B BACKBONE") - plt.xlabel("index") - plt.ylabel("difference") - plt.show(block=False) - - - # ---------------------------------------------------------------- # - # Plotting TYPE IIIs - # ---------------------------------------------------------------- # - - # WIND - fig, axes = plt.subplots(1, 1, sharex=True, figsize=(9, 13)) - waves_spec_lfr.plot(axes=axes,norm=LogNorm(vmin=waves_mm_l[0], vmax=waves_mm_l[1]), cmap="jet") - waves_spec_hfr.plot(axes=axes,norm=LogNorm(vmin=waves_mm_h[0], vmax=waves_mm_h[1]), cmap="jet") - - - # LEADING EDGE - if leadingedge ==True: - # axes.plot(waves_risetimes_l_LE, waves_testfreq_l_LE, 'k*') - # axes.plot(waves_risetimes_h_LE, waves_testfreq_h_LE, 'k*') - axes.plot(wavesrisetimes_LE, wavestestfreq_LE, 'k*') - axes.plot(fittimes_corrected_waves_LE,fitfreqs_waves_LE, "k--") - axes.plot(fittimes_corrected_waves_LE,fitfreqs_waves_LE, "y--") - # BACKBONE - if backbone ==True: - # axes.plot(waves_risetimes_l_BB, waves_testfreq_l_BB, 'k*') - # axes.plot(waves_risetimes_h_BB, waves_testfreq_h_BB, 'k*') - axes.plot(wavesrisetimes_BB, wavestestfreq_BB, 'k*') - axes.plot(fittimes_corrected_waves_BB,fitfreqs_waves_BB, "k--") - axes.plot(fittimes_corrected_waves_BB,fitfreqs_waves_BB, "y--") - - - - - axes.set_ylim(reversed(axes.get_ylim())) - axes.set_yscale('log') - axes.set_xlim(datetime(YYYY, MM, dd, 19,28), datetime(YYYY, MM, dd, HH_1,mm_1)) - axes.set_ylim([freqfitmax, freqfitmin]) - plt.subplots_adjust(hspace=0.31) - figfname = f"{figdir}/{YYYY}_{MM:02}_{dd:02}_WIND.png" - plt.savefig(figfname, dpi='figure') - - plt.show(block=False) - - - ### - # --------------------------------------------------------------------- - # FITTING TESTS - # --------------------------------------------------------------------- - if backbone == True: - testfreqs = wavestestfreq_BB - testtimes, testref = epoch2time(wavesrisetimes_BB) - - if leadingedge == True: - testfreqs = wavestestfreq_LE - testtimes,testref = epoch2time(wavesrisetimes_LE) - - popt_test, pcov_test = curve_fit(exponential_func3, testfreqs, testtimes) - times_fit = exponential_func3(np.array(testfreqs), *popt_test) - - - residuals = times_fit - testtimes - if plot_residuals == True: - fig, axs = plt.subplots(1, 2, figsize=(12, 8)) - - # Plot the data on each subplot - axs[0].plot(testtimes,testfreqs, 'k*') - axs[0].set_title('fit') - axs[0].set_xlabel('time [s]') - axs[0].set_ylabel('frequency MHz') - axs[0].plot(times_fit, testfreqs,'r--') - # axs[0].set_yscale("log") - axs[0].set_yscale("linear") - axs[0].semilogy(base=np.e) - - - - - axs[1].plot(residuals, 'r.') - axs[1].set_title('residuals') - - # Display the plot - plt.show(block=False) - - # --------------------------------------------------------------------- - # --------------------------------------------------------------------- - - ### - # STEREO A - fig, axes = plt.subplots(1, 1, sharex=True, figsize=(9, 13)) - swaves_a_spec_hfr.plot(axes=axes,norm=LogNorm(vmin=swaves_a_mm_h[0], vmax=swaves_a_mm_h[1]), cmap="jet") - # swaves_a_spec_lfr.plot(axes=axes) - - # LEADING EDGE - if leadingedge ==True: - axes.plot(swaves_a_risetimes_h_LE, swaves_a_testfreq_h_LE, 'r*') - axes.plot(fittimes_corrected_swaves_a_LE,fitfreqs_swaves_a_LE, "k--") - # BACKBONE - if backbone ==True: - axes.plot(swaves_a_risetimes_h_BB, swaves_a_testfreq_h_BB, 'r*') - axes.plot(fittimes_corrected_swaves_a_BB,fitfreqs_swaves_a_BB, "k--") - - axes.set_ylim(reversed(axes.get_ylim())) - axes.set_yscale('log') - axes.set_xlim(datetime(YYYY, MM, dd, HH_0,mm_0), datetime(YYYY, MM, dd, HH_1,mm_1)) - axes.set_ylim([freqfitmax, freqfitmin]) - plt.subplots_adjust(hspace=0.31) - figfname = f"{figdir}/{YYYY}_{MM:02}_{dd:02}_STEA.png" - plt.savefig(figfname, dpi='figure') - plt.show(block=False) - - - # STEREO B - fig, axes = plt.subplots(1, 1, sharex=True, figsize=(9, 13)) - swaves_b_spec_hfr.plot(axes=axes,norm=LogNorm(vmin=swaves_b_mm_h[0], vmax=swaves_b_mm_h[1]), cmap="jet") - # swaves_b_spec_lfr.plot(axes=axes) - - # LEADING EDGE - if leadingedge ==True: - axes.plot(swaves_b_risetimes_h_LE, swaves_b_testfreq_h_LE, 'ro',markeredgecolor="w") - axes.plot(fittimes_corrected_swaves_b_LE,fitfreqs_swaves_b_LE, "k--") - - if backbone ==True: - axes.plot(swaves_b_risetimes_h_BB, swaves_b_testfreq_h_BB, 'ro',markeredgecolor="w") - axes.plot(fittimes_corrected_swaves_b_BB,fitfreqs_swaves_b_BB, "k--") - - - axes.set_ylim(reversed(axes.get_ylim())) - axes.set_yscale('log') - axes.set_xlim(datetime(YYYY, MM, dd, HH_0,25), datetime(YYYY, MM, dd, HH_1,mm_1)) - axes.set_ylim([freqfitmax, freqfitmin]) - - plt.subplots_adjust(hspace=0.31) - # axes.set_ylim([freq4trimax,freq4trimin]) - figfname = f"{figdir}/{YYYY}_{MM:02}_{dd:02}_STEB.png" - plt.savefig(figfname, dpi='figure') - plt.show(block=False) - - - - # ---------------------------------------------------------------- # - # JOINT PLOT - # ---------------------------------------------------------------- # - # ---------------------------------------------------------------- # - # JOINT PLOT Horizontal - # ---------------------------------------------------------------- # - - fig, axes = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(20, 9)) - - waves_spec_lfr.plot(axes=axes[0], norm=LogNorm(vmin=waves_mm_l[0], vmax=waves_mm_l[1]), cmap="jet") - waves_spec_hfr.plot(axes=axes[0], norm=LogNorm(vmin=waves_mm_h[0], vmax=waves_mm_h[1]), cmap="jet") - swaves_a_spec_hfr.plot(axes=axes[1], norm=LogNorm(vmin=swaves_a_mm_h[0], vmax=swaves_a_mm_h[1]), cmap="jet") - # swaves_spec[1].plot(axes=axes[1], norm=LogNorm(vmin=swaves_mm_h[0], vmax=swaves_mm_h[1])) - swaves_b_spec_hfr.plot(axes=axes[2], norm=LogNorm(vmin=swaves_b_mm_h[0], vmax=swaves_b_mm_h[1]), cmap="jet") - - axes[0].set_title("WIND, WAVES, RAD1+RAD2") - axes[1].set_title("STEREO A, SWAVES, HFR") - axes[2].set_title("STEREO B, SWAVES, HFR") - - axes[0].set_ylabel("Frequency (MHz)") - - # # by default y-axis low to height flip so moving away from sun with time - axes[0].set_ylim(reversed(axes[0].get_ylim())) - axes[1].set_ylim(reversed(axes[1].get_ylim())) - axes[2].set_ylim(reversed(axes[2].get_ylim())) - - - # log y-axis - axes[0].set_yscale('log') - axes[1].set_yscale('log') - axes[2].set_yscale('log') - - if leadingedge == True: - times4tri_waves_LE_dt_plus30 = list(np.array(times4tri_waves_LE_dt)+timedelta(seconds=30)) - times4tri_waves_LE_dt_minus30 = list(np.array(times4tri_waves_LE_dt)-timedelta(seconds=30)) - - if backbone == True: - times4tri_waves_BB_dt_plus30 = list(np.array(times4tri_waves_BB_dt)+timedelta(seconds=30)) - times4tri_waves_BB_dt_minus30 = list(np.array(times4tri_waves_BB_dt)-timedelta(seconds=30)) - - - - # Models Leading Edge - if leadingedge ==True: - axes[0].plot(wavesrisetimes_LE, wavestestfreq_LE, c = "white", markeredgecolor="red", marker="o", ls="") - axes[1].plot(swaves_a_risetimes_h_LE, swaves_a_testfreq_h_LE, c = "white", markeredgecolor="red", marker="o", ls="") - axes[2].plot(swaves_b_risetimes_h_LE, swaves_b_testfreq_h_LE,c = "white", markeredgecolor="red", marker="o", ls="") - - - axes[0].plot(times4tri_waves_LE_dt_minus30, freqs4tri , "k--") - axes[0].plot(np.array(times4tri_waves_LE_dt), freqs4tri , "k-") - axes[0].plot(times4tri_waves_LE_dt_plus30, freqs4tri , "k--") - - axes[1].plot(np.array(times4tri_swaves_a_LE_dt)-timedelta(seconds=30), freqs4tri, "k--") - axes[1].plot(np.array(times4tri_swaves_a_LE_dt), freqs4tri, "k-") - axes[1].plot(np.array(times4tri_swaves_a_LE_dt)+timedelta(seconds=30), freqs4tri, "k--") - - axes[2].plot(np.array(times4tri_swaves_b_LE_dt)-timedelta(seconds=30), freqs4tri, "k--") - axes[2].plot(np.array(times4tri_swaves_b_LE_dt), freqs4tri, "k-") - axes[2].plot(np.array(times4tri_swaves_b_LE_dt)+timedelta(seconds=30), freqs4tri, "k--") - - # Models BACKBONE - if backbone ==True: - axes[0].plot(wavesrisetimes_BB, wavestestfreq_BB, c = "white", markeredgecolor="red", marker="o", ls="") - axes[1].plot(swaves_a_risetimes_h_BB, swaves_a_testfreq_h_BB, c = "white", markeredgecolor="red", marker="o", ls="") - axes[2].plot(swaves_b_risetimes_h_BB, swaves_b_testfreq_h_BB,c = "white", markeredgecolor="red", marker="o", ls="") - - - axes[0].plot(times4tri_waves_BB_dt_minus30, freqs4tri , "k--") - axes[0].plot(np.array(times4tri_waves_BB_dt), freqs4tri , "k-") - axes[0].plot(times4tri_waves_BB_dt_plus30, freqs4tri , "k--") - - axes[1].plot(np.array(times4tri_swaves_a_BB_dt)-timedelta(seconds=30), freqs4tri, "k--") - axes[1].plot(np.array(times4tri_swaves_a_BB_dt), freqs4tri, "k-") - axes[1].plot(np.array(times4tri_swaves_a_BB_dt)+timedelta(seconds=30), freqs4tri, "k--") - - axes[2].plot(np.array(times4tri_swaves_b_BB_dt)-timedelta(seconds=30), freqs4tri, "k--") - axes[2].plot(np.array(times4tri_swaves_b_BB_dt), freqs4tri, "k-") - axes[2].plot(np.array(times4tri_swaves_b_BB_dt)+timedelta(seconds=30), freqs4tri, "k--") - - - - - axes[0].set_ylim([freq4trimax,freq4trimin]) - # axes[1].set_ylim([freqlimmax, 0.2]) - # axes[2].set_ylim([freqlimmax, 0.2]) - - - axes[1].set_xlim(datetime(YYYY, MM, dd, HH_0,mm_0), datetime(YYYY, MM, dd, HH_1,mm_1)) - plt.subplots_adjust(left=0.041, bottom=0.096, right=0.984, top=0.93, wspace=0.132, hspace=0.31) - - plt.tick_params(axis='y', which='minor') - axes[0].yaxis.set_minor_formatter(FormatStrFormatter("%.1f")) - plt.tick_params(axis='y', which='major') - axes[0].yaxis.set_major_formatter(FormatStrFormatter("%.1f")) - figfname = f"{figdir}/{YYYY}_{MM:02}_{dd:02}_Horizontal.png" - plt.savefig(figfname, dpi='figure') - plt.show(block=False) - - - - - - if leadingedge == True: - wavesrisetimes_LE_stamp = datetime2timestamp(wavesrisetimes_LE) - spl = CubicSpline(wavestestfreq_LE, wavesrisetimes_LE_stamp) - - wavesrisetimes_LE_SWAVESfrequencies = timestamp2datetime(spl(swaves_a_testfreq_LE)) - wavesrisetimes_LE_SWAVESfrequencies = timestamp2datetime(interp(swaves_a_testfreq_LE,wavestestfreq_LE, wavesrisetimes_LE_stamp )) - - plt.figure() - plt.plot(times4tri_waves_LE_dt, freqs4tri, 'r--', label="WAVES") - plt.plot(times4tri_swaves_a_LE_dt, freqs4tri, 'b--', label="SWAVES_A") - plt.plot(times4tri_swaves_b_LE_dt, freqs4tri, 'g--', label="SWAVES_B") - - plt.plot(wavesrisetimes_LE, wavestestfreq_LE, 'r*') - plt.plot(wavesrisetimes_LE_SWAVESfrequencies, swaves_b_testfreq_LE, 'k.') # WAVES datapoints interpolated to match the SWAVES frequency channels - plt.plot(swaves_a_risetimes_LE, swaves_a_testfreq_LE, 'b*') - plt.plot(swaves_b_risetimes_LE, swaves_b_testfreq_LE, 'g*') - - plt.xlabel('epoch') - plt.ylabel('freq') - plt.yscale('log') - plt.legend(loc=1) - plt.title("Type IIIs datapoints combined") - plt.show(block=False) - - if backbone == True: - wavesrisetimes_BB_stamp = datetime2timestamp(wavesrisetimes_BB) - spl = CubicSpline(wavestestfreq_BB, wavesrisetimes_BB_stamp) - - wavesrisetimes_BB_SWAVESfrequencies = timestamp2datetime(spl(swaves_a_testfreq_BB)) - wavesrisetimes_BB_SWAVESfrequencies = timestamp2datetime( - interp(swaves_a_testfreq_BB, wavestestfreq_BB, wavesrisetimes_BB_stamp)) - - plt.figure() - plt.plot(times4tri_waves_BB_dt, freqs4tri, 'r--', label="WAVES") - plt.plot(times4tri_swaves_a_BB_dt, freqs4tri, 'b--', label="SWAVES_A") - plt.plot(times4tri_swaves_b_BB_dt, freqs4tri, 'g--', label="SWAVES_B") - - plt.plot(wavesrisetimes_BB, wavestestfreq_BB, 'r*') - plt.plot(wavesrisetimes_BB_SWAVESfrequencies, swaves_b_testfreq_BB, 'g.') - plt.plot(swaves_a_risetimes_BB, swaves_a_testfreq_BB, 'b*') - plt.plot(swaves_b_risetimes_BB, swaves_b_testfreq_BB, 'g*') - - plt.xlabel('epoch') - plt.ylabel('freq') - plt.yscale('log') - plt.legend(loc=1) - plt.show(block=False) - - typeIIIdir = mkdirectory(f"Data/TypeIII/{YYYY}_{MM:02}_{dd:02}/") - if leadingedge == True: - if (np.array_equal(fitfreqs, fitfreqs_waves_LE) and np.array_equal(fitfreqs,fitfreqs_swaves_a_LE) and np.array_equal(fitfreqs, fitfreqs_swaves_b_LE)): - typeIII_LE = ({ - 'Freqs': freqs4tri, - 'WindTime': times4tri_waves_LE_dt,#times4tri_waves_LE_dt_minus30, # - 'StereoATime': times4tri_swaves_a_LE_dt, - 'StereoBTime': times4tri_swaves_b_LE_dt, - }) - savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}_WIND_STEREO_A_STEREO_B_Freqs_{freq4trimin}_{freq4trimax}_LE_HR.pkl' - with open(savedfilepath, 'wb') as outp: - pickle.dump(typeIII_LE, outp, pickle.HIGHEST_PROTOCOL) - print(f"Saved results: {savedfilepath}") - - typeIII_LE = ({ - 'Freqs': swaves_a_testfreq_LE, - 'WindTime': wavesrisetimes_LE_SWAVESfrequencies,#times4tri_waves_LE_dt_minus30, # - 'StereoATime': swaves_a_risetimes_LE, - 'StereoBTime': swaves_b_risetimes_LE, - }) - savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}_WIND_STEREO_A_STEREO_B_Freqs_{freq4trimin}_{freq4trimax}_LE_HR_SCATTER.pkl' - with open(savedfilepath, 'wb') as outp: - pickle.dump(typeIII_LE, outp, pickle.HIGHEST_PROTOCOL) - print(f"Saved results: {savedfilepath}") - - else: - print("Missing frequencies in one of the spacecraft. Poor fit of the radio burst.") - if backbone == True: - if (np.array_equal(fitfreqs, fitfreqs_waves_BB) and np.array_equal(fitfreqs,fitfreqs_swaves_a_BB) and np.array_equal(fitfreqs, fitfreqs_swaves_b_BB)): - typeIII_BB = ({ - 'Freqs': freqs4tri, - 'WindTime': times4tri_waves_BB_dt,#times4tri_waves_BB_dt_minus30, # - 'StereoATime': times4tri_swaves_a_BB_dt, - 'StereoBTime': times4tri_swaves_b_BB_dt, - }) - savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}_WIND_STEREO_A_STEREO_B_Freqs_{freq4trimin}_{freq4trimax}_BB_HR.pkl' - with open(savedfilepath, 'wb') as outp: - pickle.dump(typeIII_BB, outp, pickle.HIGHEST_PROTOCOL) - print(f"Saved results: {savedfilepath}") - - typeIII_BB = ({ - 'Freqs': swaves_a_testfreq_BB, - 'WindTime': wavesrisetimes_BB_SWAVESfrequencies,#times4tri_waves_BB_dt_minus30, # - 'StereoATime': swaves_a_risetimes_BB, - 'StereoBTime': swaves_b_risetimes_BB, - }) - savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}_WIND_STEREO_A_STEREO_B_Freqs_{freq4trimin}_{freq4trimax}_BB_HR_SCATTER.pkl' - with open(savedfilepath, 'wb') as outp: - pickle.dump(typeIII_BB, outp, pickle.HIGHEST_PROTOCOL) - print(f"Saved results: {savedfilepath}") - - else: - print("Missing frequencies in one of the spacecraft. Poor fit of the radio burst.") diff --git a/examples/stacked_dyn_spectra_2021_12_04_LE_BB_TE_BURST_A.py b/examples/stacked_dyn_spectra_2021_12_04_LE_BB_TE_BURST_A.py new file mode 100755 index 0000000..3f9fef1 --- /dev/null +++ b/examples/stacked_dyn_spectra_2021_12_04_LE_BB_TE_BURST_A.py @@ -0,0 +1,1617 @@ +# Author: L Alberto Canizares canizares (at) cp.dias.ie +# some_file.py +import sys +sys.path.insert(1, '/Users/canizares/Library/CloudStorage/OneDrive-Personal/Work/0_PhD/Projects/BELLA_Projects/TCDSolarBELLA') + + +from datetime import datetime + +import matplotlib as mpl +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.colors import LogNorm + +import astropy.units as u +from astropy.visualization import ImageNormalize + + +from matplotlib.dates import DateFormatter +from matplotlib.ticker import FormatStrFormatter + +from astropy.constants import R_sun, au + + +import importlib + +from radiospectra.spectrogram import Spectrogram # in the process of updating old spectrogram + +from sunpy.net import attrs as a + +from bella.type_III_fitter import dynspec, openmarsis +from bella.type_III_fitter import typeIIIfitter as t3f +from bella.type_III_fitter.psp_quicklook import rfs_spec +from bella.type_III_fitter.solo_quicklook_L3_data import open_rpw_l3 +# from spacepy import pycdf +# from rpw_mono.thr.hfr.reader import read_hfr_autoch +# from maser.data import Data + +import pickle +import argparse + +importlib.reload(t3f) +importlib.reload(dynspec) + +r_sun = R_sun.value +AU=au.value + +plt.rcParams.update({'font.size': 18}) +plt.rcParams.update({'font.family': "Times New Roman"}) + +def reverse_colourmap(cmap, name='my_cmap_r'): + """ + In: + cmap, name + Out: + my_cmap_r + + Explanation: + t[0] goes from 0 to 1 + row i: x y0 y1 -> t[0] t[1] t[2] + / + / + row i+1: x y0 y1 -> t[n] t[1] t[2] + + so the inverse should do the same: + row i+1: x y1 y0 -> 1-t[0] t[2] t[1] + / + / + row i: x y1 y0 -> 1-t[n] t[2] t[1] + """ + reverse = [] + k = [] + + for key in cmap._segmentdata: + k.append(key) + channel = cmap._segmentdata[key] + data = [] + + for t in channel: + data.append((1 - t[0], t[2], t[1])) + reverse.append(sorted(data)) + + LinearL = dict(zip(k, reverse)) + my_cmap_r = mpl.colors.LinearSegmentedColormap(name, LinearL) + return my_cmap_r + + + +if __name__=="__main__": + plt.close('all') + parser = argparse.ArgumentParser(description='Process spacecraft names.') + parser.add_argument('--spacecraft', nargs='+', help='Spacecraft name(s)') + + args = parser.parse_args() + + #plt.close('all') + + # ---------------------------------------------------------------- # + " Settings " + # ---------------------------------------------------------------- # + YYYY = 2021 + MM = 12 + dd = 4 + HH_0 = 13 + mm_0 = 0 + HH_1 = 13 + mm_1 = 20 + # HH_1 = 14 + # mm_1 = 52 + # + background_subtraction = True + plot_residuals = False + savedata = True + savefigs = True + showfigs = True + + sat_hist = False # Saturated histograms + + stokesV = True + leadingedge = True + backbone = False + trailingedge = False + + + # spacecraft = ["stereo_a", "solo", "psp", "mex"] + spacecraft = [] + if not args.spacecraft: + # Manually inputing Spacecraft: + windsc = True + steasc = True + stebsc = False + solosc = True + pspsc = True + mexsc = True + + if windsc: + spacecraft.append("wind") + if steasc: + spacecraft.append("stereo_a") + if stebsc: + spacecraft.append("stereo_b") + if solosc: + spacecraft.append("solo") + if pspsc: + spacecraft.append("psp") + if mexsc: + spacecraft.append("mex") + else: + for sc in args.spacecraft: + spacecraft.append(sc.lower()) + + if "wind" in spacecraft: + windsc=True + else: + windsc=False + + if "stereo_a" in spacecraft: + steasc=True + else: + steasc=False + + if "stereo_b" in spacecraft: + stebsc=True + else: + stebsc=False + + if "solo" in spacecraft: + solosc=True + else: + solosc=False + + if "psp" in spacecraft: + pspsc=True + else: + pspsc=False + + if "mex" in spacecraft: + mexsc=True + else: + mexsc=False + + print(f"Spacecraft selected: {spacecraft}") + + # Cadence compensation or shifts + wind_cad_comp = 0 + stea_cad_comp = 0 + steb_cad_comp = 0 + solo_cad_comp = 0 + psp_cad_comp = 0 + mex_cad_comp = 0 + + local_data_dir = "/Users/canizares/Library/CloudStorage/OneDrive-Personal/Work/0_PhD/Projects/BELLA_Projects/2021_12_04" # Path to data directory + + # Wind Spacecraft LOCAL DATA ON or OFF + wind_local = True + if wind_local is True: + w_loc_file_RAD1 = f"{local_data_dir}/wind_data/waves/wi_wa_rad1_l2_60s_20211204_v01.dat" # PATH to wind data + w_loc_file_RAD2 = f"{local_data_dir}/wind_data/waves/wi_wa_rad2_l2_60s_20211204_v01.dat" # PATH to wind data + + sc_str="" + if windsc: + sc_str += '_WIND' + if steasc: + sc_str += '_STEA' + if stebsc: + sc_str += '_STEB' + if solosc: + sc_str += '_SOLO' + if pspsc: + sc_str += '_PSP' + if mexsc: + sc_str += '_MEX' + + profile = "" + if leadingedge: + profile += "_LE" + if backbone: + profile += "_BB" + if trailingedge: + profile += "_TE" + profile = profile[1:] # remove underscore from beginning + + note_str = "BURST_A" # This particular timeframe shows two type IIIs. Burst A is the earlier one. Burst B is the later one. + figdir = t3f.mkdirectory(f"Figures/{YYYY}_{MM:02}_{dd:02}/{profile}") + + my_cmap = mpl.cm.jet # Dynamic spectra colormap + + # ----------------------------------- # + " TIMING LIMIT SETTINGS " + # ----------------------------------- # + + mintime = datetime(YYYY, MM, dd, HH_0,mm_0) + maxtime = datetime(YYYY, MM, dd, HH_1,mm_1) + timelims = [mintime,maxtime] + + # TIME LIMITS SPECIFIC TO EACH SPACECRAFT + if windsc: + mintime_windsc = datetime(YYYY, MM, dd, HH_0, mm_0) + maxtime_windsc = datetime(YYYY, MM, dd, HH_1, mm_1) + timelims_windsc = [mintime_windsc, maxtime_windsc] + if steasc: + mintime_steasc = datetime(YYYY, MM, dd, 13, 0) + maxtime_steasc = datetime(YYYY, MM, dd, 15, 30) + timelims_steasc = [mintime_steasc, maxtime_steasc] + if stebsc: + mintime_stebsc = datetime(YYYY, MM, dd, HH_0, mm_0) + maxtime_stebsc = datetime(YYYY, MM, dd, HH_1, mm_1) + timelims_stebsc = [mintime_stebsc, maxtime_stebsc] + if solosc: + mintime_solosc = datetime(YYYY, MM, dd, 13, 0) + maxtime_solosc = datetime(YYYY, MM, dd, 15, 30) + timelims_solosc= [mintime_solosc, maxtime_solosc] + if pspsc: + mintime_pspsc = datetime(YYYY, MM, dd, 12, 50) + maxtime_pspsc = datetime(YYYY, MM, dd, 13, 50) + timelims_pspsc = [mintime_pspsc, maxtime_pspsc] + if mexsc: + mintime_mexsc = datetime(YYYY, MM, dd, 13, 9) # 2021-12-04 12:30:22.745000 + maxtime_mexsc = datetime(YYYY, MM, dd, 13, 14) # 2021-12-04 13:14:10.035000 + timelims_mexsc = [mintime_mexsc, maxtime_mexsc] + + mintime_mexsc_disp = datetime(YYYY, MM, dd, 13, 0) # 2021-12-04 12:30:22.745000 + maxtime_mexsc_disp = datetime(YYYY, MM, dd, 13, 30) # 2021-12-04 13:14:10.035000 + timelims_mexsc_disp = [mintime_mexsc_disp, maxtime_mexsc_disp] + + # ----------------------------------------- # + " Obtain DYNAMIC SPECTRA " + # ----------------------------------------- # + if windsc: + if wind_local is True: + print("Loading Wind Spectrogram LOCALLY") + print("----------") + # Waves RAD 1 ( low freqs) + waves_spec_lfr = dynspec.local_waves_spec_l2_60s(w_loc_file_RAD1, datatype='RAD1', kind='SMEAN', bg_subtraction=background_subtraction, lighttravelshift=0) + # Waves RAD 2 (High freqs) + waves_spec_hfr = dynspec.local_waves_spec_l2_60s(w_loc_file_RAD2, datatype='RAD2', kind='SMEAN', bg_subtraction=background_subtraction, lighttravelshift=0) + else: + print("Loading Wind Spectrogram with PYSPEDAS") + print("----------") + # Waves RAD 1 ( low freqs) + waves_spec_lfr = dynspec.waves_spec(mintime, maxtime, datatype='RAD1', bg_subtraction=background_subtraction) + # Waves RAD 2 (High freqs) + waves_spec_hfr = dynspec.waves_spec(mintime, maxtime, datatype='RAD2', bg_subtraction=background_subtraction) + + # FILL BLANK GAP + freqs_fill = [waves_spec_lfr.frequencies[-1].value, waves_spec_hfr.frequencies[0].value]* u.MHz + meta = { + 'observatory': "WIND_fill", + 'instrument': "WAVES_fill", + 'detector': "RAD2", + 'freqs': freqs_fill, + 'times': waves_spec_hfr.times, + 'wavelength': a.Wavelength(freqs_fill[0], freqs_fill[-1]), + 'start_time': waves_spec_hfr.times[0], + 'end_time': waves_spec_hfr.times[-1] + } + data_fill = np.array([waves_spec_hfr.data[0],waves_spec_hfr.data[0]]) + waves_spec_fill = Spectrogram(data_fill, meta) + print("Wind Spectrogram loaded") + print("----------\n") + if steasc: + print("Loading STEREO A Spectrogram") + print("----------") + # SWAVES A HFR + swaves_a_spec_hfr = dynspec.swaves_highres_spec(mintime, maxtime, probe='a', datatype='hfr', bg_subtraction=True) + # SWAVES A LFR + swaves_a_spec_lfr = dynspec.swaves_highres_spec(mintime, maxtime, probe='a', datatype='lfr', bg_subtraction=True) + + if stokesV: + # STOKES V + swaves_a_POL_hfr = dynspec.swaves_stokesV(mintime, maxtime, probe='a', datatype='hfr', bg_subtraction=False) + # SWAVES A LFR + swaves_a_POL_lfr = dynspec.swaves_stokesV(mintime, maxtime, probe='a', datatype='lfr', bg_subtraction=False) + + + print("STEREO A Spectrogram loaded") + print("----------\n") + if stebsc: + print("Loading STEREO B Spectrogram") + print("----------") + # SWAVES B HFR + swaves_b_spec_hfr = dynspec.swaves_highres_spec(mintime, maxtime, probe='b', datatype='hfr', bg_subtraction=True) + # SWAVES A LFR + swaves_b_spec_lfr = dynspec.swaves_highres_spec(mintime, maxtime, probe='b', datatype='lfr', bg_subtraction=True) + print("STEREO B Spectrogram loaded") + print("----------\n") + if solosc: + print("Loading SolO Spectrogram") + print("----------") + # RPW SOLO + cdf_file_path = f'{local_data_dir}/solo_data/RPW/l3/solo_L3_rpw-hfr-flux_20211204_V01.cdf' + rpw_spec_hfr = open_rpw_l3(cdf_file_path, bg_subtraction=True, lighttravelshift=solo_cad_comp) + + cdf_file_path = f"{local_data_dir}/solo_data/RPW/l3/solo_L3_rpw-tnr-flux_20211204_V01.cdf" + rpw_spec_tnr = open_rpw_l3(cdf_file_path, bg_subtraction=True, lighttravelshift=solo_cad_comp) + # rpw_mm_tnr = ImageNormalize(rpw_spec_tnr.data, interval=PercentileInterval(97.5)) + print("SolO Spectrogram loaded") + print("----------\n") + if pspsc: + print("Loading PSP Spectrogram") + print("----------") + # RFS FIELDS PSP + rfs_spec_lfr = rfs_spec(mintime, maxtime, datatype='rfs_lfr', bg_subtraction=True, lighttravelshift=0) + rfs_spec_hfr = rfs_spec(mintime, maxtime, datatype='rfs_hfr', bg_subtraction=True, lighttravelshift=0) + print("PSP Spectrogram loaded") + print("----------\n") + if mexsc: + print("Loading MEX MARSIS Spectrogram") + print("----------") + # MEX MARSIS + mex_filepath = f"{local_data_dir}/mex_data/marsis/22648/FRM_AIS_RDR_22648_ASCII_.DAT" # + mex_spec = openmarsis.marsis_spectra(mex_filepath) + # background = np.mean(mex_spec.data, axis=0) + # mex_spec.data = mex_spec.data - background + + print("MEX MARSIS Spectrogram loaded") + print("----------\n") + + + + # ----------------------------------------- # + " HISTOGRAM " + # ----------------------------------------- # + if windsc: + if sat_hist: + hist_l = np.percentile(waves_spec_lfr.data[~np.isnan(waves_spec_lfr.data)], [80,90]) + hist_h = np.percentile(waves_spec_hfr.data[~np.isnan(waves_spec_hfr.data)], [80,90]) + waves_mm_l = ImageNormalize(vmin=hist_l[0], vmax=hist_l[1]) + waves_mm_h = ImageNormalize(vmin=hist_h[0], vmax=hist_h[1]) + else: + waves_mm_l = np.percentile(waves_spec_lfr.data[~np.isnan(waves_spec_lfr.data)], [10,99.99]) + waves_mm_h = np.percentile(waves_spec_hfr.data[~np.isnan(waves_spec_hfr.data)], [10,99.99]) + if steasc: + if sat_hist: + hist_l = np.percentile(swaves_a_spec_lfr.data[~np.isnan(swaves_a_spec_lfr.data)], [90,99]) + hist_h = np.percentile(swaves_a_spec_hfr.data[~np.isnan(swaves_a_spec_hfr.data)], [90,99]) + swaves_a_mm_l = ImageNormalize(vmin=hist_l[0], vmax=hist_l[1]) + swaves_a_mm_h = ImageNormalize(vmin=hist_h[0], vmax=hist_h[1]) + else: + swaves_a_mm_l = np.percentile(swaves_a_spec_lfr.data[~np.isnan(swaves_a_spec_lfr.data)], [20, 99.2]) + swaves_a_mm_h = np.percentile(swaves_a_spec_hfr.data[~np.isnan(swaves_a_spec_hfr.data)], [20, 99.94]) + + if stokesV: + swaves_a_POL_lfr_data = swaves_a_POL_lfr.data[~np.isnan(swaves_a_POL_lfr.data)] + if len(swaves_a_POL_lfr_data)==0: + print(" WARNING : NO STOKES V DATA FOR SWAVES A LFR") + else: + swaves_a_pol_hist_l = np.percentile(swaves_a_POL_lfr_data, [30, 91]) + + swaves_a_POL_hfr_data = swaves_a_POL_hfr.data[~np.isnan(swaves_a_POL_hfr.data)] + if len(swaves_a_POL_hfr_data)==0: + print(" WARNING : NO STOKES V DATA FOR SWAVES A LFR") + else: + swaves_a_pol_hist_h = np.percentile(swaves_a_POL_hfr.data[~np.isnan(swaves_a_POL_hfr.data)], [50, 95]) + if stebsc: + if sat_hist: + hist_l = np.percentile(swaves_b_spec_lfr.data[~np.isnan(swaves_b_spec_lfr.data)], [80,90]) + hist_h = np.percentile(swaves_b_spec_hfr.data[~np.isnan(swaves_b_spec_hfr.data)], [80,90]) + swaves_b_mm_l = ImageNormalize(vmin=hist_l[0], vmax=hist_l[1]) + swaves_b_mm_h = ImageNormalize(vmin=hist_h[0], vmax=hist_h[1]) + + else: + swaves_b_mm_l = np.percentile(swaves_b_spec_lfr.data[~np.isnan(swaves_b_spec_lfr.data)], [20,99.8]) + swaves_b_mm_h = np.percentile(swaves_b_spec_hfr.data[~np.isnan(swaves_b_spec_hfr.data)], [20,99.8]) + if solosc: + if sat_hist: + hist_l = np.percentile(rpw_spec_tnr.data[~np.isnan(rpw_spec_tnr.data)], [80,90]) + hist_h = np.percentile(rpw_spec_hfr.data[~np.isnan(rpw_spec_hfr.data)], [80,90]) + rpw_mm_l = ImageNormalize(vmin=hist_l[0], vmax=hist_l[1]) + rpw_mm_h = ImageNormalize(vmin=hist_h[0], vmax=hist_h[1]) + else: + rpw_mm_l = np.percentile(rpw_spec_tnr.data[~np.isnan(rpw_spec_tnr.data)], [20,99.8]) + rpw_mm_h = np.percentile(rpw_spec_hfr.data[~np.isnan(rpw_spec_hfr.data)], [20,99.8]) + if pspsc: + if sat_hist: + rfs_mm_l = np.percentile(rfs_spec_lfr.data[~np.isnan(rfs_spec_lfr.data)], [70,99.2]) + rfs_mm_h = np.percentile(rfs_spec_hfr.data[~np.isnan(rfs_spec_hfr.data)], [70,99.2]) + else: + rfs_mm_l = np.percentile(rfs_spec_lfr.data[~np.isnan(rfs_spec_lfr.data)], [20,99.6]) + rfs_mm_h = np.percentile(rfs_spec_hfr.data[~np.isnan(rfs_spec_hfr.data)], [20,99.6]) + if mexsc: + filtered_data = mex_spec.data[~np.isnan(mex_spec.data) & (mex_spec.data > 0)] + if sat_hist: + mex_mm_h = np.percentile(filtered_data, [1,50]) + else: + mex_mm_h = np.percentile(filtered_data, [10,98]) + + + # ----------------------------------------- # + " FREQUENCY SETTINGS " + # ----------------------------------------- # + "" + # FREQUENCY RANGES + if leadingedge is True: + freqlimmin_LE = 0.5 # 0.18 + freqlimmax_LE = 10 + freqlims_LE = [freqlimmin_LE, freqlimmax_LE] + if backbone is True: + freqlimmin_BB = 0.14 # 0.18 + freqlimmax_BB = 15 + freqlims_BB = [freqlimmin_BB, freqlimmax_BB] + if trailingedge is True: + freqlimmin_TE = 0.4 # 0.18 + freqlimmax_TE = 3 + freqlims_TE = [freqlimmin_TE, freqlimmax_TE] + + freqlimmin = 0.18#0.18 + freqlimmax = 3 + freqlims = [freqlimmin, freqlimmax] + + if windsc: + freqlimmin_wind = 0.18 # 0.18 + freqlimmax_wind = 3 + freqlims_wind = [freqlimmin_wind, freqlimmax_wind] + + freqlimmin_wind_l = freqlimmin_wind + freqlimmax_wind_l = 1 + freqlims_wind_l = [freqlimmin_wind_l, freqlimmax_wind_l] + + freqlimmin_wind_h = 1 + freqlimmax_wind_h = freqlimmax_wind + freqlims_wind_h = [freqlimmin_wind_h, freqlimmax_wind_h] + if steasc: + freqlimmin_stea = 0.5 # 0.18 + freqlimmax_stea = 10 + freqlims_stea = [freqlimmin_stea, freqlimmax_stea] + + freqlimmin_stea_l = freqlimmin_stea + freqlimmax_stea_l = 1 + freqlims_stea_l = [freqlimmin_stea_l, freqlimmax_stea_l] + + freqlimmin_stea_h = 1 # 0.18 + freqlimmax_stea_h = freqlimmax_stea + freqlims_stea_h = [freqlimmin_stea_h, freqlimmax_stea_h] + if stebsc: + freqlimmin_steb = 0.5 # 0.18 + freqlimmax_steb = 10 + freqlims_steb = [freqlimmin_steb, freqlimmax_steb] + + freqlimmin_steb_l = freqlimmin_steb + freqlimmax_steb_l = 1 + freqlims_steb_l = [freqlimmin_steb_l, freqlimmax_steb_l] + + freqlimmin_steb_h = 1 # 0.18 + freqlimmax_steb_h = freqlimmax_steb + freqlims_steb_h = [freqlimmin_steb_h, freqlimmax_steb_h] + if solosc: + freqlimmin_solo = 0.5 # 0.18 + freqlimmax_solo = 5 + freqlims_solo = [freqlimmin_solo, freqlimmax_solo] + + freqlimmin_solo_l = freqlimmin_solo # 0.18 + freqlimmax_solo_l = 1 + freqlims_solo_l = [freqlimmin_solo_l, freqlimmax_solo_l] + + freqlimmin_solo_h = 1 + freqlimmax_solo_h = freqlimmax_solo + freqlims_solo_h = [freqlimmin_solo_h, freqlimmax_solo_h] + if pspsc: + freqlimmin_psp = 0.1 + freqlimmax_psp = 18 + freqlims_psp = [freqlimmin_psp, freqlimmax_psp] + + freqlimmin_psp_l = freqlimmin_psp + freqlimmax_psp_l = 1 + freqlims_psp_l = [freqlimmin_psp_l, freqlimmax_psp_l] + + freqlimmin_psp_h = 1 + freqlimmax_psp_h = freqlimmax_psp + freqlims_psp_h = [freqlimmin_psp_h, freqlimmax_psp_h] + if mexsc: + freqlimmin_mex = 0.05 + freqlimmax_mex = 10 + freqlims_mex = [freqlimmin_mex, freqlimmax_mex] + + freqlimmin_mex_l = freqlimmin_mex + freqlimmax_mex_l = 1 + freqlims_mex_l = [freqlimmin_mex_l, freqlimmax_mex_l] + + freqlimmin_mex_h = 1 + freqlimmax_mex_h = freqlimmax_mex + freqlims_mex_h = [freqlimmin_mex_h, freqlimmax_mex_h] + + # Input for fit + freqfitmin = 0.04 + freqfitmax = 50 + + # Frequencies to be extracted for triangulation + freq4trimin = 0.5 + freq4trimax = 3 + + # DISPLAY ONLY FREQUENCY LIMITS + freqlims_disp = [0.1, 15] + + + + + # ---------------------------------------------------------------- # + " FITTING Lightcurves " + # ---------------------------------------------------------------- # + "" + # MANUAL DATA EXTRACTION + # MANUAL points only because of Type III pair + print("Manual lightcurve data selection. \n ---------- \n") + freqs_manual = np.logspace(np.log10(freq4trimin), np.log10(freq4trimax), num=10) + # ------------------------- LEADING EDGE ----------------------------------------- # + if leadingedge is True: + if windsc: + # WAVES low + # waves_risetimes_l_LE =[] + # waves_riseval_l_LE = [] + # waves_testfreq_l_LE = [] + # + # # WAVES high + # waves_risetimes_h_LE = [datetime(2021, 12, 4, 13, 10, 5), + # datetime(2021, 12, 4, 13, 9, 8), + # datetime(2021, 12, 4, 13, 8, 9), + # datetime(2021, 12, 4, 13, 7, 23), + # datetime(2021, 12, 4, 13, 7, 6), + # datetime(2021, 12, 4, 13, 6, 44), + # datetime(2021, 12, 4, 13, 6, 36), + # datetime(2021, 12, 4, 13, 6, 30), + # datetime(2021, 12, 4, 13, 6, 22), + # datetime(2021, 12, 4, 13, 6, 6)] + # waves_riseval_h_LE = [] + # waves_testfreq_h_LE = np.logspace(np.log10(0.4), np.log10(3.6), num=10) + # + # waves_risetimes_h_LE.insert(10, datetime(2021, 12, 4, 13, 6, 12)) + # waves_testfreq_h_LE = np.insert(waves_testfreq_h_LE, 10, 10) + # + # # WAVES low + waves_risetimes_l_LE =[] + waves_riseval_l_LE = [] + waves_testfreq_l_LE = [] + + # waves_risetimes_l_LE.insert(0, datetime(2021, 12, 4, 13, 17, 31)) + # waves_testfreq_l_LE = np.insert(waves_testfreq_l_LE, 0, 0.14) + + # WAVES high + # waves_risetimes_h_LE = [datetime(2021, 12, 4, 13, 7, 30), + # datetime(2021, 12, 4, 13, 6, 31), + # datetime(2021, 12, 4, 13, 6, 9), + # datetime(2021, 12, 4, 13, 5, 23), + # datetime(2021, 12, 4, 13, 5, 6), + # datetime(2021, 12, 4, 13, 4, 59), + # datetime(2021, 12, 4, 13, 4, 51), + # datetime(2021, 12, 4, 13, 4, 46), + # datetime(2021, 12, 4, 13, 4, 36), + # datetime(2021, 12, 4, 13, 4, 29)] + waves_risetimes_h_LE = [datetime(2021, 12, 4, 13, 8, 48), + datetime(2021, 12, 4, 13, 7, 41), + datetime(2021, 12, 4, 13, 6, 39), + datetime(2021, 12, 4, 13, 5, 56), + datetime(2021, 12, 4, 13, 5, 34), + datetime(2021, 12, 4, 13, 5, 10), + datetime(2021, 12, 4, 13, 5, 9), + datetime(2021, 12, 4, 13, 4, 57), + datetime(2021, 12, 4, 13, 4, 51), + datetime(2021, 12, 4, 13, 4, 40)] + + waves_riseval_h_LE = [] + waves_testfreq_h_LE = np.logspace(np.log10(0.4), np.log10(3.6), num=10) + + waves_risetimes_h_LE.insert(len(waves_risetimes_h_LE), datetime(2021, 12, 4, 13, 4, 33)) + waves_testfreq_h_LE = np.insert(waves_testfreq_h_LE, len(waves_testfreq_h_LE), 10) + + # waves_risetimes_h_LE.insert(10, datetime(2021, 12, 4, 13, 4, 20)) + # waves_testfreq_h_LE = np.insert(waves_testfreq_h_LE, 10, 10) + if steasc: + # SWAVES A Low freqs + swaves_a_risetimes_l_LE = [] + swaves_a_riseval_l_LE = [] + swaves_a_testfreq_l_LE = [] + # SWAVES A High freqs + swaves_a_risetimes_h_LE = [ datetime(2021, 12, 4, 13, 14, 37), + datetime(2021, 12, 4, 13, 10, 27), + datetime(2021, 12, 4, 13, 8, 3), + datetime(2021, 12, 4, 13, 6, 54), + datetime(2021, 12, 4, 13, 5, 47), + datetime(2021, 12, 4, 13, 5, 23), + datetime(2021, 12, 4, 13, 5, 46), + datetime(2021, 12, 4, 13, 5, 15), + datetime(2021, 12, 4, 13, 5, 15), + datetime(2021, 12, 4, 13, 4, 54)] + + swaves_a_riseval_h_LE = [] + swaves_a_testfreq_h_LE = np.logspace(np.log10(0.23), np.log10(10), num=10) + if stebsc: + # SWAVES A Low freqs + swaves_b_risetimes_l_LE = [] + swaves_b_riseval_l_LE = [] + swaves_b_testfreq_l_LE = [] + # SWAVES B High freqs + swaves_b_risetimes_h_LE = [] + swaves_b_riseval_h_LE = [] + swaves_b_testfreq_h_LE = [] + if solosc: + # RPW TNR + rpw_risetimes_l_LE = [] + rpw_riseval_l_LE = [] + rpw_testfreq_l_LE = [] + + # RPW HFR + rpw_risetimes_h_LE = [datetime(2021, 12, 4, 13, 9, 20), + datetime(2021, 12, 4, 13, 8, 9), + datetime(2021, 12, 4, 13, 7, 12), + datetime(2021, 12, 4, 13, 6, 30), + datetime(2021, 12, 4, 13, 6, 8), + datetime(2021, 12, 4, 13, 5, 44), + datetime(2021, 12, 4, 13, 5, 36), + datetime(2021, 12, 4, 13, 5, 30), + datetime(2021, 12, 4, 13, 5, 22), + datetime(2021, 12, 4, 13, 5, 16)] + rpw_riseval_h_LE = [] + rpw_testfreq_h_LE = np.logspace(np.log10(0.4), np.log10(3.6), num=10) + + + + rpw_risetimes_h_LE.insert(10, datetime(2021, 12, 4, 13, 5, 6)) + rpw_testfreq_h_LE = np.insert(rpw_testfreq_h_LE, 10, 10) + + if pspsc: + # RFS LFR + rfs_risetimes_l_LE = [] + rfs_riseval_l_LE = [] + rfs_testfreq_l_LE = [] + # RFS HFR + rfs_risetimes_h_LE = [datetime(2021, 12, 4, 13, 7, 6), + datetime(2021, 12, 4, 13, 4, 51), + datetime(2021, 12, 4, 13, 3, 5), + datetime(2021, 12, 4, 13, 2, 19), + datetime(2021, 12, 4, 13, 2, 7), + datetime(2021, 12, 4, 13, 1, 42), + datetime(2021, 12, 4, 13, 1, 14), + datetime(2021, 12, 4, 13, 1, 8), + datetime(2021, 12, 4, 13, 0, 27), + datetime(2021, 12, 4, 13, 0, 16)] + rfs_riseval_h_LE = [] + rfs_testfreq_h_LE = np.logspace(np.log10(0.53), np.log10(10), num=10) + + # Correct fit at low frequency + # rfs_risetimes_h_LE.insert(0, datetime(2021, 12, 4, 13, 32, 58)) + # rfs_testfreq_h_LE = np.insert(rfs_testfreq_h_LE, 0, 0.11) + if mexsc: + # MEX MARSIS + mex_risetimes_l_LE = [] + mex_riseval_l_LE = [] + mex_testfreq_l_LE = [] + + mex_risetimes_h_LE = [datetime(2021, 12, 4, 13, 13, 47), + datetime(2021, 12, 4, 13, 12, 25), + datetime(2021, 12, 4, 13, 11, 36), + datetime(2021, 12, 4, 13, 10, 58), + datetime(2021, 12, 4, 13, 10, 38), + datetime(2021, 12, 4, 13, 10, 25), + datetime(2021, 12, 4, 13, 10, 10), + datetime(2021, 12, 4, 13, 10, 3), + datetime(2021, 12, 4, 13, 10, 0), + datetime(2021, 12, 4, 13, 9, 57)] + mex_riseval_h_LE = [] + mex_testfreq_h_LE = np.logspace(np.log10(0.5), np.log10(4.5), num=10) + # ------------------------- Backbone ----------------------------------------- # + if backbone is True: + if windsc: + # WAVES low + waves_risetimes_l_BB =[] + waves_riseval_l_BB = [] + waves_testfreq_l_BB = [] + + # WAVES high + waves_risetimes_h_BB = [] + waves_riseval_h_BB= [] + waves_testfreq_h_BB = [] + if steasc: + # SWAVES A LOW freqs + swaves_a_risetimes_l_BB = [] + swaves_a_riseval_l_BB = [] + swaves_a_testfreq_l_BB = [] + # SWAVES A High freqs + swaves_a_risetimes_h_BB = [] + swaves_a_riseval_h_BB = [] + swaves_a_testfreq_h_BB = [] + if stebsc: + # SWAVES B LOW freqs + swaves_b_risetimes_l_BB = [] + swaves_b_riseval_l_BB = [] + swaves_b_testfreq_l_BB = [] + # SWAVES B High freqs + swaves_b_risetimes_h_BB = [] + swaves_b_riseval_h_BB = [] + swaves_b_testfreq_h_BB = [] + if solosc: + # RPW TNR + rpw_risetimes_l_BB = [] + rpw_riseval_l_BB = [] + rpw_testfreq_l_BB = [] + + # RPW HFR + rpw_risetimes_h_BB = [] + rpw_riseval_h_BB = [] + rpw_testfreq_h_BB = [] + if pspsc: + # RFS LFR + rfs_risetimes_l_BB = [] + rfs_riseval_l_BB = [] + rfs_testfreq_l_BB = [] + # RFS HFR + rfs_risetimes_h_BB = [] + rfs_riseval_h_BB = [] + rfs_testfreq_h_BB = [] + if mexsc: + # MEX MARSIS + mex_risetimes_l_BB = [] + mex_riseval_l_BB = [] + mex_testfreq_l_BB = [] + + mex_risetimes_h_BB = [] + mex_riseval_h_BB = [] + mex_testfreq_h_BB = [] + # ------------------------- Trailing EDGE ----------------------------------------- # + if trailingedge is True: + if windsc: + # WAVES low + waves_risetimes_l_TE =[] + waves_riseval_l_TE = [] + waves_testfreq_l_TE = [] + + # WAVES high + waves_risetimes_h_TE = [] + waves_riseval_h_TE= [] + waves_testfreq_h_TE = [] + if steasc: + # SWAVES A LOW freqs + swaves_a_risetimes_l_TE = [] + swaves_a_riseval_l_TE = [] + swaves_a_testfreq_l_TE = [] + # SWAVES A High freqs + swaves_a_risetimes_h_TE = [] + swaves_a_riseval_h_TE = [] + swaves_a_testfreq_h_TE = [] + if stebsc: + # SWAVES B LOW freqs + swaves_b_risetimes_l_TE = [] + swaves_b_riseval_l_TE = [] + swaves_b_testfreq_l_TE = [] + # SWAVES B High freqs + swaves_b_risetimes_h_TE = [] + swaves_b_riseval_h_TE = [] + swaves_b_testfreq_h_TE = [] + if solosc: + # RPW TNR + rpw_risetimes_l_TE = [] + rpw_riseval_l_TE = [] + rpw_testfreq_l_TE = [] + + # RPW HFR + rpw_risetimes_h_TE = [] + rpw_riseval_h_TE = [] + rpw_testfreq_h_TE = [] + if pspsc: + # RFS LFR + rfs_risetimes_l_TE = [] + rfs_riseval_l_TE = [] + rfs_testfreq_l_TE = [] + # RFS HFR + rfs_risetimes_h_TE = [] + rfs_riseval_h_TE = [] + rfs_testfreq_h_TE = [] + if mexsc: + # MEX MARSIS + mex_risetimes_l_TE = [] + mex_riseval_l_TE = [] + mex_testfreq_l_TE = [] + + mex_risetimes_h_TE = [] + mex_riseval_h_TE = [] + mex_testfreq_h_TE = [] + + + # ---------------------------------------------------------------- # + " Type III FITTING " + # ---------------------------------------------------------------- # + fitfreqs = np.logspace(np.log10(freqfitmin), np.log10(freqfitmax), num=50) + freqs4tri = np.logspace(np.log10(freq4trimin), np.log10(freq4trimax), num=50) + if leadingedge is True: + print("Leading edge fitting. \n ---------- \n") + + + if windsc: + # Wind + wavesrisetimes_LE = list(waves_risetimes_l_LE) + list(waves_risetimes_h_LE) + wavesriseval_LE = list(waves_riseval_l_LE) + list(waves_riseval_h_LE) + wavestestfreq_LE = list(waves_testfreq_l_LE) + list(waves_testfreq_h_LE) + + times4tri_waves_LE_dt, fitfreqs_waves_LE, fittimes_corrected_waves_LE = t3f.typeIIIfitting(wavesrisetimes_LE, + wavestestfreq_LE, + fitfreqs, freqs4tri, + plot_residuals=False) + if steasc: + swaves_a_risetimes_LE = list(swaves_a_risetimes_h_LE) + list(swaves_a_risetimes_l_LE) + swaves_a_riseval_LE = list(swaves_a_riseval_h_LE) + list(swaves_a_riseval_l_LE) + swaves_a_testfreq_LE = list(swaves_a_testfreq_h_LE) + list(swaves_a_testfreq_l_LE) + + times4tri_swaves_a_LE_dt, fitfreqs_swaves_a_LE, fittimes_corrected_swaves_a_LE = t3f.typeIIIfitting(swaves_a_risetimes_LE, + swaves_a_testfreq_LE, + fitfreqs, freqs4tri, + plot_residuals=False) + if stebsc: + swaves_b_risetimes_LE = list(swaves_b_risetimes_h_LE) + list(swaves_b_risetimes_l_LE) + swaves_b_riseval_LE = list(swaves_b_riseval_h_LE) + list(swaves_b_riseval_l_LE) + swaves_b_testfreq_LE = list(swaves_b_testfreq_h_LE) + list(swaves_b_testfreq_l_LE) + + times4tri_swaves_b_LE_dt, fitfreqs_swaves_b_LE, fittimes_corrected_swaves_b_LE = t3f.typeIIIfitting(swaves_b_risetimes_LE, + swaves_b_testfreq_LE, + fitfreqs, freqs4tri, + plot_residuals=False) + if solosc: + rpw_risetimes_LE = list(rpw_risetimes_h_LE) + list(rpw_risetimes_l_LE) + rpw_riseval_LE = list(rpw_riseval_h_LE) + list(rpw_riseval_l_LE) + rpw_testfreq_LE = list(rpw_testfreq_h_LE) +list(rpw_testfreq_l_LE) + + times4tri_rpw_LE_dt, fitfreqs_rpw_LE, fittimes_corrected_rpw_LE = t3f.typeIIIfitting(rpw_risetimes_LE, + rpw_testfreq_LE, + fitfreqs, freqs4tri, + plot_residuals=False) + if pspsc: + rfsrisetimes_LE = list(rfs_risetimes_l_LE) + list(rfs_risetimes_h_LE) + rfsriseval_LE = list(rfs_riseval_l_LE) + list(rfs_riseval_h_LE) + rfstestfreq_LE = list(rfs_testfreq_l_LE) + list(rfs_testfreq_h_LE) + + times4tri_rfs_LE_dt, fitfreqs_rfs_LE, fittimes_corrected_rfs_LE = t3f.typeIIIfitting(rfsrisetimes_LE, + rfstestfreq_LE, + fitfreqs, freqs4tri, + plot_residuals=False) + if mexsc: + # MEX MARSIS + mex_risetimes_LE = list(mex_risetimes_l_LE) + list(mex_risetimes_h_LE) + mex_riseval_LE = list(mex_riseval_l_LE) + list(mex_riseval_h_LE) + mex_testfreq_LE = list(mex_testfreq_l_LE) + list(mex_testfreq_h_LE) + + times4tri_mex_LE_dt, fitfreqs_mex_LE, fittimes_corrected_mex_LE = t3f.typeIIIfitting(mex_risetimes_LE, + mex_testfreq_LE, + fitfreqs, freqs4tri, + plot_residuals=False) + + print("Leading edge fitting. DONE. \n ---------- \n") + + # BACKBONE + if backbone is True: + print("Leading edge fitting. \n ---------- \n") + + if windsc: + # Wind + wavesrisetimes_BB = list(waves_risetimes_l_BB) + list(waves_risetimes_h_BB) + wavesriseval_BB = list(waves_riseval_l_BB) + list(waves_riseval_h_BB) + wavestestfreq_BB = list(waves_testfreq_l_BB) + list(waves_testfreq_h_BB) + + times4tri_waves_BB_dt, fitfreqs_waves_BB, fittimes_corrected_waves_BB = t3f.typeIIIfitting( + wavesrisetimes_BB, + wavestestfreq_BB, + fitfreqs, freqs4tri, + plot_residuals=False) + if steasc: + swaves_a_risetimes_BB = list(swaves_a_risetimes_h_BB) + list(swaves_a_risetimes_l_BB) + swaves_a_riseval_BB = list(swaves_a_riseval_h_BB) + list(swaves_a_riseval_l_BB) + swaves_a_testfreq_BB = list(swaves_a_testfreq_h_BB) + list(swaves_a_testfreq_l_BB) + + times4tri_swaves_a_BB_dt, fitfreqs_swaves_a_BB, fittimes_corrected_swaves_a_BB = t3f.typeIIIfitting( + swaves_a_risetimes_BB, + swaves_a_testfreq_BB, + fitfreqs, freqs4tri, + plot_residuals=False) + if stebsc: + swaves_b_risetimes_BB = list(swaves_b_risetimes_h_BB) + list(swaves_b_risetimes_l_BB) + swaves_b_riseval_BB = list(swaves_b_riseval_h_BB) + list(swaves_b_riseval_l_BB) + swaves_b_testfreq_BB = list(swaves_b_testfreq_h_BB) + list(swaves_b_testfreq_l_BB) + + times4tri_swaves_b_BB_dt, fitfreqs_swaves_b_BB, fittimes_corrected_swaves_b_BB = t3f.typeIIIfitting( + swaves_b_risetimes_BB, + swaves_b_testfreq_BB, + fitfreqs, freqs4tri, + plot_residuals=False) + if solosc: + rpw_risetimes_BB = list(rpw_risetimes_h_BB) + list(rpw_risetimes_l_BB) + rpw_riseval_BB = list(rpw_riseval_h_BB) + list(rpw_riseval_l_BB) + rpw_testfreq_BB = list(rpw_testfreq_h_BB) + list(rpw_testfreq_l_BB) + + times4tri_rpw_BB_dt, fitfreqs_rpw_BB, fittimes_corrected_rpw_BB = t3f.typeIIIfitting(rpw_risetimes_BB, + rpw_testfreq_BB, + fitfreqs, freqs4tri, + plot_residuals=False) + if pspsc: + rfsrisetimes_BB = list(rfs_risetimes_l_BB) + list(rfs_risetimes_h_BB) + rfsriseval_BB = list(rfs_riseval_l_BB) + list(rfs_riseval_h_BB) + rfstestfreq_BB = list(rfs_testfreq_l_BB) + list(rfs_testfreq_h_BB) + + times4tri_rfs_BB_dt, fitfreqs_rfs_BB, fittimes_corrected_rfs_BB = t3f.typeIIIfitting(rfsrisetimes_BB, + rfstestfreq_BB, + fitfreqs, freqs4tri, + plot_residuals=False) + if mexsc: + # MEX MARSIS + mex_risetimes_BB = list(mex_risetimes_l_BB) + list(mex_risetimes_h_BB) + mex_riseval_BB = list(mex_riseval_l_BB) + list(mex_riseval_h_BB) + mex_testfreq_BB = list(mex_testfreq_l_BB) + list(mex_testfreq_h_BB) + + times4tri_mex_BB_dt, fitfreqs_mex_BB, fittimes_corrected_mex_BB = t3f.typeIIIfitting(mex_risetimes_BB, + mex_testfreq_BB, + fitfreqs, freqs4tri, + plot_residuals=False) + + print("Leading edge fitting. DONE. \n ---------- \n") + + # TRAILING EDGE + if trailingedge is True: + print("Leading edge fitting. \n ---------- \n") + + if windsc: + # Wind + wavesrisetimes_TE = list(waves_risetimes_l_TE) + list(waves_risetimes_h_TE) + wavesriseval_TE = list(waves_riseval_l_TE) + list(waves_riseval_h_TE) + wavestestfreq_TE = list(waves_testfreq_l_TE) + list(waves_testfreq_h_TE) + + times4tri_waves_TE_dt, fitfreqs_waves_TE, fittimes_corrected_waves_TE = t3f.typeIIIfitting( + wavesrisetimes_TE, + wavestestfreq_TE, + fitfreqs, freqs4tri, + plot_residuals=False) + if steasc: + swaves_a_risetimes_TE = list(swaves_a_risetimes_h_TE) + list(swaves_a_risetimes_l_TE) + swaves_a_riseval_TE = list(swaves_a_riseval_h_TE) + list(swaves_a_riseval_l_TE) + swaves_a_testfreq_TE = list(swaves_a_testfreq_h_TE) + list(swaves_a_testfreq_l_TE) + + times4tri_swaves_a_TE_dt, fitfreqs_swaves_a_TE, fittimes_corrected_swaves_a_TE = t3f.typeIIIfitting( + swaves_a_risetimes_TE, + swaves_a_testfreq_TE, + fitfreqs, freqs4tri, + plot_residuals=False) + if stebsc: + swaves_b_risetimes_TE = list(swaves_b_risetimes_h_TE) + list(swaves_b_risetimes_l_TE) + swaves_b_riseval_TE = list(swaves_b_riseval_h_TE) + list(swaves_b_riseval_l_TE) + swaves_b_testfreq_TE = list(swaves_b_testfreq_h_TE) + list(swaves_b_testfreq_l_TE) + + times4tri_swaves_b_TE_dt, fitfreqs_swaves_b_TE, fittimes_corrected_swaves_b_TE = t3f.typeIIIfitting( + swaves_b_risetimes_TE, + swaves_b_testfreq_TE, + fitfreqs, freqs4tri, + plot_residuals=False) + if solosc: + rpw_risetimes_TE = list(rpw_risetimes_h_TE) + list(rpw_risetimes_l_TE) + rpw_riseval_TE = list(rpw_riseval_h_TE) + list(rpw_riseval_l_TE) + rpw_testfreq_TE = list(rpw_testfreq_h_TE) + list(rpw_testfreq_l_TE) + + times4tri_rpw_TE_dt, fitfreqs_rpw_TE, fittimes_corrected_rpw_TE = t3f.typeIIIfitting(rpw_risetimes_TE, + rpw_testfreq_TE, + fitfreqs, freqs4tri, + plot_residuals=False) + if pspsc: + rfsrisetimes_TE = list(rfs_risetimes_l_TE) + list(rfs_risetimes_h_TE) + rfsriseval_TE = list(rfs_riseval_l_TE) + list(rfs_riseval_h_TE) + rfstestfreq_TE = list(rfs_testfreq_l_TE) + list(rfs_testfreq_h_TE) + + times4tri_rfs_TE_dt, fitfreqs_rfs_TE, fittimes_corrected_rfs_TE = t3f.typeIIIfitting(rfsrisetimes_TE, + rfstestfreq_TE, + fitfreqs, freqs4tri, + plot_residuals=False) + if mexsc: + # MEX MARSIS + mex_risetimes_TE = list(mex_risetimes_l_TE) + list(mex_risetimes_h_TE) + mex_riseval_TE = list(mex_riseval_l_TE) + list(mex_riseval_h_TE) + mex_testfreq_TE = list(mex_testfreq_l_TE) + list(mex_testfreq_h_TE) + + times4tri_mex_TE_dt, fitfreqs_mex_TE, fittimes_corrected_mex_TE = t3f.typeIIIfitting(mex_risetimes_TE, + mex_testfreq_TE, + fitfreqs, freqs4tri, + plot_residuals=False) + + print("Leading edge fitting. DONE. \n ---------- \n") + + # ---------------------------------------------------------------- # + " Plotting TYPE IIIs " + # ---------------------------------------------------------------- # + if windsc: + # --------------------------------------------------------- + # WIND + # ---------------------------------------------------------- + fig, axes = plt.subplots(1, 1, sharex=True, figsize=(10, 9)) + if isinstance(waves_mm_l, ImageNormalize): + waves_spec_lfr.plot(axes=axes, norm=waves_mm_l, cmap=my_cmap) + else: + waves_spec_lfr.plot(axes=axes, norm=LogNorm(vmin=waves_mm_l[0], vmax=waves_mm_l[1]), cmap=my_cmap) + if isinstance(waves_mm_h, ImageNormalize): + waves_spec_hfr.plot(axes=axes, norm=waves_mm_h, cmap=my_cmap) + else: + waves_spec_hfr.plot(axes=axes, norm=LogNorm(vmin=waves_mm_h[0], vmax=waves_mm_h[1]), cmap=my_cmap) + + # LEADING EDGE + if leadingedge is True: + # axes.plot(waves_risetimes_l_LE, waves_testfreq_l_LE, 'k*') + # axes.plot(waves_risetimes_h_LE, waves_testfreq_h_LE, 'k*') + axes.plot(wavesrisetimes_LE, wavestestfreq_LE, 'k*') + axes.plot(fittimes_corrected_waves_LE, fitfreqs_waves_LE, "k--") + axes.plot(fittimes_corrected_waves_LE, fitfreqs_waves_LE, "y--") + + # axes.plot(rpw_risetimes_LE, rpw_testfreq_LE, 'k*') + # axes.plot(fittimes_corrected_rpw_LE, fitfreqs_rpw_LE, "k--") + # axes.plot(fittimes_corrected_rpw_LE, fitfreqs_rpw_LE, "y--") + # BACKBONE + if backbone is True: + # axes.plot(waves_risetimes_l_BB, waves_testfreq_l_BB, 'k*') + # axes.plot(waves_risetimes_h_BB, waves_testfreq_h_BB, 'k*') + axes.plot(wavesrisetimes_BB, wavestestfreq_BB, 'k*') + axes.plot(fittimes_corrected_waves_BB, fitfreqs_waves_BB, "k--") + axes.plot(fittimes_corrected_waves_BB, fitfreqs_waves_BB, "y--") + + # TRAILING EDGE + if trailingedge is True: + # axes.plot(waves_risetimes_l_TE, waves_testfreq_l_TE, 'k*') + # axes.plot(waves_risetimes_h_TE, waves_testfreq_h_TE, 'k*') + axes.plot(wavesrisetimes_TE, wavestestfreq_TE, 'k*') + axes.plot(fittimes_corrected_waves_TE, fitfreqs_waves_TE, "k--") + axes.plot(fittimes_corrected_waves_TE, fitfreqs_waves_TE, "y--") + + axes.set_ylim(reversed(axes.get_ylim())) + axes.set_yscale('log') + axes.set_xlim(mintime, maxtime) + # axes.set_ylim([freqlimmax, freqlimmin]) + plt.subplots_adjust(hspace=0.31) + if showfigs: + plt.show(block=False) + else: + plt.close(fig) + if steasc: + # --------------------------------------------------------------------- + # STEREO A + # --------------------------------------------------------------------- + # + fig, axes = plt.subplots(1, 1, sharex=True, figsize=(10, 9)) + + if isinstance(swaves_a_mm_l, ImageNormalize): + swaves_a_spec_lfr.plot(axes=axes, norm=swaves_a_mm_l, cmap=my_cmap) + else: + swaves_a_spec_lfr.plot(axes=axes, norm=LogNorm(vmin=swaves_a_mm_l[0], vmax=swaves_a_mm_l[1]), cmap=my_cmap) + + if isinstance(swaves_a_mm_h, ImageNormalize): + swaves_a_spec_hfr.plot(axes=axes, norm=swaves_a_mm_h, cmap=my_cmap) + else: + swaves_a_spec_hfr.plot(axes=axes, norm=LogNorm(vmin=swaves_a_mm_h[0], vmax=swaves_a_mm_h[1]), cmap=my_cmap) + + # LEADING EDGE + if leadingedge is True: + axes.plot(swaves_a_risetimes_h_LE, swaves_a_testfreq_h_LE, 'ro', markeredgecolor="w") + axes.plot(fittimes_corrected_swaves_a_LE, fitfreqs_swaves_a_LE, "k--") + + # BACKBONE + if backbone is True: + axes.plot(swaves_a_risetimes_h_BB, swaves_a_testfreq_h_BB, 'r*') + axes.plot(fittimes_corrected_swaves_a_BB, fitfreqs_swaves_a_BB, "k--") + + # TRAILING EDGE + if trailingedge is True: + axes.plot(swaves_a_risetimes_h_TE, swaves_a_testfreq_h_TE, 'r*') + axes.plot(fittimes_corrected_swaves_a_TE, fitfreqs_swaves_a_TE, "k--") + + axes.set_yscale('log') + axes.set_xlim(mintime_steasc, maxtime_steasc) + axes.set_ylim(reversed(axes.get_ylim())) + # axes.set_ylim([freqlimmax, freqlimmin]) + plt.subplots_adjust(hspace=0.31) + if showfigs: + plt.show(block=False) + else: + plt.close(fig) + if stokesV: + # STOKES V + fig, axes = plt.subplots(1, 1, sharex=True, figsize=(10, 9)) + swaves_a_POL_hfr.plot(axes=axes, vmin=swaves_a_pol_hist_h[0], vmax=swaves_a_pol_hist_h[1], cmap=my_cmap) + # swaves_a_POL_lfr.plot(axes=axes, norm=LogNorm(vmin=swaves_a_pol_hist_l[0], vmax=swaves_a_pol_hist_l[1]), cmap=my_cmap) + + # # LEADING EDGE + # if leadingedge is True: + # axes.plot(swaves_a_risetimes_h_LE, swaves_a_testfreq_h_LE, 'ro', markeredgecolor="w") + # axes.plot(fittimes_corrected_swaves_a_LE, fitfreqs_swaves_a_LE, "k--") + # + # # BACKBONE + # if backbone is True: + # axes.plot(swaves_a_risetimes_h_BB, swaves_a_testfreq_h_BB, 'r*') + # axes.plot(fittimes_corrected_swaves_a_BB, fitfreqs_swaves_a_BB, "k--") + # + # # TRAILING EDGE + # if trailingedge is True: + # axes.plot(swaves_a_risetimes_h_TE, swaves_a_testfreq_h_TE, 'r*') + # axes.plot(fittimes_corrected_swaves_a_TE, fitfreqs_swaves_a_TE, "k--") + + axes.set_yscale('log') + axes.set_xlim(mintime_steasc, maxtime_steasc) + axes.set_ylim(reversed(axes.get_ylim())) + # axes.set_ylim([freqlimmax, freqlimmin]) + plt.subplots_adjust(hspace=0.31) + if showfigs: + plt.show(block=False) + else: + plt.close(fig) + # freqs_pol = [2, 6] + # tlim_pol = [datetime(2021,12, 4, 13,12), datetime(2021,12, 4, 13,14)] + # fpol_idx = np.where(swaves_a_POL_hfr.frequencies.value>freqs_pol[0]) and np.where(swaves_a_POL_hfr.frequencies.valuetlim_pol[0]) and np.where(swaves_a_POL_hfr.times.value 1 else axes # Adjust for when there's only one condition + ax_mapping[key] = i # Assigning an index to the axes for easy access + n_axes += 1 + titles.append(title_spec) + + # Plot spectrogram data + for j, spec_data in enumerate(data): + mm_value = mm_values[j] if len(mm_values) > j else mm_values[0] # Adjusted condition for mm_values + if isinstance(mm_value, ImageNormalize): + spec_data.plot(axes=ax, norm=mm_value, cmap=my_cmap) + else: + spec_data.plot(axes=ax, norm=LogNorm(vmin=mm_value[0], vmax=mm_value[1]), cmap=my_cmap) + + ax.set_title(title_spec) + + # Plot model lines for the current condition if applicable + for model_type in models: + if key in models[model_type]: + fittimes, fitfreqs, times4tri, freqs4tri = models[model_type][key] + ax.plot(fittimes, fitfreqs, "r-", label=model_type + " fit") # Model fit line + ax.plot(times4tri, freqs4tri, "w-", lw=2, label=model_type + " data") # Model data points/triangles + + # # It might be useful to add a legend if you're plotting model lines + # if key in models.get('LE', {}) or key in models.get('BB', {}) or key in models.get('TE', {}): + # ax.legend(loc='upper right') + + for each in axes: + each.set_ylim(reversed(each.get_ylim())) + each.set_yscale('log') + + + axes[0].set_ylim([freqlims_disp[1],freqlims_disp[0]]) + # axes[1].set_ylim([freqlimmax, 0.2]) + # axes[2].set_ylim([freqlimmax, 0.2]) + date_format = DateFormatter('%H:%M') + axes[0].xaxis.set_major_formatter(date_format) + + axes[0].set_xlim(datetime(YYYY, MM, dd, HH_0,mm_0), datetime(YYYY, MM, dd, HH_1,mm_1)) + plt.subplots_adjust(left=0.05, bottom=0.096, right=0.984, top=0.93, wspace=0.1, hspace=0.31) + plt.tick_params(axis='y', which='minor') + # axes[0].yaxis.set_minor_formatter(FormatStrFormatter("%.1f")) + plt.tick_params(axis='y', which='major') + axes[0].yaxis.set_major_formatter(FormatStrFormatter("%.1f")) + axes[0].set_ylabel("Frequency (MHz)") + # plt.tight_layout() + if savefigs: + if note_str!="": + note = f"_{note_str}" + figfname = f"{figdir}/{YYYY}_{MM:02}_{dd:02}{sc_str}_Horizontal{note}.png" + plt.savefig(figfname, dpi='figure') + print(f"\nFigure saved: {figfname}") + + if showfigs: + plt.show(block=False) + else: + plt.close(fig) + # ---------------------------------------------------------------- # + "END JOINT PLOT" + # ---------------------------------------------------------------- # + + + "" + # ---------------------------------------------------------------- # + "Save the data" + # ---------------------------------------------------------------- # + if savedata is True: + typeIIIdir = t3f.mkdirectory(f"Data/TypeIII/{YYYY}_{MM:02}_{dd:02}/") + if leadingedge is True: + typeIII_LE = {'Freqs': freqs4tri} # Start with 'Freqs' key-value pair + # Conditional assignments based on the presence of spacecraft data + if windsc: + typeIII_LE['WindTime'] = times4tri_waves_LE_dt + if steasc: + typeIII_LE['StereoATime'] = times4tri_swaves_a_LE_dt + if stebsc: + typeIII_LE['StereoBTime'] = times4tri_swaves_b_LE_dt + if solosc: + typeIII_LE['SoloTime'] = times4tri_rpw_LE_dt + if pspsc: + typeIII_LE['PSPTime'] = times4tri_rfs_LE_dt + if mexsc: + typeIII_LE['MEXTime'] = times4tri_mex_LE_dt + if note_str=="": + savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}{sc_str}_Freqs_{freq4trimin}_{freq4trimax}_LE.pkl' + else: + savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}{sc_str}_Freqs_{freq4trimin}_{freq4trimax}_LE_{note_str}.pkl' + with open(savedfilepath, 'wb') as outp: + pickle.dump(typeIII_LE, outp, pickle.HIGHEST_PROTOCOL) + print(f"\nSaved results: {savedfilepath}") + if backbone is True: + typeIII_BB = {'Freqs': freqs4tri} # Start with 'Freqs' key-value pair + # Conditional assignments based on the presence of spacecraft data + if windsc: + typeIII_BB['WindTime'] = times4tri_waves_BB_dt + if steasc: + typeIII_BB['StereoATime'] = times4tri_swaves_a_BB_dt + if stebsc: + typeIII_BB['StereoBTime'] = times4tri_swaves_b_BB_dt + if solosc: + typeIII_BB['SoloTime'] = times4tri_rpw_BB_dt + if pspsc: + typeIII_BB['PSPTime'] = times4tri_rfs_BB_dt + if mexsc: + typeIII_BB['MEXTime'] = times4tri_mex_BB_dt + + if note_str=="": + savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}{sc_str}_Freqs_{freq4trimin}_{freq4trimax}_BB.pkl' + else: + savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}{sc_str}_Freqs_{freq4trimin}_{freq4trimax}_BB_{note_str}.pkl' + with open(savedfilepath, 'wb') as outp: + pickle.dump(typeIII_BB, outp, pickle.HIGHEST_PROTOCOL) + print(f"\nSaved results: {savedfilepath}") + if trailingedge is True: + typeIII_TE = {'Freqs': freqs4tri} # Start with 'Freqs' key-value pair + # Conditional assignments based on the presence of spacecraft data + if windsc: + typeIII_TE['WindTime'] = times4tri_waves_TE_dt + if steasc: + typeIII_TE['StereoATime'] = times4tri_swaves_a_TE_dt + if stebsc: + typeIII_TE['StereoBTime'] = times4tri_swaves_b_TE_dt + if solosc: + typeIII_TE['SoloTime'] = times4tri_rpw_TE_dt + if pspsc: + typeIII_TE['PSPTime'] = times4tri_rfs_TE_dt + if mexsc: + typeIII_TE['MEXTime'] = times4tri_mex_TE_dt + if note_str=="": + savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}{sc_str}_Freqs_{freq4trimin}_{freq4trimax}_TE.pkl' + else: + savedfilepath = f'{typeIIIdir}typeIII_{YYYY}{MM:02}{dd:02}_{HH_0:02}{mm_0:02}{sc_str}_Freqs_{freq4trimin}_{freq4trimax}_TE_{note_str}.pkl' + with open(savedfilepath, 'wb') as outp: + pickle.dump(typeIII_TE, outp, pickle.HIGHEST_PROTOCOL) + print(f"\nSaved results: {savedfilepath}") + + else: + print("Warning. savedata is False, data was not saved.") + + + + """ --------------------------------------------------------------- """ + # VERTICAL PLOT + """ --------------------------------------------------------------- """ + + # Create subplots based on the dynamically determined number of axes + fig, axes = plt.subplots(len(conditions),1, sharex=True, sharey=True, figsize=(13, 13)) + + # Set face color for each subplot + cmap_lowest = my_cmap(0) + for ax in axes: + ax.set_facecolor(cmap_lowest) + + # Plot each spectrogram on the allocated axes + for i, (key, title_spec, data, mm_values) in enumerate(conditions): + ax = axes[i] if len(conditions) > 1 else axes # Adjust for when there's only one condition + ax_mapping[key] = i # Assigning an index to the axes for easy access + n_axes += 1 + titles.append(title_spec) + + # Plot spectrogram data + for j, spec_data in enumerate(data): + mm_value = mm_values[j] if len(mm_values) > j else mm_values[0] # Adjusted condition for mm_values + if isinstance(mm_value, ImageNormalize): + spec_data.plot(axes=ax, norm=mm_value, cmap=my_cmap) + else: + spec_data.plot(axes=ax, norm=LogNorm(vmin=mm_value[0], vmax=mm_value[1]), cmap=my_cmap) + + ax.set_title(title_spec) + + # Plot model lines for the current condition if applicable + for model_type in models: + if key in models[model_type]: + fittimes, fitfreqs, times4tri, freqs4tri = models[model_type][key] + ax.plot(fittimes, fitfreqs, "r-", label=model_type + " fit") # Model fit line + ax.plot(times4tri, freqs4tri, "w-", lw=2, label=model_type + " data") # Model data points/triangles + + # # It might be useful to add a legend if you're plotting model lines + # if key in models.get('LE', {}) or key in models.get('BB', {}) or key in models.get('TE', {}): + # ax.legend(loc='upper right') + + for each in axes: + each.set_ylim(reversed(each.get_ylim())) + each.set_yscale('log') + + + axes[0].set_ylim([freqlims_disp[1],freqlims_disp[0]]) + # axes[1].set_ylim([freqlimmax, 0.2]) + # axes[2].set_ylim([freqlimmax, 0.2]) + + + axes[0].set_xlim(datetime(YYYY, MM, dd, HH_0,mm_0), datetime(YYYY, MM, dd, HH_1,mm_1)) + plt.subplots_adjust(left=0.041, bottom=0.096, right=0.984, top=0.93, wspace=0.132, hspace=0.31) + plt.tick_params(axis='y', which='minor') + axes[0].yaxis.set_minor_formatter(FormatStrFormatter("%.1f")) + plt.tick_params(axis='y', which='major') + axes[0].yaxis.set_major_formatter(FormatStrFormatter("%.1f")) + axes[0].set_ylabel("Frequency (MHz)") + plt.tight_layout() + if savefigs: + if note_str!="": + note = f"_{note_str}" + figfname = f"{figdir}/{YYYY}_{MM:02}_{dd:02}{sc_str}_Vertical{note}.png" + plt.savefig(figfname, dpi='figure') + print(f"\nFigure saved: {figfname}") + + if showfigs: + plt.show(block=False) + else: + plt.close(fig) + + + report_cadences = f""" + Report CADENCES + Wind: + lfr:{dynspec.check_spectro_cadence(waves_spec_lfr).total_seconds()}s + hfr:{dynspec.check_spectro_cadence(waves_spec_lfr).total_seconds()}s + STEREO: + + """ + + for each in conditions: + print(f"{each[1]}:") + for spec in each[2]: + print(f"{dynspec.check_spectro_cadence(spec).total_seconds()}s") diff --git a/licenses/MIT.rst b/licenses/MIT.rst index 067dd37..ed12a52 100644 --- a/licenses/MIT.rst +++ b/licenses/MIT.rst @@ -18,4 +18,4 @@ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. \ No newline at end of file +SOFTWARE. diff --git a/pytest.ini b/pytest.ini index 8fb845d..aad2fa4 100644 --- a/pytest.ini +++ b/pytest.ini @@ -11,6 +11,7 @@ norecursedirs = *.egg-info examples bella/_dev + bella/type_III_fitter .history bella/extern doctest_plus = enabled @@ -31,3 +32,6 @@ filterwarnings = # A list of warnings to ignore follows. If you add to this list, you MUST # add a comment or ideally a link to an issue that explains why the warning # is being ignored + ignore::RuntimeWarning + ignore::DeprecationWarning + ignore::FutureWarning