From 9384a4bb7a636d70a8cdbac920ff45c768493252 Mon Sep 17 00:00:00 2001 From: buduam Date: Wed, 11 Jun 2025 17:40:03 +0200 Subject: [PATCH 1/6] Fixed J4 perturbation --- cascade.py/dynamics/_simple_earth.py | 160 +++++++++++++++++---------- 1 file changed, 102 insertions(+), 58 deletions(-) diff --git a/cascade.py/dynamics/_simple_earth.py b/cascade.py/dynamics/_simple_earth.py index ca123eb3..bbff31f5 100644 --- a/cascade.py/dynamics/_simple_earth.py +++ b/cascade.py/dynamics/_simple_earth.py @@ -9,16 +9,16 @@ # This little helper returns the heyoka expression for the density using # an exponential fit import typing +import numpy as np import heyoka as hy - +from cascade.dynamics import kepler def _compute_atmospheric_density(h): """ Returns the heyoka expression for the atmosheric density in kg.m^3. Input is the altitude in m. """ - import numpy as np - import heyoka as hy + # This array is produced by fitting best_x = np.array( @@ -45,6 +45,43 @@ def _compute_atmospheric_density(h): retval += alpha * hy.exp(-(h - gamma) * beta) return retval +def ECI2ECEF(days_since_J2000): + """ + This function returns the rotation matrix from ECI to ECEF at a given number of days elapsed since J2000 + Args: + - days_since_J2000 (`float`): days elapsed since J2000 + + Returns: + - `list`: Earth rotation matrix + """ + era = ( + 2 + * np.pi + * (0.7790572732640 + 1.00273781191135448 * days_since_J2000) # Earth rotation angle + ) + R = [[hy.cos(era), hy.sin(era), 0], [-hy.sin(era), hy.cos(era), 0], [0, 0, 1]] + return R + + +def _compute_density_thermonets(r, f107, f107a, ap): + """ + Returns the heyoka expression for the atmosheric density in kg.m^3, computed through ThermoNets + """ + + # days elapsed since the reference epoch J2000 (1st Jan 2000 12:00) + days_since_J2000 = hy.time / 86400 + + doy = days_since_J2000 + + + xyz_ecef = np.matmul(ECI2ECEF(days_since_J2000), r) # matrix multiplication r_ECEF = R(ECI2ECEF) @ r_ECI + + h, lat, lon = hy.model.cart2geo([xyz_ecef[0], xyz_ecef[1], xyz_ecef[2]]) # compute geodetic cooordinates [h is in meters] + + density_nn = hy.model.nrlmsise00_tn(geodetic=[h / 1000, lat, lon], f107=f107, f107a=f107a, ap=ap, time_expr=doy) # compute density [h must be in kilometers] + + return density_nn + def simple_earth( J2: bool = True, @@ -55,6 +92,7 @@ def simple_earth( moon: bool = False, SRP: bool = False, drag: bool = True, + thermonets: bool = False, ) -> typing.List[typing.Tuple[hy.expression, hy.expression]]: """Perturbed dynamics around the Earth. @@ -92,33 +130,23 @@ def simple_earth( Returns: The dynamics in SI units. Can be used to instantiate a :class:`~cascade.sim`. """ - from cascade.dynamics import kepler - import heyoka as hy - import numpy as np + # constants (final underscore reminds us its not SI) GMe_ = 3.986004407799724e5 # [km^3/sec^2] GMo_ = 1.32712440018e11 # [km^3/sec^2] GMm_ = 4.9028e3 # [km^3/sec^2] Re_ = 6378.1363 # [km] - C20 = ( - -4.84165371736e-4 - ) # (J2 in m^5/s^2 is 1.75553E25, C20 is the Stokes coefficient) - C22 = 2.43914352398e-6 - S22 = -1.40016683654e-6 - J3_dim_value = -2.61913e29 # (m^6/s^2) is # name is to differentiate from kwarg - J4_adim_value = -1.61989759991697e-6 # [-] - theta_g = ( - np.pi / 180 - ) * 280.4606 # [rad] # This value defines the rotation of the Earth fixed system at t0 - nu_e = (np.pi / 180) * ( - 4.178074622024230e-3 - ) # [rad/sec] # This value represents the Earth spin angular velocity. + + + + theta_g = (np.pi / 180) * 280.4606 # [rad] # This value defines the rotation of the Earth fixed system at t0 + nu_e = (np.pi / 180) * (4.178074622024230e-3) # [rad/sec] # This value represents the Earth spin angular velocity. nu_o = (np.pi / 180) * (1.1407410259335311e-5) # [rad/sec] nu_ma = (np.pi / 180) * (1.512151961904581e-4) # [rad/sec] nu_mp = (np.pi / 180) * (1.2893925235125941e-6) # [rad/sec] nu_ms = (np.pi / 180) * (6.128913003523574e-7) # [rad/sec] - alpha_o_ = 1.49619e8 # [km] + alpha_o_ = 1.49619e8 # 1 AU in [km] epsilon = (np.pi / 180) * 23.4392911 # [rad] phi_o = (np.pi / 180) * 357.5256 # [rad] Omega_plus_w = (np.pi / 180) * 282.94 # [rad] @@ -128,64 +156,67 @@ def simple_earth( x, y, z, vx, vy, vz = hy.make_vars("x", "y", "z", "vx", "vy", "vz") # Create Keplerian dynamics in SI units. + Re_SI = Re_ * 1000 GMe_SI = GMe_ * 1e9 dyn = kepler(mu=GMe_SI) # Define the radius squared - magr2 = x**2 + y**2 + z**2 + r2 = x**2 + y**2 + z**2 + # Define the radius + r = hy.sqrt(r2) - Re_SI = Re_ * 1000 + if J2: - J2term1 = GMe_SI * (Re_SI**2) * np.sqrt(5) * C20 / (2 * magr2 ** (1 / 2)) - J2term2 = 3 / (magr2**2) - J2term3 = 15 * (z**2) / (magr2**3) - fJ2x = J2term1 * x * (J2term2 - J2term3) - fJ2y = J2term1 * y * (J2term2 - J2term3) - fJ2z = J2term1 * z * (3 * J2term2 - J2term3) + J2_adim = -0.1082635854 * 1e-2 # from 'https://en.wikipedia.org/wiki/Geopotential_spherical_harmonic_model' + J2 = - J2_adim * GMe_SI * Re_SI**2 + u_J2 = J2/(2*r**5)*(3*z**2 - r2) + fJ2x = - hy.diff(u_J2,x) + fJ2y = - hy.diff(u_J2,y) + fJ2z = - hy.diff(u_J2,z) dyn[3] = (dyn[3][0], dyn[3][1] + fJ2x) dyn[4] = (dyn[4][0], dyn[4][1] + fJ2y) dyn[5] = (dyn[5][0], dyn[5][1] + fJ2z) if J3: - magr9 = magr2 ** (1 / 2) * magr2**4 # r**9 - fJ3x = J3_dim_value * x * y / magr9 * (10 * z**2 - 15 / 2 * (x**2 + y**2)) - fJ3y = J3_dim_value * z * y / magr9 * (10 * z**2 - 15 / 2 * (x**2 + y**2)) - fJ3z = ( - J3_dim_value - * 1 - / magr9 - * ( - 4 * z**2 * (z**2 - 3 * (x**2 + y**2)) - + 3 / 2 * (x**2 + y**2) ** 2 - ) - ) + J3_adim = 0.2532435346 * 1e-5 # from 'https://en.wikipedia.org/wiki/Geopotential_spherical_harmonic_model' + J3 = - J3_adim * GMe_SI * Re_SI**3 + u_J3 = J3*z/(2*r**7)*(5*z**2-3*r2) + fJ3x = - hy.diff(u_J3,x) + fJ3y = - hy.diff(u_J3,y) + fJ3z = - hy.diff(u_J3,z) dyn[3] = (dyn[3][0], dyn[3][1] + fJ3x) dyn[4] = (dyn[4][0], dyn[4][1] + fJ3y) dyn[5] = (dyn[5][0], dyn[5][1] + fJ3z) if J4: - fJ4x = ((15*J4_adim_value*GMe_SI*(Re_SI**4)*x)/(8*((x**2 + y**2 + z**2)**3.5))) * (1 - ((14*(z**2))/((x**2 + y**2 + z**2))) + (21*(z**4)/((x**2 + y**2 + z**2)**2))) - fJ4y = ((15*J4_adim_value*GMe_SI*(Re_SI**4)*y)/(8*((x**2 + y**2 + z**2)**3.5))) * (1 - ((14*(z**2))/((x**2 + y**2 + z**2))) + (21*(z**4)/((x**2 + y**2 + z**2)**2))) - fJ4z = ((15*J4_adim_value*GMe_SI*(Re_SI**4)*z)/(8*((x**2 + y**2 + z**2)**3.5))) * (5 - ((70*(z**2))/(3*(x**2 + y**2 + z**2))) + (21*(z**4)/((x**2 + y**2 + z**2)**2))) + J4_adim = 0.1619331205*1e-5 # from 'https://en.wikipedia.org/wiki/Geopotential_spherical_harmonic_model' + J4 = - J4_adim * GMe_SI * Re_SI**4 + u_J4 = J4/8*(35*z**4 - 30*r2*z**2 + 3*r2**2)/r**9 + fJ4x = - hy.diff(u_J4,x) + fJ4y = - hy.diff(u_J4,y) + fJ4z = - hy.diff(u_J4,z) dyn[3] = (dyn[3][0], dyn[3][1] + fJ4x) dyn[4] = (dyn[4][0], dyn[4][1] + fJ4y) dyn[5] = (dyn[5][0], dyn[5][1] + fJ4z) if C22S22: + C22 = 2.43914352398e-6 + S22 = -1.40016683654e-6 + X = x * hy.cos(theta_g + nu_e * hy.time) + y * hy.sin(theta_g + nu_e * hy.time) Y = -x * hy.sin(theta_g + nu_e * hy.time) + y * hy.cos(theta_g + nu_e * hy.time) Z = z C22term1 = ( - 5 * GMe_SI * (Re_SI**2) * np.sqrt(15) * C22 / (2 * magr2 ** (7 / 2)) + 5 * GMe_SI * (Re_SI**2) * np.sqrt(15) * C22 / (2 * r2 ** (7 / 2)) ) - C22term2 = GMe_SI * (Re_SI**2) * np.sqrt(15) * C22 / (magr2 ** (5 / 2)) + C22term2 = GMe_SI * (Re_SI**2) * np.sqrt(15) * C22 / (r2 ** (5 / 2)) fC22X = C22term1 * X * (Y**2 - X**2) + C22term2 * X fC22Y = C22term1 * Y * (Y**2 - X**2) - C22term2 * Y fC22Z = C22term1 * Z * (Y**2 - X**2) - S22term1 = 5 * GMe_SI * (Re_SI**2) * np.sqrt(15) * S22 / (magr2 ** (7.0 / 2)) - S22term2 = GMe_SI * (Re_SI**2) * np.sqrt(15) * S22 / (magr2 ** (5.0 / 2)) + S22term1 = 5 * GMe_SI * (Re_SI**2) * np.sqrt(15) * S22 / (r2 ** (7.0 / 2)) + S22term2 = GMe_SI * (Re_SI**2) * np.sqrt(15) * S22 / (r2 ** (5.0 / 2)) fS22X = -S22term1 * (X**2) * Y + S22term2 * Y fS22Y = -S22term1 * X * (Y**2) + S22term2 * X fS22Z = -S22term1 * X * Y * Z @@ -331,16 +362,23 @@ def simple_earth( dyn[4] = (dyn[4][0], dyn[4][1] + fMoonY) dyn[5] = (dyn[5][0], dyn[5][1] + fMoonZ) - drag_par_idx = 0 if drag: # Adds the drag force. + BSTAR = hy.par[0] magv2 = vx**2 + vy**2 + vz**2 magv = hy.sqrt(magv2) - # Here we consider a spherical Earth ... would be easy to account for the oblateness effect. - altitude = hy.sqrt(magr2) - Re_SI - density = _compute_atmospheric_density(altitude) + + if thermonets: + f107 = hy.par[1] + f107a = hy.par[2] + ap = hy.par[3] + density = _compute_density_thermonets(r = [x, y, z], f107 = f107, f107a = f107a, ap = ap) # time must be seconds elapsed since J2000 + else: + altitude = r - Re_SI # Here we consider a spherical Earth ... would be easy to account for the oblateness effect. + density = _compute_atmospheric_density(altitude) + ref_density = 0.1570 / Re_SI - fdrag = density / ref_density * hy.par[drag_par_idx] * magv + fdrag = density / ref_density * BSTAR * magv fdragx = -fdrag * vx fdragy = -fdrag * vy fdragz = -fdrag * vz @@ -348,14 +386,20 @@ def simple_earth( dyn[4] = (dyn[4][0], dyn[4][1] + fdragy) dyn[5] = (dyn[5][0], dyn[5][1] + fdragz) - srp_par_idx = 0 + if SRP: + if drag and thermonets: + k = hy.par[4] # c_r * A / m where c_r is the reflectivity, A is the area facing the Sun and m is the S/C mass (D.A. Vallado - Fundamentals of Astrodynamics and Applications - 4th Edition - Section 8.6.4) + elif drag and not thermonets: + k = hy.par[1] + elif not drag: + k = hy.par[0] + PSRP_SI = PSRP_ / 1000.0 # [kg/(m*sec^2)] - alpha_o_SI = alpha_o_ * 1000.0 # [m] - if drag: - srp_par_idx = 1 + alpha_o_SI = alpha_o_ * 1000.0 # 1 AU in [m] + SRPterm = ( - hy.par[srp_par_idx] * PSRP_SI * (alpha_o_SI**2) / (magRRo2 ** (3.0 / 2)) + k * PSRP_SI * (alpha_o_SI**2) / (magRRo2 ** (3.0 / 2)) ) fSRPX = SRPterm * (x - Xo) fSRPY = SRPterm * (y - Yo) @@ -364,4 +408,4 @@ def simple_earth( dyn[4] = (dyn[4][0], dyn[4][1] + fSRPY) dyn[5] = (dyn[5][0], dyn[5][1] + fSRPZ) - return dyn + return dyn \ No newline at end of file From 802c1e76685dbd405bd929d9a647e4dfff68d811 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Thu, 12 Jun 2025 09:48:01 +0200 Subject: [PATCH 2/6] Update gha_ci.yml update osx CI machine --- .github/workflows/gha_ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/gha_ci.yml b/.github/workflows/gha_ci.yml index 0f65b609..2d1bce83 100644 --- a/.github/workflows/gha_ci.yml +++ b/.github/workflows/gha_ci.yml @@ -97,7 +97,7 @@ jobs: cd c:\ python -c "from cascade import test; test.run_test_suite()" osx_11: - runs-on: macos-11 + runs-on: macos-15 steps: - uses: actions/checkout@v4 - name: Build From 49d03b2bad8fa0133a54d7502eadec49daced474 Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Thu, 12 Jun 2025 09:51:21 +0200 Subject: [PATCH 3/6] Update gha_conda_osx.sh --- tools/gha_conda_osx.sh | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/tools/gha_conda_osx.sh b/tools/gha_conda_osx.sh index f42006a7..bfbc3326 100644 --- a/tools/gha_conda_osx.sh +++ b/tools/gha_conda_osx.sh @@ -6,16 +6,13 @@ set -x # Exit on error. set -e -# Core deps. -sudo apt-get install wget # Install conda+deps. -wget https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh -O miniforge.sh +wget https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-MacOSX-x86_64.sh -O miniforge3.sh export deps_dir=$HOME/local -export PATH="$HOME/miniforge/bin:$PATH" -bash miniforge.sh -b -p $HOME/miniforge -conda env create -f cascade_devel.yml -q -p $deps_dir -source activate $deps_dir +export PATH="$HOME/miniforge3/bin:$PATH" +bash miniforge3.sh -b -p $HOME/miniforge3 +mamba env create --file=kep3_devel.yml -q -p $deps_dir source activate $deps_dir mkdir build From 88ddc4d1660dfb2ee18b84d5c6b52462135aeb8c Mon Sep 17 00:00:00 2001 From: Dario Izzo Date: Thu, 12 Jun 2025 09:53:45 +0200 Subject: [PATCH 4/6] Update gha_conda_osx.sh typo --- tools/gha_conda_osx.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/gha_conda_osx.sh b/tools/gha_conda_osx.sh index bfbc3326..eeeac13f 100644 --- a/tools/gha_conda_osx.sh +++ b/tools/gha_conda_osx.sh @@ -12,7 +12,7 @@ wget https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge export deps_dir=$HOME/local export PATH="$HOME/miniforge3/bin:$PATH" bash miniforge3.sh -b -p $HOME/miniforge3 -mamba env create --file=kep3_devel.yml -q -p $deps_dir +mamba env create --file=cascade_devel.yml -q -p $deps_dir source activate $deps_dir mkdir build From 406764a7539f8e884b616d760af928afeeeaed31 Mon Sep 17 00:00:00 2001 From: buduam Date: Thu, 12 Jun 2025 11:36:52 +0200 Subject: [PATCH 5/6] Added dynamics tests --- cascade.py/dynamics/__init__.py | 2 +- cascade.py/test.py | 65 +++++++++++++++++++++++++-------- 2 files changed, 50 insertions(+), 17 deletions(-) diff --git a/cascade.py/dynamics/__init__.py b/cascade.py/dynamics/__init__.py index 76e9f73b..37855fa5 100644 --- a/cascade.py/dynamics/__init__.py +++ b/cascade.py/dynamics/__init__.py @@ -26,5 +26,5 @@ from ..core import _kepler as kepler # Import the pure python symbols -from ._simple_earth import simple_earth +from ._simple_earth import simple_earth, _compute_density_thermonets from ._mascon_asteroid import mascon_asteroid, mascon_asteroid_energy diff --git a/cascade.py/test.py b/cascade.py/test.py index 1198fc5a..07ceabbc 100644 --- a/cascade.py/test.py +++ b/cascade.py/test.py @@ -13,39 +13,37 @@ class dynamics_test_case(_ut.TestCase): def runTest(self): self.test_simple_earth_api() self.test_kepler_equivalence() + self.test_perturbation_magnitudes() def test_simple_earth_api(self): from .dynamics import simple_earth simple_earth( - J2=True, - J3=False, - C22S22=False, - sun=False, - moon=False, - SRP=False, - drag=False, + J2=True, J3=False, J4=False, C22S22=False, sun=False, moon=False, SRP=False, drag=False, thermonets=False + ) + simple_earth( + J2=True, J3=True, J4=False, C22S22=False, sun=False, moon=False, SRP=False, drag=False, thermonets=False ) simple_earth( - J2=True, J3=True, C22S22=False, sun=False, moon=False, SRP=False, drag=False + J2=True, J3=True, J4=True, C22S22=False, sun=False, moon=False, SRP=False, drag=False, thermonets=False ) simple_earth( - J2=True, J3=True, C22S22=True, sun=False, moon=False, SRP=False, drag=False + J2=True, J3=True, J4=True, C22S22=True, sun=False, moon=False, SRP=False, drag=False, thermonets=False ) simple_earth( - J2=True, J3=True, C22S22=True, sun=True, moon=False, SRP=False, drag=False + J2=True, J3=True, J4=True, C22S22=True, sun=True, moon=False, SRP=False, drag=False, thermonets=False ) simple_earth( - J2=True, J3=True, C22S22=True, sun=True, moon=True, SRP=False, drag=False + J2=True, J3=True, J4=True, C22S22=True, sun=True, moon=True, SRP=False, drag=False, thermonets=False ) simple_earth( - J2=True, J3=True, C22S22=True, sun=True, moon=True, SRP=True, drag=False + J2=True, J3=True, J4=True, C22S22=True, sun=True, moon=True, SRP=True, drag=False, thermonets=False ) simple_earth( - J2=True, J3=True, C22S22=True, sun=True, moon=True, SRP=False, drag=True + J2=True, J3=True, J4=True, C22S22=True, sun=True, moon=True, SRP=True, drag=True, thermonets=False ) simple_earth( - J2=True, J3=True, C22S22=True, sun=True, moon=True, SRP=True, drag=True + J2=True, J3=True, J4=True, C22S22=True, sun=True, moon=True, SRP=True, drag=True, thermonets=True ) def test_kepler_equivalence(self): @@ -58,57 +56,92 @@ def test_kepler_equivalence(self): self.assertEqual(dyn1, dyn2) def test_perturbation_magnitudes(self): - from .dynamics import simple_earth, kepler + from .dynamics import simple_earth, kepler, _compute_density_thermonets import numpy as np from heyoka import cfunc, make_vars dynkep = simple_earth( J2=False, J3=False, + J4=False, C22S22=False, sun=False, moon=False, SRP=False, drag=False, + thermonets=False ) dynJ2 = simple_earth( J2=True, J3=False, + J4=False, C22S22=False, sun=False, moon=False, SRP=False, drag=False, + thermonets=False ) dynJ3 = simple_earth( J2=False, J3=True, + J4=False, C22S22=False, sun=False, moon=False, SRP=False, drag=False, + thermonets=False + ) + dynJ4 = simple_earth( + J2=False, + J3=False, + J4=True, + C22S22=False, + sun=False, + moon=False, + SRP=False, + drag=False, + thermonets=False ) # Dynamical variables. - x, y, z, vx, vy, vz = make_vars("x", "y", "z", "vx", "vy", "vz") + x, y, z, f107, f107a, ap = make_vars("x", "y", "z", "f107","f107a","ap") dynkep_c = cfunc([dynkep[i][1] for i in [3, 4, 5]], vars=[x,y,z]) dynJ2_c = cfunc([dynJ2[i][1] for i in [3, 4, 5]], vars=[x,y,z]) dynJ3_c = cfunc([dynJ3[i][1] for i in [3, 4, 5]], vars=[x,y,z]) + dynJ4_c = cfunc([dynJ4[i][1] for i in [3, 4, 5]], vars=[x,y,z]) # We compute the various acceleration magnitudes at 7000 km pos = np.array([7000000.0, 0.0, 0.0]) acckep = dynkep_c(pos) accJ2 = dynJ2_c(pos) - dynkep_c(pos) accJ3 = dynJ3_c(pos) - dynkep_c(pos) + accJ4 = dynJ4_c(pos) - dynkep_c(pos) # And check magnitudes self.assertTrue(np.linalg.norm(acckep) > 8) self.assertTrue(np.linalg.norm(accJ2) > 0.01) self.assertTrue(np.linalg.norm(accJ3) > 0.00001) + self.assertTrue(np.linalg.norm(accJ4) > 0.00001) self.assertTrue(np.linalg.norm(acckep) < 10) self.assertTrue(np.linalg.norm(accJ2) < 0.1) self.assertTrue(np.linalg.norm(accJ3) < 0.0001) + self.assertTrue(np.linalg.norm(accJ4) < 0.0001) + + # Check atmospheric density from ThermoNets + # Space weather indices at J2000 (1st January 2000, 12:00:00 TT) + f107 = 129.9 + f107a = 166.2 + ap = 30 + + density_func = _compute_density_thermonets(r=[x,y,z],f107=f107,f107a=f107a,ap=ap) + density_c = cfunc([density_func], vars=[x,y,z,f107,f107a,ap]) + density = density_c([pos[0],pos[1],pos[2],f107,f107a,ap],time=0.) + + # Check magnitude + self.assertTrue(density > 5e-14) + self.assertTrue(density < 1e-13) class sim_test_case(_ut.TestCase): From e8197379853ea6ac38d4f7d2810fc8a85c3cd349 Mon Sep 17 00:00:00 2001 From: buduam Date: Thu, 12 Jun 2025 13:56:24 +0200 Subject: [PATCH 6/6] modified test --- cascade.py/test.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cascade.py/test.py b/cascade.py/test.py index 07ceabbc..c8eca43e 100644 --- a/cascade.py/test.py +++ b/cascade.py/test.py @@ -131,13 +131,13 @@ def test_perturbation_magnitudes(self): # Check atmospheric density from ThermoNets # Space weather indices at J2000 (1st January 2000, 12:00:00 TT) - f107 = 129.9 - f107a = 166.2 - ap = 30 + f107_val = 129.9 + f107a_val = 166.2 + ap_val = 30 density_func = _compute_density_thermonets(r=[x,y,z],f107=f107,f107a=f107a,ap=ap) density_c = cfunc([density_func], vars=[x,y,z,f107,f107a,ap]) - density = density_c([pos[0],pos[1],pos[2],f107,f107a,ap],time=0.) + density = density_c([pos[0],pos[1],pos[2],f107_val,f107a_val,ap_val],time=0.) # Check magnitude self.assertTrue(density > 5e-14)