From aeb25b49c9ed512e92952c6679ab7df570372dbf Mon Sep 17 00:00:00 2001 From: RomanG Date: Fri, 4 Feb 2022 19:14:51 +0100 Subject: [PATCH 01/39] Add general test + output check + Add config parameter: full-auto (which skips user interactions and blocks plot showing for testing) + radiocc.core.run() now returns list of Scenario that has been processed. + Apply full auto-mode in reconstruct (user interaction), R8 and R14 (to block plot showing) --- .gitignore | 2 + README.md | 2 +- pyproject.toml | 8 +- radiocc/__init__.py | 2 +- radiocc/config.py | 21 +++++- radiocc/core.py | 15 +++- radiocc/old/R14_Plot_check.py | 124 +++++++++++++++++++------------ radiocc/old/R8_Plot_Doppler_1.py | 2 +- radiocc/reconstruct.py | 19 +++-- setup.cfg | 2 +- tests/input.md | 3 + tests/test_radiocc.py | 68 ++++++++++++++++- 12 files changed, 199 insertions(+), 69 deletions(-) create mode 100644 tests/input.md diff --git a/.gitignore b/.gitignore index 3c52565..ce0565e 100644 --- a/.gitignore +++ b/.gitignore @@ -14,5 +14,7 @@ radiocc/assets/interface/icon.ico to_process/ to_process_bkp/ results/ +tests/to_process/ +tests/results radiocc.log radiocc.yaml diff --git a/README.md b/README.md index 4b4fb1b..776ca09 100644 --- a/README.md +++ b/README.md @@ -94,7 +94,7 @@ Licensed under the [Apache 2.0 license][license file]. [license badge]: https://img.shields.io/badge/License-Apache%202.0-blue.svg [coverage badge]: https://img.shields.io/badge/coverage-0%25-red [coverage url]: https://github.com/pytest-dev/pytest-cov -[version badge]: https://img.shields.io/badge/version-0.6.15-blue +[version badge]: https://img.shields.io/badge/version-0.6.16-blue [python url]: https://www.python.org/ [python badge]: https://img.shields.io/badge/python->=3.8,<3.11-blue [pre-commit url]: https://pre-commit.com diff --git a/pyproject.toml b/pyproject.toml index e5c962a..464e9f4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "radiocc" -version = "0.6.15" +version = "0.6.16" description = "Radio occultations" license = "Apache-2.0" authors = ["Ananya Krishnan ", "Greg Henry "] @@ -85,7 +85,7 @@ testpaths = "tests" [tool.mypy] files = ["radiocc", "tests"] -exclude = "radiocc/old" +exclude = "radiocc/old|personal" follow_imports = "silent" strict = true ignore_missing_imports = true @@ -94,8 +94,8 @@ show_error_codes = true [tool.black] line_length = 88 -exclude = "radiocc/old" +exclude = "radiocc/old|personal" [tool.isort] profile = "black" -skip = "radiocc/old" +skip = ["radiocc/old", "personal"] diff --git a/radiocc/__init__.py b/radiocc/__init__.py index 7d8e83e..cea6526 100644 --- a/radiocc/__init__.py +++ b/radiocc/__init__.py @@ -24,7 +24,7 @@ __all__ = [ # Project variables. NAME = "radiocc" -VERSION = "0.6.15" +VERSION = "0.6.16" # Logging settings. LOG_LEVEL = logging.DEBUG diff --git a/radiocc/config.py b/radiocc/config.py index 073d472..7464c91 100644 --- a/radiocc/config.py +++ b/radiocc/config.py @@ -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 + __full_auto = __DEFAULT_FULL_AUTO + def __init__(self) -> None: pass @@ -66,7 +70,7 @@ class Cfg: CFG_FILE = radiocc.utils.read_yaml(radiocc.CFG_PATH) self.__parse_config_file(CFG_FILE) - def __parse_config_file(self, CFG_FILE: DotMap) -> None: + def __parse_config_file(self, CFG_FILE: DotMap) -> None: # noqa:C901 """Apply variables from the config to the actual config.""" if CFG_FILE["to_process"] not in (None, DotMap()): self.to_process = Path(CFG_FILE["to_process"]) @@ -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.""" diff --git a/radiocc/core.py b/radiocc/core.py index 7fa45d4..7293482 100644 --- a/radiocc/core.py +++ b/radiocc/core.py @@ -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,10 +26,15 @@ 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: - enable_gui() + if ( + radiocc.cfg.graphical_interface + and radiocc.gui is None + and not radiocc.cfg.full_auto + ): + enable_gui() check_cfg() @@ -48,6 +54,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 +70,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.""" diff --git a/radiocc/old/R14_Plot_check.py b/radiocc/old/R14_Plot_check.py index aa9b971..a9dd5eb 100644 --- a/radiocc/old/R14_Plot_check.py +++ b/radiocc/old/R14_Plot_check.py @@ -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 diff --git a/radiocc/old/R8_Plot_Doppler_1.py b/radiocc/old/R8_Plot_Doppler_1.py index 77338ae..dd77815 100644 --- a/radiocc/old/R8_Plot_Doppler_1.py +++ b/radiocc/old/R8_Plot_Doppler_1.py @@ -121,7 +121,7 @@ def PLOT_Dop1(distance, Doppler, Doppler_debias, Doppler_biasfit, ET, PLOT_DIR, figure.savefig(PLOT_DIR + "/" + "Doppler_1.svg", dpi=100) if radiocc.gui is not None: radiocc.gui.draw() - else: + elif not radiocc.cfg.full_auto: plt.draw() plt.pause(1.0) input("Press [enter] to continue.") diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index 06cb497..bc4bea0 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -258,7 +258,12 @@ def run( # noqa: C901 radiocc.gui.go_next = False stay_while = not radiocc.gui.go_next else: - HAPPY_or_NOT = str(input("Enter 'YES' if the debias is good and 'NO' if bad: ")) + if radiocc.cfg.full_auto: + HAPPY_or_NOT = "YES" + else: + HAPPY_or_NOT = str( + input("Enter 'YES' if the debias is good and 'NO' if bad: ") + ) stay_while = HAPPY_or_NOT == "NO" while stay_while: # Updating the regression order or altitude value triggers an update @@ -461,11 +466,15 @@ def run( # noqa: C901 radiocc.gui.go_next = False stay_while = not radiocc.gui.go_next else: - RefCorr_GD_or_NOT = str( - input( - "Enter 'YES' if the refractivity correction is good and 'NO' if bad: " + if radiocc.cfg.full_auto: + RefCorr_GD_or_NOT = "YES" + else: + RefCorr_GD_or_NOT = str( + input( + "Enter 'YES' if the refractivity correction is good" + " and 'NO' if bad: " + ) ) - ) stay_while = RefCorr_GD_or_NOT == "NO" while stay_while: if radiocc.gui is not None: diff --git a/setup.cfg b/setup.cfg index 8f6875c..ad39be0 100644 --- a/setup.cfg +++ b/setup.cfg @@ -4,4 +4,4 @@ ignore = W503 max-complexity = 10 per-file-ignores = radiocc/__init__.py:F401 files = radiocc,tests -exclude = radiocc/old +exclude = radiocc/old,personal diff --git a/tests/input.md b/tests/input.md new file mode 100644 index 0000000..8dec92c --- /dev/null +++ b/tests/input.md @@ -0,0 +1,3 @@ +# Cases tested + ++ MEX-M-MRS-1-2-3-PRM-0235-V1.0 diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index 91a9d3e..7f67f5f 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -4,9 +4,69 @@ radiocc """ -# from radiocc.core import start # noqa: F401 +from pathlib import Path +import radiocc -def test() -> None: - """Test start""" - assert True +# Constants for tests +NUMBER_OF_TEXT_OUTPUTS = 2 +NUMBER_OF_FIGURE_OUTPUTS = 15 + + +def test_code_is_working() -> None: + """Test if the code runs without error""" + radiocc.cfg.to_process = Path("tests/to_process") + radiocc.cfg.results = Path("tests/results") + radiocc.cfg.full_auto = True + radiocc.core.run() + + +def test_number_of_text_outputs() -> None: + # Run the code a new type + # Lire avec Pathlib le nombre de fichiers dans le dossier results/DATA + radiocc.cfg.to_process = Path("tests/to_process") + radiocc.cfg.results = Path("tests/results") + radiocc.cfg.full_auto = True + SCENARIOS = iter(radiocc.core.run()) + SCENARIO = next(SCENARIOS) + TEXT_OUTPUTS = list( + radiocc.utils.directories( + SCENARIO.results(radiocc.cfg.results) / "DATA", "*.txt" + ) + ) + assert len(TEXT_OUTPUTS) == NUMBER_OF_TEXT_OUTPUTS + + +def test_number_of_figure_outputs() -> None: + radiocc.cfg.to_process = Path("tests/to_process") + radiocc.cfg.results = Path("tests/results") + radiocc.cfg.full_auto = True + SCENARIOS = iter(radiocc.core.run()) + SCENARIO = next(SCENARIOS) + FIGURE_OUTPUTS = list( + radiocc.utils.directories( + SCENARIO.results(radiocc.cfg.results) / "PLOTS", ["*.png", "*.svg"] + ) + ) + assert len(FIGURE_OUTPUTS) == NUMBER_OF_FIGURE_OUTPUTS + + +def test_number_of_lines_and_columns() -> None: + radiocc.cfg.to_process = Path("tests/to_process") + radiocc.cfg.results = Path("tests/results") + radiocc.cfg.full_auto = True + SCENARIOS = iter(radiocc.core.run()) + SCENARIO = next(SCENARIOS) + INPUTS = radiocc.utils.directories( + SCENARIO.TO_PROCESS / "DATA/LEVEL02/OPEN_LOOP/DSN/DPX", "*.TAB" + ) + INPUT = next(INPUTS) + NUMBER_OF_LINES_INPUT = len(open(INPUT).readlines()) + + TEXT_OUTPUTS = list( + radiocc.utils.directories( + SCENARIO.results(radiocc.cfg.results) / "DATA", "*.txt" + ) + ) + for TEXT_OUTPUT in TEXT_OUTPUTS: + assert len(open(TEXT_OUTPUT).readlines()) == NUMBER_OF_LINES_INPUT -- GitLab From e4e833da5293c1789d8b97900de2c4d5b8446478 Mon Sep 17 00:00:00 2001 From: RomanG Date: Sat, 5 Feb 2022 01:54:31 +0100 Subject: [PATCH 02/39] update test number of lines and columns in output + The test nows actually tests for the number of columns in the X-ATMO-1.txt and X-IONO-1.txt output files. --- tests/test_radiocc.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index 7f67f5f..bf97e08 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -6,11 +6,14 @@ radiocc from pathlib import Path +import numpy as np + import radiocc # Constants for tests -NUMBER_OF_TEXT_OUTPUTS = 2 -NUMBER_OF_FIGURE_OUTPUTS = 15 +NUMBER_OF_TEXT_OUTPUT = 2 +NUMBER_OF_FIGURE_OUTPUT = 15 +NUMBER_OF_COLUMNS_OUTPUT = {"ATMO": 16, "IONO": 11} def test_code_is_working() -> None: @@ -34,7 +37,7 @@ def test_number_of_text_outputs() -> None: SCENARIO.results(radiocc.cfg.results) / "DATA", "*.txt" ) ) - assert len(TEXT_OUTPUTS) == NUMBER_OF_TEXT_OUTPUTS + assert len(TEXT_OUTPUTS) == NUMBER_OF_TEXT_OUTPUT def test_number_of_figure_outputs() -> None: @@ -48,7 +51,7 @@ def test_number_of_figure_outputs() -> None: SCENARIO.results(radiocc.cfg.results) / "PLOTS", ["*.png", "*.svg"] ) ) - assert len(FIGURE_OUTPUTS) == NUMBER_OF_FIGURE_OUTPUTS + assert len(FIGURE_OUTPUTS) == NUMBER_OF_FIGURE_OUTPUT def test_number_of_lines_and_columns() -> None: @@ -63,10 +66,14 @@ def test_number_of_lines_and_columns() -> None: INPUT = next(INPUTS) NUMBER_OF_LINES_INPUT = len(open(INPUT).readlines()) - TEXT_OUTPUTS = list( - radiocc.utils.directories( - SCENARIO.results(radiocc.cfg.results) / "DATA", "*.txt" - ) - ) - for TEXT_OUTPUT in TEXT_OUTPUTS: + PATH_RESULTS = SCENARIO.results(radiocc.cfg.results) + TEXT_OUTPUTS = { + "ATMO": PATH_RESULTS / "DATA/X-ATMO-1.txt", + "IONO": PATH_RESULTS / "DATA/X-IONO-1.txt", + } + + for LAYER, TEXT_OUTPUT in TEXT_OUTPUTS.items(): + # Test if the number of input lines matches the number of output lines assert len(open(TEXT_OUTPUT).readlines()) == NUMBER_OF_LINES_INPUT + # Test if the X-ATMO and X-IONO output have the correct number of columns + assert np.loadtxt(TEXT_OUTPUT).shape[1] == NUMBER_OF_COLUMNS_OUTPUT[LAYER] -- GitLab From a749cc7ccb27fd7f90471e57594bf6d2b6fc44a1 Mon Sep 17 00:00:00 2001 From: RomanG Date: Sat, 5 Feb 2022 02:32:26 +0100 Subject: [PATCH 03/39] Add test for output consistency + The test compares all values of the produced text output to the expected text results stored in tests/expected_results. + It compares NaN's as equal. --- tests/test_radiocc.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index bf97e08..faf6989 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -77,3 +77,29 @@ def test_number_of_lines_and_columns() -> None: assert len(open(TEXT_OUTPUT).readlines()) == NUMBER_OF_LINES_INPUT # Test if the X-ATMO and X-IONO output have the correct number of columns assert np.loadtxt(TEXT_OUTPUT).shape[1] == NUMBER_OF_COLUMNS_OUTPUT[LAYER] + + +def test_output_consistency() -> None: + radiocc.cfg.to_process = Path("tests/to_process") + radiocc.cfg.results = Path("tests/results") + radiocc.cfg.full_auto = True + SCENARIOS = iter(radiocc.core.run()) + SCENARIO = next(SCENARIOS) + + PATH_RESULTS = SCENARIO.results(radiocc.cfg.results) + PATH_EXPECTED_RESULTS = Path("tests/expected_results") / SCENARIO.TO_PROCESS.name + TEXT_OUTPUTS = { + "ATMO": PATH_RESULTS / "DATA/X-ATMO-1.txt", + "IONO": PATH_RESULTS / "DATA/X-IONO-1.txt", + } + EXPECTED_TEXT_OUTPUTS = { + "ATMO": PATH_EXPECTED_RESULTS / "DATA/X-ATMO-1.txt", + "IONO": PATH_EXPECTED_RESULTS / "DATA/X-IONO-1.txt", + } + for TEXT_OUTPUT, EXPECTED_TEXT_OUTPUT in zip( + TEXT_OUTPUTS.values(), EXPECTED_TEXT_OUTPUTS.values() + ): + # Test if all values of the produced and expected output are equal + assert np.array_equal( + np.loadtxt(TEXT_OUTPUT), np.loadtxt(EXPECTED_TEXT_OUTPUT), equal_nan=True + ) -- GitLab From 2eceaa02cda8a08f625b1c18495b869db977724e Mon Sep 17 00:00:00 2001 From: RomanG Date: Sat, 5 Feb 2022 13:23:06 +0100 Subject: [PATCH 04/39] Refactor and make unit tests for R4. + R4 has been refactored using unit functions and has now typing annotations for mypy. + Unit tests for each corresponding unit function if R4 have been made + reconstruct has been adapted to work with the refactored R4 + Constants for the dimension of state vectors have been added in radiocc.constants --- radiocc/constants.py | 8 + radiocc/old/R4_State_Vectors_twoway_v2.py | 247 +++++----------------- radiocc/reconstruct.py | 51 +++-- tests/unit_tests_R4.py | 48 +++++ 4 files changed, 142 insertions(+), 212 deletions(-) create mode 100644 tests/unit_tests_R4.py diff --git a/radiocc/constants.py b/radiocc/constants.py index 285dee4..612fdfa 100644 --- a/radiocc/constants.py +++ b/radiocc/constants.py @@ -37,3 +37,11 @@ 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 +VELOCITY_INDEX_START = 3 +POSITION_INDICES = slice(POSITION_INDEX_START, POSITION_INDEX_START + DIMENSION_VECTOR) +VELOCITY_INDICES = slice(VELOCITY_INDEX_START, VELOCITY_INDEX_START + DIMENSION_VECTOR) diff --git a/radiocc/old/R4_State_Vectors_twoway_v2.py b/radiocc/old/R4_State_Vectors_twoway_v2.py index 5497878..91629e4 100644 --- a/radiocc/old/R4_State_Vectors_twoway_v2.py +++ b/radiocc/old/R4_State_Vectors_twoway_v2.py @@ -6,203 +6,66 @@ Created on Tue Nov 5 15:26:40 2019 @author: ananya """ -import os + +from typing import Tuple import numpy as np import spiceypy as spy +from nptyping import NDArray from pudb import set_trace as bp # noqa -from radiocc.model import ProcessType +import radiocc # noqa + + +def get_epoch(etobs: NDArray, obs: int, direct: str, targ: int) -> NDArray: + epoch_targ = np.full_like(etobs, np.nan) + for i, et in enumerate(etobs): + epoch_targ[i], _ = spy.ltime(et, obs, direct, targ) + return epoch_targ + + +def get_state(etobs: NDArray, obs: int, targ: int, ref: str, abcorr: str) -> NDArray: + state_targ = np.array(spy.spkezr(str(targ), etobs, ref, abcorr, str(obs))[0]) * 1e3 + return state_targ def ephemerides( - N_data, ET, SC, GS, Mars, DATA_DIR, Ref_frame, FOLDER_TYPE: ProcessType -): - # os.chdir(DATA_DIR) - # tB = np.zeros((N_data,)) - tMEX = np.zeros((N_data,)) - tA = np.zeros((N_data,)) - t0_up = np.zeros((N_data,)) - t0_dn = np.zeros((N_data,)) - t0_upp = np.zeros((N_data,)) - - # - # - # - # - # - # - # - # - # GS(to_up) ---> t_Mex ----> G(to_dn) - - ET_MEX_Mars = np.zeros((N_data,)) - - pos_GS_dn = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - vel_GS_dn = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - pos_MEX_dn = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - vel_MEX_dn = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - pos_GS_up = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - vel_GS_up = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - pos_MEX_up = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - vel_MEX_up = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - pos_GS_upp = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - vel_GS_upp = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - pos_MEX_upp = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - vel_MEX_upp = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - pos_MEX = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], float) - vel_MEX = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], float) - - # This routine computes the transmit (or receive) time of a signal at a specified target, given the receive (or transmit) time at a specified observer. The elapsed time between transmit and receive is also returned. - # Direction the signal travels ( "->" or "<-" ) - # etobs I Epoch of a signal at some observer - # obs I NAIF ID of some observer - # dir I Direction the signal travels ( "->" or "<-" ) - # targ I NAIF ID of the target object - # ettarg O Epoch of the signal at the target - # elapsd O Time between transmit and receipt of the signal - - # void ltime_c ( etobs, obs, ConstSpiceChar * dir, SpiceInt targ, - # output * ettarg, - # output * elapsd ) - - # spkezr Return the state (position and velocity) of a target body relative to an observing body, optionally corrected for light time (planetary aberration) and stellar aberration. - # void spkezr_c (targ, et, ref, abcorr, obs) output starg[6],*lt ) - - # SUBROUTINE LTIME ( Epoch of a signal at some observer, NAIF-id of some observer, direction the signal travel, AIF-id of the target object, Time between transmit and receipt of the signal ) - - for i in range(N_data): - # Downlink - # tB = ET[i] - [ET_MEX_GS, t_0] = spy.ltime( - ET[i], int(GS), "<-", int(SC) - ) # Epoch of the signal at MEX tdown - tMEX[ - i - ] = ET_MEX_GS # time at which MEX send signal to mars via downlink = time at which MEX received signal from Mars via uplink - - # [ET_Mars_GS, t_0] = spy.ltime(tMEX[i], int(Mars), '<-', int(SC)) # Epoch of the signal at Mars t down - - [ET_Mars_GS, t_0] = spy.ltime( - tMEX[i], int(SC), "->", int(Mars) - ) # Epoch of the signal at Mars t down - t0_dn[i] = ET_Mars_GS # The signal was transmitted at: t0_dn from Mars - - # [ET_Mars_GS, t_0] = spy.ltime(tMEX[i], int(Mars), '->', int(SC)) # Epoch of the signal at MARS tup - [ET_Mars_GS, t_0] = spy.ltime( - tMEX[i], int(SC), "<-", int(Mars) - ) # Epoch of the signal at Mars tup - t0_up[i] = ET_Mars_GS # time at which Mars send signal to MEX uplink - - # [ states1, lt1 ] = spy.spkezr( 'DSS-65', t0_dn[i] , 'marsiau', "NONE", 'mars' ) - [states1, lt1] = spy.spkezr( - GS, t0_dn[i], Ref_frame, "NONE", "mars" - ) # position of GS at T_Mars wrt Mars - pos_GS_dn[i] = states1[:3] * 1e3 - vel_GS_dn[i] = states1[3:] * 1e3 - - [states2, lt2] = spy.spkezr( - FOLDER_TYPE.name, t0_dn[i], Ref_frame, "NONE", "mars" - ) # position of MEX at T_Mars wrt Mars - pos_MEX_dn[i] = states2[:3] * 1e3 - vel_MEX_dn[i] = states2[3:] * 1e3 - - # uplink - - [states3, lt3] = spy.spkezr( - GS, t0_up[i], Ref_frame, "NONE", "mars" - ) # position of GS tup - pos_GS_up[i] = states3[:3] * 1e3 - vel_GS_up[i] = states3[3:] * 1e3 - - [states4, lt2] = spy.spkezr( - FOLDER_TYPE.name, t0_up[i], Ref_frame, "NONE", "mars" - ) # position of MEX tup - pos_MEX_up[i] = states4[:3] * 1e3 - vel_MEX_up[i] = states4[3:] * 1e3 - - [ET_GS_SC, t_0] = spy.ltime(tMEX[i], int(SC), "<-", int(GS)) - tA[i] = ET_GS_SC - - # Used for routine R6, not in R9 ! - - [ET_MEX_Marss, one_way_t_span] = spy.ltime(tMEX[i], int(Mars), "<-", int(SC)) - ET_MEX_Mars[i] = ET_MEX_Marss - - [states5, lt] = spy.spkezr( - FOLDER_TYPE.name, ET_MEX_Mars[i], Ref_frame, "NONE", "mars" - ) - pos_MEX[i] = states5[:3] * 1e3 - vel_MEX[i] = states5[3:] * 1e3 - - [ET_GS_Mars, t_0] = spy.ltime(tA[i], int(Mars), "<-", int(GS)) - t0_upp[i] = ET_GS_Mars - - [states6, lt6] = spy.spkezr(GS, t0_upp[i], Ref_frame, "NONE", "mars") - pos_GS_upp[i] = states6[:3] * 1e3 - vel_GS_upp[i] = states6[3:] * 1e3 - - [states7, lt7] = spy.spkezr( - FOLDER_TYPE.name, t0_upp[i], Ref_frame, "NONE", "mars" - ) - pos_MEX_upp[i] = states7[:3] * 1e3 - vel_MEX_upp[i] = states7[3:] * 1e3 - - # print('\tR5: Done') - # print('') - # print('') - # print('') - # print('') - - return ( - t0_up, - pos_GS_up, - vel_GS_up, - pos_MEX_up, - vel_MEX_up, - tMEX, - t0_dn, - pos_GS_dn, - vel_GS_dn, - pos_MEX_dn, - vel_MEX_dn, - t0_upp, - pos_GS_upp, - vel_GS_upp, - pos_MEX_upp, - vel_MEX_upp, - pos_MEX, - vel_MEX, - ET_MEX_Mars, - ) + ET: NDArray, ID_SC: int, ID_GS: int, ID_P: int, Ref_frame: str, abcorr: str +) -> Tuple[NDArray, NDArray, NDArray, NDArray, NDArray, NDArray]: + + # ================================ + # Signal emission/reception epochs + # ================================ + + # Epoch at which SC emits downlink signal to GS = Epoch at which SC receives uplink signal from GS + tSC = get_epoch(ET, ID_GS, "<-", ID_SC) + + # Epoch at which the downlink signal is received by Planet from SC + t0_dn = get_epoch(tSC, ID_SC, "->", ID_P) + + # Epoch at which the uplink signal is received by SC from Planet + t0_up = get_epoch(tSC, ID_SC, "<-", ID_P) + + # Didn't really undestand what it was yet. + ET_SC_P = get_epoch(tSC, ID_P, "<-", ID_SC) + + # ==================================================== + # State vectors of GS and SC relative to Planet center + # ==================================================== + + # State vector of GS wrt Planet at t0_dn + state_GS_dn = get_state(t0_dn, ID_P, ID_GS, Ref_frame, abcorr) + + # State vector of SC wrt Planet at t0_dn + state_SC_dn = get_state(t0_dn, ID_P, ID_SC, Ref_frame, abcorr) + + # State vector of GS wrt Planet at t0_up + state_GS_up = get_state(t0_up, ID_P, ID_GS, Ref_frame, abcorr) + + # State vector of SC wrt Planet at t0_up + state_SC_up = get_state(t0_up, ID_P, ID_SC, Ref_frame, abcorr) + + # Used for routine R6, not in R9 ! + state_SC = get_state(ET_SC_P, ID_P, ID_SC, Ref_frame, abcorr) + + return state_GS_up, state_SC_up, state_GS_dn, state_SC_dn, state_SC, ET_SC_P diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index bc4bea0..62312ae 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -100,33 +100,44 @@ def run( # noqa: C901 T_BC_low = constants.T_BC_low T_BC_med = constants.T_BC_med T_BC_upp = constants.T_BC_upp + + POSITION_INDICES = constants.POSITION_INDICES + VELOCITY_INDICES = constants.VELOCITY_INDICES + # ================================================================================= - # R5: State_Vectors_twoway + # R4: State_Vectors_twoway ( - t0_up, - pos_GS_up, - vel_GS_up, - pos_MEX_up, - vel_MEX_up, - tMEX, - t0_dn, - pos_GS_dn, - vel_GS_dn, - pos_MEX_dn, - vel_MEX_dn, - t0_upp, - pos_GS_upp, - vel_GS_upp, - pos_MEX_upp, - vel_MEX_upp, - pos_MEX, - vel_MEX, + state_GS_up, + state_MEX_up, + state_GS_dn, + state_MEX_dn, + state_MEX, ET_MEX_Mars, ) = R4_State_Vectors_twoway_v2.ephemerides( - N_data, ET, NAIF_SC, NAIF_GS, NAIF_P, DATA_PRO, Ref_frame, l2_data.FOLDER_TYPE + ET, + l2_data.SPACECRAFT_NAIF_CODE, + l2_data.dsn_station_naif_code, + l2_data.PLANET_NAIF_CODE, + Ref_frame, + "NONE", ) + pos_GS_up = state_GS_up[:, POSITION_INDICES] + vel_GS_up = state_GS_up[:, VELOCITY_INDICES] + + pos_MEX_up = state_MEX_up[:, POSITION_INDICES] + vel_MEX_up = state_MEX_up[:, VELOCITY_INDICES] + + pos_GS_dn = state_GS_dn[:, POSITION_INDICES] + vel_GS_dn = state_GS_dn[:, VELOCITY_INDICES] + + pos_MEX_dn = state_MEX_dn[:, POSITION_INDICES] + vel_MEX_dn = state_MEX_dn[:, VELOCITY_INDICES] + + pos_MEX = state_MEX[:, POSITION_INDICES] + vel_MEX = state_MEX[:, VELOCITY_INDICES] + # R6: footprint d_lat, d_lon, dist_Mars_Sun, spclon_SC, spclat_SC, re, f = R5_Foot_Print.long_lat( N_data, pos_MEX, vel_MEX, ET, ET_MEX_Mars, Ref_frame, l2_data.FOLDER_TYPE diff --git a/tests/unit_tests_R4.py b/tests/unit_tests_R4.py new file mode 100644 index 0000000..e4b6bf7 --- /dev/null +++ b/tests/unit_tests_R4.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python3 + +""" +radiocc +""" + +from pathlib import Path + +import numpy as np +import spiceypy as spy + +import radiocc +from radiocc.old import R4_State_Vectors_twoway_v2 + +MARS_NAIF_ID = 499 +EARTH_NAIF_ID = 399 +TESTED_CASE = "MEX-M-MRS-1-2-3-PRM-0235-V1.0" +LENGTH_TEST_ARRAY = 10 + + +def test_get_epoch() -> None: + SPK_PATH = Path("tests/to_process") / TESTED_CASE / "EXTRAS/ANCILLARY/SPICE/SPK" + FILES = list(radiocc.utils.directories(SPK_PATH, "*.BSP")) + spy.kclear() + for FILE in FILES: + spy.furnsh(str(FILE)) + ET = np.linspace(0, 3600, LENGTH_TEST_ARRAY) + result = R4_State_Vectors_twoway_v2.get_epoch(ET, EARTH_NAIF_ID, "->", MARS_NAIF_ID) + # Test if output has the right dimension + assert result.shape == (LENGTH_TEST_ARRAY,) + # Test if output doesn't contain NaN's + assert not np.isnan(result).all() + + +def test_get_state() -> None: + SPK_PATH = Path("tests/to_process") / TESTED_CASE / "EXTRAS/ANCILLARY/SPICE/SPK" + FILES = list(radiocc.utils.directories(SPK_PATH, "*.BSP")) + spy.kclear() + for FILE in FILES: + spy.furnsh(str(FILE)) + ET = np.linspace(0, 3600, LENGTH_TEST_ARRAY) + result = R4_State_Vectors_twoway_v2.get_state( + ET, EARTH_NAIF_ID, MARS_NAIF_ID, "J2000", "NONE" + ) + # Test if output has the right dimension + assert result.shape == (LENGTH_TEST_ARRAY, radiocc.constants.DIMENSION_STATE) + # Test if output doesn't contain NaN's + assert not np.isnan(result).all() -- GitLab From f42019d040289aab2ee3bbf53f0a72e87cada077 Mon Sep 17 00:00:00 2001 From: RomanG Date: Sat, 5 Feb 2022 13:31:29 +0100 Subject: [PATCH 05/39] Change R4 unit functions names to be slightly more explicit --- radiocc/old/R4_State_Vectors_twoway_v2.py | 24 ++++++++++++----------- tests/unit_tests_R4.py | 6 ++++-- 2 files changed, 17 insertions(+), 13 deletions(-) diff --git a/radiocc/old/R4_State_Vectors_twoway_v2.py b/radiocc/old/R4_State_Vectors_twoway_v2.py index 91629e4..c336d19 100644 --- a/radiocc/old/R4_State_Vectors_twoway_v2.py +++ b/radiocc/old/R4_State_Vectors_twoway_v2.py @@ -17,14 +17,16 @@ from pudb import set_trace as bp # noqa import radiocc # noqa -def get_epoch(etobs: NDArray, obs: int, direct: str, targ: int) -> NDArray: +def get_epoch_at_target(etobs: NDArray, obs: int, direct: str, targ: int) -> NDArray: epoch_targ = np.full_like(etobs, np.nan) for i, et in enumerate(etobs): epoch_targ[i], _ = spy.ltime(et, obs, direct, targ) return epoch_targ -def get_state(etobs: NDArray, obs: int, targ: int, ref: str, abcorr: str) -> NDArray: +def get_state_at_target( + etobs: NDArray, obs: int, targ: int, ref: str, abcorr: str +) -> NDArray: state_targ = np.array(spy.spkezr(str(targ), etobs, ref, abcorr, str(obs))[0]) * 1e3 return state_targ @@ -38,34 +40,34 @@ def ephemerides( # ================================ # Epoch at which SC emits downlink signal to GS = Epoch at which SC receives uplink signal from GS - tSC = get_epoch(ET, ID_GS, "<-", ID_SC) + tSC = get_epoch_at_target(ET, ID_GS, "<-", ID_SC) # Epoch at which the downlink signal is received by Planet from SC - t0_dn = get_epoch(tSC, ID_SC, "->", ID_P) + t0_dn = get_epoch_at_target(tSC, ID_SC, "->", ID_P) # Epoch at which the uplink signal is received by SC from Planet - t0_up = get_epoch(tSC, ID_SC, "<-", ID_P) + t0_up = get_epoch_at_target(tSC, ID_SC, "<-", ID_P) # Didn't really undestand what it was yet. - ET_SC_P = get_epoch(tSC, ID_P, "<-", ID_SC) + ET_SC_P = get_epoch_at_target(tSC, ID_P, "<-", ID_SC) # ==================================================== # State vectors of GS and SC relative to Planet center # ==================================================== # State vector of GS wrt Planet at t0_dn - state_GS_dn = get_state(t0_dn, ID_P, ID_GS, Ref_frame, abcorr) + state_GS_dn = get_state_at_target(t0_dn, ID_P, ID_GS, Ref_frame, abcorr) # State vector of SC wrt Planet at t0_dn - state_SC_dn = get_state(t0_dn, ID_P, ID_SC, Ref_frame, abcorr) + state_SC_dn = get_state_at_target(t0_dn, ID_P, ID_SC, Ref_frame, abcorr) # State vector of GS wrt Planet at t0_up - state_GS_up = get_state(t0_up, ID_P, ID_GS, Ref_frame, abcorr) + state_GS_up = get_state_at_target(t0_up, ID_P, ID_GS, Ref_frame, abcorr) # State vector of SC wrt Planet at t0_up - state_SC_up = get_state(t0_up, ID_P, ID_SC, Ref_frame, abcorr) + state_SC_up = get_state_at_target(t0_up, ID_P, ID_SC, Ref_frame, abcorr) # Used for routine R6, not in R9 ! - state_SC = get_state(ET_SC_P, ID_P, ID_SC, Ref_frame, abcorr) + state_SC = get_state_at_target(ET_SC_P, ID_P, ID_SC, Ref_frame, abcorr) return state_GS_up, state_SC_up, state_GS_dn, state_SC_dn, state_SC, ET_SC_P diff --git a/tests/unit_tests_R4.py b/tests/unit_tests_R4.py index e4b6bf7..e59fd4a 100644 --- a/tests/unit_tests_R4.py +++ b/tests/unit_tests_R4.py @@ -25,7 +25,9 @@ def test_get_epoch() -> None: for FILE in FILES: spy.furnsh(str(FILE)) ET = np.linspace(0, 3600, LENGTH_TEST_ARRAY) - result = R4_State_Vectors_twoway_v2.get_epoch(ET, EARTH_NAIF_ID, "->", MARS_NAIF_ID) + result = R4_State_Vectors_twoway_v2.get_epoch_at_target( + ET, EARTH_NAIF_ID, "->", MARS_NAIF_ID + ) # Test if output has the right dimension assert result.shape == (LENGTH_TEST_ARRAY,) # Test if output doesn't contain NaN's @@ -39,7 +41,7 @@ def test_get_state() -> None: for FILE in FILES: spy.furnsh(str(FILE)) ET = np.linspace(0, 3600, LENGTH_TEST_ARRAY) - result = R4_State_Vectors_twoway_v2.get_state( + result = R4_State_Vectors_twoway_v2.get_state_at_target( ET, EARTH_NAIF_ID, MARS_NAIF_ID, "J2000", "NONE" ) # Test if output has the right dimension -- GitLab From c28bf70e1837037191d5800abd0a2948f337d984 Mon Sep 17 00:00:00 2001 From: RomanG Date: Tue, 15 Feb 2022 12:39:47 +0100 Subject: [PATCH 06/39] Edit global consistency test to compare within tol --- tests/test_radiocc.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index faf6989..9431f1e 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -79,7 +79,7 @@ def test_number_of_lines_and_columns() -> None: assert np.loadtxt(TEXT_OUTPUT).shape[1] == NUMBER_OF_COLUMNS_OUTPUT[LAYER] -def test_output_consistency() -> None: +def test_output_global_consistency() -> None: radiocc.cfg.to_process = Path("tests/to_process") radiocc.cfg.results = Path("tests/results") radiocc.cfg.full_auto = True @@ -100,6 +100,10 @@ def test_output_consistency() -> None: TEXT_OUTPUTS.values(), EXPECTED_TEXT_OUTPUTS.values() ): # Test if all values of the produced and expected output are equal - assert np.array_equal( - np.loadtxt(TEXT_OUTPUT), np.loadtxt(EXPECTED_TEXT_OUTPUT), equal_nan=True + assert np.allclose( + np.loadtxt(TEXT_OUTPUT), + np.loadtxt(EXPECTED_TEXT_OUTPUT), + equal_nan=True, + atol=0, + rtol=0, ) -- GitLab From 4614c88038165bb1d1099e29bfc7741c7f9a9e60 Mon Sep 17 00:00:00 2001 From: RomanG Date: Tue, 15 Feb 2022 12:42:23 +0100 Subject: [PATCH 07/39] Refactor R6. Make unit tests for R6 + Refactor R6 by breaking the routine into small subroutines + Make unit tests for each subroutine + Implemented 2 functions in utils.py to compute the length of number and to round a float or an array to a number of significant figures. --- radiocc/old/R6_Framework_Conversion.py | 350 ++++++++----------------- radiocc/reconstruct.py | 87 ++++-- radiocc/utils.py | 14 + tests/unit_tests_R6.py | 127 +++++++++ 4 files changed, 311 insertions(+), 267 deletions(-) create mode 100644 tests/unit_tests_R6.py diff --git a/radiocc/old/R6_Framework_Conversion.py b/radiocc/old/R6_Framework_Conversion.py index ede4932..ebf2c4e 100644 --- a/radiocc/old/R6_Framework_Conversion.py +++ b/radiocc/old/R6_Framework_Conversion.py @@ -1,247 +1,113 @@ +from typing import Tuple + import numpy as np -import scipy as sc -import spiceypy +import spiceypy as spy +from nptyping import NDArray from pudb import set_trace as bp - -def Occ_Plane_Conversion(N_data, pos_MEX, vel_MEX, pos_GS, vel_GS): - - gamma = np.zeros((N_data,)) - - for i in range(N_data): - - vsep = spiceypy.spiceypy.vsep(pos_MEX[i], pos_GS[i]) - - gamma[i] = vsep - (sc.pi / 2.0) - - r_MEX = np.zeros((N_data,)) - - z_MEX = np.zeros((N_data,)) - - z_GS = np.zeros((N_data,)) - - for i in range(N_data): - - r_MEX[i] = np.linalg.norm(pos_MEX[i]) * np.cos(gamma[i]) - - z_MEX[i] = np.linalg.norm(pos_MEX[i]) * np.sin(gamma[i]) - - z_GS[i] = -np.linalg.norm(pos_GS[i]) - - minus_pos_GS = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - for i in range(N_data): - - minus_pos_GS[i] = -pos_GS[i] - - vr_MEX = np.zeros((N_data,)) - - vz_MEX = np.zeros((N_data,)) - - vr_GS = np.zeros((N_data,)) - - vz_GS = np.zeros((N_data,)) - - for i in range(N_data): - - n_vec = np.cross(pos_GS[i], pos_MEX[i]) - - r_vec = np.cross(minus_pos_GS[i], n_vec) - - M = 3 - - B = np.array( - [ - [n_vec[0], r_vec[0], minus_pos_GS[i][0], 1, 0, 0], - [n_vec[1], r_vec[1], minus_pos_GS[i][1], 0, 1, 0], - [n_vec[2], r_vec[2], minus_pos_GS[i][2], 0, 0, 1], - ], - dtype=float, - ) - - G = -1 - - while G < M - 1: - - G += 1 - - B[G, :] = B[G, :] / B[G, G] - - for J in range(M): - - if J == G: - - K = G - - else: - - if K == M - 1: - - K = 0 - - else: - - K += 1 - - B[K, :] = B[K, :] - B[K, G] * B[G, :] - - P = B[:, 3:6] - - beta_MEX = ( - vel_MEX[i][0] * P[1, 0] + vel_MEX[i][1] * P[1, 1] + vel_MEX[i][2] * P[1, 2] - ) - - gamma_MEX = ( - vel_MEX[i][0] * P[2, 0] + vel_MEX[i][1] * P[2, 1] + vel_MEX[i][2] * P[2, 2] - ) - - vr_MEX[i] = -np.sqrt( - (beta_MEX * r_vec[0]) ** 2 - + (beta_MEX * r_vec[1]) ** 2 - + (beta_MEX * r_vec[2]) ** 2 - ) - - vz_MEX[i] = np.sqrt( - (gamma_MEX * minus_pos_GS[i][0]) ** 2 - + (gamma_MEX * minus_pos_GS[i][1]) ** 2 - + (gamma_MEX * minus_pos_GS[i][2]) ** 2 - ) - - beta_GS = ( - vel_GS[i][0] * P[1, 0] + vel_GS[i][1] * P[1, 1] + vel_GS[i][2] * P[1, 2] - ) - - gamma_GS = ( - vel_GS[i][0] * P[2, 0] + vel_GS[i][1] * P[2, 1] + vel_GS[i][2] * P[2, 2] - ) - - vr_GS[i] = np.sqrt( - (beta_GS * r_vec[0]) ** 2 - + (beta_GS * r_vec[1]) ** 2 - + (beta_GS * r_vec[2]) ** 2 - ) - - vz_GS[i] = np.sqrt( - (gamma_GS * minus_pos_GS[i][0]) ** 2 - + (gamma_GS * minus_pos_GS[i][1]) ** 2 - + (gamma_GS * minus_pos_GS[i][2]) ** 2 - ) - - delta_s = np.zeros((N_data,)) - - beta_e = np.zeros((N_data,)) - - for i in range(N_data): - - a = np.array([r_MEX[i], z_MEX[i]], dtype=float) - - b = np.array([0, z_GS[i]], dtype=float) - - diff_vec = a - b - +from radiocc.constants import POSITION_INDICES, VELOCITY_INDICES +from radiocc.utils import round_to_N_significant_figures + + +def compute_SC_angle(pos_SC: NDArray, pos_GS: NDArray) -> NDArray: + """Compute the angle between the r-axis of the occultation plane + and the direction vector of Mars to the spacecraft. This angle is + defined as gamma in Fig. 20 of Fjeldbo 1971.""" + gamma = np.full(pos_SC.shape[0], np.nan) + for i, (p_SC, p_GS) in enumerate(zip(pos_SC, pos_GS)): + gamma[i] = spy.vsep(p_SC, p_GS) - 0.5 * np.pi + return gamma + + +def position_in_occ_plane( + pos_SC: NDArray, pos_GS: NDArray, gamma: NDArray +) -> Tuple[NDArray, NDArray, NDArray]: + """Transform the direction vectors of Mars to the spacecraft and + Mars to the ground station into the occultation coordinate system (r,z).""" + r_SC = np.linalg.norm(pos_SC, axis=1) * np.cos(gamma) + z_SC = np.linalg.norm(pos_SC, axis=1) * np.sin(gamma) + z_GS = -np.linalg.norm(pos_GS, axis=1) + return r_SC, z_SC, z_GS + + +def transition_matrix_occ_plane(p_SC: NDArray, p_GS: NDArray) -> NDArray: + """Compute the transition matrix from the equatorial geocentric basis (x,y,z) + to the occultation plane basis (n, r, z).""" + n_vec = np.cross(p_GS, p_SC) # vector normal to occ. plane + n_vec /= np.linalg.norm(n_vec) + r_vec = -np.cross(p_GS, n_vec) # r-axis of occ. plane (r = z x n) + r_vec /= np.linalg.norm(r_vec) + z_vec = -p_GS / np.linalg.norm(p_GS) + + base_matrix = np.stack( + (n_vec, r_vec, z_vec), axis=1 + ) # Transition matrix from equatorial frame to occ. plane + return np.linalg.inv(base_matrix) + + +def velocity_in_occ_plane( + pos_SC: NDArray, pos_GS: NDArray, vel_SC: NDArray, vel_GS: NDArray +) -> Tuple[NDArray, NDArray, NDArray, NDArray]: + """Using a transition matrix, transform the velocity vectors of the spacecraft + and the ground station relative to Mars center into the occultation coordinate system.""" + vr_SC = np.full(vel_SC.shape[0], np.nan) + vz_SC = np.full_like(vr_SC, np.nan) + vr_GS = np.full_like(vr_SC, np.nan) + vz_GS = np.full_like(vr_SC, np.nan) + + for i, (p_SC, p_GS, v_SC, v_GS) in enumerate(zip(pos_SC, pos_GS, vel_SC, vel_GS)): + trans_mat = transition_matrix_occ_plane(p_SC, p_GS) + _, vr_SC[i], vz_SC[i] = trans_mat.dot(v_SC) + _, vr_GS[i], vz_GS[i] = trans_mat.dot(v_GS) + vr_SC[i] = -np.abs(vr_SC[i]) + vz_SC[i] = np.abs(vz_SC[i]) + vr_GS[i] = np.abs(vr_GS[i]) + vz_GS[i] = np.abs(vz_GS[i]) + return vr_SC, vz_SC, vr_GS, vz_GS + + +def compute_delta_s(r_SC: NDArray, z_SC: NDArray, z_GS: NDArray) -> NDArray: + """Compute the angle between direction vector of the ground station to Mars + (z-axis of the occultation plane) and the direction vector of the ground station + to the spacecraft. This angle is defined as delta_s in Fig. 20 of Fjeldbo 1971.""" + delta_s = np.full(r_SC.shape, np.nan) + + for i, (r_SC_i, z_SC_i, z_GS_i) in enumerate(zip(r_SC, z_SC, z_GS)): + # delta_s[i] = spy.vsepg( + # np.array([0, -z_GS_i]), np.array([r_SC_i, z_SC_i - z_GS_i]), 2 + # ) # Defined in Fjeldbo 1971 delta_s[i] = np.arccos( - np.dot(-b, diff_vec) / np.linalg.norm(-z_GS[i]) / np.linalg.norm(diff_vec) - ) - - beta_e[i] = (np.pi / 2.0) - np.arccos( - np.dot(-b, diff_vec) / np.linalg.norm(-z_GS[i]) / np.linalg.norm(diff_vec) - ) - - # print('\tR7: Done') - # print('gamma[3000]:',gamma[3000], 'beta_e[3000]:', beta_e[3000], 'delta_s[3000]:', delta_s[3000]) - return r_MEX, z_MEX, z_GS, vr_MEX, vz_MEX, vr_GS, vz_GS, gamma, beta_e, delta_s - - -def Occ_Plane_Conversion_2(N_data, pos_MEX, vel_MEX, pos_GS, vel_GS): - - z_vec = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - n_vec = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - r_vec = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - pos_MEX_new = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - vel_MEX_new = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - pos_GS_new = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - vel_GS_new = np.full( - (N_data, 3), np.nan - ) # np.array([[np.nan for j in range(3)] for i in range(N_data)], dtype = float) - - r_MEX = np.zeros((N_data,)) - - z_MEX = np.zeros((N_data,)) - - z_GS = np.zeros((N_data,)) - - vr_MEX = np.zeros((N_data,)) - - vz_MEX = np.zeros((N_data,)) - - vr_GS = np.zeros((N_data,)) - - vz_GS = np.zeros((N_data,)) - - gamma = np.zeros((N_data,)) - - for i in range(N_data): - - vsep = spiceypy.spiceypy.vsep(pos_MEX[i], pos_GS[i]) - - gamma[i] = vsep - (sc.pi / 2.0) - - z_vec[i] = -pos_GS[i] / np.linalg.norm(pos_GS[i]) - - n_vec[i] = np.cross(pos_MEX[i] / np.linalg.norm(pos_MEX[i]), z_vec[i]) / np.sin( - (sc.pi / 2.0) - gamma[i] + np.dot(np.array([0, -z_GS_i]), np.array([r_SC_i, z_SC_i - z_GS_i])) + / np.abs(z_GS_i) + / np.linalg.norm(np.array([r_SC_i, z_SC_i - z_GS_i])) ) - - r_vec[i] = np.cross(z_vec[i], n_vec[i]) - - H_old_to_new = np.array([z_vec[i][:], n_vec[i][:], r_vec[i][:]], dtype=float) - - pos_MEX_new[i] = np.dot(H_old_to_new, pos_MEX[i]) - - pos_GS_new[i] = np.dot(H_old_to_new, pos_GS[i]) - - vel_MEX_new[i] = np.dot(H_old_to_new, vel_MEX[i]) - - vel_GS_new[i] = np.dot(H_old_to_new, vel_GS[i]) - - r_MEX[i] = pos_MEX_new[i][2] - - z_MEX[i] = pos_MEX_new[i][0] - - z_GS[i] = pos_GS_new[i][0] - - vr_MEX[i] = vel_MEX_new[i][2] - - vz_MEX[i] = vel_MEX_new[i][0] - - vr_GS[i] = vel_GS_new[i][2] - - vz_GS[i] = vel_GS_new[i][0] - - # print('') - # print('') - # print('\tR7: Done') - # print('') - # print('') - - return r_MEX, z_MEX, z_GS, vr_MEX, vz_MEX, vr_GS, vz_GS + return delta_s + + +def Occ_Plane_Conversion( + state_SC: NDArray, state_GS: NDArray +) -> Tuple[NDArray, NDArray, NDArray, NDArray]: + """ + + Transform the state vectors of the spacecraft and the ground station + relative to Mars center expressed in the equatorial geocentric frame into + the occultation plane coordinate system (r,z). + + Compute the angle between the r-axis of the occultation plane and the + Mars-SC direction vector -> gamma in Fig. 20 of Fjeldbo 1971. + + Compute the angle between the z-axis of the occultation plane and the + GS-SC direction vector -> delta_s in Fig. 20 of Fjeldbo 1971. + """ + N_data = state_SC.shape[0] + pos_SC = state_SC[:, POSITION_INDICES] + vel_SC = state_SC[:, VELOCITY_INDICES] + pos_GS = state_GS[:, POSITION_INDICES] + vel_GS = state_GS[:, VELOCITY_INDICES] + + gamma = compute_SC_angle(pos_SC, pos_GS) + r_SC, z_SC, z_GS = position_in_occ_plane(pos_SC, pos_GS, gamma) + vr_SC, vz_SC, vr_GS, vz_GS = velocity_in_occ_plane(pos_SC, pos_GS, vel_SC, vel_GS) + + state_occ_SC = np.stack((r_SC, z_SC, vr_SC, vz_SC), axis=1) + state_occ_GS = np.stack((np.zeros(N_data), z_GS, vr_GS, vz_GS), axis=1) + + delta_s = compute_delta_s(r_SC, z_SC, z_GS) + return state_occ_SC, state_occ_GS, gamma, delta_s diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index 62312ae..8bf9efe 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -138,42 +138,79 @@ def run( # noqa: C901 pos_MEX = state_MEX[:, POSITION_INDICES] vel_MEX = state_MEX[:, VELOCITY_INDICES] - # R6: footprint + # R5: footprint d_lat, d_lon, dist_Mars_Sun, spclon_SC, spclat_SC, re, f = R5_Foot_Print.long_lat( N_data, pos_MEX, vel_MEX, ET, ET_MEX_Mars, Ref_frame, l2_data.FOLDER_TYPE ) - # R7 coordinate frames + # # R6 coordinate frames + # ( + # 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, + # ) = R6_Framework_Conversion.Occ_Plane_Conversion( + # N_data, pos_MEX_up, vel_MEX_up, pos_GS_up, vel_GS_up + # ) + + # # R6 coordinate frames + # ( + # 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, + # ) = R6_Framework_Conversion.Occ_Plane_Conversion( + # N_data, pos_MEX_dn, vel_MEX_dn, pos_GS_dn, vel_GS_dn + # ) + + # R6 framework conversion (uplink) ( - 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, gamma_up, - beta_e_up, delta_s_up, - ) = R6_Framework_Conversion.Occ_Plane_Conversion( - N_data, pos_MEX_up, vel_MEX_up, pos_GS_up, vel_GS_up - ) + ) = R6_Framework_Conversion.Occ_Plane_Conversion(state_MEX_up, state_GS_up) + + 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] - # R7 coordinate frames + z_GS_up = state_occ_GS_up[:, 1] + vr_GS_up = state_occ_GS_up[:, 2] + vz_GS_up = state_occ_GS_up[:, 3] + + # R6 framework conversion (downlink) ( - r_MEX_dn, - z_MEX_dn, - z_GS_dn, - vr_MEX_dn, - vz_MEX_dn, - vr_GS_dn, - vz_GS_dn, + state_occ_SC_dn, + state_occ_GS_dn, gamma_dn, - beta_e_dn, delta_s_dn, - ) = R6_Framework_Conversion.Occ_Plane_Conversion( - N_data, pos_MEX_dn, vel_MEX_dn, pos_GS_dn, vel_GS_dn - ) + ) = R6_Framework_Conversion.Occ_Plane_Conversion(state_MEX_dn, state_GS_dn) + + 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 * numpy.pi - delta_s_up + beta_e_dn = 0.5 * numpy.pi - delta_s_dn if l2_data.FOLDER_TYPE == ProcessType.MAVEN: # R18 : Occ_info : diff --git a/radiocc/utils.py b/radiocc/utils.py index da898fa..b721a7b 100644 --- a/radiocc/utils.py +++ b/radiocc/utils.py @@ -9,10 +9,12 @@ import itertools from pathlib import Path from typing import Any, Dict, Iterable, Iterator, List, Optional, Union +import numpy as np import ruamel.yaml import yaml from colored import attr, fg from dotmap import DotMap +from nptyping import NDArray from pudb import set_trace as bp # noqa DOTMAP_NONE = (None, DotMap()) @@ -158,3 +160,15 @@ def read_yaml(PATH: Path) -> DotMap: """Load the yaml file.""" CFGF_DICT, IND, BSI = ruamel.yaml.util.load_yaml_guess_indent(PATH.open()) return DotMap(CFGF_DICT) + + +def number_length(value: Union[float, NDArray], base: int) -> Union[int, NDArray]: + return np.floor(np.log(np.abs(value)) / np.log(base)).astype(int) + 1 + + +def round_to_N_significant_figures(value: Union[float, NDArray], figures: int) -> float: + n = figures - number_length(value, 10) + if isinstance(value, float): + return round(value * 10 ** n) / 10 ** n + else: + return (value * 10 ** n).round() / 10 ** n diff --git a/tests/unit_tests_R6.py b/tests/unit_tests_R6.py new file mode 100644 index 0000000..835446c --- /dev/null +++ b/tests/unit_tests_R6.py @@ -0,0 +1,127 @@ +#!/usr/bin/env python3 + +""" +radiocc +""" + +from pathlib import Path +from zlib import Z_SYNC_FLUSH + +import numpy as np +import spiceypy as spy + +import radiocc +from radiocc.old import R6_Framework_Conversion +from radiocc.utils import number_length, round_to_N_significant_figures + +LENGTH_TEST_ARRAY = 2 +DIMENSION_VECTOR = radiocc.constants.DIMENSION_VECTOR +PATH_EXPECTED_OUTPUT = Path("tests/expected_unit_tests_R6") + +# Test input (first 2 lines from MEX 235) +pos_SC_test = np.array( + [ + [-7726575.13322412, -7666711.32894175, -7958119.5346874], + [-7726451.93234389, -7667213.36854496, -7957816.58827798], + ] +) +pos_GS_test = np.array( + [ + [192537483124.12, 265959242098.085, 113482080420.543], + [192537459864.319, 265959250204.846, 113482084397.897], + ] +) + +gamma_test = np.array([1.2397252315812128, 1.2397669744753683]) + +r_SC_test = np.array([4382958.64228272, 4382438.07915611]) +z_SC_test = np.array([12751461.32641642, 12751678.38508523]) +z_GS_test = np.array([-347394852353.5334, -347394846967.8551]) + + +def test_compute_SC_angle() -> None: + + output = R6_Framework_Conversion.compute_SC_angle(pos_SC_test, pos_GS_test) + # Test if output has the right dimension + assert output.shape == (LENGTH_TEST_ARRAY,) + # Test if output doesn't contain NaN's + assert not np.isnan(output).all() + # Test if output of test is element-wise equal to the expected output within + # an absolute tolerance corresponding to the number of significant figures + assert np.allclose( + output, + gamma_test, + atol=0, + rtol=1e-6, + equal_nan=False, + ) + + +def test_position_in_occ_plane() -> None: + gamma_rounded = round_to_N_significant_figures(gamma_test, 15) + output = R6_Framework_Conversion.position_in_occ_plane( + pos_SC_test, pos_GS_test, gamma_rounded + ) + # Test if output has the right dimension + assert np.array(output).shape == (3, LENGTH_TEST_ARRAY) + # Test if output doesn't contain NaN's + assert not np.isnan(output).all() + # Test if output of test is element-wise equal to the expected output within + # an absolute tolerance corresponding to the number of significant figures + + expected_output = (r_SC_test, z_SC_test, z_GS_test) + assert np.allclose( + np.array(output), + np.array(expected_output), + atol=0, + rtol=1e-6, + equal_nan=False, + ) + + +def test_transition_matrix_occ_plane() -> None: + output = R6_Framework_Conversion.transition_matrix_occ_plane( + pos_SC_test, pos_GS_test + ) + assert not np.isnan(output).all() + + +def test_velocity() -> None: + vel_SC_test = np.array( + [ + [240.58413463, -980.54769211, 591.63687863], + [240.65325798, -980.4791224, 591.70807734], + ] + ) + vel_GS_test = np.array( + [ + [-45427.78743875, 15832.98573322, 7768.01058944], + [-45427.7801956, 15832.99253347, 7768.00905843], + ] + ) + + vr_SC = np.array([-1016.962309512641, -1016.9834579503646]) + output = R6_Framework_Conversion.velocity_in_occ_plane( + pos_SC_test, pos_GS_test, vel_SC_test, vel_GS_test + ) + assert np.allclose( + output[0], + vr_SC, + atol=0, + rtol=1e-10, + equal_nan=False, + ) + + +def test_compute_delta_s() -> None: + + delta_s = R6_Framework_Conversion.compute_delta_s(r_SC_test, z_SC_test, z_GS_test) + expected_delta_s = np.array([1.2616189367258302e-05, 1.2614684479727514e-05]) + + assert np.allclose( + delta_s, + expected_delta_s, + atol=0, + rtol=1e-5, + equal_nan=False, + ) -- GitLab From 3dbf2e08b22c1f616f4e167f71521d411881d7e6 Mon Sep 17 00:00:00 2001 From: RomanG Date: Tue, 15 Feb 2022 12:48:47 +0100 Subject: [PATCH 08/39] Removed call to unused R5 routine. --- radiocc/reconstruct.py | 37 ------------------------------------- 1 file changed, 37 deletions(-) diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index 8bf9efe..1b3e7b2 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -138,43 +138,6 @@ def run( # noqa: C901 pos_MEX = state_MEX[:, POSITION_INDICES] vel_MEX = state_MEX[:, VELOCITY_INDICES] - # R5: footprint - d_lat, d_lon, dist_Mars_Sun, spclon_SC, spclat_SC, re, f = R5_Foot_Print.long_lat( - N_data, pos_MEX, vel_MEX, ET, ET_MEX_Mars, Ref_frame, l2_data.FOLDER_TYPE - ) - - # # R6 coordinate frames - # ( - # 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, - # ) = R6_Framework_Conversion.Occ_Plane_Conversion( - # N_data, pos_MEX_up, vel_MEX_up, pos_GS_up, vel_GS_up - # ) - - # # R6 coordinate frames - # ( - # 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, - # ) = R6_Framework_Conversion.Occ_Plane_Conversion( - # N_data, pos_MEX_dn, vel_MEX_dn, pos_GS_dn, vel_GS_dn - # ) - # R6 framework conversion (uplink) ( state_occ_SC_up, -- GitLab From c744a38702e024d9e92e9c79f4c69090fbdd2bd2 Mon Sep 17 00:00:00 2001 From: RomanG Date: Tue, 15 Feb 2022 12:52:20 +0100 Subject: [PATCH 09/39] Remove now useless temporary variables. --- radiocc/reconstruct.py | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index 1b3e7b2..6b1d904 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -123,21 +123,6 @@ def run( # noqa: C901 "NONE", ) - pos_GS_up = state_GS_up[:, POSITION_INDICES] - vel_GS_up = state_GS_up[:, VELOCITY_INDICES] - - pos_MEX_up = state_MEX_up[:, POSITION_INDICES] - vel_MEX_up = state_MEX_up[:, VELOCITY_INDICES] - - pos_GS_dn = state_GS_dn[:, POSITION_INDICES] - vel_GS_dn = state_GS_dn[:, VELOCITY_INDICES] - - pos_MEX_dn = state_MEX_dn[:, POSITION_INDICES] - vel_MEX_dn = state_MEX_dn[:, VELOCITY_INDICES] - - pos_MEX = state_MEX[:, POSITION_INDICES] - vel_MEX = state_MEX[:, VELOCITY_INDICES] - # R6 framework conversion (uplink) ( state_occ_SC_up, -- GitLab From 8662a4092dd5681d68f3b73b15205f5710c6314f Mon Sep 17 00:00:00 2001 From: RomanG Date: Tue, 15 Feb 2022 15:27:21 +0100 Subject: [PATCH 10/39] Divide R7 into unit functions. Make tests for each function. + Refactor R7 by improving code formatting and legibility + Adapt reconstruct to the new R7 --- radiocc/old/R7_Offset_Correction.py | 137 ++++++++++++--------------- radiocc/reconstruct.py | 138 ++++++++-------------------- tests/unit_tests_R7.py | 42 +++++++++ 3 files changed, 142 insertions(+), 175 deletions(-) create mode 100644 tests/unit_tests_R7.py diff --git a/radiocc/old/R7_Offset_Correction.py b/radiocc/old/R7_Offset_Correction.py index 55e9813..56509b5 100644 --- a/radiocc/old/R7_Offset_Correction.py +++ b/radiocc/old/R7_Offset_Correction.py @@ -1,19 +1,53 @@ -import os -from typing import Optional +from typing import Optional, Tuple import numpy as np +from nptyping import NDArray from pudb import set_trace as bp # noqa -from scipy import signal as sg import radiocc -#Type = "Quadratic" or "Linear", threshold [m] -def Off_Cor(Type, threshold, Doppler, Diff_doppler, ET, N_data, distance, CODE_DIR,i_integral, DATA_TYPE: Optional[radiocc.model.RadioDataType]): - - Doppler_debias = np.full(N_data, np.nan) - Doppler_biasfit = np.full(N_data, np.nan) - delET = np.full(N_data, np.nan) +def find_corr_threshold_index(distance: NDArray, threshold: float) -> int: + """Find index of the altitude corresponding to the bias correction threshold.""" + i_corr = 0 + while distance[i_corr] > threshold: + i_corr += 1 + if i_corr == len(distance): + i_corr -= 1 + return i_corr + + +def fit_Doppler_bias( + delET: NDArray, Doppler: NDArray, i_integral: int, i_corr: int, Type: str +) -> Tuple[NDArray, NDArray]: + """Fit the raw Doppler residuals above the bias correction threshold + by a polynomial of degree 1 or 2. + Output the coefficients of the fitting polynomial and the values of + the polynomial.""" + if Type == "Linear": + deg = 1 + elif Type == "Quadratic": + deg = 2 + else: + print("Type is not recognized") + + poly_coefs = np.polyfit(delET[i_integral:i_corr], Doppler[i_integral:i_corr], deg) + Doppler_biasfit = np.polyval(poly_coefs, delET) + return poly_coefs, Doppler_biasfit + + +def Off_Cor( + Type: str, + threshold: int, + Doppler: NDArray, + ET: NDArray, + distance: NDArray, + i_integral: int, + DATA_TYPE: Optional[radiocc.model.RadioDataType], +) -> Tuple[NDArray, NDArray, NDArray, NDArray, int, NDArray]: + + """Corrects the offset bias in the raw Doppler residuals by fitting the residuals + above the upper atmosphere by a Linear or Quadratic function.""" # If egress, we reverse the matrices to have them from atmosphere to surface. if DATA_TYPE is not None and DATA_TYPE == radiocc.model.RadioDataType.EGRESS: @@ -21,73 +55,24 @@ def Off_Cor(Type, threshold, Doppler, Diff_doppler, ET, N_data, distance, CODE_D Doppler = Doppler[::-1] ET = ET[::-1] - # altitude threshold for fit - Corr_ind = 0 - for i in range(len(distance)): - if distance[i] > threshold: - Corr_ind += 1 - if Corr_ind == len(distance): - Corr_ind = Corr_ind-1 - # Search for maximum - - if len(Doppler)%2==1 : - taille = len(Doppler) - else: - taille = len(Doppler)-1 - - Smoothed_Doppler = sg.savgol_filter(Doppler,taille,40) - deriv = np.gradient(Smoothed_Doppler) - - good=np.logical_and(deriv[np.array(np.where(deriv[:-1]<1E-5))-1]<0,\ - deriv[np.array(np.where(deriv[:-1]<1E-5))+1]>0) - - - targ=np.where(deriv[:-1]<1E-5) - index = targ[0][good[0]] - good_one = index[-1] - - # FIT - if Type == 'Fit': - p = np.polyfit(np.append(ET[:Corr_ind],ET[good_one])\ - ,np.append(Doppler[:Corr_ind],Doppler[good_one]),2) - - for i in range(i_integral,N_data): - Doppler_biasfit[i] = p[0] * ET[i]**2 + p[1]*ET[i] + p[2] - -########################################################################################### - for i in range(N_data): - delET[i]= ET[i]-ET[0] - - # bias fit - if Type == 'Linear': - p = np.polyfit(delET[i_integral:Corr_ind],Doppler[i_integral:Corr_ind],1) - p0 = p[0] - p1 = p[1] - for i in range(N_data): - Doppler_biasfit[i] = p[0] * (delET[i]) + p[1] - - elif Type == 'Quadratic': - p = np.polyfit(delET[i_integral:Corr_ind],Doppler[i_integral:Corr_ind],2) - p0 = p[0] - p1 = p[1] - p2 = p[2] - for i in range(N_data): - Doppler_biasfit[i] = p[0] * (delET[i])**2 + p[1]*(delET[i]) + p[2] - - else: - print('Type is not recognized') - - # subtract bias and scale residual - for i in range(N_data): - Doppler_debias[i] = (Doppler[i]- Doppler_biasfit[i]) - #print(ET[13542],Doppler[13542],Doppler_biasfit[13542], Doppler_debias[13542]) - # back to code directory (why in this subroutine?) - # os.chdir(CODE_DIR) - if Type == 'Linear': - return ET, delET, Doppler_debias, Doppler_biasfit, Corr_ind, p0, p1 - - if Type == 'Quadraticuadratic': - return ET, delET, Doppler_debias, Doppler_biasfit, Corr_ind, p0, p1,p2 + delET = ET - ET[0] + + # Find index of the altitude threshold for bias correction + i_corr = find_corr_threshold_index(distance, threshold) + # Fit raw Doppler residuals above the threshold by a polynomial + poly_coefs, Doppler_biasfit = fit_Doppler_bias( + delET, Doppler, i_integral, i_corr, Type + ) + # Subtract the fitted bias from the raw residuals + Doppler_debias = Doppler - Doppler_biasfit + return ( + ET, + delET, + Doppler_debias, + Doppler_biasfit, + i_corr, + poly_coefs, + ) diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index 6b1d904..40363a3 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -187,60 +187,27 @@ def run( # noqa: C901 else: DATA_TYPE = None + ( + ET, + delET, + Doppler_debias, + Doppler_biasfit, + Corr_ind, + poly_coefs, + ) = R7_Offset_Correction.Off_Cor( + B_Cor_Type, + Threshold_Cor, + Doppler, + ET, + distance, + i_integral, + DATA_TYPE, + ) + if B_Cor_Type == "Linear": - ( - ET, - delET, - Doppler_debias, - Doppler_biasfit, - Corr_ind, - p0, - p1, - ) = R7_Offset_Correction.Off_Cor( - B_Cor_Type, - Threshold_Cor, - Doppler, - Diff_doppler, - ET, - N_data, - distance, - CODE_DIR, - i_integral, - DATA_TYPE, - # Threshold_Surface, - # Threshold_int, - # l2_data.FOLDER_TYPE, - ) + p0, p1 = poly_coefs if B_Cor_Type == "Quadratic": - ( - i_Surface, - ET, - delET, - Doppler_debias, - Doppler_biasfit, - Corr_ind, - p0, - p1, - p2, - ) = R7_Offset_Correction.Off_Cor( - i_integral, - i_Surface, - ET, - delET, - B_Cor_Type, - Threshold_Cor, - Doppler, - Diff_doppler, - ET, - N_data, - distance, - CODE_DIR, - i_integral, - DATA_TYPE, - # Threshold_Surface, - # Threshold_int, - # l2_data.FOLDER_TYPE, - ) + p0, p1, p2 = poly_coefs # R5: Plot doppler after bias correction R8_Plot_Doppler_1.PLOT_Dop1( @@ -277,55 +244,28 @@ def run( # noqa: C901 else: NEW_THRESHOLD = float(input("Enter new threshold altitude = ")) NEW_TYPE = str(input("Enter new type, Linear or Quadratic: ")) + + ( + ET, + delET, + Doppler_debias, + Doppler_biasfit, + Corr_ind, + poly_coefs, + ) = R7_Offset_Correction.Off_Cor( + NEW_TYPE, + NEW_THRESHOLD, + Doppler, + ET, + distance, + i_integral, + DATA_TYPE, + ) + if NEW_TYPE == "Linear": - ( - ET, - delET, - Doppler_debias, - Doppler_biasfit, - Corr_ind, - p0, - p1, - ) = R7_Offset_Correction.Off_Cor( - NEW_TYPE, - NEW_THRESHOLD, - Doppler, - Diff_doppler, - ET, - N_data, - distance, - CODE_DIR, - i_integral, - DATA_TYPE, - # Threshold_Surface, - # Threshold_int, - # l2_data.FOLDER_TYPE, - ) + p0, p1 = poly_coefs if NEW_TYPE == "Quadratic": - ( - ET, - delET, - Doppler_debias, - Doppler_biasfit, - Corr_ind, - p0, - p1, - p2, - ) = R7_Offset_Correction.Off_Cor( - NEW_TYPE, - NEW_THRESHOLD, - Doppler, - Diff_doppler, - ET, - N_data, - distance, - CODE_DIR, - i_integral, - DATA_TYPE, - # Threshold_Surface, - # Threshold_int, - # l2_data.FOLDER_TYPE, - ) + p0, p1, p2 = poly_coefs R8_Plot_Doppler_1.PLOT_Dop1( distance, Doppler, Doppler_debias, Doppler_biasfit, ET, PLOT_DIR, N_data diff --git a/tests/unit_tests_R7.py b/tests/unit_tests_R7.py new file mode 100644 index 0000000..6d938be --- /dev/null +++ b/tests/unit_tests_R7.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 + +""" +radiocc +""" + +from pathlib import Path +from turtle import distance + +import numpy as np + +import radiocc +from radiocc.old import R7_Offset_Correction +from radiocc.utils import number_length, round_to_N_significant_figures + +PATH_EXPECTED_OUTPUT = Path("tests/expected_unit_tests_R7") +Doppler_test = np.loadtxt(PATH_EXPECTED_OUTPUT / "Doppler.txt") +delET_test = np.loadtxt(PATH_EXPECTED_OUTPUT / "delET.txt") +i_corr_test = 921 + + +def test_find_corr_threshold_index() -> None: + expected_output = i_corr_test + threshold = 3.9e6 + distance = np.loadtxt(PATH_EXPECTED_OUTPUT / "distance.txt") + output = R7_Offset_Correction.find_corr_threshold_index(distance, threshold) + assert output == expected_output + + +def test_fit_Doppler_bias() -> None: + poly_coefs = np.array([6.182063507272804e-06, 0.05037674850270424]) + Doppler_biasfit = np.loadtxt(PATH_EXPECTED_OUTPUT / "Doppler_biasfit.txt") + + i_integral = 0 + + output = R7_Offset_Correction.fit_Doppler_bias( + delET_test, Doppler_test, i_integral, i_corr_test, Type="Linear" + ) + # The polynomial fit coefficients should compare exactly + assert np.array_equal(output[0], poly_coefs) + + assert np.array_equal(output[1], Doppler_biasfit) -- GitLab From e4be282ce7791e847c525f5f71c7743c4e3a2ffd Mon Sep 17 00:00:00 2001 From: RomanG Date: Tue, 15 Feb 2022 15:30:20 +0100 Subject: [PATCH 11/39] Remove useless imports --- tests/unit_tests_R7.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/unit_tests_R7.py b/tests/unit_tests_R7.py index 6d938be..e4a04bc 100644 --- a/tests/unit_tests_R7.py +++ b/tests/unit_tests_R7.py @@ -5,13 +5,11 @@ radiocc """ from pathlib import Path -from turtle import distance import numpy as np import radiocc from radiocc.old import R7_Offset_Correction -from radiocc.utils import number_length, round_to_N_significant_figures PATH_EXPECTED_OUTPUT = Path("tests/expected_unit_tests_R7") Doppler_test = np.loadtxt(PATH_EXPECTED_OUTPUT / "Doppler.txt") -- GitLab From ac7e20c7af8bd53f780c1c3e2f76cb85709a4f2e Mon Sep 17 00:00:00 2001 From: RomanG Date: Tue, 15 Feb 2022 23:19:39 +0100 Subject: [PATCH 12/39] Refactor R17. Removed unused part of the code. --- radiocc/old/R17_Run_Routines1.py | 175 +++++++++++++++++++++++++++++++ radiocc/reconstruct.py | 70 +------------ 2 files changed, 179 insertions(+), 66 deletions(-) diff --git a/radiocc/old/R17_Run_Routines1.py b/radiocc/old/R17_Run_Routines1.py index 5f3bd67..7e44565 100644 --- a/radiocc/old/R17_Run_Routines1.py +++ b/radiocc/old/R17_Run_Routines1.py @@ -8,6 +8,25 @@ Created on Tue Sep 1 12:14:39 2020 from pudb import set_trace as bp # noqa: F401 +from radiocc.constants import ( + CC, + PI, + G, + M_Planet, + R_Planet, + T_BC_low, + T_BC_med, + T_BC_upp, + Threshold_neut_upp, + c, + e, + eps0, + kB, + ks, + kx, + m_bar, + me, +) from radiocc.old import ( R9_BendAng_ImpParam_dn, R9_BendAng_ImpParam_up, @@ -21,6 +40,162 @@ from radiocc.old import ( def run( + Doppler, + Doppler_debias, + distance, + Diff_doppler, + item, + i_integral, + i_Surface, + PLOT_DIR, + state_occ_SC_dn, + state_occ_GS_dn, + gamma_dn, + delta_s_dn, + Bande, + fsup, +): + + Doppler_debias_dn_iono_x = Doppler_debias / (1.0 + (kx ** 2)) + Doppler_debias_dn_iono_s = Doppler_debias / (1.0 + (ks ** 2)) + Doppler_debias_dn_atmo_x = Doppler_debias / 2.0 + Doppler_debias_dn_atmo_s = Doppler_debias / 2.0 + Diff_Doppler_debias_dn_iono_s = Diff_doppler / (1.0 - (9.0 / 121.0)) + + N_data = Doppler.shape[0] + + 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_dn = 0.5 * PI - delta_s_dn + + if Bande == "X": + + if item == "IONO": + Doppler_debias_dn = Doppler_debias_dn_iono_x + if item == "ATMO": + Doppler_debias_dn = Doppler_debias_dn_atmo_x + + if Bande == "S": + + if item == "IONO": + Doppler_debias_dn = Doppler_debias_dn_iono_s + if item == "ATMO": + Doppler_debias_dn = Doppler_debias_dn_atmo_s + + if Bande == "Diff": + + if item == "IONO": + Doppler_debias_dn = Diff_Doppler_debias_dn_iono_s + + if item == "ATMO": + Doppler_debias_dn = Doppler_debias_dn_atmo_s + + # R8: bending angle & impact parameter down + (imp_param, bend_ang, 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, + Bande, + fsup, + c, + ) + + # R11 + ( + ref_index, + refractivity, + bend_radius, + Sum_tot, + ) = R11_Refractivity_and_Bending_Radius_v2.Abel_Analytical( + N_data, i_integral, i_Surface, imp_param, bend_ang, R_Planet + ) + + if item == "IONO": + + Ne = R12_Electron_Density.Elec( + N_data, refractivity, e, eps0, me, Bande, fsup, c + ) + TEC = R12_Electron_Density.TEC_Calc( + N_data, Ne, i_integral, i_Surface, bend_radius, R_Planet + ) + + R14_Plot_check.PLOT1( + distance, Doppler, Doppler_debias, refractivity, Ne, PLOT_DIR, item + ) + + return refractivity, Ne, TEC + + if item == "ATMO": + + # R13: neutral density (Level 04) + Nn, i_neut_upp = R13_Neutral_Number_Density.Neut( + N_data, bend_radius, refractivity, Threshold_neut_upp, i_Surface, CC, kB + ) + # R13a: pressure & temperature (Level 04) + P_low, T_low = R13a_Pressure_and_Temperature.Temp_Pre( + N_data, + T_BC_low, + Nn, + bend_radius, + i_neut_upp, + i_Surface, + kB, + G, + m_bar, + M_Planet, + refractivity, + ) + P_med, T_med = R13a_Pressure_and_Temperature.Temp_Pre( + N_data, + T_BC_med, + Nn, + bend_radius, + i_neut_upp, + i_Surface, + kB, + G, + m_bar, + M_Planet, + refractivity, + ) + P_upp, T_upp = R13a_Pressure_and_Temperature.Temp_Pre( + N_data, + T_BC_upp, + Nn, + bend_radius, + i_neut_upp, + i_Surface, + kB, + G, + m_bar, + M_Planet, + refractivity, + ) + + R14_Plot_check.PLOT2( + distance, Doppler, Doppler_debias, refractivity, Nn, PLOT_DIR, item + ) + + return refractivity, Nn, i_neut_upp, P_low, P_med, P_upp, T_low, T_med, T_upp + + +def run_old( Doppler, Doppler_debias, distance, diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index 40363a3..e997765 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -288,44 +288,13 @@ def run( # noqa: C901 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, - 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, + state_occ_SC_dn, + state_occ_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, ) if item == "ATMO": ( @@ -346,44 +315,13 @@ def run( # noqa: C901 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, - 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, + state_occ_SC_dn, + state_occ_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, ) if radiocc.gui is not None: -- GitLab From c87b9c66eb5af88d202f5509866864ce7470b215 Mon Sep 17 00:00:00 2001 From: RomanG Date: Wed, 16 Feb 2022 00:27:33 +0100 Subject: [PATCH 13/39] Refactor R12. Adapt R17 --- radiocc/old/R12_Electron_Density.py | 43 ++++++++++++----------------- radiocc/old/R17_Run_Routines1.py | 9 ++---- 2 files changed, 21 insertions(+), 31 deletions(-) diff --git a/radiocc/old/R12_Electron_Density.py b/radiocc/old/R12_Electron_Density.py index 6fad98f..7b18970 100644 --- a/radiocc/old/R12_Electron_Density.py +++ b/radiocc/old/R12_Electron_Density.py @@ -8,45 +8,38 @@ Created on Fri Nov 8 10:46:27 2019 import numpy as np -def Elec(N_data, refractivity, e, eps0, me,Bande,fsup,c): +from radiocc.constants import c, e, eps0, ks, kx, me - if Bande == 'X': - f_TX=(fsup)*(880./749.) # downward X-band - if Bande == 'S': - f_TX=(fsup)*(240./749.) # downward S-band - if Bande == 'Diff': - f_TX=(fsup)*(240./749.) # downward S-band +def Elec(refractivity, Bande, fsup): - Ne = np.full(N_data,np.nan)#np.array([np.nan for i in range( N_data )], float) + if Bande == "X": + f_TX = fsup * kx # downward X-band + if Bande == "S": + f_TX = fsup * ks # downward S-band + if Bande == "Diff": + f_TX = fsup * ks # downward S-band - for i in range(N_data): + Ne = -refractivity * 1e-6 / 40.31 * f_TX ** 2 - Ne[i] = - (refractivity[i]*1E-6)/40.31*f_TX**2 - #Ne_dn[i] = -2*np.pi*(refractivity[i]*1E-6)/(2.818E-15*(c/f_TX)**2) + # Ne_dn[i] = -2*np.pi*(refractivity[i]*1E-6)/(2.818E-15*(c/f_TX)**2) + + # Ne_dn[i] = - refractivity[i]*(8*np.pi**2*me*eps0*f_TX**2)/(e**2)*1E-6 - #Ne_dn[i] = - refractivity[i]*(8*np.pi**2*me*eps0*f_TX**2)/(e**2)*1E-6 - #print(f_TX) - #print(Ne[1600]) return Ne -def TEC_Calc(N_data, Ne, Corr_ind, i_Surface, bend_radius, R_Mars): +def TEC_Calc(Ne, Corr_ind, i_Surface, bend_radius): - TEC = np.full(N_data,np.nan)#np.array([np.nan for i in range( N_data )], float) - TEC_loc = 0. - for i in range(Corr_ind+1, i_Surface): + TEC = np.full_like(Ne, np.nan) + TEC_loc = 0.0 + for i in range(Corr_ind + 1, i_Surface): if Ne[i] > 0: - TEC_loc += Ne[i] * ((bend_radius[i-1]-bend_radius[i])) / 1.0e16 + TEC_loc += Ne[i] * (bend_radius[i - 1] - bend_radius[i]) / 1.0e16 TEC[i] = TEC_loc - #print np.shape(TEC) - #print('') - #print('') - print('\t ELECTRON DENSITY DONE') - #print('') - #print('') + print("\t ELECTRON DENSITY DONE") return TEC diff --git a/radiocc/old/R17_Run_Routines1.py b/radiocc/old/R17_Run_Routines1.py index 7e44565..57d381b 100644 --- a/radiocc/old/R17_Run_Routines1.py +++ b/radiocc/old/R17_Run_Routines1.py @@ -128,12 +128,9 @@ def run( if item == "IONO": - Ne = R12_Electron_Density.Elec( - N_data, refractivity, e, eps0, me, Bande, fsup, c - ) - TEC = R12_Electron_Density.TEC_Calc( - N_data, Ne, i_integral, i_Surface, bend_radius, R_Planet - ) + Ne = R12_Electron_Density.Elec(refractivity, Bande, fsup) + + TEC = R12_Electron_Density.TEC_Calc(Ne, i_integral, i_Surface, bend_radius) R14_Plot_check.PLOT1( distance, Doppler, Doppler_debias, refractivity, Ne, PLOT_DIR, item -- GitLab From 59bb15881fd2e9a7cbede2cf8e713600b7510338 Mon Sep 17 00:00:00 2001 From: RomanG Date: Wed, 16 Feb 2022 12:10:11 +0100 Subject: [PATCH 14/39] Add docstring to functions --- radiocc/old/R4_State_Vectors_twoway_v2.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/radiocc/old/R4_State_Vectors_twoway_v2.py b/radiocc/old/R4_State_Vectors_twoway_v2.py index c336d19..388566e 100644 --- a/radiocc/old/R4_State_Vectors_twoway_v2.py +++ b/radiocc/old/R4_State_Vectors_twoway_v2.py @@ -18,6 +18,8 @@ import radiocc # noqa def get_epoch_at_target(etobs: NDArray, obs: int, direct: str, targ: int) -> NDArray: + """Compute the transmit (or receive) epoch of a signal at a target specified by its + NAIF ID given the receive (or transmit) epoch at an observer specified by its NAIF ID.""" epoch_targ = np.full_like(etobs, np.nan) for i, et in enumerate(etobs): epoch_targ[i], _ = spy.ltime(et, obs, direct, targ) @@ -27,6 +29,8 @@ def get_epoch_at_target(etobs: NDArray, obs: int, direct: str, targ: int) -> NDA def get_state_at_target( etobs: NDArray, obs: int, targ: int, ref: str, abcorr: str ) -> NDArray: + """Compute the state vector of a target body specified by its NAIF ID relative + to an observing body specified by its NAIF ID at a given epoch.""" state_targ = np.array(spy.spkezr(str(targ), etobs, ref, abcorr, str(obs))[0]) * 1e3 return state_targ @@ -34,6 +38,10 @@ def get_state_at_target( def ephemerides( ET: NDArray, ID_SC: int, ID_GS: int, ID_P: int, Ref_frame: str, abcorr: str ) -> Tuple[NDArray, NDArray, NDArray, NDArray, NDArray, NDArray]: + """Compute the state vectors of the spacecraft and the ground station relative + to the Planet center at the epoch at which the uplink signal is received by + the SC and the epoch at which the downlink signal is received by the Planet + from the SC""" # ================================ # Signal emission/reception epochs -- GitLab From a39ff93c730886f67a3e0736cdc80bb5427c19b5 Mon Sep 17 00:00:00 2001 From: RomanG Date: Wed, 16 Feb 2022 12:30:11 +0100 Subject: [PATCH 15/39] Add the frequency ratios kx and ks --- radiocc/constants.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/radiocc/constants.py b/radiocc/constants.py index 612fdfa..011237f 100644 --- a/radiocc/constants.py +++ b/radiocc/constants.py @@ -4,6 +4,7 @@ Physical constants. """ +PI = 3.141592653589793 c = 299792458.0 # m s^-1 kB = 1.38064852e-23 # m^2 kg s^-2 K^-1 G = 6.673e-11 # m^3 kg^-1 s^-2 @@ -28,6 +29,10 @@ T_BC_low = 130 # Ref: Patzol data processing files T_BC_med = 165 # Ref : Patzol data processing files T_BC_upp = 200 # Ref : Patzol data processing files +# transponder ratios for downlink signal +kx = 880 / 749 # X-band +ks = 240 / 749 # S-band + SC = "MEX" # '-41' Planet = "MARS" # '499' -- GitLab From 634e159c9e3b473bbe9b1abf8b3efd251c5a903b Mon Sep 17 00:00:00 2001 From: RomanG Date: Thu, 17 Feb 2022 14:38:44 +0100 Subject: [PATCH 16/39] Add ephemerides.py (previously R4) --- radiocc/ephemerides.py | 80 ++++++++++++++++++++++++++++++++++++++++++ radiocc/reconstruct.py | 5 ++- tests/unit_tests_R4.py | 10 +++--- 3 files changed, 87 insertions(+), 8 deletions(-) create mode 100644 radiocc/ephemerides.py diff --git a/radiocc/ephemerides.py b/radiocc/ephemerides.py new file mode 100644 index 0000000..a72ece5 --- /dev/null +++ b/radiocc/ephemerides.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 + +""" +Compute ephemerides (epoch and state vector) of the ground station and the spacecraft +relative to the planet center at uplink/downlink reception/transmission time. +""" + + +from typing import Tuple + +import numpy as np +import spiceypy as spy +from nptyping import NDArray +from pudb import set_trace as bp # noqa + +import radiocc # noqa + + +def get_epoch_at_target(etobs: NDArray, obs: int, direct: str, targ: int) -> NDArray: + """Compute the transmit (or receive) epoch of a signal at a target specified by its + NAIF ID given the receive (or transmit) epoch at an observer specified by its NAIF ID.""" + epoch_targ = np.full_like(etobs, np.nan) + for i, et in enumerate(etobs): + epoch_targ[i], _ = spy.ltime(et, obs, direct, targ) + return epoch_targ + + +def get_state_at_target( + etobs: NDArray, obs: int, targ: int, ref: str, abcorr: str +) -> NDArray: + """Compute the state vector of a target body specified by its NAIF ID relative + to an observing body specified by its NAIF ID at a given epoch.""" + state_targ = np.array(spy.spkezr(str(targ), etobs, ref, abcorr, str(obs))[0]) * 1e3 + return state_targ + + +def get_ephemerides( + ET: NDArray, ID_SC: int, ID_GS: int, ID_P: int, Ref_frame: str, abcorr: str +) -> Tuple[NDArray, NDArray, NDArray, NDArray, NDArray, NDArray]: + """Compute the state vectors of the spacecraft and the ground station relative + to the Planet center at the epoch at which the uplink signal is received by + the SC and the epoch at which the downlink signal is received by the Planet + from the SC""" + + # ================================ + # Signal emission/reception epochs + # ================================ + + # Epoch at which SC emits downlink signal to GS = Epoch at which SC receives uplink signal from GS + tSC = get_epoch_at_target(ET, ID_GS, "<-", ID_SC) + + # Epoch at which the downlink signal is received by Planet from SC + t0_dn = get_epoch_at_target(tSC, ID_SC, "->", ID_P) + + # Epoch at which the uplink signal is received by SC from Planet + t0_up = get_epoch_at_target(tSC, ID_SC, "<-", ID_P) + + # Didn't really undestand what it was yet. + ET_SC_P = get_epoch_at_target(tSC, ID_P, "<-", ID_SC) + + # ==================================================== + # State vectors of GS and SC relative to Planet center + # ==================================================== + + # State vector of GS wrt Planet at t0_dn + state_GS_dn = get_state_at_target(t0_dn, ID_P, ID_GS, Ref_frame, abcorr) + + # State vector of SC wrt Planet at t0_dn + state_SC_dn = get_state_at_target(t0_dn, ID_P, ID_SC, Ref_frame, abcorr) + + # State vector of GS wrt Planet at t0_up + state_GS_up = get_state_at_target(t0_up, ID_P, ID_GS, Ref_frame, abcorr) + + # State vector of SC wrt Planet at t0_up + state_SC_up = get_state_at_target(t0_up, ID_P, ID_SC, Ref_frame, abcorr) + + # Used for routine R6, not in R9 ! + state_SC = get_state_at_target(ET_SC_P, ID_P, ID_SC, Ref_frame, abcorr) + + return state_GS_up, state_SC_up, state_GS_dn, state_SC_dn, state_SC, ET_SC_P diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index e997765..aca1fae 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -12,7 +12,7 @@ import numpy from pudb import set_trace as bp # noqa import radiocc -from radiocc import constants, old +from radiocc import constants, ephemerides, old from radiocc.model import ( Export, L2Data, @@ -22,7 +22,6 @@ from radiocc.model import ( Scenario, ) from radiocc.old import ( - R4_State_Vectors_twoway_v2, R5_Foot_Print, R6_Framework_Conversion, R7_Offset_Correction, @@ -114,7 +113,7 @@ def run( # noqa: C901 state_MEX_dn, state_MEX, ET_MEX_Mars, - ) = R4_State_Vectors_twoway_v2.ephemerides( + ) = ephemerides.get_ephemerides( ET, l2_data.SPACECRAFT_NAIF_CODE, l2_data.dsn_station_naif_code, diff --git a/tests/unit_tests_R4.py b/tests/unit_tests_R4.py index e59fd4a..bc25774 100644 --- a/tests/unit_tests_R4.py +++ b/tests/unit_tests_R4.py @@ -10,7 +10,7 @@ import numpy as np import spiceypy as spy import radiocc -from radiocc.old import R4_State_Vectors_twoway_v2 +from radiocc import ephemerides MARS_NAIF_ID = 499 EARTH_NAIF_ID = 399 @@ -25,9 +25,8 @@ def test_get_epoch() -> None: for FILE in FILES: spy.furnsh(str(FILE)) ET = np.linspace(0, 3600, LENGTH_TEST_ARRAY) - result = R4_State_Vectors_twoway_v2.get_epoch_at_target( - ET, EARTH_NAIF_ID, "->", MARS_NAIF_ID - ) + result = ephemerides.get_epoch_at_target(ET, EARTH_NAIF_ID, "->", MARS_NAIF_ID) + spy.kclear() # Test if output has the right dimension assert result.shape == (LENGTH_TEST_ARRAY,) # Test if output doesn't contain NaN's @@ -41,9 +40,10 @@ def test_get_state() -> None: for FILE in FILES: spy.furnsh(str(FILE)) ET = np.linspace(0, 3600, LENGTH_TEST_ARRAY) - result = R4_State_Vectors_twoway_v2.get_state_at_target( + result = ephemerides.get_state_at_target( ET, EARTH_NAIF_ID, MARS_NAIF_ID, "J2000", "NONE" ) + spy.kclear() # Test if output has the right dimension assert result.shape == (LENGTH_TEST_ARRAY, radiocc.constants.DIMENSION_STATE) # Test if output doesn't contain NaN's -- GitLab From 1050645f79d600cb407ae77585069f5bdae6926d Mon Sep 17 00:00:00 2001 From: RomanG Date: Thu, 17 Feb 2022 14:40:53 +0100 Subject: [PATCH 17/39] Minor refactor R12, R13, R13a. Major refactor R17 --- radiocc/old/R12_Electron_Density.py | 17 ++-- radiocc/old/R13_Neutral_Number_Density.py | 20 ++++- radiocc/old/R13a_Pressure_and_Temperature.py | 81 ++++++++++++++------ radiocc/old/R17_Run_Routines1.py | 46 +---------- 4 files changed, 89 insertions(+), 75 deletions(-) diff --git a/radiocc/old/R12_Electron_Density.py b/radiocc/old/R12_Electron_Density.py index 7b18970..87b73de 100644 --- a/radiocc/old/R12_Electron_Density.py +++ b/radiocc/old/R12_Electron_Density.py @@ -7,12 +7,13 @@ Created on Fri Nov 8 10:46:27 2019 """ import numpy as np +from nptyping import NDArray -from radiocc.constants import c, e, eps0, ks, kx, me +from radiocc.constants import ks, kx -def Elec(refractivity, Bande, fsup): - +def Elec(refractivity: NDArray, Bande: NDArray, fsup: float) -> NDArray: + """Compute the electron density in the ionosphere from the refractivity.""" if Bande == "X": f_TX = fsup * kx # downward X-band if Bande == "S": @@ -20,6 +21,7 @@ def Elec(refractivity, Bande, fsup): if Bande == "Diff": f_TX = fsup * ks # downward S-band + # Wang et al. 2018 Eq. 12 Ne = -refractivity * 1e-6 / 40.31 * f_TX ** 2 # Ne_dn[i] = -2*np.pi*(refractivity[i]*1E-6)/(2.818E-15*(c/f_TX)**2) @@ -29,14 +31,15 @@ def Elec(refractivity, Bande, fsup): return Ne -def TEC_Calc(Ne, Corr_ind, i_Surface, bend_radius): - +def TEC_Calc( + Ne: NDArray, Corr_ind: int, i_Surface: int, bend_radius: NDArray +) -> NDArray: + """Compute the Total Electron Content along the ray in the ionosphere""" TEC = np.full_like(Ne, np.nan) TEC_loc = 0.0 - for i in range(Corr_ind + 1, i_Surface): + for i in range(Corr_ind + 1, i_Surface): if Ne[i] > 0: - TEC_loc += Ne[i] * (bend_radius[i - 1] - bend_radius[i]) / 1.0e16 TEC[i] = TEC_loc diff --git a/radiocc/old/R13_Neutral_Number_Density.py b/radiocc/old/R13_Neutral_Number_Density.py index 52a88fb..3daf757 100644 --- a/radiocc/old/R13_Neutral_Number_Density.py +++ b/radiocc/old/R13_Neutral_Number_Density.py @@ -1,9 +1,24 @@ import numpy as np from pudb import set_trace as bp # noqa:F401 +from radiocc.constants import CC -# import matplotlib.pyplot as plt -def Neut(N_data, bend_radius, refractivity, d_neut_upp, i_Surface, C1, kB): + +def Neut(bend_radius, refractivity, d_neut_upp, i_Surface): + + i_neut_upp = np.where(bend_radius >= d_neut_upp)[0][-2] + + Nn = np.full_like(refractivity, np.nan) + + for i in range(i_neut_upp, i_Surface): + Nn[i] = refractivity[i] / CC * 1e-6 + + print("\t NEUTRAL DENSITY DONE") + + return Nn, i_neut_upp + + +def Neut_old(N_data, bend_radius, refractivity, d_neut_upp, i_Surface, C1, kB): i_neut_upp = np.where(bend_radius >= d_neut_upp)[0][-2] Nn = np.array([np.nan for i in range(N_data)], float) @@ -16,5 +31,6 @@ def Neut(N_data, bend_radius, refractivity, d_neut_upp, i_Surface, C1, kB): print("\t NEUTRAL DENSITY DONE") # print('') # print('') + bp() return Nn, i_neut_upp diff --git a/radiocc/old/R13a_Pressure_and_Temperature.py b/radiocc/old/R13a_Pressure_and_Temperature.py index 52a41b3..a7a8fbf 100644 --- a/radiocc/old/R13a_Pressure_and_Temperature.py +++ b/radiocc/old/R13a_Pressure_and_Temperature.py @@ -1,36 +1,69 @@ import numpy as np -#import scipy as sc -def Temp_Pre(N_data, T_BC, Nn, bend_radius, i_neut_upp, i_Surface, kB, G, m_bar, M_Mars,refractivity): +from radiocc.constants import G, M_Planet, kB, m_bar + +# import scipy as sc + + +def Temp_Pre(T_BC, Nn, bend_radius, i_neut_upp, i_Surface): + P = np.full_like(Nn, np.nan) + T = np.full_like(Nn, np.nan) - P = np.array([np.nan for i in range( N_data )], float) - T = np.array([np.nan for i in range( N_data )], float) - - #T[i_neut_upp] = T_BC P[i_neut_upp] = T_BC * Nn[i_neut_upp] * kB - for i in range(i_neut_upp+1, i_Surface): - - gg = G * M_Mars / (bend_radius[i])**2 - #gg1 = G * M_Mars / (bend_radius[i-1])**2 + for i in range(i_neut_upp + 1, i_Surface): - dz = (bend_radius[i-1] - bend_radius[i]) + gg = G * M_Planet / (bend_radius[i]) ** 2 + dz = bend_radius[i - 1] - bend_radius[i] - P[i] = P[i-1] + 0.5 * (Nn[i] + Nn[i-1]) * gg * m_bar * dz - # P[i] = P[i-1] + 0.5 * (Nn[i]*gg + Nn[i-1]*gg1) * m_bar * dz - # T[i] = T[i-1]*refractivity[i_neut_upp]/refractivity[i-1] + T[i-1] + 0.5 * (Nn[i]*gg + Nn[i-1]*gg1) * m_bar * dz / ( Nn[i] * kB) - # T[i] = T[i-1] + 0.5 * (Nn[i]*gg + Nn[i-1]*gg1) * m_bar * dz / ( Nn[i] * kB) + P[i] = P[i - 1] + 0.5 * (Nn[i] + Nn[i - 1]) * gg * m_bar * dz - for i in range(N_data): + T = P / (Nn * kB) - #P[i] = T[i] * ( Nn[i] * kB ) - T[i] = P[i] / ( Nn[i] * kB ) - return P, T - #print('') - #print('') - print('\tR14: Done') - #print('') - #print('') +def Temp_Pre_old( + N_data, + T_BC, + Nn, + bend_radius, + i_neut_upp, + i_Surface, + kB, + G, + m_bar, + M_Mars, + refractivity, +): + + P = np.array([np.nan for i in range(N_data)], float) + T = np.array([np.nan for i in range(N_data)], float) + + # T[i_neut_upp] = T_BC + P[i_neut_upp] = T_BC * Nn[i_neut_upp] * kB + + for i in range(i_neut_upp + 1, i_Surface): + + gg = G * M_Mars / (bend_radius[i]) ** 2 + # gg1 = G * M_Mars / (bend_radius[i-1])**2 + + dz = bend_radius[i - 1] - bend_radius[i] + + P[i] = P[i - 1] + 0.5 * (Nn[i] + Nn[i - 1]) * gg * m_bar * dz + # P[i] = P[i-1] + 0.5 * (Nn[i]*gg + Nn[i-1]*gg1) * m_bar * dz + # T[i] = T[i-1]*refractivity[i_neut_upp]/refractivity[i-1] + T[i-1] + 0.5 * (Nn[i]*gg + Nn[i-1]*gg1) * m_bar * dz / ( Nn[i] * kB) + # T[i] = T[i-1] + 0.5 * (Nn[i]*gg + Nn[i-1]*gg1) * m_bar * dz / ( Nn[i] * kB) + + for i in range(N_data): + + # P[i] = T[i] * ( Nn[i] * kB ) + T[i] = P[i] / (Nn[i] * kB) + + return P, T + + # print('') + # print('') + print("\tR14: Done") + # print('') + # print('') diff --git a/radiocc/old/R17_Run_Routines1.py b/radiocc/old/R17_Run_Routines1.py index 57d381b..8f0db4e 100644 --- a/radiocc/old/R17_Run_Routines1.py +++ b/radiocc/old/R17_Run_Routines1.py @@ -9,23 +9,15 @@ Created on Tue Sep 1 12:14:39 2020 from pudb import set_trace as bp # noqa: F401 from radiocc.constants import ( - CC, PI, - G, - M_Planet, R_Planet, T_BC_low, T_BC_med, T_BC_upp, Threshold_neut_upp, c, - e, - eps0, - kB, ks, kx, - m_bar, - me, ) from radiocc.old import ( R9_BendAng_ImpParam_dn, @@ -142,47 +134,17 @@ def run( # R13: neutral density (Level 04) Nn, i_neut_upp = R13_Neutral_Number_Density.Neut( - N_data, bend_radius, refractivity, Threshold_neut_upp, i_Surface, CC, kB + bend_radius, refractivity, Threshold_neut_upp, i_Surface ) # R13a: pressure & temperature (Level 04) P_low, T_low = R13a_Pressure_and_Temperature.Temp_Pre( - N_data, - T_BC_low, - Nn, - bend_radius, - i_neut_upp, - i_Surface, - kB, - G, - m_bar, - M_Planet, - refractivity, + T_BC_low, Nn, bend_radius, i_neut_upp, i_Surface ) P_med, T_med = R13a_Pressure_and_Temperature.Temp_Pre( - N_data, - T_BC_med, - Nn, - bend_radius, - i_neut_upp, - i_Surface, - kB, - G, - m_bar, - M_Planet, - refractivity, + T_BC_med, Nn, bend_radius, i_neut_upp, i_Surface ) P_upp, T_upp = R13a_Pressure_and_Temperature.Temp_Pre( - N_data, - T_BC_upp, - Nn, - bend_radius, - i_neut_upp, - i_Surface, - kB, - G, - m_bar, - M_Planet, - refractivity, + T_BC_upp, Nn, bend_radius, i_neut_upp, i_Surface ) R14_Plot_check.PLOT2( -- GitLab From 0e8f8acaa980be66c4b6d209f219367bd7f1aca4 Mon Sep 17 00:00:00 2001 From: RomanG Date: Thu, 17 Feb 2022 21:25:19 +0100 Subject: [PATCH 18/39] Finish refactor of R6 and R7 with unit tests. Renaming. + Replaced the R6 routine by radiocc/occultation_plane.py + Replaced the R7 routine by radiocc/offset_correction.py + Adapted reconstruct.py and unit tests --- radiocc/occultation_plane.py | 114 ++++++++++++++++++ radiocc/offset_correction.py | 77 ++++++++++++ radiocc/reconstruct.py | 12 +- tests/unit_tests_occultation_plane.py | 164 ++++++++++++++++++++++++++ tests/unit_tests_offset_correction.py | 76 ++++++++++++ 5 files changed, 436 insertions(+), 7 deletions(-) create mode 100644 radiocc/occultation_plane.py create mode 100644 radiocc/offset_correction.py create mode 100644 tests/unit_tests_occultation_plane.py create mode 100644 tests/unit_tests_offset_correction.py diff --git a/radiocc/occultation_plane.py b/radiocc/occultation_plane.py new file mode 100644 index 0000000..0b681f6 --- /dev/null +++ b/radiocc/occultation_plane.py @@ -0,0 +1,114 @@ +#!/usr/bin/env python3 + +from typing import Tuple + +import numpy as np +import spiceypy as spy +from nptyping import NDArray +from pudb import set_trace as bp + +from radiocc.constants import POSITION_INDICES, VELOCITY_INDICES + + +def compute_SC_angle(pos_SC: NDArray, pos_GS: NDArray) -> NDArray: + """Compute the angle between the r-axis of the occultation plane + and the direction vector of Mars to the spacecraft. This angle is + defined as gamma in Fig. 20 of Fjeldbo 1971.""" + gamma = np.full(pos_SC.shape[0], np.nan) + for i, (p_SC, p_GS) in enumerate(zip(pos_SC, pos_GS)): + gamma[i] = spy.vsep(p_SC, p_GS) - 0.5 * np.pi + return gamma + + +def position_in_occ_plane( + pos_SC: NDArray, pos_GS: NDArray, gamma: NDArray +) -> Tuple[NDArray, NDArray, NDArray]: + """Transform the direction vectors of Mars to the spacecraft and + Mars to the ground station into the occultation coordinate system (r,z).""" + r_SC = np.linalg.norm(pos_SC, axis=1) * np.cos(gamma) + z_SC = np.linalg.norm(pos_SC, axis=1) * np.sin(gamma) + z_GS = -np.linalg.norm(pos_GS, axis=1) + return r_SC, z_SC, z_GS + + +def transition_matrix_occ_plane(p_SC: NDArray, p_GS: NDArray) -> NDArray: + """Compute the transition matrix from the equatorial geocentric basis (x,y,z) + to the occultation plane basis (n, r, z).""" + n_vec = np.cross(p_GS, p_SC) # vector normal to occ. plane + n_vec /= np.linalg.norm(n_vec) + r_vec = -np.cross(p_GS, n_vec) # r-axis of occ. plane (r = z x n) + r_vec /= np.linalg.norm(r_vec) + z_vec = -p_GS / np.linalg.norm(p_GS) + + base_matrix = np.stack( + (n_vec, r_vec, z_vec), axis=1 + ) # Transition matrix from equatorial frame to occ. plane + return np.linalg.inv(base_matrix) + + +def velocity_in_occ_plane( + pos_SC: NDArray, pos_GS: NDArray, vel_SC: NDArray, vel_GS: NDArray +) -> Tuple[NDArray, NDArray, NDArray, NDArray]: + """Using a transition matrix, transform the velocity vectors of the spacecraft + and the ground station relative to Mars center into the occultation coordinate system.""" + vr_SC = np.full(vel_SC.shape[0], np.nan) + vz_SC = np.full_like(vr_SC, np.nan) + vr_GS = np.full_like(vr_SC, np.nan) + vz_GS = np.full_like(vr_SC, np.nan) + + for i, (p_SC, p_GS, v_SC, v_GS) in enumerate(zip(pos_SC, pos_GS, vel_SC, vel_GS)): + trans_mat = transition_matrix_occ_plane(p_SC, p_GS) + _, vr_SC[i], vz_SC[i] = trans_mat.dot(v_SC) + _, vr_GS[i], vz_GS[i] = trans_mat.dot(v_GS) + vr_SC[i] = -np.abs(vr_SC[i]) + vz_SC[i] = np.abs(vz_SC[i]) + vr_GS[i] = np.abs(vr_GS[i]) + vz_GS[i] = np.abs(vz_GS[i]) + return vr_SC, vz_SC, vr_GS, vz_GS + + +def compute_delta_s(r_SC: NDArray, z_SC: NDArray, z_GS: NDArray) -> NDArray: + """Compute the angle between direction vector of the ground station to Mars + (z-axis of the occultation plane) and the direction vector of the ground station + to the spacecraft. This angle is defined as delta_s in Fig. 20 of Fjeldbo 1971.""" + delta_s = np.full(r_SC.shape, np.nan) + + for i, (r_SC_i, z_SC_i, z_GS_i) in enumerate(zip(r_SC, z_SC, z_GS)): + # delta_s[i] = spy.vsepg( + # np.array([0, -z_GS_i]), np.array([r_SC_i, z_SC_i - z_GS_i]), 2 + # ) # Defined in Fjeldbo 1971 + delta_s[i] = np.arccos( + np.dot(np.array([0, -z_GS_i]), np.array([r_SC_i, z_SC_i - z_GS_i])) + / np.abs(z_GS_i) + / np.linalg.norm(np.array([r_SC_i, z_SC_i - z_GS_i])) + ) + return delta_s + + +def Occ_Plane_Conversion( + state_SC: NDArray, state_GS: NDArray +) -> Tuple[NDArray, NDArray, NDArray, NDArray]: + """ + + Transform the state vectors of the spacecraft and the ground station + relative to Mars center expressed in the equatorial geocentric frame into + the occultation plane coordinate system (r,z). + + Compute the angle between the r-axis of the occultation plane and the + Mars-SC direction vector -> gamma in Fig. 20 of Fjeldbo 1971. + + Compute the angle between the z-axis of the occultation plane and the + GS-SC direction vector -> delta_s in Fig. 20 of Fjeldbo 1971. + """ + N_data = state_SC.shape[0] + pos_SC = state_SC[:, POSITION_INDICES] + vel_SC = state_SC[:, VELOCITY_INDICES] + pos_GS = state_GS[:, POSITION_INDICES] + vel_GS = state_GS[:, VELOCITY_INDICES] + + gamma = compute_SC_angle(pos_SC, pos_GS) + r_SC, z_SC, z_GS = position_in_occ_plane(pos_SC, pos_GS, gamma) + vr_SC, vz_SC, vr_GS, vz_GS = velocity_in_occ_plane(pos_SC, pos_GS, vel_SC, vel_GS) + + state_occ_SC = np.stack((r_SC, z_SC, vr_SC, vz_SC), axis=1) + state_occ_GS = np.stack((np.zeros(N_data), z_GS, vr_GS, vz_GS), axis=1) + + delta_s = compute_delta_s(r_SC, z_SC, z_GS) + return state_occ_SC, state_occ_GS, gamma, delta_s diff --git a/radiocc/offset_correction.py b/radiocc/offset_correction.py new file mode 100644 index 0000000..ff6fc74 --- /dev/null +++ b/radiocc/offset_correction.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 + +from typing import Optional, Tuple + +import numpy as np +from nptyping import NDArray +from pudb import set_trace as bp # noqa + +import radiocc + + +def find_threshold_index(distance: NDArray, threshold: float) -> int: + """Find index of the altitude closest to the threshold.""" + for i_corr, dist in enumerate(distance): + if dist <= threshold: + break + return i_corr + + +def fit_Doppler_bias( + delET: NDArray, Doppler: NDArray, i_integral: int, i_corr: int, Type: str +) -> Tuple[NDArray, NDArray]: + """Fit the raw Doppler residuals above the bias correction threshold + by a polynomial of degree 1 or 2. + Output the coefficients of the fitting polynomial and the values of + the polynomial.""" + if Type == "Linear": + deg = 1 + elif Type == "Quadratic": + deg = 2 + else: + print("Type is not recognized") + + poly_coefs = np.polyfit(delET[i_integral:i_corr], Doppler[i_integral:i_corr], deg) + Doppler_biasfit = np.polyval(poly_coefs, delET) + return poly_coefs, Doppler_biasfit + + +def apply_correction( + Type: str, + threshold: int, + Doppler: NDArray, + ET: NDArray, + distance: NDArray, + i_integral: int, + DATA_TYPE: Optional[radiocc.model.RadioDataType], +) -> Tuple[NDArray, NDArray, NDArray, NDArray, int, NDArray]: + """Corrects the offset bias in the raw Doppler residuals by fitting the residuals + above the upper atmosphere by a Linear or Quadratic function.""" + + # If egress, we reverse the matrices to have them from atmosphere to surface. + if DATA_TYPE is not None and DATA_TYPE == radiocc.model.RadioDataType.EGRESS: + distance = distance[::-1] + Doppler = Doppler[::-1] + ET = ET[::-1] + + delET = ET - ET[0] + + # Find index of the altitude threshold for bias correction + i_corr = find_threshold_index(distance, threshold) + + # Fit raw Doppler residuals above the threshold by a polynomial + poly_coefs, Doppler_biasfit = fit_Doppler_bias( + delET, Doppler, i_integral, i_corr, Type + ) + + # Subtract the fitted bias from the raw residuals + Doppler_debias = Doppler - Doppler_biasfit + + return ( + ET, + delET, + Doppler_debias, + Doppler_biasfit, + i_corr, + poly_coefs, + ) diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index aca1fae..cd38e19 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -12,7 +12,7 @@ import numpy from pudb import set_trace as bp # noqa import radiocc -from radiocc import constants, ephemerides, old +from radiocc import constants, ephemerides, occultation_plane, offset_correction, old from radiocc.model import ( Export, L2Data, @@ -23,8 +23,6 @@ from radiocc.model import ( ) from radiocc.old import ( R5_Foot_Print, - R6_Framework_Conversion, - R7_Offset_Correction, R8_Plot_Doppler_1, R15_Errors, R16_Biasfit_Correction, @@ -128,7 +126,7 @@ def run( # noqa: C901 state_occ_GS_up, gamma_up, delta_s_up, - ) = R6_Framework_Conversion.Occ_Plane_Conversion(state_MEX_up, state_GS_up) + ) = occultation_plane.Occ_Plane_Conversion(state_MEX_up, state_GS_up) r_MEX_up = state_occ_SC_up[:, 0] z_MEX_up = state_occ_SC_up[:, 1] @@ -145,7 +143,7 @@ def run( # noqa: C901 state_occ_GS_dn, gamma_dn, delta_s_dn, - ) = R6_Framework_Conversion.Occ_Plane_Conversion(state_MEX_dn, state_GS_dn) + ) = occultation_plane.Occ_Plane_Conversion(state_MEX_dn, state_GS_dn) r_MEX_dn = state_occ_SC_dn[:, 0] z_MEX_dn = state_occ_SC_dn[:, 1] @@ -193,7 +191,7 @@ def run( # noqa: C901 Doppler_biasfit, Corr_ind, poly_coefs, - ) = R7_Offset_Correction.Off_Cor( + ) = offset_correction.apply_correction( B_Cor_Type, Threshold_Cor, Doppler, @@ -251,7 +249,7 @@ def run( # noqa: C901 Doppler_biasfit, Corr_ind, poly_coefs, - ) = R7_Offset_Correction.Off_Cor( + ) = offset_correction.apply_correction( NEW_TYPE, NEW_THRESHOLD, Doppler, diff --git a/tests/unit_tests_occultation_plane.py b/tests/unit_tests_occultation_plane.py new file mode 100644 index 0000000..ac899ec --- /dev/null +++ b/tests/unit_tests_occultation_plane.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python3 + +""" +Unit tests for the functions of occultation_plane.py +""" + +import numpy as np +import pytest + +import radiocc +from radiocc import occultation_plane +from radiocc.utils import round_to_N_significant_figures + +LENGTH_TEST_ARRAY = 2 +DIMENSION_VECTOR = radiocc.constants.DIMENSION_VECTOR + +# Test input (first 2 lines from MEX 235) +pos_SC_sample = np.array( + [ + [-7726575.13322412, -7666711.32894175, -7958119.5346874], + [-7726451.93234389, -7667213.36854496, -7957816.58827798], + ] +) +pos_GS_sample = np.array( + [ + [192537483124.12, 265959242098.085, 113482080420.543], + [192537459864.319, 265959250204.846, 113482084397.897], + ] +) + +vel_SC_sample = np.array( + [ + [240.58413463, -980.54769211, 591.63687863], + [240.65325798, -980.4791224, 591.70807734], + ] +) +vel_GS_sample = np.array( + [ + [-45427.78743875, 15832.98573322, 7768.01058944], + [-45427.7801956, 15832.99253347, 7768.00905843], + ] +) + +gamma_sample = np.array([1.2397252315812128, 1.2397669744753683]) + +r_SC_sample = np.array([4382958.64228272, 4382438.07915611]) +z_SC_sample = np.array([12751461.32641642, 12751678.38508523]) +z_GS_sample = np.array([-347394852353.5334, -347394846967.8551]) + +vr_SC_sample = np.array([-1016.962309512641, -1016.9834579503646]) +vz_SC_sample = np.array([424.0824883669944, 423.9684628262237]) +vr_GS_sample = np.array([7681.6916741349305, 7679.542472103281]) +vz_GS_sample = np.array([10518.558189401045, 10518.546132036772]) + + +@pytest.mark.parametrize( + "pos_SC_input, pos_GS_input, expected_gamma", + [(pos_SC_sample, pos_GS_sample, gamma_sample)], +) +def test_compute_SC_angle(pos_SC_input, pos_GS_input, expected_gamma) -> None: + + output = occultation_plane.compute_SC_angle(pos_SC_input, pos_GS_input) + # Test if output has the right dimension + assert output.shape == (LENGTH_TEST_ARRAY,) + # Test if output doesn't contain NaN's + assert not np.isnan(output).all() + # Test if output of test is element-wise equal to the expected output within + # a relative tolerance + assert np.allclose( + output, + expected_gamma, + atol=0, + rtol=1e-8, + equal_nan=False, + ) + + +@pytest.mark.parametrize( + "pos_SC_input, pos_GS_input, gamma_input, expected_output", + [ + ( + pos_SC_sample, + pos_GS_sample, + round_to_N_significant_figures(gamma_sample, 15), + (r_SC_sample, z_SC_sample, z_GS_sample), + ) + ], +) +def test_position_in_occ_plane( + pos_SC_input, pos_GS_input, gamma_input, expected_output +) -> None: + output = occultation_plane.position_in_occ_plane( + pos_SC_input, pos_GS_input, gamma_input + ) + # Test if output has the right dimension + assert np.array(output).shape == (3, LENGTH_TEST_ARRAY) + # Test if output doesn't contain NaN's + assert not np.isnan(output).all() + # Test if output of test is element-wise equal to the expected output within + # a relative tolerance + + assert np.allclose( + np.array(output), + np.array(expected_output), + atol=0, + rtol=1e-8, + equal_nan=False, + ) + + +def test_transition_matrix_occ_plane() -> None: + output = occultation_plane.transition_matrix_occ_plane(pos_SC_sample, pos_GS_sample) + assert not np.isnan(output).all() + + +@pytest.mark.parametrize( + "pos_SC_input, pos_GS_input, vel_SC_input, vel_GS_input, expected_output", + [ + ( + pos_SC_sample, + pos_GS_sample, + vel_SC_sample, + vel_GS_sample, + (vr_SC_sample, vz_SC_sample, vr_GS_sample, vz_GS_sample), + ) + ], +) +def test_velocity_in_occ_plane( + pos_SC_input, pos_GS_input, vel_SC_input, vel_GS_input, expected_output +) -> None: + + output = occultation_plane.velocity_in_occ_plane( + pos_SC_input, pos_GS_input, vel_SC_input, vel_GS_input + ) + assert np.allclose( + np.array(output), + np.array(expected_output), + atol=0, + rtol=1e-8, + equal_nan=False, + ) + + +@pytest.mark.parametrize( + "r_SC_input, z_SC_input, z_GS_input, expected_delta_s", + [ + ( + r_SC_sample, + z_SC_sample, + z_GS_sample, + np.array([1.2616189367258302e-05, 1.2614684479727514e-05]), + ) + ], +) +def test_compute_delta_s(r_SC_input, z_SC_input, z_GS_input, expected_delta_s) -> None: + delta_s = occultation_plane.compute_delta_s(r_SC_input, z_SC_input, z_GS_input) + + assert np.allclose( + delta_s, + expected_delta_s, + atol=0, + rtol=1e-8, + equal_nan=False, + ) diff --git a/tests/unit_tests_offset_correction.py b/tests/unit_tests_offset_correction.py new file mode 100644 index 0000000..662ee03 --- /dev/null +++ b/tests/unit_tests_offset_correction.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 + +""" +Unit tests for the functions of offset_correction.py +""" + +from pathlib import Path + +import numpy as np +import pytest + +import radiocc +from radiocc import offset_correction + +PATH_EXPECTED_OUTPUT = Path("tests/expected_unit_tests_R7") + + +@pytest.mark.parametrize( + "distance_input, threshold_input, expected_index", + [ + (np.arange(5, 0, -1, dtype=float), 3.2, 2), # threshold closer to floor + (np.arange(5, 0, -1, dtype=float), 3.8, 2), # closer to ceiling + (np.arange(1, 0, -0.1, dtype=float), 0.8, 2), # equal to one value + (np.arange(1, 0, -0.1, dtype=float), 1.5, 0), # larger than first value + (np.arange(1, 0, -0.1, dtype=float), -0.5, 9), # lower than last value + ( + pytest.param( + np.arange(1, 0, -0.2, dtype=float), + 0.9, + 0, + marks=pytest.mark.xfail, # test fail: wrong index + ) + ), + ( + np.loadtxt(PATH_EXPECTED_OUTPUT / "distance.txt"), + 3.9e6, + 921, + ), # test with values from MEX case + ], +) +def test_find_threshold_index(distance_input, threshold_input, expected_index) -> None: + assert ( + offset_correction.find_threshold_index(distance_input, threshold_input) + == expected_index + ) + + +@pytest.mark.parametrize( + "delET_input, Doppler_input, i_integral_input, i_corr_input, Type, expected_poly_coefs, expected_Doppler_biasfit", + [ + ( + np.loadtxt(PATH_EXPECTED_OUTPUT / "delET.txt"), + np.loadtxt(PATH_EXPECTED_OUTPUT / "Doppler.txt"), + 0, + 921, + "Linear", + np.array([6.182063507272804e-06, 0.05037674850270424]), + np.loadtxt(PATH_EXPECTED_OUTPUT / "Doppler_biasfit.txt"), + ) + ], +) +def test_fit_Doppler_bias( + delET_input, + Doppler_input, + i_integral_input, + i_corr_input, + Type, + expected_poly_coefs, + expected_Doppler_biasfit, +) -> None: + output = offset_correction.fit_Doppler_bias( + delET_input, Doppler_input, i_integral_input, i_corr_input, Type + ) + + assert np.array_equal(output[0], expected_poly_coefs) + assert np.array_equal(output[1], expected_Doppler_biasfit) -- GitLab From df3aab65f68ce2a975695139ba261535e01f3cbb Mon Sep 17 00:00:00 2001 From: RomanG Date: Thu, 17 Feb 2022 21:28:32 +0100 Subject: [PATCH 19/39] Rename unit tests R4. Delete outdated unit tests --- tests/unit_tests_R6.py | 127 ------------------ tests/unit_tests_R7.py | 40 ------ ..._tests_R4.py => unit_tests_ephemerides.py} | 2 +- 3 files changed, 1 insertion(+), 168 deletions(-) delete mode 100644 tests/unit_tests_R6.py delete mode 100644 tests/unit_tests_R7.py rename tests/{unit_tests_R4.py => unit_tests_ephemerides.py} (96%) diff --git a/tests/unit_tests_R6.py b/tests/unit_tests_R6.py deleted file mode 100644 index 835446c..0000000 --- a/tests/unit_tests_R6.py +++ /dev/null @@ -1,127 +0,0 @@ -#!/usr/bin/env python3 - -""" -radiocc -""" - -from pathlib import Path -from zlib import Z_SYNC_FLUSH - -import numpy as np -import spiceypy as spy - -import radiocc -from radiocc.old import R6_Framework_Conversion -from radiocc.utils import number_length, round_to_N_significant_figures - -LENGTH_TEST_ARRAY = 2 -DIMENSION_VECTOR = radiocc.constants.DIMENSION_VECTOR -PATH_EXPECTED_OUTPUT = Path("tests/expected_unit_tests_R6") - -# Test input (first 2 lines from MEX 235) -pos_SC_test = np.array( - [ - [-7726575.13322412, -7666711.32894175, -7958119.5346874], - [-7726451.93234389, -7667213.36854496, -7957816.58827798], - ] -) -pos_GS_test = np.array( - [ - [192537483124.12, 265959242098.085, 113482080420.543], - [192537459864.319, 265959250204.846, 113482084397.897], - ] -) - -gamma_test = np.array([1.2397252315812128, 1.2397669744753683]) - -r_SC_test = np.array([4382958.64228272, 4382438.07915611]) -z_SC_test = np.array([12751461.32641642, 12751678.38508523]) -z_GS_test = np.array([-347394852353.5334, -347394846967.8551]) - - -def test_compute_SC_angle() -> None: - - output = R6_Framework_Conversion.compute_SC_angle(pos_SC_test, pos_GS_test) - # Test if output has the right dimension - assert output.shape == (LENGTH_TEST_ARRAY,) - # Test if output doesn't contain NaN's - assert not np.isnan(output).all() - # Test if output of test is element-wise equal to the expected output within - # an absolute tolerance corresponding to the number of significant figures - assert np.allclose( - output, - gamma_test, - atol=0, - rtol=1e-6, - equal_nan=False, - ) - - -def test_position_in_occ_plane() -> None: - gamma_rounded = round_to_N_significant_figures(gamma_test, 15) - output = R6_Framework_Conversion.position_in_occ_plane( - pos_SC_test, pos_GS_test, gamma_rounded - ) - # Test if output has the right dimension - assert np.array(output).shape == (3, LENGTH_TEST_ARRAY) - # Test if output doesn't contain NaN's - assert not np.isnan(output).all() - # Test if output of test is element-wise equal to the expected output within - # an absolute tolerance corresponding to the number of significant figures - - expected_output = (r_SC_test, z_SC_test, z_GS_test) - assert np.allclose( - np.array(output), - np.array(expected_output), - atol=0, - rtol=1e-6, - equal_nan=False, - ) - - -def test_transition_matrix_occ_plane() -> None: - output = R6_Framework_Conversion.transition_matrix_occ_plane( - pos_SC_test, pos_GS_test - ) - assert not np.isnan(output).all() - - -def test_velocity() -> None: - vel_SC_test = np.array( - [ - [240.58413463, -980.54769211, 591.63687863], - [240.65325798, -980.4791224, 591.70807734], - ] - ) - vel_GS_test = np.array( - [ - [-45427.78743875, 15832.98573322, 7768.01058944], - [-45427.7801956, 15832.99253347, 7768.00905843], - ] - ) - - vr_SC = np.array([-1016.962309512641, -1016.9834579503646]) - output = R6_Framework_Conversion.velocity_in_occ_plane( - pos_SC_test, pos_GS_test, vel_SC_test, vel_GS_test - ) - assert np.allclose( - output[0], - vr_SC, - atol=0, - rtol=1e-10, - equal_nan=False, - ) - - -def test_compute_delta_s() -> None: - - delta_s = R6_Framework_Conversion.compute_delta_s(r_SC_test, z_SC_test, z_GS_test) - expected_delta_s = np.array([1.2616189367258302e-05, 1.2614684479727514e-05]) - - assert np.allclose( - delta_s, - expected_delta_s, - atol=0, - rtol=1e-5, - equal_nan=False, - ) diff --git a/tests/unit_tests_R7.py b/tests/unit_tests_R7.py deleted file mode 100644 index e4a04bc..0000000 --- a/tests/unit_tests_R7.py +++ /dev/null @@ -1,40 +0,0 @@ -#!/usr/bin/env python3 - -""" -radiocc -""" - -from pathlib import Path - -import numpy as np - -import radiocc -from radiocc.old import R7_Offset_Correction - -PATH_EXPECTED_OUTPUT = Path("tests/expected_unit_tests_R7") -Doppler_test = np.loadtxt(PATH_EXPECTED_OUTPUT / "Doppler.txt") -delET_test = np.loadtxt(PATH_EXPECTED_OUTPUT / "delET.txt") -i_corr_test = 921 - - -def test_find_corr_threshold_index() -> None: - expected_output = i_corr_test - threshold = 3.9e6 - distance = np.loadtxt(PATH_EXPECTED_OUTPUT / "distance.txt") - output = R7_Offset_Correction.find_corr_threshold_index(distance, threshold) - assert output == expected_output - - -def test_fit_Doppler_bias() -> None: - poly_coefs = np.array([6.182063507272804e-06, 0.05037674850270424]) - Doppler_biasfit = np.loadtxt(PATH_EXPECTED_OUTPUT / "Doppler_biasfit.txt") - - i_integral = 0 - - output = R7_Offset_Correction.fit_Doppler_bias( - delET_test, Doppler_test, i_integral, i_corr_test, Type="Linear" - ) - # The polynomial fit coefficients should compare exactly - assert np.array_equal(output[0], poly_coefs) - - assert np.array_equal(output[1], Doppler_biasfit) diff --git a/tests/unit_tests_R4.py b/tests/unit_tests_ephemerides.py similarity index 96% rename from tests/unit_tests_R4.py rename to tests/unit_tests_ephemerides.py index bc25774..bc92be6 100644 --- a/tests/unit_tests_R4.py +++ b/tests/unit_tests_ephemerides.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 """ -radiocc +Unit tests for the functions of ephemerides.py """ from pathlib import Path -- GitLab From b6b2c0a95062e5c8eebf38e70cced20eeda1fdeb Mon Sep 17 00:00:00 2001 From: RomanG Date: Fri, 18 Feb 2022 12:27:23 +0100 Subject: [PATCH 20/39] Refactor R9. Replace by ray_parameter.py --- radiocc/old/R17_Run_Routines1.py | 23 ++-- radiocc/ray_parameters.py | 180 +++++++++++++++++++++++++++++++ tests/unit_tests_R9.py | 145 +++++++++++++++++++++++++ 3 files changed, 336 insertions(+), 12 deletions(-) create mode 100644 radiocc/ray_parameters.py create mode 100644 tests/unit_tests_R9.py diff --git a/radiocc/old/R17_Run_Routines1.py b/radiocc/old/R17_Run_Routines1.py index 8f0db4e..1516223 100644 --- a/radiocc/old/R17_Run_Routines1.py +++ b/radiocc/old/R17_Run_Routines1.py @@ -8,6 +8,7 @@ Created on Tue Sep 1 12:14:39 2020 from pudb import set_trace as bp # noqa: F401 +from radiocc import ray_parameters from radiocc.constants import ( PI, R_Planet, @@ -90,22 +91,20 @@ def run( Doppler_debias_dn = Doppler_debias_dn_atmo_s # R8: bending angle & impact parameter down - (imp_param, bend_ang, 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, + + ( + imp_param, + bend_ang, + delta_r_dn, + beta_r_dn, + ) = ray_parameters.compute_bending_angle_impact_param( + state_occ_SC_dn, + state_occ_GS_dn, gamma_dn, - beta_e_dn, delta_s_dn, Doppler_debias_dn, - Bande, fsup, - c, + Bande, ) # R11 diff --git a/radiocc/ray_parameters.py b/radiocc/ray_parameters.py new file mode 100644 index 0000000..c471906 --- /dev/null +++ b/radiocc/ray_parameters.py @@ -0,0 +1,180 @@ +#!/usr/bin/env python3 + +from typing import Tuple, Union + +import numpy as np +from nptyping import NDArray +from pudb import set_trace as bp + +from radiocc.constants import c, ks, kx + +"""Compute of the two geometrical ray parameters: bending angle and impact parameter +from the occultation measurements and ephemerides of the spacecraft and ground station.""" + + +def coeffs_equation( + beta_r: float, + delta_r: float, + vrs: float, + vzs: float, + vrt: float, + vzt: float, + rs: float, + zs: float, + zt: float, + ga: float, + b_e: float, + d_s: float, +) -> Tuple[float, float, float, float]: + """Compute the coefficients of the system of 2 linearized equations + expressed in Eq. A3 and A4 of Fjeldbo 1971.""" + b11 = -vrs * np.sin(b_e - beta_r) + vzs * np.cos(b_e - beta_r) + b12 = -vrt * np.cos(d_s - delta_r) + vzt * np.sin(d_s - delta_r) + b21 = np.sqrt(rs ** 2 + zs ** 2) * np.cos(b_e - ga - beta_r) + b22 = zt * np.cos(d_s - delta_r) + + return b11, b12, b21, b22 + + +def equation_rhs( + Doppler_debias: float, + fs: float, + beta_r: float, + delta_r: float, + vrs: float, + vzs: float, + vrt: float, + vzt: float, + rs: float, + zs: float, + zt: float, + ga: float, + b_e: float, + d_s: float, +) -> Tuple[float, float]: + """Compute the right-hand side (k1 and k2) of the system of 2 linearized equations + expressed in Eq. A3 and A4 of Fjeldbo 1971.""" + k1 = ( + c * Doppler_debias / fs + + vrs * (np.cos(b_e - beta_r) - np.cos(b_e)) + + vzs * (np.sin(b_e - beta_r) - np.sin(b_e)) + - vrt * (np.sin(d_s - delta_r) - np.sin(d_s)) + - vzt * (np.cos(d_s - delta_r) - np.cos(d_s)) + ) + + k2 = zt * np.sin(d_s - delta_r) + np.sqrt(rs ** 2 + zs ** 2) * np.sin( + b_e - ga - beta_r + ) + + return k1, k2 + + +def compute_angles_correction( + b11: float, b12: float, b21: float, b22: float, k1: float, k2: float +) -> Tuple[float, float]: + """Solve Eqs A3 and A4 (Fjeldbo 1971) by substitution.""" + beta_r_corr = (k1 - b12 * k2 / b22) / (b11 - b21 / b22 * b12) + delta_r_corr = (k2 - b21 * beta_r_corr) / b22 + return beta_r_corr, delta_r_corr + + +def compute_impact_parameter( + r_SC: Union[float, NDArray], + z_SC: Union[float, NDArray], + beta_e: Union[float, NDArray], + beta_r: Union[float, NDArray], + gamma: Union[float, NDArray], +) -> Union[float, NDArray]: + """Compute the impact parameter (distance from the planetary center of mass + to the ray-path asymptote) following Eq A6 of Fjeldbo 1971.""" + return np.sqrt(r_SC ** 2 + z_SC ** 2) * np.sin(beta_e - beta_r - gamma) + + +def compute_bending_angle_impact_param( + state_occ_SC: NDArray, + state_occ_GS: NDArray, + gamma: NDArray, + delta_s: NDArray, + Doppler_debias: NDArray, + fsup: float, + Bande: str, +) -> Tuple[NDArray, NDArray, NDArray, NDArray]: + """Compute bending ange and impact parameter.""" + N_data = state_occ_SC.shape[0] + + r_SC = state_occ_SC[:, 0] + z_SC = state_occ_SC[:, 1] + vr_SC = state_occ_SC[:, 2] + vz_SC = state_occ_SC[:, 3] + + z_GS = state_occ_GS[:, 1] + vr_GS = state_occ_GS[:, 2] + vz_GS = state_occ_GS[:, 3] + + beta_e = 0.5 * np.pi - delta_s + + delta_r_array = np.full(N_data, np.nan) + beta_r_array = np.full(N_data, np.nan) + + if Bande == "X": + fs = fsup * kx # downward X-band frequency + else: + fs = fsup * ks # downward S-band frequency + + # We follow convention from Fjeldbo 1971 for coordinates and velocities in occ. plane + for i, (vrs, vzs, vrt, vzt, rs, zs, zt, ga, b_e, d_s) in enumerate( + zip(vr_SC, vz_SC, vr_GS, vz_GS, r_SC, z_SC, z_GS, gamma, beta_e, delta_s) + ): + + # We solve a system of 2 equations to find the delta_r and beta_r angles (see Fjeldbo 1971) + # We start with rays passing outside the atmosphere (delta_r = beta_r = 0) + delta_r = 0.0 + beta_r = 0.0 + + epsilon = 999 + stock = 0 + while np.abs(epsilon) > 1e-3: + # Computing coefficients of eq. A3 and A4 (Fjeldbo 1971) + b11, b12, b21, b22 = coeffs_equation( + beta_r, delta_r, vrs, vzs, vrt, vzt, rs, zs, zt, ga, b_e, d_s + ) + + k1, k2 = equation_rhs( + Doppler_debias[i], + fs, + beta_r, + delta_r, + vrs, + vzs, + vrt, + vzt, + rs, + zs, + zt, + ga, + b_e, + d_s, + ) + + # We solve the corrections for delta_r and beta_r from eq. A3 and A4 (Fjeldbo 1971) + beta_r_corr, delta_r_corr = compute_angles_correction( + b11, b12, b21, b22, k1, k2 + ) + + delta_r += delta_r_corr + beta_r += beta_r_corr + + epsilon = delta_r_corr - stock + stock = delta_r_corr + + delta_r_array[i] = delta_r + beta_r_array[i] = beta_r + + # Compute bending angle (Fjeldbo 1971, eq. A5) + bend_ang = delta_r_array + beta_r_array + # Compute impact parameter (Fjeldbo 1971, eq. A6) + imp_param = compute_impact_parameter(r_SC, z_SC, beta_e, beta_r_array, gamma) + + print("\t BENDING ANGLE UP DONE") + + return imp_param, bend_ang, delta_r_array, beta_r_array diff --git a/tests/unit_tests_R9.py b/tests/unit_tests_R9.py new file mode 100644 index 0000000..440c721 --- /dev/null +++ b/tests/unit_tests_R9.py @@ -0,0 +1,145 @@ +#!/usr/bin/env python3 + +""" +Unit tests for the functions of R9 +""" + +import numpy as np +import pytest + +import radiocc +from radiocc import ray_parameters + +test_input = { + "Doppler_debias": -0.005117963887076302, + "fs": 8420122709.105474, + "beta_r": 0, + "delta_r": 0, + "vrs": -1016.9660257949688, + "vzs": 424.0624558972604, + "vrt": 7681.31413033953, + "vzt": 10518.556071108187, + "rs": 4382867.188086733, + "zs": 12751499.46453215, + "zt": -347394851407.35156, + "ga": 1.2397325652323983, + "b_e": 1.5707837108695317, + "d_s": 1.2615925364885982e-05, +} + +test_output = { + "b11": 1016.9713756543314, + "b12": -7681.181428409908, + "b21": 12751554.757442705, + "b22": -347394851379.7056, + "k1": -0.00018222145052620506, + "k2": -1.2016121335327625, +} + + +@pytest.mark.parametrize( + "input, expected_output", + [ + ( + ( + test_input["beta_r"], + test_input["delta_r"], + test_input["vrs"], + test_input["vzs"], + test_input["vrt"], + test_input["vzt"], + test_input["rs"], + test_input["zs"], + test_input["zt"], + test_input["ga"], + test_input["b_e"], + test_input["d_s"], + ), + ( + test_output["b11"], + test_output["b12"], + test_output["b21"], + test_output["b22"], + ), + ) + ], +) +def test_coeffs_equation(input, expected_output) -> None: + output = ray_parameters.coeffs_equation(*input) + assert np.array_equal(np.array(output), np.array(expected_output)) + + +@pytest.mark.parametrize( + "input, expected_output", + [ + ( + ( + test_input["Doppler_debias"], + test_input["fs"], + test_input["beta_r"], + test_input["delta_r"], + test_input["vrs"], + test_input["vzs"], + test_input["vrt"], + test_input["vzt"], + test_input["rs"], + test_input["zs"], + test_input["zt"], + test_input["ga"], + test_input["b_e"], + test_input["d_s"], + ), + ( + test_output["k1"], + test_output["k2"], + ), + ) + ], +) +def test_equation_rhs(input, expected_output) -> None: + output = ray_parameters.equation_rhs(*input) + assert np.array_equal(np.array(output), np.array(expected_output)) + + +@pytest.mark.parametrize( + "input, expected_output", + [ + ( + ( + test_output["b11"], + test_output["b12"], + test_output["b21"], + test_output["b22"], + test_output["k1"], + test_output["k2"], + ), + ( + -1.7920406843024083e-07, + -3.1189821999626236e-12, + ), + ) + ], +) +def test_compute_angles_correction(input, expected_output) -> None: + output = ray_parameters.compute_angles_correction(*input) + assert np.array_equal(np.array(output), np.array(expected_output)) + + +@pytest.mark.parametrize( + "input, expected_output", + [ + ( + ( + test_input["rs"], + test_input["zs"], + test_input["b_e"], + -1.7920406843024083e-07, + test_input["ga"], + ), + 4382708.600902833, + ) + ], +) +def test_compute_impact_parameter(input, expected_output) -> None: + output = ray_parameters.compute_impact_parameter(*input) + assert np.array_equal(np.array(output), np.array(expected_output)) -- GitLab From ce42684a70cae43854afe97fb9f3db96fcc2f9c2 Mon Sep 17 00:00:00 2001 From: RomanG Date: Fri, 18 Feb 2022 12:28:43 +0100 Subject: [PATCH 21/39] Rename unit test R9. --- tests/{unit_tests_R9.py => unit_tests_ray_parameters.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tests/{unit_tests_R9.py => unit_tests_ray_parameters.py} (100%) diff --git a/tests/unit_tests_R9.py b/tests/unit_tests_ray_parameters.py similarity index 100% rename from tests/unit_tests_R9.py rename to tests/unit_tests_ray_parameters.py -- GitLab From 433d0807517b9d3e4845d10cd3f024e0956cb657 Mon Sep 17 00:00:00 2001 From: RomanG Date: Fri, 18 Feb 2022 18:44:26 +0100 Subject: [PATCH 22/39] Rename unit tests to be launched by pytest --- ...sts_ephemerides.py => test_ephemerides.py} | 0 ...ion_plane.py => test_occultation_plane.py} | 53 +++++++++++++------ ...orrection.py => test_offset_correction.py} | 23 ++++---- ...y_parameters.py => test_ray_parameters.py} | 11 ++-- 4 files changed, 56 insertions(+), 31 deletions(-) rename tests/{unit_tests_ephemerides.py => test_ephemerides.py} (100%) rename tests/{unit_tests_occultation_plane.py => test_occultation_plane.py} (76%) rename tests/{unit_tests_offset_correction.py => test_offset_correction.py} (83%) rename tests/{unit_tests_ray_parameters.py => test_ray_parameters.py} (91%) diff --git a/tests/unit_tests_ephemerides.py b/tests/test_ephemerides.py similarity index 100% rename from tests/unit_tests_ephemerides.py rename to tests/test_ephemerides.py diff --git a/tests/unit_tests_occultation_plane.py b/tests/test_occultation_plane.py similarity index 76% rename from tests/unit_tests_occultation_plane.py rename to tests/test_occultation_plane.py index ac899ec..67a4f50 100644 --- a/tests/unit_tests_occultation_plane.py +++ b/tests/test_occultation_plane.py @@ -4,8 +4,11 @@ Unit tests for the functions of occultation_plane.py """ +from typing import Tuple + import numpy as np import pytest +from nptyping import NDArray import radiocc from radiocc import occultation_plane @@ -19,45 +22,51 @@ pos_SC_sample = np.array( [ [-7726575.13322412, -7666711.32894175, -7958119.5346874], [-7726451.93234389, -7667213.36854496, -7957816.58827798], - ] + ], + dtype=np.float64, ) pos_GS_sample = np.array( [ [192537483124.12, 265959242098.085, 113482080420.543], [192537459864.319, 265959250204.846, 113482084397.897], - ] + ], + dtype=np.float64, ) vel_SC_sample = np.array( [ [240.58413463, -980.54769211, 591.63687863], [240.65325798, -980.4791224, 591.70807734], - ] + ], + dtype=np.float64, ) vel_GS_sample = np.array( [ [-45427.78743875, 15832.98573322, 7768.01058944], [-45427.7801956, 15832.99253347, 7768.00905843], - ] + ], + dtype=np.float64, ) -gamma_sample = np.array([1.2397252315812128, 1.2397669744753683]) +gamma_sample = np.array([1.2397252315812128, 1.2397669744753683], dtype=np.float64) -r_SC_sample = np.array([4382958.64228272, 4382438.07915611]) -z_SC_sample = np.array([12751461.32641642, 12751678.38508523]) -z_GS_sample = np.array([-347394852353.5334, -347394846967.8551]) +r_SC_sample = np.array([4382958.64228272, 4382438.07915611], dtype=np.float64) +z_SC_sample = np.array([12751461.32641642, 12751678.38508523], dtype=np.float64) +z_GS_sample = np.array([-347394852353.5334, -347394846967.8551], dtype=np.float64) -vr_SC_sample = np.array([-1016.962309512641, -1016.9834579503646]) -vz_SC_sample = np.array([424.0824883669944, 423.9684628262237]) -vr_GS_sample = np.array([7681.6916741349305, 7679.542472103281]) -vz_GS_sample = np.array([10518.558189401045, 10518.546132036772]) +vr_SC_sample = np.array([-1016.962309512641, -1016.9834579503646], dtype=np.float64) +vz_SC_sample = np.array([424.0824883669944, 423.9684628262237], dtype=np.float64) +vr_GS_sample = np.array([7681.6916741349305, 7679.542472103281], dtype=np.float64) +vz_GS_sample = np.array([10518.558189401045, 10518.546132036772], dtype=np.float64) @pytest.mark.parametrize( "pos_SC_input, pos_GS_input, expected_gamma", [(pos_SC_sample, pos_GS_sample, gamma_sample)], ) -def test_compute_SC_angle(pos_SC_input, pos_GS_input, expected_gamma) -> None: +def test_compute_SC_angle( + pos_SC_input: NDArray, pos_GS_input: NDArray, expected_gamma: NDArray +) -> None: output = occultation_plane.compute_SC_angle(pos_SC_input, pos_GS_input) # Test if output has the right dimension @@ -87,7 +96,10 @@ def test_compute_SC_angle(pos_SC_input, pos_GS_input, expected_gamma) -> None: ], ) def test_position_in_occ_plane( - pos_SC_input, pos_GS_input, gamma_input, expected_output + pos_SC_input: NDArray, + pos_GS_input: NDArray, + gamma_input: NDArray, + expected_output: Tuple[NDArray, NDArray, NDArray], ) -> None: output = occultation_plane.position_in_occ_plane( pos_SC_input, pos_GS_input, gamma_input @@ -126,7 +138,11 @@ def test_transition_matrix_occ_plane() -> None: ], ) def test_velocity_in_occ_plane( - pos_SC_input, pos_GS_input, vel_SC_input, vel_GS_input, expected_output + pos_SC_input: NDArray, + pos_GS_input: NDArray, + vel_SC_input: NDArray, + vel_GS_input: NDArray, + expected_output: Tuple[NDArray, NDArray, NDArray, NDArray], ) -> None: output = occultation_plane.velocity_in_occ_plane( @@ -152,7 +168,12 @@ def test_velocity_in_occ_plane( ) ], ) -def test_compute_delta_s(r_SC_input, z_SC_input, z_GS_input, expected_delta_s) -> None: +def test_compute_delta_s( + r_SC_input: NDArray, + z_SC_input: NDArray, + z_GS_input: NDArray, + expected_delta_s: NDArray, +) -> None: delta_s = occultation_plane.compute_delta_s(r_SC_input, z_SC_input, z_GS_input) assert np.allclose( diff --git a/tests/unit_tests_offset_correction.py b/tests/test_offset_correction.py similarity index 83% rename from tests/unit_tests_offset_correction.py rename to tests/test_offset_correction.py index 662ee03..3b21913 100644 --- a/tests/unit_tests_offset_correction.py +++ b/tests/test_offset_correction.py @@ -8,8 +8,8 @@ from pathlib import Path import numpy as np import pytest +from nptyping import NDArray -import radiocc from radiocc import offset_correction PATH_EXPECTED_OUTPUT = Path("tests/expected_unit_tests_R7") @@ -38,7 +38,9 @@ PATH_EXPECTED_OUTPUT = Path("tests/expected_unit_tests_R7") ), # test with values from MEX case ], ) -def test_find_threshold_index(distance_input, threshold_input, expected_index) -> None: +def test_find_threshold_index( + distance_input: NDArray, threshold_input: float, expected_index: int +) -> None: assert ( offset_correction.find_threshold_index(distance_input, threshold_input) == expected_index @@ -46,7 +48,8 @@ def test_find_threshold_index(distance_input, threshold_input, expected_index) - @pytest.mark.parametrize( - "delET_input, Doppler_input, i_integral_input, i_corr_input, Type, expected_poly_coefs, expected_Doppler_biasfit", + "delET_input, Doppler_input, i_integral_input, i_corr_input, Type, " + "expected_poly_coefs, expected_Doppler_biasfit", [ ( np.loadtxt(PATH_EXPECTED_OUTPUT / "delET.txt"), @@ -60,13 +63,13 @@ def test_find_threshold_index(distance_input, threshold_input, expected_index) - ], ) def test_fit_Doppler_bias( - delET_input, - Doppler_input, - i_integral_input, - i_corr_input, - Type, - expected_poly_coefs, - expected_Doppler_biasfit, + delET_input: NDArray, + Doppler_input: NDArray, + i_integral_input: int, + i_corr_input: int, + Type: str, + expected_poly_coefs: NDArray, + expected_Doppler_biasfit: NDArray, ) -> None: output = offset_correction.fit_Doppler_bias( delET_input, Doppler_input, i_integral_input, i_corr_input, Type diff --git a/tests/unit_tests_ray_parameters.py b/tests/test_ray_parameters.py similarity index 91% rename from tests/unit_tests_ray_parameters.py rename to tests/test_ray_parameters.py index 440c721..35d6cf9 100644 --- a/tests/unit_tests_ray_parameters.py +++ b/tests/test_ray_parameters.py @@ -4,10 +4,11 @@ Unit tests for the functions of R9 """ +from typing import Any, Tuple # noqa + import numpy as np import pytest -import radiocc from radiocc import ray_parameters test_input = { @@ -64,7 +65,7 @@ test_output = { ) ], ) -def test_coeffs_equation(input, expected_output) -> None: +def test_coeffs_equation(input: Any, expected_output: Any) -> None: output = ray_parameters.coeffs_equation(*input) assert np.array_equal(np.array(output), np.array(expected_output)) @@ -96,7 +97,7 @@ def test_coeffs_equation(input, expected_output) -> None: ) ], ) -def test_equation_rhs(input, expected_output) -> None: +def test_equation_rhs(input: Any, expected_output: Any) -> None: output = ray_parameters.equation_rhs(*input) assert np.array_equal(np.array(output), np.array(expected_output)) @@ -120,7 +121,7 @@ def test_equation_rhs(input, expected_output) -> None: ) ], ) -def test_compute_angles_correction(input, expected_output) -> None: +def test_compute_angles_correction(input: Any, expected_output: Any) -> None: output = ray_parameters.compute_angles_correction(*input) assert np.array_equal(np.array(output), np.array(expected_output)) @@ -140,6 +141,6 @@ def test_compute_angles_correction(input, expected_output) -> None: ) ], ) -def test_compute_impact_parameter(input, expected_output) -> None: +def test_compute_impact_parameter(input: Any, expected_output: float) -> None: output = ray_parameters.compute_impact_parameter(*input) assert np.array_equal(np.array(output), np.array(expected_output)) -- GitLab From d74523097e0c83a3b9badb00da6e326f4524b173 Mon Sep 17 00:00:00 2001 From: RomanG Date: Fri, 18 Feb 2022 18:47:46 +0100 Subject: [PATCH 23/39] Correct flake8 and mypy errors --- radiocc/ephemerides.py | 15 +++++++++++---- radiocc/occultation_plane.py | 18 ++++++++++-------- radiocc/offset_correction.py | 6 +++--- radiocc/ray_parameters.py | 15 +++++++-------- radiocc/reconstruct.py | 7 +------ radiocc/utils.py | 4 +++- 6 files changed, 35 insertions(+), 30 deletions(-) diff --git a/radiocc/ephemerides.py b/radiocc/ephemerides.py index a72ece5..cef8188 100644 --- a/radiocc/ephemerides.py +++ b/radiocc/ephemerides.py @@ -18,7 +18,7 @@ import radiocc # noqa def get_epoch_at_target(etobs: NDArray, obs: int, direct: str, targ: int) -> NDArray: """Compute the transmit (or receive) epoch of a signal at a target specified by its - NAIF ID given the receive (or transmit) epoch at an observer specified by its NAIF ID.""" + NAIF ID given the receive (or transmit) epoch at an observer (NAIF ID).""" epoch_targ = np.full_like(etobs, np.nan) for i, et in enumerate(etobs): epoch_targ[i], _ = spy.ltime(et, obs, direct, targ) @@ -27,7 +27,7 @@ def get_epoch_at_target(etobs: NDArray, obs: int, direct: str, targ: int) -> NDA def get_state_at_target( etobs: NDArray, obs: int, targ: int, ref: str, abcorr: str -) -> NDArray: +) -> NDArray[float]: """Compute the state vector of a target body specified by its NAIF ID relative to an observing body specified by its NAIF ID at a given epoch.""" state_targ = np.array(spy.spkezr(str(targ), etobs, ref, abcorr, str(obs))[0]) * 1e3 @@ -36,7 +36,14 @@ def get_state_at_target( def get_ephemerides( ET: NDArray, ID_SC: int, ID_GS: int, ID_P: int, Ref_frame: str, abcorr: str -) -> Tuple[NDArray, NDArray, NDArray, NDArray, NDArray, NDArray]: +) -> Tuple[ + NDArray[float], + NDArray[float], + NDArray[float], + NDArray[float], + NDArray[float], + NDArray[float], +]: """Compute the state vectors of the spacecraft and the ground station relative to the Planet center at the epoch at which the uplink signal is received by the SC and the epoch at which the downlink signal is received by the Planet @@ -46,7 +53,7 @@ def get_ephemerides( # Signal emission/reception epochs # ================================ - # Epoch at which SC emits downlink signal to GS = Epoch at which SC receives uplink signal from GS + # Epoch at which SC emits downlink signal to GS tSC = get_epoch_at_target(ET, ID_GS, "<-", ID_SC) # Epoch at which the downlink signal is received by Planet from SC diff --git a/radiocc/occultation_plane.py b/radiocc/occultation_plane.py index 0b681f6..c850e5a 100644 --- a/radiocc/occultation_plane.py +++ b/radiocc/occultation_plane.py @@ -5,7 +5,6 @@ from typing import Tuple import numpy as np import spiceypy as spy from nptyping import NDArray -from pudb import set_trace as bp from radiocc.constants import POSITION_INDICES, VELOCITY_INDICES @@ -31,10 +30,10 @@ def position_in_occ_plane( return r_SC, z_SC, z_GS -def transition_matrix_occ_plane(p_SC: NDArray, p_GS: NDArray) -> NDArray: +def transition_matrix_occ_plane(p_SC: NDArray[float], p_GS: NDArray[float]) -> NDArray: """Compute the transition matrix from the equatorial geocentric basis (x,y,z) to the occultation plane basis (n, r, z).""" - n_vec = np.cross(p_GS, p_SC) # vector normal to occ. plane + n_vec: NDArray = np.cross(p_GS, p_SC) # vector normal to occ. plane n_vec /= np.linalg.norm(n_vec) r_vec = -np.cross(p_GS, n_vec) # r-axis of occ. plane (r = z x n) r_vec /= np.linalg.norm(r_vec) @@ -47,10 +46,13 @@ def transition_matrix_occ_plane(p_SC: NDArray, p_GS: NDArray) -> NDArray: def velocity_in_occ_plane( - pos_SC: NDArray, pos_GS: NDArray, vel_SC: NDArray, vel_GS: NDArray + pos_SC: NDArray[float], + pos_GS: NDArray[float], + vel_SC: NDArray[float], + vel_GS: NDArray[float], ) -> Tuple[NDArray, NDArray, NDArray, NDArray]: - """Using a transition matrix, transform the velocity vectors of the spacecraft - and the ground station relative to Mars center into the occultation coordinate system.""" + """Using a transition matrix, transform the velocity vectors of the spacecraft and + ground station relative to Mars center into the occultation coordinate system.""" vr_SC = np.full(vel_SC.shape[0], np.nan) vz_SC = np.full_like(vr_SC, np.nan) vr_GS = np.full_like(vr_SC, np.nan) @@ -107,8 +109,8 @@ def Occ_Plane_Conversion( r_SC, z_SC, z_GS = position_in_occ_plane(pos_SC, pos_GS, gamma) vr_SC, vz_SC, vr_GS, vz_GS = velocity_in_occ_plane(pos_SC, pos_GS, vel_SC, vel_GS) - state_occ_SC = np.stack((r_SC, z_SC, vr_SC, vz_SC), axis=1) - state_occ_GS = np.stack((np.zeros(N_data), z_GS, vr_GS, vz_GS), axis=1) + state_occ_SC: NDArray = np.stack((r_SC, z_SC, vr_SC, vz_SC), axis=1) + state_occ_GS: NDArray = np.stack((np.zeros(N_data), z_GS, vr_GS, vz_GS), axis=1) delta_s = compute_delta_s(r_SC, z_SC, z_GS) return state_occ_SC, state_occ_GS, gamma, delta_s diff --git a/radiocc/offset_correction.py b/radiocc/offset_correction.py index ff6fc74..1e95170 100644 --- a/radiocc/offset_correction.py +++ b/radiocc/offset_correction.py @@ -10,11 +10,11 @@ import radiocc def find_threshold_index(distance: NDArray, threshold: float) -> int: - """Find index of the altitude closest to the threshold.""" + """Find index in a sorted list of distances closest to the threshold.""" for i_corr, dist in enumerate(distance): if dist <= threshold: - break - return i_corr + return i_corr + return len(distance) - 1 def fit_Doppler_bias( diff --git a/radiocc/ray_parameters.py b/radiocc/ray_parameters.py index c471906..9ba4cf4 100644 --- a/radiocc/ray_parameters.py +++ b/radiocc/ray_parameters.py @@ -4,12 +4,11 @@ from typing import Tuple, Union import numpy as np from nptyping import NDArray -from pudb import set_trace as bp from radiocc.constants import c, ks, kx -"""Compute of the two geometrical ray parameters: bending angle and impact parameter -from the occultation measurements and ephemerides of the spacecraft and ground station.""" +"""Compute the two geometrical ray parameters (bending angle and impact parameter) from +the occultation measurements and ephemerides of the spacecraft and ground station.""" def coeffs_equation( @@ -121,18 +120,18 @@ def compute_bending_angle_impact_param( else: fs = fsup * ks # downward S-band frequency - # We follow convention from Fjeldbo 1971 for coordinates and velocities in occ. plane + # Naming convention from Fjeldbo 1971 for coordinates and velocities in occ. plane for i, (vrs, vzs, vrt, vzt, rs, zs, zt, ga, b_e, d_s) in enumerate( zip(vr_SC, vz_SC, vr_GS, vz_GS, r_SC, z_SC, z_GS, gamma, beta_e, delta_s) ): - # We solve a system of 2 equations to find the delta_r and beta_r angles (see Fjeldbo 1971) + # We solve a system of 2 equations to find the delta_r and beta_r angles # We start with rays passing outside the atmosphere (delta_r = beta_r = 0) delta_r = 0.0 beta_r = 0.0 - epsilon = 999 - stock = 0 + epsilon = 999.0 + stock = 0.0 while np.abs(epsilon) > 1e-3: # Computing coefficients of eq. A3 and A4 (Fjeldbo 1971) b11, b12, b21, b22 = coeffs_equation( @@ -156,7 +155,7 @@ def compute_bending_angle_impact_param( d_s, ) - # We solve the corrections for delta_r and beta_r from eq. A3 and A4 (Fjeldbo 1971) + # Solve correction terms of delta_r and beta_r beta_r_corr, delta_r_corr = compute_angles_correction( b11, b12, b21, b22, k1, k2 ) diff --git a/radiocc/reconstruct.py b/radiocc/reconstruct.py index cd38e19..4ace1c8 100644 --- a/radiocc/reconstruct.py +++ b/radiocc/reconstruct.py @@ -12,7 +12,7 @@ import numpy from pudb import set_trace as bp # noqa import radiocc -from radiocc import constants, ephemerides, occultation_plane, offset_correction, old +from radiocc import constants, ephemerides, occultation_plane, offset_correction from radiocc.model import ( Export, L2Data, @@ -22,7 +22,6 @@ from radiocc.model import ( Scenario, ) from radiocc.old import ( - R5_Foot_Print, R8_Plot_Doppler_1, R15_Errors, R16_Biasfit_Correction, @@ -71,7 +70,6 @@ def run( # noqa: C901 RESULTS = scenario.results(radiocc.cfg.results) DATA_PRO = str(scenario.TO_PROCESS.parent) DATA_ID_i_profile = str(scenario.TO_PROCESS.resolve()) - CODE_DIR = str(Path(old.__file__).parent.resolve()) # DATA_DIR = str(scenario.TO_PROCESS.parent.parent.resolve()) PLOT_DIR = str((RESULTS / ResultsFolders.PLOTS.name).resolve()) EPHE_DIR = str((RESULTS / ResultsFolders.EPHEMERIDS.name).resolve()) @@ -98,9 +96,6 @@ def run( # noqa: C901 T_BC_med = constants.T_BC_med T_BC_upp = constants.T_BC_upp - POSITION_INDICES = constants.POSITION_INDICES - VELOCITY_INDICES = constants.VELOCITY_INDICES - # ================================================================================= # R4: State_Vectors_twoway diff --git a/radiocc/utils.py b/radiocc/utils.py index b721a7b..98f34a6 100644 --- a/radiocc/utils.py +++ b/radiocc/utils.py @@ -166,7 +166,9 @@ def number_length(value: Union[float, NDArray], base: int) -> Union[int, NDArray return np.floor(np.log(np.abs(value)) / np.log(base)).astype(int) + 1 -def round_to_N_significant_figures(value: Union[float, NDArray], figures: int) -> float: +def round_to_N_significant_figures( + value: Union[float, NDArray], figures: int +) -> Union[float, NDArray]: n = figures - number_length(value, 10) if isinstance(value, float): return round(value * 10 ** n) / 10 ** n -- GitLab From 1aac0d2738a3ac3f8db91cf230568b7cc2503d4c Mon Sep 17 00:00:00 2001 From: RomanG Date: Fri, 18 Feb 2022 18:48:52 +0100 Subject: [PATCH 24/39] Add pytest-pudb to dependencies --- poetry.lock | 38 +++++++++++++++++++++++++++++++++++++- pyproject.toml | 1 + 2 files changed, 38 insertions(+), 1 deletion(-) diff --git a/poetry.lock b/poetry.lock index 4cd09b2..4794439 100644 --- a/poetry.lock +++ b/poetry.lock @@ -863,6 +863,21 @@ toml = "*" [package.extras] testing = ["fields", "hunter", "process-tests", "six", "pytest-xdist", "virtualenv"] +[[package]] +name = "pytest-pudb" +version = "0.7.0" +description = "Pytest PuDB debugger integration" +category = "dev" +optional = false +python-versions = "*" + +[package.dependencies] +pudb = "*" +pytest = ">=2.0" + +[package.extras] +dev = ["pexpect", "tox", "flake8"] + [[package]] name = "pytest-sugar" version = "0.9.4" @@ -1463,7 +1478,7 @@ test = [] [metadata] lock-version = "1.1" python-versions = ">=3.8,<3.11" -content-hash = "4afce8ef32df2423154cc22e29e4ddd34155d56c3edb6745d0fb003975ef048a" +content-hash = "692f58fb478315f64ed18c10ac88be0e552f081019faad9c181f9fa5f6a9a9b9" [metadata.files] alabaster = [ @@ -1702,6 +1717,7 @@ greenlet = [ {file = "greenlet-1.1.2-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:97e5306482182170ade15c4b0d8386ded995a07d7cc2ca8f27958d34d6736497"}, {file = "greenlet-1.1.2-cp310-cp310-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:e6a36bb9474218c7a5b27ae476035497a6990e21d04c279884eb10d9b290f1b1"}, {file = "greenlet-1.1.2-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:abb7a75ed8b968f3061327c433a0fbd17b729947b400747c334a9c29a9af6c58"}, + {file = "greenlet-1.1.2-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:b336501a05e13b616ef81ce329c0e09ac5ed8c732d9ba7e3e983fcc1a9e86965"}, {file = "greenlet-1.1.2-cp310-cp310-win_amd64.whl", hash = "sha256:14d4f3cd4e8b524ae9b8aa567858beed70c392fdec26dbdb0a8a418392e71708"}, {file = "greenlet-1.1.2-cp35-cp35m-macosx_10_14_x86_64.whl", hash = "sha256:17ff94e7a83aa8671a25bf5b59326ec26da379ace2ebc4411d690d80a7fbcf23"}, {file = "greenlet-1.1.2-cp35-cp35m-manylinux1_x86_64.whl", hash = "sha256:9f3cba480d3deb69f6ee2c1825060177a22c7826431458c697df88e6aeb3caee"}, @@ -1714,6 +1730,7 @@ greenlet = [ {file = "greenlet-1.1.2-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:f9d29ca8a77117315101425ec7ec2a47a22ccf59f5593378fc4077ac5b754fce"}, {file = "greenlet-1.1.2-cp36-cp36m-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:21915eb821a6b3d9d8eefdaf57d6c345b970ad722f856cd71739493ce003ad08"}, {file = "greenlet-1.1.2-cp36-cp36m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:eff9d20417ff9dcb0d25e2defc2574d10b491bf2e693b4e491914738b7908168"}, + {file = "greenlet-1.1.2-cp36-cp36m-musllinux_1_1_x86_64.whl", hash = "sha256:b8c008de9d0daba7b6666aa5bbfdc23dcd78cafc33997c9b7741ff6353bafb7f"}, {file = "greenlet-1.1.2-cp36-cp36m-win32.whl", hash = "sha256:32ca72bbc673adbcfecb935bb3fb1b74e663d10a4b241aaa2f5a75fe1d1f90aa"}, {file = "greenlet-1.1.2-cp36-cp36m-win_amd64.whl", hash = "sha256:f0214eb2a23b85528310dad848ad2ac58e735612929c8072f6093f3585fd342d"}, {file = "greenlet-1.1.2-cp37-cp37m-macosx_10_14_x86_64.whl", hash = "sha256:b92e29e58bef6d9cfd340c72b04d74c4b4e9f70c9fa7c78b674d1fec18896dc4"}, @@ -1722,6 +1739,7 @@ greenlet = [ {file = "greenlet-1.1.2-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:1e12bdc622676ce47ae9abbf455c189e442afdde8818d9da983085df6312e7a1"}, {file = "greenlet-1.1.2-cp37-cp37m-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:8c790abda465726cfb8bb08bd4ca9a5d0a7bd77c7ac1ca1b839ad823b948ea28"}, {file = "greenlet-1.1.2-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:f276df9830dba7a333544bd41070e8175762a7ac20350786b322b714b0e654f5"}, + {file = "greenlet-1.1.2-cp37-cp37m-musllinux_1_1_x86_64.whl", hash = "sha256:8c5d5b35f789a030ebb95bff352f1d27a93d81069f2adb3182d99882e095cefe"}, {file = "greenlet-1.1.2-cp37-cp37m-win32.whl", hash = "sha256:64e6175c2e53195278d7388c454e0b30997573f3f4bd63697f88d855f7a6a1fc"}, {file = "greenlet-1.1.2-cp37-cp37m-win_amd64.whl", hash = "sha256:b11548073a2213d950c3f671aa88e6f83cda6e2fb97a8b6317b1b5b33d850e06"}, {file = "greenlet-1.1.2-cp38-cp38-macosx_10_14_x86_64.whl", hash = "sha256:9633b3034d3d901f0a46b7939f8c4d64427dfba6bbc5a36b1a67364cf148a1b0"}, @@ -1730,6 +1748,7 @@ greenlet = [ {file = "greenlet-1.1.2-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:e859fcb4cbe93504ea18008d1df98dee4f7766db66c435e4882ab35cf70cac43"}, {file = "greenlet-1.1.2-cp38-cp38-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:00e44c8afdbe5467e4f7b5851be223be68adb4272f44696ee71fe46b7036a711"}, {file = "greenlet-1.1.2-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ec8c433b3ab0419100bd45b47c9c8551248a5aee30ca5e9d399a0b57ac04651b"}, + {file = "greenlet-1.1.2-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:2bde6792f313f4e918caabc46532aa64aa27a0db05d75b20edfc5c6f46479de2"}, {file = "greenlet-1.1.2-cp38-cp38-win32.whl", hash = "sha256:288c6a76705dc54fba69fbcb59904ae4ad768b4c768839b8ca5fdadec6dd8cfd"}, {file = "greenlet-1.1.2-cp38-cp38-win_amd64.whl", hash = "sha256:8d2f1fb53a421b410751887eb4ff21386d119ef9cde3797bf5e7ed49fb51a3b3"}, {file = "greenlet-1.1.2-cp39-cp39-macosx_10_14_x86_64.whl", hash = "sha256:166eac03e48784a6a6e0e5f041cfebb1ab400b394db188c48b3a84737f505b67"}, @@ -1738,6 +1757,7 @@ greenlet = [ {file = "greenlet-1.1.2-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:b1692f7d6bc45e3200844be0dba153612103db241691088626a33ff1f24a0d88"}, {file = "greenlet-1.1.2-cp39-cp39-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:7227b47e73dedaa513cdebb98469705ef0d66eb5a1250144468e9c3097d6b59b"}, {file = "greenlet-1.1.2-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:7ff61ff178250f9bb3cd89752df0f1dd0e27316a8bd1465351652b1b4a4cdfd3"}, + {file = "greenlet-1.1.2-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:0051c6f1f27cb756ffc0ffbac7d2cd48cb0362ac1736871399a739b2885134d3"}, {file = "greenlet-1.1.2-cp39-cp39-win32.whl", hash = "sha256:f70a9e237bb792c7cc7e44c531fd48f5897961701cdaa06cf22fc14965c496cf"}, {file = "greenlet-1.1.2-cp39-cp39-win_amd64.whl", hash = "sha256:013d61294b6cd8fe3242932c1c5e36e5d1db2c8afb58606c5a67efce62c1f5fd"}, {file = "greenlet-1.1.2.tar.gz", hash = "sha256:e30f5ea4ae2346e62cedde8794a56858a67b878dd79f7df76a0767e356b1744a"}, @@ -2225,6 +2245,10 @@ pytest-cov = [ {file = "pytest-cov-2.12.1.tar.gz", hash = "sha256:261ceeb8c227b726249b376b8526b600f38667ee314f910353fa318caa01f4d7"}, {file = "pytest_cov-2.12.1-py2.py3-none-any.whl", hash = "sha256:261bb9e47e65bd099c89c3edf92972865210c36813f80ede5277dceb77a4a62a"}, ] +pytest-pudb = [ + {file = "pytest-pudb-0.7.0.tar.gz", hash = "sha256:0ea87316d39c82163d340c28a168e08a163b8d3f484e60a53c9fd5eefe432c63"}, + {file = "pytest_pudb-0.7.0-py2.py3-none-any.whl", hash = "sha256:21e96fc16f313a7bd75e1df1b151de8a721144318b0ae8350208d6554222005a"}, +] pytest-sugar = [ {file = "pytest-sugar-0.9.4.tar.gz", hash = "sha256:b1b2186b0a72aada6859bea2a5764145e3aaa2c1cfbb23c3a19b5f7b697563d3"}, ] @@ -2254,18 +2278,26 @@ pyyaml = [ {file = "PyYAML-5.4.1-cp27-cp27mu-manylinux1_x86_64.whl", hash = "sha256:bb4191dfc9306777bc594117aee052446b3fa88737cd13b7188d0e7aa8162185"}, {file = "PyYAML-5.4.1-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:6c78645d400265a062508ae399b60b8c167bf003db364ecb26dcab2bda048253"}, {file = "PyYAML-5.4.1-cp36-cp36m-manylinux1_x86_64.whl", hash = "sha256:4e0583d24c881e14342eaf4ec5fbc97f934b999a6828693a99157fde912540cc"}, + {file = "PyYAML-5.4.1-cp36-cp36m-manylinux2014_aarch64.whl", hash = "sha256:72a01f726a9c7851ca9bfad6fd09ca4e090a023c00945ea05ba1638c09dc3347"}, + {file = "PyYAML-5.4.1-cp36-cp36m-manylinux2014_s390x.whl", hash = "sha256:895f61ef02e8fed38159bb70f7e100e00f471eae2bc838cd0f4ebb21e28f8541"}, {file = "PyYAML-5.4.1-cp36-cp36m-win32.whl", hash = "sha256:3bd0e463264cf257d1ffd2e40223b197271046d09dadf73a0fe82b9c1fc385a5"}, {file = "PyYAML-5.4.1-cp36-cp36m-win_amd64.whl", hash = "sha256:e4fac90784481d221a8e4b1162afa7c47ed953be40d31ab4629ae917510051df"}, {file = "PyYAML-5.4.1-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:5accb17103e43963b80e6f837831f38d314a0495500067cb25afab2e8d7a4018"}, {file = "PyYAML-5.4.1-cp37-cp37m-manylinux1_x86_64.whl", hash = "sha256:e1d4970ea66be07ae37a3c2e48b5ec63f7ba6804bdddfdbd3cfd954d25a82e63"}, + {file = "PyYAML-5.4.1-cp37-cp37m-manylinux2014_aarch64.whl", hash = "sha256:cb333c16912324fd5f769fff6bc5de372e9e7a202247b48870bc251ed40239aa"}, + {file = "PyYAML-5.4.1-cp37-cp37m-manylinux2014_s390x.whl", hash = "sha256:fe69978f3f768926cfa37b867e3843918e012cf83f680806599ddce33c2c68b0"}, {file = "PyYAML-5.4.1-cp37-cp37m-win32.whl", hash = "sha256:dd5de0646207f053eb0d6c74ae45ba98c3395a571a2891858e87df7c9b9bd51b"}, {file = "PyYAML-5.4.1-cp37-cp37m-win_amd64.whl", hash = "sha256:08682f6b72c722394747bddaf0aa62277e02557c0fd1c42cb853016a38f8dedf"}, {file = "PyYAML-5.4.1-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:d2d9808ea7b4af864f35ea216be506ecec180628aced0704e34aca0b040ffe46"}, {file = "PyYAML-5.4.1-cp38-cp38-manylinux1_x86_64.whl", hash = "sha256:8c1be557ee92a20f184922c7b6424e8ab6691788e6d86137c5d93c1a6ec1b8fb"}, + {file = "PyYAML-5.4.1-cp38-cp38-manylinux2014_aarch64.whl", hash = "sha256:fd7f6999a8070df521b6384004ef42833b9bd62cfee11a09bda1079b4b704247"}, + {file = "PyYAML-5.4.1-cp38-cp38-manylinux2014_s390x.whl", hash = "sha256:bfb51918d4ff3d77c1c856a9699f8492c612cde32fd3bcd344af9be34999bfdc"}, {file = "PyYAML-5.4.1-cp38-cp38-win32.whl", hash = "sha256:fa5ae20527d8e831e8230cbffd9f8fe952815b2b7dae6ffec25318803a7528fc"}, {file = "PyYAML-5.4.1-cp38-cp38-win_amd64.whl", hash = "sha256:0f5f5786c0e09baddcd8b4b45f20a7b5d61a7e7e99846e3c799b05c7c53fa696"}, {file = "PyYAML-5.4.1-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:294db365efa064d00b8d1ef65d8ea2c3426ac366c0c4368d930bf1c5fb497f77"}, {file = "PyYAML-5.4.1-cp39-cp39-manylinux1_x86_64.whl", hash = "sha256:74c1485f7707cf707a7aef42ef6322b8f97921bd89be2ab6317fd782c2d53183"}, + {file = "PyYAML-5.4.1-cp39-cp39-manylinux2014_aarch64.whl", hash = "sha256:d483ad4e639292c90170eb6f7783ad19490e7a8defb3e46f97dfe4bacae89122"}, + {file = "PyYAML-5.4.1-cp39-cp39-manylinux2014_s390x.whl", hash = "sha256:fdc842473cd33f45ff6bce46aea678a54e3d21f1b61a7750ce3c498eedfe25d6"}, {file = "PyYAML-5.4.1-cp39-cp39-win32.whl", hash = "sha256:49d4cdd9065b9b6e206d0595fee27a96b5dd22618e7520c33204a4a3239d5b10"}, {file = "PyYAML-5.4.1-cp39-cp39-win_amd64.whl", hash = "sha256:c20cfa2d49991c8b4147af39859b167664f2ad4561704ee74c1de03318e898db"}, {file = "PyYAML-5.4.1.tar.gz", hash = "sha256:607774cbba28732bfa802b54baa7484215f530991055bb562efbed5b2f20a45e"}, @@ -2291,6 +2323,10 @@ rfc3986 = [ {file = "ruamel.yaml-0.17.20.tar.gz", hash = "sha256:4b8a33c1efb2b443a93fcaafcfa4d2e445f8e8c29c528d9f5cdafb7cc9e4004c"}, ] "ruamel.yaml.clib" = [ + {file = "ruamel.yaml.clib-0.2.6-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:6e7be2c5bcb297f5b82fee9c665eb2eb7001d1050deaba8471842979293a80b0"}, + {file = "ruamel.yaml.clib-0.2.6-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl", hash = "sha256:221eca6f35076c6ae472a531afa1c223b9c29377e62936f61bc8e6e8bdc5f9e7"}, + {file = "ruamel.yaml.clib-0.2.6-cp310-cp310-win32.whl", hash = "sha256:1070ba9dd7f9370d0513d649420c3b362ac2d687fe78c6e888f5b12bf8bc7bee"}, + {file = "ruamel.yaml.clib-0.2.6-cp310-cp310-win_amd64.whl", hash = "sha256:77df077d32921ad46f34816a9a16e6356d8100374579bc35e15bab5d4e9377de"}, {file = "ruamel.yaml.clib-0.2.6-cp35-cp35m-macosx_10_6_intel.whl", hash = "sha256:cfdb9389d888c5b74af297e51ce357b800dd844898af9d4a547ffc143fa56751"}, {file = "ruamel.yaml.clib-0.2.6-cp35-cp35m-manylinux1_x86_64.whl", hash = "sha256:7b2927e92feb51d830f531de4ccb11b320255ee95e791022555971c466af4527"}, {file = "ruamel.yaml.clib-0.2.6-cp35-cp35m-win32.whl", hash = "sha256:ada3f400d9923a190ea8b59c8f60680c4ef8a4b0dfae134d2f2ff68429adfab5"}, diff --git a/pyproject.toml b/pyproject.toml index 464e9f4..7bb75d7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -68,6 +68,7 @@ types-PyYAML = "^5.4.10" Sphinx = "^4.4.0" pydata-sphinx-theme = "^0.8.0" sphinx-autoapi = "^1.8.4" +pytest-pudb = "^0.7.0" [build-system] requires = ["poetry-core>=1.0.0"] -- GitLab From 41dc2d26c06c5391cdda02a1e7ed7d615a25f5e2 Mon Sep 17 00:00:00 2001 From: RomanG Date: Fri, 18 Feb 2022 19:15:26 +0100 Subject: [PATCH 25/39] Add README in /tests --- tests/README.md | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 tests/README.md diff --git a/tests/README.md b/tests/README.md new file mode 100644 index 0000000..5cc8d4d --- /dev/null +++ b/tests/README.md @@ -0,0 +1,8 @@ +# tests radiocc + +## Prerequisites + +Some tests require additional external data to work. + +## Usage + -- GitLab From 8b3d95b5c615acdffccc6a9448e4c31ae4c76c9f Mon Sep 17 00:00:00 2001 From: RomanG Date: Mon, 21 Feb 2022 12:06:30 +0100 Subject: [PATCH 26/39] Run the code only one time for all tests --- tests/test_radiocc.py | 39 +++++++++------------------------------ 1 file changed, 9 insertions(+), 30 deletions(-) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index 9431f1e..aabc6db 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -4,6 +4,7 @@ radiocc """ +from calendar import c from pathlib import Path import numpy as np @@ -15,23 +16,17 @@ NUMBER_OF_TEXT_OUTPUT = 2 NUMBER_OF_FIGURE_OUTPUT = 15 NUMBER_OF_COLUMNS_OUTPUT = {"ATMO": 16, "IONO": 11} - -def test_code_is_working() -> None: - """Test if the code runs without error""" - radiocc.cfg.to_process = Path("tests/to_process") - radiocc.cfg.results = Path("tests/results") - radiocc.cfg.full_auto = True - radiocc.core.run() +# We don't put it in the fixture because it would be too +# long to rerun the code for each test. So we put it in +# a constant for now. +radiocc.cfg.to_process = Path("tests/to_process") +radiocc.cfg.results = Path("tests/results") +radiocc.cfg.full_auto = True +SCENARIOS = iter(radiocc.core.run()) +SCENARIO = next(SCENARIOS) def test_number_of_text_outputs() -> None: - # Run the code a new type - # Lire avec Pathlib le nombre de fichiers dans le dossier results/DATA - radiocc.cfg.to_process = Path("tests/to_process") - radiocc.cfg.results = Path("tests/results") - radiocc.cfg.full_auto = True - SCENARIOS = iter(radiocc.core.run()) - SCENARIO = next(SCENARIOS) TEXT_OUTPUTS = list( radiocc.utils.directories( SCENARIO.results(radiocc.cfg.results) / "DATA", "*.txt" @@ -41,11 +36,6 @@ def test_number_of_text_outputs() -> None: def test_number_of_figure_outputs() -> None: - radiocc.cfg.to_process = Path("tests/to_process") - radiocc.cfg.results = Path("tests/results") - radiocc.cfg.full_auto = True - SCENARIOS = iter(radiocc.core.run()) - SCENARIO = next(SCENARIOS) FIGURE_OUTPUTS = list( radiocc.utils.directories( SCENARIO.results(radiocc.cfg.results) / "PLOTS", ["*.png", "*.svg"] @@ -55,11 +45,6 @@ def test_number_of_figure_outputs() -> None: def test_number_of_lines_and_columns() -> None: - radiocc.cfg.to_process = Path("tests/to_process") - radiocc.cfg.results = Path("tests/results") - radiocc.cfg.full_auto = True - SCENARIOS = iter(radiocc.core.run()) - SCENARIO = next(SCENARIOS) INPUTS = radiocc.utils.directories( SCENARIO.TO_PROCESS / "DATA/LEVEL02/OPEN_LOOP/DSN/DPX", "*.TAB" ) @@ -80,12 +65,6 @@ def test_number_of_lines_and_columns() -> None: def test_output_global_consistency() -> None: - radiocc.cfg.to_process = Path("tests/to_process") - radiocc.cfg.results = Path("tests/results") - radiocc.cfg.full_auto = True - SCENARIOS = iter(radiocc.core.run()) - SCENARIO = next(SCENARIOS) - PATH_RESULTS = SCENARIO.results(radiocc.cfg.results) PATH_EXPECTED_RESULTS = Path("tests/expected_results") / SCENARIO.TO_PROCESS.name TEXT_OUTPUTS = { -- GitLab From c9d248cc59fa408dc540257aefb70a6aa9e9a280 Mon Sep 17 00:00:00 2001 From: RomanG Date: Mon, 21 Feb 2022 12:59:00 +0100 Subject: [PATCH 27/39] Add script to fetch test data. Run script with ci. --- .gitlab-ci.yml | 2 ++ tests/fetch_test_data.sh | 8 ++++++++ tests/test_offset_correction.py | 2 +- tests/test_radiocc.py | 9 +++++---- 4 files changed, 16 insertions(+), 5 deletions(-) create mode 100755 tests/fetch_test_data.sh diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 0c66cba..784dbf6 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -17,6 +17,8 @@ linter: - pip install --user -U pip poetry - poetry export --without-hashes --dev -f requirements.txt -o requirements.txt - pip install -r requirements.txt + - cd tests && ./fetch_test_data.sh + - pwd script: - pre-commit run --all-files diff --git a/tests/fetch_test_data.sh b/tests/fetch_test_data.sh new file mode 100755 index 0000000..89a60e4 --- /dev/null +++ b/tests/fetch_test_data.sh @@ -0,0 +1,8 @@ +#!/usr/bin/env bash + +# Download test data +wget https://cloud-as.oma.be/index.php/s/akpCYM4giGzkZyT/download +# Unzip archive and delete archive +unzip download +rm download +# Now, test data are located in radiocc/tests/test_data diff --git a/tests/test_offset_correction.py b/tests/test_offset_correction.py index 3b21913..4f2d1ce 100644 --- a/tests/test_offset_correction.py +++ b/tests/test_offset_correction.py @@ -12,7 +12,7 @@ from nptyping import NDArray from radiocc import offset_correction -PATH_EXPECTED_OUTPUT = Path("tests/expected_unit_tests_R7") +PATH_EXPECTED_OUTPUT = Path("tests/test_data/expected_offset_correction") @pytest.mark.parametrize( diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index aabc6db..d1a1a7d 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -4,7 +4,6 @@ radiocc """ -from calendar import c from pathlib import Path import numpy as np @@ -19,8 +18,8 @@ NUMBER_OF_COLUMNS_OUTPUT = {"ATMO": 16, "IONO": 11} # We don't put it in the fixture because it would be too # long to rerun the code for each test. So we put it in # a constant for now. -radiocc.cfg.to_process = Path("tests/to_process") -radiocc.cfg.results = Path("tests/results") +radiocc.cfg.to_process = Path("tests/test_data/to_process") +radiocc.cfg.results = Path("tests/test_data/results") radiocc.cfg.full_auto = True SCENARIOS = iter(radiocc.core.run()) SCENARIO = next(SCENARIOS) @@ -66,7 +65,9 @@ def test_number_of_lines_and_columns() -> None: def test_output_global_consistency() -> None: PATH_RESULTS = SCENARIO.results(radiocc.cfg.results) - PATH_EXPECTED_RESULTS = Path("tests/expected_results") / SCENARIO.TO_PROCESS.name + PATH_EXPECTED_RESULTS = ( + Path("tests/test_data/expected_results") / SCENARIO.TO_PROCESS.name + ) TEXT_OUTPUTS = { "ATMO": PATH_RESULTS / "DATA/X-ATMO-1.txt", "IONO": PATH_RESULTS / "DATA/X-IONO-1.txt", -- GitLab From 3b630c9b9b264a556f36aa0e8c322b53993edc0c Mon Sep 17 00:00:00 2001 From: RomanG Date: Mon, 21 Feb 2022 13:10:15 +0100 Subject: [PATCH 28/39] Remove blank line --- tests/README.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/README.md b/tests/README.md index 5cc8d4d..b1b5a1a 100644 --- a/tests/README.md +++ b/tests/README.md @@ -2,7 +2,8 @@ ## Prerequisites -Some tests require additional external data to work. - -## Usage - +Several tests require additional data which is not provided with the code to work. +The tests in `test_radiocc.py` use the MEX case `MEX-M-MRS-1-2-3-PRM-0235-V1.0` +as a test case to check the basic functions of the code and the consistency of +its output. +This particular MEX case must be put inside the `radiocc/tests/to_process` folder -- GitLab From 8420cd8fc3aead379b50f4ef2501da03671d4c0f Mon Sep 17 00:00:00 2001 From: RomanG Date: Mon, 21 Feb 2022 13:15:46 +0100 Subject: [PATCH 29/39] Edit ci command --- .gitlab-ci.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 784dbf6..0c32fd5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -17,8 +17,7 @@ linter: - pip install --user -U pip poetry - poetry export --without-hashes --dev -f requirements.txt -o requirements.txt - pip install -r requirements.txt - - cd tests && ./fetch_test_data.sh - - pwd + - cd tests && ./fetch_test_data.sh && cd .. script: - pre-commit run --all-files -- GitLab From ca65cfeb85a9437f94112cc291515369c3898294 Mon Sep 17 00:00:00 2001 From: RomanG Date: Mon, 21 Feb 2022 14:44:44 +0100 Subject: [PATCH 30/39] Fix full_auto support --- radiocc/old/R8_Plot_Doppler_1.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/radiocc/old/R8_Plot_Doppler_1.py b/radiocc/old/R8_Plot_Doppler_1.py index dd77815..58fb96c 100644 --- a/radiocc/old/R8_Plot_Doppler_1.py +++ b/radiocc/old/R8_Plot_Doppler_1.py @@ -5,7 +5,7 @@ Created on Sun Aug 30 19:52:33 2020 @author: ananya """ - +import matplotlib import matplotlib.pyplot as plt import numpy as np import scipy.stats as sc @@ -15,14 +15,14 @@ import radiocc def PLOT_Dop1(distance, Doppler, Doppler_debias, Doppler_biasfit, ET, PLOT_DIR, N_data): - import matplotlib - if radiocc.gui is not None: figure = radiocc.gui.figure figure.clear() else: - matplotlib.use("TkAgg") figure = plt.figure(6) + + if not radiocc.cfg.full_auto: + matplotlib.use("TkAgg") plt.ion() axis_left = figure.add_subplot(121) -- GitLab From 755d1edc2901fa8066112a0061728f93cce7e2f9 Mon Sep 17 00:00:00 2001 From: RomanG Date: Mon, 21 Feb 2022 14:59:25 +0100 Subject: [PATCH 31/39] Add tolerance in test of output consistency --- tests/test_ephemerides.py | 8 ++++++-- tests/test_radiocc.py | 2 +- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/tests/test_ephemerides.py b/tests/test_ephemerides.py index bc92be6..fa15f4b 100644 --- a/tests/test_ephemerides.py +++ b/tests/test_ephemerides.py @@ -19,7 +19,9 @@ LENGTH_TEST_ARRAY = 10 def test_get_epoch() -> None: - SPK_PATH = Path("tests/to_process") / TESTED_CASE / "EXTRAS/ANCILLARY/SPICE/SPK" + SPK_PATH = ( + Path("tests/test_data/to_process") / TESTED_CASE / "EXTRAS/ANCILLARY/SPICE/SPK" + ) FILES = list(radiocc.utils.directories(SPK_PATH, "*.BSP")) spy.kclear() for FILE in FILES: @@ -34,7 +36,9 @@ def test_get_epoch() -> None: def test_get_state() -> None: - SPK_PATH = Path("tests/to_process") / TESTED_CASE / "EXTRAS/ANCILLARY/SPICE/SPK" + SPK_PATH = ( + Path("tests/test_data/to_process") / TESTED_CASE / "EXTRAS/ANCILLARY/SPICE/SPK" + ) FILES = list(radiocc.utils.directories(SPK_PATH, "*.BSP")) spy.kclear() for FILE in FILES: diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index d1a1a7d..70fd902 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -85,5 +85,5 @@ def test_output_global_consistency() -> None: np.loadtxt(EXPECTED_TEXT_OUTPUT), equal_nan=True, atol=0, - rtol=0, + rtol=1e-8, ) -- GitLab From 3431b46f9264516c3dabd80dbc93b749ededc562 Mon Sep 17 00:00:00 2001 From: RomanG Date: Mon, 21 Feb 2022 17:32:05 +0100 Subject: [PATCH 32/39] Change tolerance --- tests/test_radiocc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index 70fd902..6aee244 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -85,5 +85,5 @@ def test_output_global_consistency() -> None: np.loadtxt(EXPECTED_TEXT_OUTPUT), equal_nan=True, atol=0, - rtol=1e-8, + rtol=1e-5, ) -- GitLab From 93ed6148b324352be66c16ab5df3c35654d834aa Mon Sep 17 00:00:00 2001 From: RomanG Date: Mon, 21 Feb 2022 17:45:34 +0100 Subject: [PATCH 33/39] Lower tolerance --- tests/test_radiocc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index 6aee244..ad8a328 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -85,5 +85,5 @@ def test_output_global_consistency() -> None: np.loadtxt(EXPECTED_TEXT_OUTPUT), equal_nan=True, atol=0, - rtol=1e-5, + rtol=1e-4, ) -- GitLab From ca96d6d51eceb838cb490f79329902dbc2cd5063 Mon Sep 17 00:00:00 2001 From: RomanG Date: Mon, 21 Feb 2022 17:58:04 +0100 Subject: [PATCH 34/39] Lower tolerance --- tests/test_radiocc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index ad8a328..3d06280 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -85,5 +85,5 @@ def test_output_global_consistency() -> None: np.loadtxt(EXPECTED_TEXT_OUTPUT), equal_nan=True, atol=0, - rtol=1e-4, + rtol=1e-2, ) -- GitLab From c8e0998989e45951603d69c888b1a3c58fc0a8b6 Mon Sep 17 00:00:00 2001 From: RomanG Date: Tue, 22 Feb 2022 11:56:30 +0100 Subject: [PATCH 35/39] Change tolerance to 0 --- tests/test_occultation_plane.py | 72 ++++++++++++++++++++++++--------- 1 file changed, 53 insertions(+), 19 deletions(-) diff --git a/tests/test_occultation_plane.py b/tests/test_occultation_plane.py index 67a4f50..d0a8c15 100644 --- a/tests/test_occultation_plane.py +++ b/tests/test_occultation_plane.py @@ -18,46 +18,80 @@ LENGTH_TEST_ARRAY = 2 DIMENSION_VECTOR = radiocc.constants.DIMENSION_VECTOR # Test input (first 2 lines from MEX 235) +# pos_SC_sample = np.array( +# [ +# [-7726575.13322412, -7666711.32894175, -7958119.5346874], +# [-7726451.93234389, -7667213.36854496, -7957816.58827798], +# ], +# dtype=np.float64, +# ) + pos_SC_sample = np.array( [ - [-7726575.13322412, -7666711.32894175, -7958119.5346874], - [-7726451.93234389, -7667213.36854496, -7957816.58827798], + [-7726575.133224118, -7666711.328941751, -7958119.534687403], + [-7726451.932343891, -7667213.368544964, -7957816.588277978], ], dtype=np.float64, ) + +# pos_GS_sample = np.array( +# [ +# [192537483124.12, 265959242098.085, 113482080420.543], +# [192537459864.319, 265959250204.846, 113482084397.897], +# ], +# dtype=np.float64, +# ) + pos_GS_sample = np.array( [ - [192537483124.12, 265959242098.085, 113482080420.543], - [192537459864.319, 265959250204.846, 113482084397.897], + [192537483124.12003, 265959242098.08502, 113482080420.54332], + [192537459864.31918, 265959250204.84573, 113482084397.89697], ], dtype=np.float64, ) +# vel_SC_sample = np.array( +# [ +# [240.58413463, -980.54769211, 591.63687863], +# [240.65325798, -980.4791224, 591.70807734], +# ], +# dtype=np.float64, +# ) + vel_SC_sample = np.array( [ - [240.58413463, -980.54769211, 591.63687863], - [240.65325798, -980.4791224, 591.70807734], + [240.58413463353858, -980.5476921106875, 591.63687863378], + [240.65325797954853, -980.4791223961633, 591.7080773418134], ], dtype=np.float64, ) + +# vel_GS_sample = np.array( +# [ +# [-45427.78743875, 15832.98573322, 7768.01058944], +# [-45427.7801956, 15832.99253347, 7768.00905843], +# ], +# dtype=np.float64, +# ) + vel_GS_sample = np.array( [ - [-45427.78743875, 15832.98573322, 7768.01058944], - [-45427.7801956, 15832.99253347, 7768.00905843], + [-45427.78743875351, 15832.985733222018, 7768.010589440946], + [-45427.78019560141, 15832.992533468696, 7768.009058428624], ], dtype=np.float64, ) gamma_sample = np.array([1.2397252315812128, 1.2397669744753683], dtype=np.float64) -r_SC_sample = np.array([4382958.64228272, 4382438.07915611], dtype=np.float64) -z_SC_sample = np.array([12751461.32641642, 12751678.38508523], dtype=np.float64) +r_SC_sample = np.array([4382958.642282683, 4382438.079156135], dtype=np.float64) +z_SC_sample = np.array([12751461.326416431, 12751678.385085225], dtype=np.float64) z_GS_sample = np.array([-347394852353.5334, -347394846967.8551], dtype=np.float64) -vr_SC_sample = np.array([-1016.962309512641, -1016.9834579503646], dtype=np.float64) -vz_SC_sample = np.array([424.0824883669944, 423.9684628262237], dtype=np.float64) -vr_GS_sample = np.array([7681.6916741349305, 7679.542472103281], dtype=np.float64) -vz_GS_sample = np.array([10518.558189401045, 10518.546132036772], dtype=np.float64) +vr_SC_sample = np.array([-1016.962309512641, -1016.9834579503647], dtype=np.float64) +vz_SC_sample = np.array([424.0824883669945, 423.9684628262238], dtype=np.float64) +vr_GS_sample = np.array([7681.69167413493, 7679.542472103282], dtype=np.float64) +vz_GS_sample = np.array([10518.558189401043, 10518.546132036772], dtype=np.float64) @pytest.mark.parametrize( @@ -79,7 +113,7 @@ def test_compute_SC_angle( output, expected_gamma, atol=0, - rtol=1e-8, + rtol=0, equal_nan=False, ) @@ -90,7 +124,7 @@ def test_compute_SC_angle( ( pos_SC_sample, pos_GS_sample, - round_to_N_significant_figures(gamma_sample, 15), + gamma_sample, (r_SC_sample, z_SC_sample, z_GS_sample), ) ], @@ -115,7 +149,7 @@ def test_position_in_occ_plane( np.array(output), np.array(expected_output), atol=0, - rtol=1e-8, + rtol=0, equal_nan=False, ) @@ -152,7 +186,7 @@ def test_velocity_in_occ_plane( np.array(output), np.array(expected_output), atol=0, - rtol=1e-8, + rtol=0, equal_nan=False, ) @@ -180,6 +214,6 @@ def test_compute_delta_s( delta_s, expected_delta_s, atol=0, - rtol=1e-8, + rtol=0, equal_nan=False, ) -- GitLab From d342feef91783f1191f6b584ba0c8175073025f0 Mon Sep 17 00:00:00 2001 From: RomanG Date: Tue, 22 Feb 2022 12:03:17 +0100 Subject: [PATCH 36/39] Fix unused import --- tests/test_occultation_plane.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_occultation_plane.py b/tests/test_occultation_plane.py index d0a8c15..544157d 100644 --- a/tests/test_occultation_plane.py +++ b/tests/test_occultation_plane.py @@ -12,7 +12,6 @@ from nptyping import NDArray import radiocc from radiocc import occultation_plane -from radiocc.utils import round_to_N_significant_figures LENGTH_TEST_ARRAY = 2 DIMENSION_VECTOR = radiocc.constants.DIMENSION_VECTOR -- GitLab From b1f030d88e1dd5a08f91e24a1cea037276695ac8 Mon Sep 17 00:00:00 2001 From: RomanG Date: Wed, 23 Feb 2022 15:50:44 +0100 Subject: [PATCH 37/39] Lower tolerance to 0 --- tests/test_radiocc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index 3d06280..d1a1a7d 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -85,5 +85,5 @@ def test_output_global_consistency() -> None: np.loadtxt(EXPECTED_TEXT_OUTPUT), equal_nan=True, atol=0, - rtol=1e-2, + rtol=0, ) -- GitLab From 0396b822d8ddb5b57da729533495488a5cca9f19 Mon Sep 17 00:00:00 2001 From: RomanG Date: Wed, 23 Feb 2022 16:42:49 +0100 Subject: [PATCH 38/39] Add tests/test_data to .gitignore. Put non 0 rtol --- .gitignore | 1 + tests/test_radiocc.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index ce0565e..61de3a0 100644 --- a/.gitignore +++ b/.gitignore @@ -16,5 +16,6 @@ to_process_bkp/ results/ tests/to_process/ tests/results +tests/test_data/ radiocc.log radiocc.yaml diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index d1a1a7d..3b77812 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -85,5 +85,5 @@ def test_output_global_consistency() -> None: np.loadtxt(EXPECTED_TEXT_OUTPUT), equal_nan=True, atol=0, - rtol=0, + rtol=1e-10, ) -- GitLab From a6eda6d7d4512ae486e409278c0258bc217e57fb Mon Sep 17 00:00:00 2001 From: RomanG Date: Wed, 23 Feb 2022 17:03:10 +0100 Subject: [PATCH 39/39] change rtol --- tests/test_radiocc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index 3b77812..6aee244 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -85,5 +85,5 @@ def test_output_global_consistency() -> None: np.loadtxt(EXPECTED_TEXT_OUTPUT), equal_nan=True, atol=0, - rtol=1e-10, + rtol=1e-5, ) -- GitLab