Commit da2590c1 authored by RomanG's avatar RomanG
Browse files

Refactor R4,6,9,17. Create full_auto mode. Make tests.

parent 60b32b08
......@@ -13,5 +13,6 @@ rsc/img/*
to_process/
to_process_bkp/
results/
tests/to_process/
radiocc.log
radiocc.yaml
......@@ -47,6 +47,7 @@ logging.basicConfig(
datefmt=LOG_DATEFMT,
)
logging.getLogger("matplotlib").setLevel(logging.WARNING)
logging.getLogger("scipy").setLevel(logging.WARNING)
LOGGER = logging.getLogger(NAME)
LOGGER.info("Logging enabled.")
......
......@@ -29,6 +29,7 @@ class Cfg:
__DEFAULT_KERNEL_VERSION = KernelVersion.LATEST
__DEFAULT_RADIO_DATA_TYPE = RadioDataType.INGRESS
__DEFAULT_GRAPHICAL_INTERFACE = False
__DEFAULT_FULL_AUTO = False
# Path to the folder containing the data to be processed.
__to_process = __DEFAULT_TO_PROCESS
......@@ -57,6 +58,9 @@ class Cfg:
# Whether to use the graphical interface.
__graphical_interface = __DEFAULT_GRAPHICAL_INTERFACE
# Whether to run the code in full auto mode (no user interactions and no shown plots)
__full_auto = __DEFAULT_FULL_AUTO
def __init__(self) -> None:
pass
......@@ -95,6 +99,9 @@ class Cfg:
if CFG_FILE["graphical_interface"] not in (None, DotMap()):
self.graphical_interface = CFG_FILE["graphical_interface"]
if CFG_FILE["full_auto"] not in (None, DotMap()):
self.full_auto = CFG_FILE["full_auto"]
def to_dict(self, DEFAULT_ARE_NONE: bool = False) -> Dict[str, Optional[Any]]:
"""Convert config to hashmap."""
if DEFAULT_ARE_NONE:
......@@ -108,6 +115,7 @@ class Cfg:
kernel_version=None,
radio_data_type=None,
graphical_interface=None,
full_auto=None,
)
else:
return dict(
......@@ -120,6 +128,7 @@ class Cfg:
kernel_version=self.kernel_version.name,
radio_data_type=self.radio_data_type.name,
graphical_interface=self.graphical_interface,
full_auto=self.full_auto,
)
@property
......@@ -212,6 +221,16 @@ class Cfg:
"""Set whether to use the graphical interface."""
self.__graphical_interface = VALUE
@property
def full_auto(self) -> bool:
"""Get whether to run the code in full auto mode"""
return self.__full_auto
@full_auto.setter
def full_auto(self, VALUE: bool) -> None:
"""Set whether to run the code in full auto mode"""
self.__full_auto = VALUE
def generate_config(FORCE_OVERWRITE: bool = False) -> None:
"""Generate a config file `radiocc.yaml` in the current directory."""
......
......@@ -38,7 +38,7 @@ Ref_frame = "J2000"
Ref_frame_FP = "marsiau"
# Ref_frame = 'IAU_MARS'
# Constants for state vectors
DIMENSION_VECTOR = 3
DIMENSION_STATE = 2 * DIMENSION_VECTOR
POSITION_INDEX_START = 0
......
......@@ -6,6 +6,7 @@ Radio occultation.
import warnings
from pathlib import Path
from typing import Iterable, List # noqa: F401
from pudb import set_trace as bp # noqa: F401
......@@ -25,9 +26,13 @@ def enable_gui() -> None:
radiocc.gui.set_altitude(constants.Threshold_Cor)
def run() -> None:
def run() -> List[Scenario]:
"""Run radiocc."""
if radiocc.cfg.graphical_interface and radiocc.gui is None:
if (
radiocc.cfg.graphical_interface
and radiocc.gui is None
and not radiocc.cfg.full_auto
):
enable_gui()
check_cfg()
......@@ -48,6 +53,7 @@ def run() -> None:
if not len(TO_PROCESS_DIRS):
warnings.warn(f"{radiocc.cfg.to_process} directory is empty.", Warning)
SCENARIOS = []
for INDEX_PROCESS, PROCESS_PATH in enumerate(TO_PROCESS_DIRS):
print(f"Reading {PROCESS_PATH.name}..")
......@@ -63,10 +69,14 @@ def run() -> None:
band_export(SCENARIO)
SCENARIOS.append(SCENARIO)
if radiocc.gui is not None:
radiocc.gui.label_data.set_label("FINISHED!")
radiocc.gui.mainloop()
return SCENARIOS
def run_scenario(SCENARIO: Scenario) -> None:
"""Run a scenario."""
......
......@@ -6,33 +6,38 @@ Created on Sun Aug 30 20:17:48 2020
@author: ananya
"""
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as sc
from pudb import set_trace as bp # noqa: F401
import radiocc
from pudb import set_trace as bp #noqa: F401
def PLOT1( distance,Doppler,Doppler_debias,refractivity, Ne,PLOT_DIR,item):
def PLOT1(distance, Doppler, Doppler_debias, refractivity, Ne, PLOT_DIR, item):
if radiocc.gui is not None:
figure = radiocc.gui.figure
figure.clear()
else:
figure = plt.figure(6)
plt.ion()
plt.show()
if not radiocc.cfg.graphical_interface:
plt.ion()
plt.show()
Doppler=np.array(Doppler)
Doppler_debias=np.array(Doppler_debias)
cond_dop = Doppler<1
cond_deb = Doppler_debias<1
Doppler = np.array(Doppler)
Doppler_debias = np.array(Doppler_debias)
cond_dop = Doppler < 1
cond_deb = Doppler_debias < 1
axis_left = figure.add_subplot(131)
axis_middle = figure.add_subplot(132)
axis_right = figure.add_subplot(133)
figure.tight_layout(pad=3.0)
figure.suptitle("Radio Occultation")
figure.text(0.04, 0.5, "Altitude (km)", va="center", ha="center", rotation="vertical")
figure.text(
0.04, 0.5, "Altitude (km)", va="center", ha="center", rotation="vertical"
)
axis_left.plot(
Doppler, distance / 1000, ":", color="black", label="Raw L2 Data"
) # (1+880./749.)
......@@ -40,16 +45,27 @@ def PLOT1( distance,Doppler,Doppler_debias,refractivity, Ne,PLOT_DIR,item):
Doppler_debias, distance / 1000, ":", color="blue", label="Debias Raw L2 Data"
) # (1+880./749.)
Doppler = Doppler[cond_dop]
Doppler_debias = Doppler_debias[cond_deb]
min_x1 = np.nanmedian(Doppler_debias)- sc.median_absolute_deviation(Doppler_debias,nan_policy='omit')*5
min_x2 = np.nanmedian(Doppler)- sc.median_absolute_deviation(Doppler,nan_policy='omit')*5
max_x1 = np.nanmedian(Doppler_debias)+ sc.median_absolute_deviation(Doppler_debias,nan_policy='omit')*5
max_x2 = np.nanmedian(Doppler)+ sc.median_absolute_deviation(Doppler,nan_policy='omit')*5
min_x = min([min_x1,min_x2])
max_x = max([max_x1,max_x2])
min_x1 = (
np.nanmedian(Doppler_debias)
- sc.median_absolute_deviation(Doppler_debias, nan_policy="omit") * 5
)
min_x2 = (
np.nanmedian(Doppler)
- sc.median_absolute_deviation(Doppler, nan_policy="omit") * 5
)
max_x1 = (
np.nanmedian(Doppler_debias)
+ sc.median_absolute_deviation(Doppler_debias, nan_policy="omit") * 5
)
max_x2 = (
np.nanmedian(Doppler)
+ sc.median_absolute_deviation(Doppler, nan_policy="omit") * 5
)
min_x = min([min_x1, min_x2])
max_x = max([max_x1, max_x2])
axis_left.set_xlabel("Doppler Frequency Residual (Hz)")
axis_left.set_ylim(3400, 4000)
......@@ -71,7 +87,11 @@ def PLOT1( distance,Doppler,Doppler_debias,refractivity, Ne,PLOT_DIR,item):
axis_middle.grid(True)
axis_right.plot(
np.array(Ne)* 1e-6 , np.array(distance / 1000), ":", color="blue", label="after correction"
np.array(Ne) * 1e-6,
np.array(distance / 1000),
":",
color="blue",
label="after correction",
)
axis_right.set_xlabel("Electron density (el/cm$^{3}$)")
axis_right.set_xlim(-5.0e4, 5.0e4)
......@@ -82,41 +102,41 @@ def PLOT1( distance,Doppler,Doppler_debias,refractivity, Ne,PLOT_DIR,item):
if radiocc.gui is not None:
radiocc.gui.draw()
else:
elif not radiocc.cfg.full_auto:
plt.pause(1.0)
input("Press [enter] to continue.")
plt.close()
print("\t PLOT CHECK DONE")
return
print('\t PLOT CHECK DONE')
return
def PLOT2( distance,Doppler,Doppler_debias,refractivity, Nn,PLOT_DIR,item):
def PLOT2(distance, Doppler, Doppler_debias, refractivity, Nn, PLOT_DIR, item):
"""
Doppler=np.array(Doppler)
Doppler_debias=np.array(Doppler_debias)
cond_dop = Doppler<1
cond_dop = Doppler<1
cond_deb = Doppler_debias<1
fig6 = plt.figure(6)
fig6 = plt.figure(6)
plt.ion()
plt.show()
plt.plot(Doppler,distance/1000,':', color='black',label = 'Raw L2 Data')#(1+880./749.)
plt.plot(Doppler_debias,distance/1000,':', color='blue', label = 'Debias Raw L2 Data')#(1+880./749.)
Doppler = Doppler[cond_dop]
Doppler_debias = Doppler_debias[cond_deb]
min_x1 = np.nanmedian(Doppler_debias)- sc.median_absolute_deviation(Doppler_debias,nan_policy='omit')*5
min_x2 = np.nanmedian(Doppler)- sc.median_absolute_deviation(Doppler,nan_policy='omit')*5
max_x1 = np.nanmedian(Doppler_debias)+ sc.median_absolute_deviation(Doppler_debias,nan_policy='omit')*5
max_x2 = np.nanmedian(Doppler)+ sc.median_absolute_deviation(Doppler,nan_policy='omit')*5
max_x2 = np.nanmedian(Doppler)+ sc.median_absolute_deviation(Doppler,nan_policy='omit')*5
min_x = min([min_x1,min_x2])
max_x = max([max_x1,max_x2])
plt.xlabel('Doppler Frequency Residuals (Hz)')
plt.ylabel('Altitude (km)')
plt.ylim(3400, 4000)
plt.xlim(min_x, max_x)
plt.grid(True)
......@@ -126,24 +146,24 @@ def PLOT2( distance,Doppler,Doppler_debias,refractivity, Nn,PLOT_DIR,item):
plt.pause(1.0)
input("Press [enter] to continue.")
#plt.close()
fig10 = plt.figure(10)
fig10 = plt.figure(10)
plt.ion()
plt.show()
plt.show()
plt.plot(np.array(refractivity),np.array(distance/1000),':', color='blue', label ='after correction')
plt.xlabel('Refractivity')
plt.ylabel('Altitude (km)')
plt.ylim(3400,4000)
plt.xlim(-0.05,0.05)# S band
plt.grid(True)
fig10.set_size_inches(6, 10) #Inch
plt.savefig(PLOT_DIR+'/'+'refractivity-before-corr.png',dpi=600)
plt.pause(1.0)
plt.pause(1.0)
input("Press [enter] to continue.")
#plt.close()
"""
......@@ -153,29 +173,35 @@ def PLOT2( distance,Doppler,Doppler_debias,refractivity, Nn,PLOT_DIR,item):
figure.clear()
else:
figure = plt.figure(4)
plt.ion()
plt.show()
if not radiocc.cfg.full_auto:
plt.ion()
plt.show()
axis = figure.add_subplot(111)
figure.suptitle("Radio Occultation")
axis.plot(np.array(Nn), np.array(distance / 1000), ":", color="blue", label="after correction")
axis.plot(
np.array(Nn),
np.array(distance / 1000),
":",
color="blue",
label="after correction",
)
axis.set_xlabel("Neutral density (el/cm$^{3}$)")
axis.set_ylabel("Altitude (km)")
axis.grid(True)
#axis.set_xlim(xmin=-1e22, xmax=3e23)
# axis.set_xlim(xmin=-1e22, xmax=3e23)
# axis.ylim(3400, 4000)
figure.savefig(PLOT_DIR + "/" + "Nn_before-Corr.svg", dpi=100)
if radiocc.gui is not None:
radiocc.gui.draw()
else:
elif not radiocc.cfg.full_auto:
plt.pause(1.0)
input("Press [enter] to continue.")
plt.close()
print('\t PLOT CHECK DONE')
return
print("\t PLOT CHECK DONE")
return
......@@ -6,8 +6,25 @@ Created on Tue Sep 1 12:14:39 2020
@author: ananya
"""
import numpy as np
from pudb import set_trace as bp # noqa: F401
from radiocc.constants import (
CC,
G,
M_Planet,
R_Planet,
T_BC_low,
T_BC_med,
T_BC_upp,
Threshold_neut_upp,
c,
e,
eps0,
kB,
m_bar,
me,
)
from radiocc.old import (
R9_BendAng_ImpParam_dn,
R9_BendAng_ImpParam_up,
......@@ -28,45 +45,44 @@ def run(
item,
i_integral,
i_Surface,
R_Planet,
e,
eps0,
me,
PLOT_DIR,
Threshold_neut_upp,
CC,
kB,
N_data,
r_MEX_up,
z_MEX_up,
z_GS_up,
vr_MEX_up,
vz_MEX_up,
vr_GS_up,
vz_GS_up,
state_occ_SC_up,
state_occ_GS_up,
state_occ_SC_dn,
state_occ_GS_dn,
gamma_up,
beta_e_up,
delta_s_up,
r_MEX_dn,
z_MEX_dn,
z_GS_dn,
vr_MEX_dn,
vz_MEX_dn,
vr_GS_dn,
vz_GS_dn,
gamma_dn,
beta_e_dn,
delta_s_dn,
Bande,
fsup,
c,
G,
m_bar,
M_Mars,
T_BC_low,
T_BC_med,
T_BC_upp,
):
N_data = Doppler.shape[0]
r_MEX_up = state_occ_SC_up[:, 0]
z_MEX_up = state_occ_SC_up[:, 1]
vr_MEX_up = state_occ_SC_up[:, 2]
vz_MEX_up = state_occ_SC_up[:, 3]
z_GS_up = state_occ_GS_up[:, 1]
vr_GS_up = state_occ_GS_up[:, 2]
vz_GS_up = state_occ_GS_up[:, 3]
r_MEX_dn = state_occ_SC_dn[:, 0]
z_MEX_dn = state_occ_SC_dn[:, 1]
vr_MEX_dn = state_occ_SC_dn[:, 2]
vz_MEX_dn = state_occ_SC_dn[:, 3]
z_GS_dn = state_occ_GS_dn[:, 1]
vr_GS_dn = state_occ_GS_dn[:, 2]
vz_GS_dn = state_occ_GS_dn[:, 3]
beta_e_up = 0.5 * np.pi - delta_s_up
beta_e_dn = 0.5 * np.pi - delta_s_dn
M_Mars = M_Planet
kx = 880.0 / 749.0
ks = 240.0 / 749.0
Doppler_debias_dn_iono_x = Doppler_debias / (1.0 + (kx * kx))
......@@ -83,270 +99,97 @@ def run(
Diff_Doppler_debias_dn_iono_s = Diff_doppler / (1.0 - (9.0 / 121.0))
if Bande == "X":
if Bande != "Diff":
if item == "IONO":
if Bande == "X":
# R8: bending angle & impact parameter up
(
imp_param_up,
bend_ang_up,
delta_r_up,
beta_r_up,
) = R9_BendAng_ImpParam_up.Bend_up(
N_data,
r_MEX_up,
z_MEX_up,
z_GS_up,
vr_MEX_up,
vz_MEX_up,
vr_GS_up,
vz_GS_up,
gamma_up,
beta_e_up,
delta_s_up,
Doppler_debias_up_iono_x,
fsup,
c,
)
# R8: bending angle & impact parameter down
(
imp_param_dn,
bend_ang_dn,
delta_r_dn,
beta_r_dn,
) = R9_BendAng_ImpParam_dn.Bend_dn(
N_data,
r_MEX_dn,
z_MEX_dn,
z_GS_dn,
vr_MEX_dn,
vz_MEX_dn,
vr_GS_dn,
vz_GS_dn,
gamma_dn,
beta_e_dn,
delta_s_dn,
Doppler_debias_dn_iono_x,
Bande,
fsup,
c,
)
# R9: Average bending angle & impact parameter
imp_param, bend_ang = R10_Avg_BendAng_ImpParam.Avg(
imp_param_up, bend_ang_up, imp_param_dn, bend_ang_dn, Bande
)
if item == "IONO":
if item == "ATMO":
# R8: bending angle & impact parameter up
(
imp_param_up,
bend_ang_up,
delta_r_up,
beta_r_up,
) = R9_BendAng_ImpParam_up.Bend_up(
N_data,
r_MEX_up,
z_MEX_up,
z_GS_up,
vr_MEX_up,
vz_MEX_up,
vr_GS_up,
vz_GS_up,
gamma_up,
beta_e_up,
delta_s_up,
Doppler_debias_up_atmo_x,
fsup,
c,
)
# R8: bending angle & impact parameter down
(
imp_param_dn,
bend_ang_dn,
delta_r_dn,
beta_r_dn,
) = R9_BendAng_ImpParam_dn.Bend_dn(
N_data,
r_MEX_dn,
z_MEX_dn,
z_GS_dn,
vr_MEX_dn,
vz_MEX_dn,
vr_GS_dn,
vz_GS_dn,
gamma_dn,
beta_e_dn,
delta_s_dn,
Doppler_debias_dn_atmo_x,
Bande,
fsup,
c,
)
# R9: Average bending angle & impact parameter
imp_param, bend_ang = R10_Avg_BendAng_ImpParam.Avg(
imp_param_up, bend_ang_up, imp_param_dn, bend_ang_dn, Bande
)
if Bande == "S":
if item == "IONO":
# R8: bending angle & impact parameter up
(
imp_param_up,
bend_ang_up,
delta_r_up,
beta_r_up,
) = R9_BendAng_ImpParam_up.Bend_up(
N_data,
r_MEX_up,
z_MEX_up,
z_GS_up,
vr_MEX_up,
vz_MEX_up,
vr_GS_up,
vz_GS_up,
gamma_up,
beta_e_up,
delta_s_up,
Doppler_debias_up_iono_s,
fsup,
c,
)
# R8: bending angle & impact parameter down
(
imp_param_dn,
bend_ang_dn,
delta_r_dn,
beta_r_dn,
) = R9_BendAng_ImpParam_dn.Bend_dn(
N_data,
r_MEX_dn,
z_MEX_dn,
z_GS_dn,
vr_MEX_dn,
vz_MEX_dn,
vr_GS_dn,
vz_GS_dn,
gamma_dn,
beta_e_dn,
delta_s_dn,
Doppler_debias_dn_iono_s,
Bande,
fsup,
c,