diff --git a/.gitignore b/.gitignore index 3c52565cc9bcb6ffecccb480f6cb9f4035ec4487..61de3a02e9f0b13dfaf6ff6cb792447fb6084730 100644 --- a/.gitignore +++ b/.gitignore @@ -14,5 +14,8 @@ radiocc/assets/interface/icon.ico to_process/ to_process_bkp/ results/ +tests/to_process/ +tests/results +tests/test_data/ radiocc.log radiocc.yaml diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 0c66cba796254bd91e8a48f571bd6ce1f3d20fb2..0c32fd5c18209a79a076c1952901334cdd3ef36f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -17,6 +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 && cd .. script: - pre-commit run --all-files diff --git a/README.md b/README.md index 4b4fb1bbb99c76bb7419c699b51e449b9d887d1b..776ca0929b28bf2ad9b405f680c57beaaef185d3 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/poetry.lock b/poetry.lock index 4cd09b2ee49861142dbc05f3cdbcdd0c8df145cf..479443980632ccccbc9dc0370545b722b283b78c 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 e5c962a8d7f5c8891d15c5e1cbf84ae050f12514..7bb75d72d4073db8c5c035086aceb8a512ec1554 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 "] @@ -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"] @@ -85,7 +86,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 +95,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 7d8e83ec0585e3f1630110b23a222357db935f01..cea652615a141e6ded17a1c2c0245f6598026e80 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 073d4727c83946e87f8ed1052b03fda1808b4758..7464c9144519b73c67fa531e4dbddbeab2c39c78 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/constants.py b/radiocc/constants.py index 285dee4b4e336a1d9bf472fb8503de52924511d1..011237f18a1cba1776db57fc93ce56f6c5c98ac3 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' @@ -37,3 +42,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/core.py b/radiocc/core.py index 7fa45d473f9ab8c1162d9005b69bd003344ece4d..729348200b82f9abc1178c7c32b9c003e52afd74 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/ephemerides.py b/radiocc/ephemerides.py new file mode 100644 index 0000000000000000000000000000000000000000..cef81880abd0d4890a5bc71a76c0b48c8cbb5e55 --- /dev/null +++ b/radiocc/ephemerides.py @@ -0,0 +1,87 @@ +#!/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 (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[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 + return state_targ + + +def get_ephemerides( + ET: NDArray, ID_SC: int, ID_GS: int, ID_P: int, Ref_frame: str, abcorr: str +) -> 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 + from the SC""" + + # ================================ + # Signal emission/reception epochs + # ================================ + + # 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 + 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/occultation_plane.py b/radiocc/occultation_plane.py new file mode 100644 index 0000000000000000000000000000000000000000..c850e5a0bb29ca50b57203846d49d647e0cd655e --- /dev/null +++ b/radiocc/occultation_plane.py @@ -0,0 +1,116 @@ +#!/usr/bin/env python3 + +from typing import Tuple + +import numpy as np +import spiceypy as spy +from nptyping import NDArray + +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[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: 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) + 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[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 + 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: 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 new file mode 100644 index 0000000000000000000000000000000000000000..1e95170df61dbda556f5913b7f1b9846b7e3bd3c --- /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 in a sorted list of distances closest to the threshold.""" + for i_corr, dist in enumerate(distance): + if dist <= threshold: + return i_corr + return len(distance) - 1 + + +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/old/R12_Electron_Density.py b/radiocc/old/R12_Electron_Density.py index 6fad98fb89d1798ecf15f1329e53d68dc6d43635..87b73deeadcab8ae1bc1eb3619ba4b01df6b65db 100644 --- a/radiocc/old/R12_Electron_Density.py +++ b/radiocc/old/R12_Electron_Density.py @@ -7,46 +7,42 @@ Created on Fri Nov 8 10:46:27 2019 """ import numpy as np +from nptyping import NDArray -def Elec(N_data, refractivity, e, eps0, me,Bande,fsup,c): +from radiocc.constants import ks, kx - 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: 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": + f_TX = fsup * ks # downward S-band + if Bande == "Diff": + f_TX = fsup * ks # downward S-band - Ne = np.full(N_data,np.nan)#np.array([np.nan for i in range( N_data )], float) + # Wang et al. 2018 Eq. 12 + Ne = -refractivity * 1e-6 / 40.31 * f_TX ** 2 - for i in range(N_data): + # Ne_dn[i] = -2*np.pi*(refractivity[i]*1E-6)/(2.818E-15*(c/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] = - 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): - - 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): +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): 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/R13_Neutral_Number_Density.py b/radiocc/old/R13_Neutral_Number_Density.py index 52a88fb4e459d8e5d0949f2660a1920a7c7aedce..3daf7578d67cc321239f56049e1d538eac326ec5 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 52a41b3042e19f1cb86d78cf337279cd119e1aec..a7a8fbfdfec18e9fceb295ce3137ad96e451f63c 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/R14_Plot_check.py b/radiocc/old/R14_Plot_check.py index aa9b971fb8da0a2f85712f61b7afa99481848cac..a9dd5eb315400cd82793542ae318587282f8565b 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/R17_Run_Routines1.py b/radiocc/old/R17_Run_Routines1.py index 5f3bd67d261677d6d065ac1f48e8e92fb99b719a..1516223d234dd22356e74171224923d7fb35c305 100644 --- a/radiocc/old/R17_Run_Routines1.py +++ b/radiocc/old/R17_Run_Routines1.py @@ -8,6 +8,18 @@ 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, + T_BC_low, + T_BC_med, + T_BC_upp, + Threshold_neut_upp, + c, + ks, + kx, +) from radiocc.old import ( R9_BendAng_ImpParam_dn, R9_BendAng_ImpParam_up, @@ -21,6 +33,127 @@ 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, + ) = ray_parameters.compute_bending_angle_impact_param( + state_occ_SC_dn, + state_occ_GS_dn, + gamma_dn, + delta_s_dn, + Doppler_debias_dn, + fsup, + Bande, + ) + + # 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(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 + ) + + return refractivity, Ne, TEC + + if item == "ATMO": + + # R13: neutral density (Level 04) + Nn, i_neut_upp = R13_Neutral_Number_Density.Neut( + bend_radius, refractivity, Threshold_neut_upp, i_Surface + ) + # R13a: pressure & temperature (Level 04) + P_low, T_low = R13a_Pressure_and_Temperature.Temp_Pre( + T_BC_low, Nn, bend_radius, i_neut_upp, i_Surface + ) + P_med, T_med = R13a_Pressure_and_Temperature.Temp_Pre( + T_BC_med, Nn, bend_radius, i_neut_upp, i_Surface + ) + P_upp, T_upp = R13a_Pressure_and_Temperature.Temp_Pre( + T_BC_upp, Nn, bend_radius, i_neut_upp, i_Surface + ) + + 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/old/R4_State_Vectors_twoway_v2.py b/radiocc/old/R4_State_Vectors_twoway_v2.py index 5497878f1447c4ff9fc5f7f8a522cc903e5c01f3..388566e2518cc7a2c627ae9a26f25f47cce9ffb7 100644 --- a/radiocc/old/R4_State_Vectors_twoway_v2.py +++ b/radiocc/old/R4_State_Vectors_twoway_v2.py @@ -6,203 +6,76 @@ 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_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 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]: + """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/old/R6_Framework_Conversion.py b/radiocc/old/R6_Framework_Conversion.py index ede4932b04e6967054e7600d798f6a9e0ca34944..ebf2c4e4235af2b50a5979bf58525a9d21a60dde 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/old/R7_Offset_Correction.py b/radiocc/old/R7_Offset_Correction.py index 55e98136b9b995702eae954cab92cd556a49e3d4..56509b5551b5f8d0104cb43c9295b00138026a99 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/old/R8_Plot_Doppler_1.py b/radiocc/old/R8_Plot_Doppler_1.py index 77338aeec5e51c688fca092bc113586dad4ee293..58fb96c3a4681d8153806a035fcfd3d3e601d7bb 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) @@ -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/ray_parameters.py b/radiocc/ray_parameters.py new file mode 100644 index 0000000000000000000000000000000000000000..9ba4cf4af3fa1d234e391e65744299e35c50b833 --- /dev/null +++ b/radiocc/ray_parameters.py @@ -0,0 +1,179 @@ +#!/usr/bin/env python3 + +from typing import Tuple, Union + +import numpy as np +from nptyping import NDArray + +from radiocc.constants import c, ks, kx + +"""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( + 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 + + # 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 + # We start with rays passing outside the atmosphere (delta_r = beta_r = 0) + delta_r = 0.0 + beta_r = 0.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( + 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, + ) + + # Solve correction terms of delta_r and beta_r + 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/radiocc/reconstruct.py b/radiocc/reconstruct.py index 06cb4970814b8dd3df524d3bc103e4faa5a835c6..4ace1c8418041d5467d6987384826d3867e8d984 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, occultation_plane, offset_correction from radiocc.model import ( Export, L2Data, @@ -22,10 +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, R8_Plot_Doppler_1, R15_Errors, R16_Biasfit_Correction, @@ -74,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()) @@ -100,69 +95,62 @@ 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 + # ================================================================================= - # 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 - ) - - # 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 + ) = ephemerides.get_ephemerides( + ET, + l2_data.SPACECRAFT_NAIF_CODE, + l2_data.dsn_station_naif_code, + l2_data.PLANET_NAIF_CODE, + Ref_frame, + "NONE", ) - # R7 coordinate frames + # 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 - ) + ) = 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] + vr_MEX_up = state_occ_SC_up[:, 2] + vz_MEX_up = state_occ_SC_up[:, 3] + + z_GS_up = state_occ_GS_up[:, 1] + vr_GS_up = state_occ_GS_up[:, 2] + vz_GS_up = state_occ_GS_up[:, 3] - # R7 coordinate frames + # 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 - ) + ) = 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] + 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 : @@ -191,60 +179,27 @@ def run( # noqa: C901 else: DATA_TYPE = None + ( + ET, + delET, + Doppler_debias, + Doppler_biasfit, + Corr_ind, + poly_coefs, + ) = offset_correction.apply_correction( + 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( @@ -258,7 +213,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 @@ -276,55 +236,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, + ) = offset_correction.apply_correction( + 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 @@ -347,44 +280,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": ( @@ -405,44 +307,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: @@ -461,11 +332,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/radiocc/utils.py b/radiocc/utils.py index da898fa21c232449756bb54680e5d8210d796992..98f34a6da3e618033502b91a884d04cc924f00cb 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,17 @@ 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 +) -> Union[float, NDArray]: + 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/setup.cfg b/setup.cfg index 8f6875ca9d1a8ab2b47ef97b83c2a2832a75a28c..ad39be01ddd58e5916ebcb484e010a8c7950252f 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/README.md b/tests/README.md new file mode 100644 index 0000000000000000000000000000000000000000..b1b5a1a07a333dd03aede2471d4d42859d1cc41d --- /dev/null +++ b/tests/README.md @@ -0,0 +1,9 @@ +# tests radiocc + +## Prerequisites + +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 diff --git a/tests/fetch_test_data.sh b/tests/fetch_test_data.sh new file mode 100755 index 0000000000000000000000000000000000000000..89a60e40424e4ff0b3bb3c2a8f7e37d4618f9359 --- /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/input.md b/tests/input.md new file mode 100644 index 0000000000000000000000000000000000000000..8dec92c9d28703f3ca04f6c7b888b0628d57fb92 --- /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_ephemerides.py b/tests/test_ephemerides.py new file mode 100644 index 0000000000000000000000000000000000000000..fa15f4b7b033a49803ba487ea61b8c29e850d300 --- /dev/null +++ b/tests/test_ephemerides.py @@ -0,0 +1,54 @@ +#!/usr/bin/env python3 + +""" +Unit tests for the functions of ephemerides.py +""" + +from pathlib import Path + +import numpy as np +import spiceypy as spy + +import radiocc +from radiocc import ephemerides + +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/test_data/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 = 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 + assert not np.isnan(result).all() + + +def test_get_state() -> None: + 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: + spy.furnsh(str(FILE)) + ET = np.linspace(0, 3600, LENGTH_TEST_ARRAY) + 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 + assert not np.isnan(result).all() diff --git a/tests/test_occultation_plane.py b/tests/test_occultation_plane.py new file mode 100644 index 0000000000000000000000000000000000000000..544157d34cf64d2c205dd38655c279600cc808fe --- /dev/null +++ b/tests/test_occultation_plane.py @@ -0,0 +1,218 @@ +#!/usr/bin/env python3 + +""" +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 + +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.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.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.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.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.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.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( + "pos_SC_input, pos_GS_input, expected_gamma", + [(pos_SC_sample, pos_GS_sample, gamma_sample)], +) +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 + 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=0, + equal_nan=False, + ) + + +@pytest.mark.parametrize( + "pos_SC_input, pos_GS_input, gamma_input, expected_output", + [ + ( + pos_SC_sample, + pos_GS_sample, + gamma_sample, + (r_SC_sample, z_SC_sample, z_GS_sample), + ) + ], +) +def test_position_in_occ_plane( + 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 + ) + # 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=0, + 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: 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( + 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=0, + 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: 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( + delta_s, + expected_delta_s, + atol=0, + rtol=0, + equal_nan=False, + ) diff --git a/tests/test_offset_correction.py b/tests/test_offset_correction.py new file mode 100644 index 0000000000000000000000000000000000000000..4f2d1cee92183d29af84609b489169ac3f5f4709 --- /dev/null +++ b/tests/test_offset_correction.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 + +""" +Unit tests for the functions of offset_correction.py +""" + +from pathlib import Path + +import numpy as np +import pytest +from nptyping import NDArray + +from radiocc import offset_correction + +PATH_EXPECTED_OUTPUT = Path("tests/test_data/expected_offset_correction") + + +@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: NDArray, threshold_input: float, expected_index: int +) -> 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: 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 + ) + + assert np.array_equal(output[0], expected_poly_coefs) + assert np.array_equal(output[1], expected_Doppler_biasfit) diff --git a/tests/test_radiocc.py b/tests/test_radiocc.py index 91a9d3ef675c74f69117b66aefbdf0bd3d5f0f14..6aee24481ec0e9d9eb66be0ed0ffdfca9a7aa3c4 100644 --- a/tests/test_radiocc.py +++ b/tests/test_radiocc.py @@ -4,9 +4,86 @@ radiocc """ -# from radiocc.core import start # noqa: F401 +from pathlib import Path +import numpy as np -def test() -> None: - """Test start""" - assert True +import radiocc + +# Constants for tests +NUMBER_OF_TEXT_OUTPUT = 2 +NUMBER_OF_FIGURE_OUTPUT = 15 +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/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) + + +def test_number_of_text_outputs() -> None: + TEXT_OUTPUTS = list( + radiocc.utils.directories( + SCENARIO.results(radiocc.cfg.results) / "DATA", "*.txt" + ) + ) + assert len(TEXT_OUTPUTS) == NUMBER_OF_TEXT_OUTPUT + + +def test_number_of_figure_outputs() -> None: + FIGURE_OUTPUTS = list( + radiocc.utils.directories( + SCENARIO.results(radiocc.cfg.results) / "PLOTS", ["*.png", "*.svg"] + ) + ) + assert len(FIGURE_OUTPUTS) == NUMBER_OF_FIGURE_OUTPUT + + +def test_number_of_lines_and_columns() -> None: + 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()) + + 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] + + +def test_output_global_consistency() -> None: + PATH_RESULTS = SCENARIO.results(radiocc.cfg.results) + 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", + } + 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.allclose( + np.loadtxt(TEXT_OUTPUT), + np.loadtxt(EXPECTED_TEXT_OUTPUT), + equal_nan=True, + atol=0, + rtol=1e-5, + ) diff --git a/tests/test_ray_parameters.py b/tests/test_ray_parameters.py new file mode 100644 index 0000000000000000000000000000000000000000..35d6cf976a7309398f9e00f85ed5721cb7e2c5da --- /dev/null +++ b/tests/test_ray_parameters.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python3 + +""" +Unit tests for the functions of R9 +""" + +from typing import Any, Tuple # noqa + +import numpy as np +import pytest + +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: Any, expected_output: Any) -> 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: Any, expected_output: Any) -> 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: Any, expected_output: Any) -> 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: Any, expected_output: float) -> None: + output = ray_parameters.compute_impact_parameter(*input) + assert np.array_equal(np.array(output), np.array(expected_output))