Potential units mismatch in `atmospheric_drag`

[Nosudrum]

🐞 Problem

Per its documentation, atmospheric_drag in poliastro.core.perturbations is computed as follows :

$$ \vec{p} = -\dfrac{1}{2}\rho v\left ( \dfrac{C_{d}A}{m} \right )\vec{v} $$

with :

• $\rho$ in kg/m3
• $v$ in km/s
• $C_{d}$ unitless
• $A/m$ in km2/kg

Mixing m (from $\rho$) and km (from $\vec{v}$) obviously cannot lead to correct units for the perturbation acceleration $\vec{p}$, unless the quantities are all associated to a unit and conversion happens automatically.

This is as I understand it not the case, as poliastro.core does not work with units anymore (per #1256 ?). Indeed, the cartesian state used by atmospheric_drag contains values in km and km/s but without using said units linked in quantity objects.

Therefore, when mixed with $\rho$ in kg/m3 following the documentation there will be a units mismatch with no automatic conversion.

You can see that behaviour for yourself with the following code :

import numpy as np
from astropy import units as u
from numpy.linalg import norm
from poliastro.bodies import Earth
from poliastro.core.perturbations import atmospheric_drag
from poliastro.earth.atmosphere import jacchia

atm_jacchia = jacchia.Jacchia77(1050 * u.K)
alt = 500 * u.km
rho = atm_jacchia.density(alt)

t0 = 0 * u.s
state = np.array([7200, 0, 0, 0, -1.03551869, 7.36809836])  # in km and km/s
k = Earth.k.to(u.km ** 3 / u.s ** 2)  # does not matter here
C_D = 1.2 * u.one
A = 20 * u.m ** 2
m = 2000 * u.kg
A_over_m = (A / m).to(u.km ** 2 / u.kg)  # conversion to km2/kg as per documentation

# Compute the drag acceleration using the poliastro function
acc_drag = atmospheric_drag(t0, state, k, C_D, A_over_m, rho)
print(f'Drag acceleration (poliastro): {acc_drag}')

# Compute the drag acceleration using the same code as the poliastro function
v_vec = state[3:]
v = norm(v_vec)
B = C_D * A_over_m

acc_drag_manual = -(1.0 / 2.0) * rho * B * v * v_vec
print(f'Drag acceleration (manual): {acc_drag_manual}')

Which yields the following output :

Drag acceleration (poliastro): [-0.00000000e+00  3.39640926e-20 -2.41667077e-19]
Drag acceleration (manual): [-0.00000000e+00  3.39640926e-20 -2.41667077e-19] km2 / m3

Same values, proving that there is no conversion in atmospheric_drag, but different units : atmospheric_drag, for a reason I do not know, only returns the raw value, without units, while the same computation performed "locally" returns the same value in km2/m3. This weird unit, which is equivalent to 1/km after conversion, proves that the km2/s2 coming from $\vec{v}$ is missing.

Now, providing poliastro.core with unitless floats (as I've understood is best practice) only hides the problem as the real issue lays in the documentation, which asks for mismatching units.

🖥 Please paste the output of following commands

conda info -a

     active environment : poliastro-main
    active env location : C:\Users\a.muller\Anaconda3\envs\poliastro-main
            shell level : 2
       user config file : C:\Users\a.muller\.condarc
 populated config files : C:\Users\a.muller\.condarc
          conda version : 4.12.0
    conda-build version : 3.21.6
         python version : 3.9.7.final.0
       virtual packages : __cuda=11.3=0
                          __win=0=0
                          __archspec=1=x86_64
       base environment : C:\Users\a.muller\Anaconda3  (writable)
      conda av data dir : C:\Users\a.muller\Anaconda3\etc\conda
  conda av metadata url : None
           channel URLs : https://conda.anaconda.org/tudat-team/win-64
                          https://conda.anaconda.org/tudat-team/noarch
                          https://repo.anaconda.com/pkgs/main/win-64
                          https://repo.anaconda.com/pkgs/main/noarch
                          https://repo.anaconda.com/pkgs/r/win-64
                          https://repo.anaconda.com/pkgs/r/noarch
                          https://repo.anaconda.com/pkgs/msys2/win-64
                          https://repo.anaconda.com/pkgs/msys2/noarch
          package cache : C:\Users\a.muller\Anaconda3\pkgs
                          C:\Users\a.muller\.conda\pkgs
                          C:\Users\a.muller\AppData\Local\conda\conda\pkgs
       envs directories : C:\Users\a.muller\Anaconda3\envs
                          C:\Users\a.muller\.conda\envs
                          C:\Users\a.muller\AppData\Local\conda\conda\envs
               platform : win-64
             user-agent : conda/4.12.0 requests/2.26.0 CPython/3.9.7 Windows/10 Windows/10.0.19041
          administrator : False
             netrc file : None
           offline mode : False

# conda environments:
#
base                     C:\Users\a.muller\Anaconda3
2022_muller_arnaud       C:\Users\a.muller\Anaconda3\envs\2022_muller_arnaud
orekit-python            C:\Users\a.muller\Anaconda3\envs\orekit-python
poliastro                C:\Users\a.muller\Anaconda3\envs\poliastro
poliastro-main        *  C:\Users\a.muller\Anaconda3\envs\poliastro-main
poliastro-test           C:\Users\a.muller\Anaconda3\envs\poliastro-test
tudat-space              C:\Users\a.muller\Anaconda3\envs\tudat-space

sys.version: 3.9.7 (default, Sep 16 2021, 16:59:28) [...
sys.prefix: C:\Users\a.muller\Anaconda3
sys.executable: C:\Users\a.muller\Anaconda3\python.exe
conda location: C:\Users\a.muller\Anaconda3\lib\site-packages\conda
conda-build: C:\Users\a.muller\Anaconda3\Scripts\conda-build.exe
conda-content-trust: C:\Users\a.muller\Anaconda3\Scripts\conda-content-trust.exe
conda-convert: C:\Users\a.muller\Anaconda3\Scripts\conda-convert.exe
conda-debug: C:\Users\a.muller\Anaconda3\Scripts\conda-debug.exe
conda-develop: C:\Users\a.muller\Anaconda3\Scripts\conda-develop.exe
conda-env: C:\Users\a.muller\Anaconda3\Scripts\conda-env.exe
conda-index: C:\Users\a.muller\Anaconda3\Scripts\conda-index.exe
conda-inspect: C:\Users\a.muller\Anaconda3\Scripts\conda-inspect.exe
conda-metapackage: C:\Users\a.muller\Anaconda3\Scripts\conda-metapackage.exe
conda-pack: C:\Users\a.muller\Anaconda3\Scripts\conda-pack.exe
conda-render: C:\Users\a.muller\Anaconda3\Scripts\conda-render.exe
conda-repo: C:\Users\a.muller\Anaconda3\Scripts\conda-repo.exe
conda-server: C:\Users\a.muller\Anaconda3\Scripts\conda-server.exe
conda-skeleton: C:\Users\a.muller\Anaconda3\Scripts\conda-skeleton.exe
conda-token: C:\Users\a.muller\Anaconda3\Scripts\conda-token.exe
conda-verify: C:\Users\a.muller\Anaconda3\Scripts\conda-verify.exe
user site dirs:

CIO_TEST: <not set>
CONDA_DEFAULT_ENV: poliastro-main
CONDA_EXE: C:\Users\a.muller\Anaconda3\Scripts\conda.exe
CONDA_PREFIX: C:\Users\a.muller\Anaconda3\envs\poliastro-main
CONDA_PREFIX_1: C:\Users\a.muller\Anaconda3
CONDA_PROMPT_MODIFIER: (poliastro-main)
CONDA_PYTHON_EXE: C:\Users\a.muller\Anaconda3\python.exe
CONDA_ROOT: C:\Users\a.muller\Anaconda3
CONDA_SHLVL: 2
CURL_CA_BUNDLE: <not set>
HOMEPATH: \
PATH: C:\Users\a.muller\Anaconda3;C:\Users\a.muller\Anaconda3\Library\mingw-w64\bin;C:\Users\a.muller\Anaconda3\Library\usr\bin;C:\Users\a.muller\Anaconda3\Library\bin;C:\Users\a.muller\Anaconda3\Scripts;C:\Users\a.muller\Anaconda3\bin;C:\Users\a.muller\Anaconda3\envs\poliastro-main;C:\Users\a.muller\Anaconda3\envs\poliastro-main\Library\mingw-w64\bin;C:\Users\a.muller\Anaconda3\envs\poliastro-main\Library\usr\bin;C:\Users\a.muller\Anaconda3\envs\poliastro-main\Library\bin;C:\Users\a.muller\Anaconda3\envs\poliastro-main\Scripts;C:\Users\a.muller\Anaconda3\envs\poliastro-main\bin;C:\Users\a.muller\Anaconda3\condabin;C:\Python365\Scripts;C:\Python365;C:\Program Files (x86)\Common Files\Oracle\Java\javapath;C:\windows\system32;C:\windows;C:\windows\System32\Wbem;C:\windows\System32\WindowsPowerShell\v1.0;C:\windows\System32\OpenSSH;C:\Users\a.muller\AppData\Local\Microsoft\WindowsApps;C:\Users\a.muller\AppData\Local\JetBrains\PyCharm 2022.1\bin;.;C:\Users\a.muller\AppData\Local\Programs\Git\cmd;C:\Users\a.muller\AppData\Local\GitHubDesktop\bin
PSMODULEPATH: D:\a.muller\Documents\WindowsPowerShell\Modules;C:\Program Files\WindowsPowerShell\Modules;C:\windows\system32\WindowsPowerShell\v1.0\Modules
REQUESTS_CA_BUNDLE: <not set>
SSL_CERT_FILE: <not set>

conda list

# packages in environment at C:\Users\a.muller\Anaconda3\envs\poliastro-main:
#
# Name                    Version                   Build  Channel
astropy                   5.0.4                    pypi_0    pypi
astroquery                0.4.6                    pypi_0    pypi
beautifulsoup4            4.11.1                   pypi_0    pypi
ca-certificates           2022.4.26            haa95532_0
certifi                   2022.5.18.1      py39haa95532_0
charset-normalizer        2.0.12                   pypi_0    pypi
cycler                    0.11.0                   pypi_0    pypi
fonttools                 4.33.3                   pypi_0    pypi
git                       2.34.1               haa95532_0
html5lib                  1.1                      pypi_0    pypi
idna                      3.3                      pypi_0    pypi
importlib-metadata        4.11.3                   pypi_0    pypi
jpeg                      9e                   h2bbff1b_0
jplephem                  2.17                     pypi_0    pypi
keyring                   23.5.0                   pypi_0    pypi
kiwisolver                1.4.2                    pypi_0    pypi
libcurl                   7.82.0               h86230a5_0
libssh2                   1.10.0               hcd4344a_0
libtiff                   4.2.0                hd0e1b90_0
llvmlite                  0.38.0                   pypi_0    pypi
lz4-c                     1.9.3                h2bbff1b_1
matplotlib                3.5.2                    pypi_0    pypi
numba                     0.55.1                   pypi_0    pypi
numpy                     1.21.6                   pypi_0    pypi
openssl                   1.1.1o               h2bbff1b_0
packaging                 21.3                     pypi_0    pypi
pandas                    1.4.2                    pypi_0    pypi
pillow                    9.1.1                    pypi_0    pypi
pip                       21.2.4           py39haa95532_0
plotly                    5.8.0                    pypi_0    pypi
poliastro                 0.17.dev0                pypi_0    pypi
proj                      8.2.1                h5ed7ab8_0
pyerfa                    2.0.0.1                  pypi_0    pypi
pyparsing                 3.0.9                    pypi_0    pypi
pyproj                    3.3.1                    pypi_0    pypi
python                    3.9.12               h6244533_0
python-dateutil           2.8.2                    pypi_0    pypi
pytz                      2022.1                   pypi_0    pypi
pyvo                      1.3                      pypi_0    pypi
pywin32-ctypes            0.2.0                    pypi_0    pypi
pyyaml                    6.0                      pypi_0    pypi
requests                  2.27.1                   pypi_0    pypi
scipy                     1.8.1                    pypi_0    pypi
setuptools                61.2.0           py39haa95532_0
six                       1.16.0             pyhd3eb1b0_1
soupsieve                 2.3.2.post1              pypi_0    pypi
sqlite                    3.38.3               h2bbff1b_0
tenacity                  8.0.1            py39haa95532_0
tzdata                    2022a                hda174b7_0
urllib3                   1.26.9                   pypi_0    pypi
vc                        14.2                 h21ff451_1
vs2015_runtime            14.27.29016          h5e58377_2
webencodings              0.5.1                    pypi_0    pypi
wheel                     0.37.1             pyhd3eb1b0_0
wincertstore              0.2              py39haa95532_2
xz                        5.2.5                h8cc25b3_1
zipp                      3.8.0                    pypi_0    pypi
zlib                      1.2.12               h8cc25b3_2
zstd                      1.4.9                h19a0ad4_0

pip freeze

astropy==5.0.4
astroquery==0.4.6
beautifulsoup4==4.11.1
certifi==2022.5.18.1
charset-normalizer==2.0.12
cycler==0.11.0
fonttools==4.33.3
html5lib==1.1
idna==3.3
importlib-metadata==4.11.3
jplephem==2.17
keyring==23.5.0
kiwisolver==1.4.2
llvmlite==0.38.0
matplotlib==3.5.2
numba==0.55.1
numpy==1.21.6
packaging==21.3
pandas==1.4.2
Pillow==9.1.1
plotly==5.8.0
poliastro @ git+https://github.com/poliastro/poliastro@c15b75ba3ee02f14b6ccdad01f1a704607439d3a
pyerfa==2.0.0.1
pyparsing==3.0.9
pyproj==3.3.1
python-dateutil==2.8.2
pytz==2022.1
pyvo==1.3
pywin32-ctypes==0.2.0
PyYAML==6.0
requests==2.27.1
scipy==1.8.1
six @ file:///tmp/build/80754af9/six_1644875935023/work
soupsieve==2.3.2.post1
tenacity @ file:///C:/ci/tenacity_1626248381338/work
urllib3==1.26.9
webencodings==0.5.1
wincertstore==0.2
zipp==3.8.0

💡 Possible solutions

Simply update the documentation to ask for parameters in matching units (= ask for $\rho$ in kg/km3).

[astrojuanlu]

atmospheric_drag, for a reason I do not know, only returns the raw value

Because this:

poliastro/src/poliastro/core/perturbations.py

Line 131 in c15b75b

@jit

Jitted functions "swallow" the units and produce always raw values, that's why we say that

This is as I understand it not the case, as poliastro.core does not work with units anymore

The idea is that the user should always use the High Level ("safe") API, or otherwise take care of unit handling themselves if using the Core API (poliastro.core).

Simply update the documentation to ask for parameters in matching units (= ask for $\rho$ in kg/km3).

Where exactly do you suggest we do this?

[Nosudrum]

Understood. It seems like I'm a bit of an edge case user as I'm mostly interacting with the Core API ^^

Where exactly do you suggest we do this?

here, change kg/m^3 to kg/km^3

poliastro/src/poliastro/core/perturbations.py

Lines 153 to 154 in c15b75b

	rho : float
	Air density at corresponding state (kg/m^3)

such that the documentation gets updated with matching units

[astrojuanlu]

Oh, I understand now! Do you want to propose the change yourself? https://github.com/poliastro/poliastro/edit/main/src/poliastro/core/perturbations.py

[astrojuanlu]

I see we were doing it correctly internally:

poliastro/tests/tests_twobody/test_perturbations.py

Lines 352 to 360 in 0bf54d3

	rho = coesa76.density((H - R) * u.km).to_value(u.kg / u.km**3)
	ax, ay, az = atmospheric_drag(
	t0,
	u_,
	k,
	C_D=C_D,
	A_over_m=A_over_m,
	rho=rho,

[Nosudrum]

Documentation change in PR #1514

[Nosudrum]

I've noticed a similar unit mismatch in the documentation of radiation_pressure. PR incoming.

[Nosudrum]

Another documentation change in PR #1515

[astrojuanlu]

Thanks a lot! Since this issue was closed in #1514, I'm closing. If you find other small annoyances like these, there's no need to open a new one, just send the PR directly.