Skip to content

Commit

Permalink
use new core api
Browse files Browse the repository at this point in the history
  • Loading branch information
s-m-e committed Jan 8, 2024
1 parent 5efd002 commit a7feb50
Showing 1 changed file with 26 additions and 22 deletions.
48 changes: 26 additions & 22 deletions docs/source/examples/revisiting-lamberts-problem-in-python.myst.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,9 @@ jupytext:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.14.0
jupytext_version: 1.16.0
kernelspec:
display_name: Python 3
display_name: Python 3 (ipykernel)
language: python
name: python3
---
Expand All @@ -15,25 +15,25 @@ kernelspec:

The Izzo algorithm to solve the Lambert problem is available in hapsira and was implemented from [this paper](https://arxiv.org/abs/1403.2705).

```{code-cell}
```{code-cell} ipython3
from cycler import cycler
from matplotlib import pyplot as plt
import numpy as np
from hapsira.core import iod
from hapsira.core.iod import compute_y_vf, tof_equation_y_vf, compute_T_min_gf, find_xy_gf
from hapsira.iod import izzo
```

## Part 1: Reproducing the original figure

```{code-cell}
```{code-cell} ipython3
x = np.linspace(-1, 2, num=1000)
M_list = 0, 1, 2, 3
ll_list = 1, 0.9, 0.7, 0, -0.7, -0.9, -1
```

```{code-cell}
```{code-cell} ipython3
fig, ax = plt.subplots(figsize=(10, 8))
ax.set_prop_cycle(
cycler("linestyle", ["-", "--"])
Expand All @@ -43,8 +43,8 @@ for M in M_list:
for ll in ll_list:
T_x0 = np.zeros_like(x)
for ii in range(len(x)):
y = iod._compute_y(x[ii], ll)
T_x0[ii] = iod._tof_equation_y(x[ii], y, 0.0, ll, M)
y = compute_y_vf(x[ii], ll)
T_x0[ii] = tof_equation_y_vf(x[ii], y, 0.0, ll, M)
if M == 0 and ll == 1:
T_x0[x > 0] = np.nan
elif M > 0:
Expand Down Expand Up @@ -86,46 +86,50 @@ ax.set_ylabel("$T$")

## Part 2: Locating $T_{min}$

```{code-cell}
```{code-cell} ipython3
:tags: [nbsphinx-thumbnail]
for M in M_list:
for ll in ll_list:
x_T_min, T_min = iod._compute_T_min(ll, M, 10, 1e-8)
x_T_min, T_min = compute_T_min_gf(ll, M, 10, 1e-8)
ax.plot(x_T_min, T_min, "kx", mew=2)
fig
```

## Part 3: Try out solution

```{code-cell}
```{code-cell} ipython3
T_ref = 1
ll_ref = 0
x_ref, _ = iod._find_xy(
ll_ref, T_ref, M=0, numiter=10, lowpath=True, rtol=1e-8
x_ref, _ = find_xy_gf(
ll_ref, T_ref,
0, # M
10, # numiter
True, # lowpath
1e-8, # rtol
)
x_ref
```

```{code-cell}
```{code-cell} ipython3
ax.plot(x_ref, T_ref, "o", mew=2, mec="red", mfc="none")
fig
```

## Part 4: Run some examples

```{code-cell}
```{code-cell} ipython3
from astropy import units as u
from hapsira.bodies import Earth
```

### Single revolution

```{code-cell}
```{code-cell} ipython3
k = Earth.k
r0 = [15945.34, 0.0, 0.0] * u.km
r = [12214.83399, 10249.46731, 0.0] * u.km
Expand All @@ -138,7 +142,7 @@ v0, v = izzo.lambert(k, r0, r, tof)
v
```

```{code-cell}
```{code-cell} ipython3
k = Earth.k
r0 = [5000.0, 10000.0, 2100.0] * u.km
r = [-14600.0, 2500.0, 7000.0] * u.km
Expand All @@ -153,7 +157,7 @@ v

### Multiple revolutions

```{code-cell}
```{code-cell} ipython3
k = Earth.k
r0 = [22592.145603, -1599.915239, -19783.950506] * u.km
r = [1922.067697, 4054.157051, -8925.727465] * u.km
Expand All @@ -169,20 +173,20 @@ expected_va_r = [-2.45759553, 1.16945801, 0.43161258] * u.km / u.s
expected_vb_r = [-5.53841370, 0.01822220, 5.49641054] * u.km / u.s
```

```{code-cell}
```{code-cell} ipython3
v0, v = izzo.lambert(k, r0, r, tof, M=0)
v
```

```{code-cell}
```{code-cell} ipython3
_, v_l = izzo.lambert(k, r0, r, tof, M=1, lowpath=True)
_, v_r = izzo.lambert(k, r0, r, tof, M=1, lowpath=False)
```

```{code-cell}
```{code-cell} ipython3
v_l
```

```{code-cell}
```{code-cell} ipython3
v_r
```

0 comments on commit a7feb50

Please sign in to comment.