```python
import sympy as sp
from sympy import Rational
from sympy import pprint

```


```python
Mc, Mm, Rp, rh_c, rh_m, Vc, Vm, Mp = sp.symbols('Mc, Mm, Rp, rh_c, rh_m, Vc, Vm, Mp')
```


```python
MOI_c = Rational(2,5) * Mc * (Rational(1,2) * Rp) ** 2
MOI_m = Rational(2,5) * Mm * (Rp **2 - (Rational(1,2) * Rp) ** 2)
```


```python
MOI = MOI_c + MOI_m

pprint(MOI) # Verified
```

         2          2
    Mc⋅Rp    3⋅Mm⋅Rp 
    ────── + ────────
      10        10   



```python
# Substitute mass by rho*V
MOI = MOI.subs([
    (Mm, rh_m * Vm),
    (Mc, rh_c * Vc)
])

# Mantle volume according to V_sphere = (4/3) pi * r ** 3
Vm_val = Rational(4,3) * sp.pi * ( Rp ** 3 - (Rational(1,2) * Rp) ** 3)

# Core volume according to V_sphere = (4/3) pi * r ** 3
Vc_val = Rational(4,3) * sp.pi * (Rational(1,2) * Rp) ** 3

# Substitute the volumes by their symbolic expressions.
MOI = MOI.subs([
    (Vm, Vm_val),
    (Vc, Vc_val)
])

# Substitute rho_c by 4* rho_m for MOI in terms of rho_m
MOI_rhm = MOI.subs([
    (rh_c, 4 * rh_m)
])

# MOI in terms of rh_c
MOI_rhc = MOI.subs([(rh_m, Rational(1,4) * rh_c)])

# Expression for MOI
pprint(MOI_rhm)
pprint(MOI_rhc)
```

        5    
    π⋅Rp ⋅rhₘ
    ─────────
        15   
        5     
    π⋅Rp ⋅rh_c
    ──────────
        60    



```python
# Summation of both mass contributions.
M = Mm + Mc

# Replace Mass by rho*V
M = M.subs([
    (Mm, rh_m * Vm),
    (Mc, rh_c * Vc)
])

# Replace volume by expressions.
M = M.subs([
    (Vm, Vm_val),
    (Vc, Vc_val)
])

# Put Mass in terms of only rho_c
M = M.subs([(rh_m, Rational(1,4) * rh_c)])

# Put rho_c in temrs of total mass
rh_c_sub = sp.solve(sp.Eq(Mp, M), rh_c)[0]

# Replace MOI expression with rho_c by rho_c found by mass.
Solved = MOI_rhc.subs([
    (rh_c, rh_c_sub)
])

MOI_p = sp.simplify(Solved)

pprint(MOI_p / (Rational(2,5) * Mp * (Rp**2)))
```

    1/11


## Problem 5


```python
import numpy as np
```


```python
mass = np.array([
    59736,
    893.3,
    0.65,
    214.7
]) * 10 ** 20

radius_surface = np.array([
    6371.0,
    1821.3,
    248.0,
    1352.6
]) * 10 ** 3

G = 6.674 * 10 ** (-11)
```


```python
# print(np)
# V_esc = np.sqrt(np.divide(G*mass,radius_surface))
x = 2* G*mass/ radius_surface
# print(x)
import sympy as sp
vesc = sp.Matrix(np.sqrt(x.astype(float)))


```


```python
g = G * mass / radius_surface / radius_surface

print(g)
```

    [9.822162851846851 1.7972999454236644 0.07053362382934444
     0.7832122663239973]



```python
vesc
```




    Matrix([
    [11187.2248148606],
    [2558.68028116063],
    [187.041913536391],
    [1455.59122794131]])




```python
vesc

```




    Matrix([
    [11187.2248148606],
    [2558.68028116063],
    [187.041913536391],
    [1455.59122794131]])




```python
def v0(h, g):
    inc  = 2 * h * g
    inc = inc.astype(float)
    return np.sqrt(inc)
```


```python
h = np.array([60, 65000, 10E12, 10000])
```


```python
from sympy import pprint
v0 = v0(h, g)
```


```python
vesc = np.array(vesc)
```


```python
vesc = vesc.flatten()
```


```python
v0/vesc
print(v0)
```

    [3.43316114e+01 4.83372520e+02 1.18771734e+06 1.25156883e+02]



```python
alpha = (2 * h / radius_surface + (h**2/radius_surface/radius_surface))
```


```python
print(alpha)
```

    [1.88354364e-05 7.26512774e-02 1.62591059e+15 1.48409964e-02]



```python
np.set_printoptions(precision=4)
```


```python
print(alpha)
```

    [1.8835e-05 7.2651e-02 1.6259e+15 1.4841e-02]



```python

```

d3f5164e5a93be7dd1c7396a680d0361ddecbcd9
ipynb
Jupyter Notebook
planetary_sciences_assignment_2.ipynb
MattTurnock/PlanetarySciencesMatt
81954d2182d9577bd7327a98c45963ae42968df4
MIT
planetary_sciences_assignment_2.ipynb
MattTurnock/PlanetarySciencesMatt
81954d2182d9577bd7327a98c45963ae42968df4
MIT
planetary_sciences_assignment_2.ipynb
MattTurnock/PlanetarySciencesMatt
81954d2182d9577bd7327a98c45963ae42968df4
MIT
Qwen/Qwen-72B
1. YES
2. YES