forked from TheAlgorithms/Python
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbasic_orbital_capture.py
181 lines (132 loc) · 5.43 KB
/
basic_orbital_capture.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
from math import pow, sqrt
from scipy.constants import G, c, pi
"""
These two functions will return the radii of impact for a target object
of mass M and radius R as well as it's effective cross sectional area σ(sigma).
That is to say any projectile with velocity v passing within σ, will impact the
target object with mass M. The derivation of which is given at the bottom
of this file.
The derivation shows that a projectile does not need to aim directly at the target
body in order to hit it, as R_capture>R_target. Astronomers refer to the effective
cross section for capture as σ=π*R_capture**2.
This algorithm does not account for an N-body problem.
"""
def capture_radii(
target_body_radius: float, target_body_mass: float, projectile_velocity: float
) -> float:
"""
Input Params:
-------------
target_body_radius: Radius of the central body SI units: meters | m
target_body_mass: Mass of the central body SI units: kilograms | kg
projectile_velocity: Velocity of object moving toward central body
SI units: meters/second | m/s
Returns:
--------
>>> capture_radii(6.957e8, 1.99e30, 25000.0)
17209590691.0
>>> capture_radii(-6.957e8, 1.99e30, 25000.0)
Traceback (most recent call last):
...
ValueError: Radius cannot be less than 0
>>> capture_radii(6.957e8, -1.99e30, 25000.0)
Traceback (most recent call last):
...
ValueError: Mass cannot be less than 0
>>> capture_radii(6.957e8, 1.99e30, c+1)
Traceback (most recent call last):
...
ValueError: Cannot go beyond speed of light
Returned SI units:
------------------
meters | m
"""
if target_body_mass < 0:
raise ValueError("Mass cannot be less than 0")
elif target_body_radius < 0:
raise ValueError("Radius cannot be less than 0")
elif projectile_velocity > c:
raise ValueError("Cannot go beyond speed of light")
else:
escape_velocity_squared = (2 * G * target_body_mass) / target_body_radius
capture_radius = target_body_radius * sqrt(
1 + escape_velocity_squared / pow(projectile_velocity, 2)
)
return round(capture_radius, 0)
def capture_area(capture_radius: float) -> float:
"""
Input Param:
------------
capture_radius: The radius of orbital capture and impact for a central body of
mass M and a projectile moving towards it with velocity v
SI units: meters | m
Returns:
--------
>>> capture_area(17209590691)
9.304455331329126e+20
>>> capture_area(-1)
Traceback (most recent call last):
...
ValueError: Cannot have a capture radius less than 0
Returned SI units:
------------------
meters*meters | m**2
"""
if capture_radius < 0:
raise ValueError("Cannot have a capture radius less than 0")
else:
sigma = pi * pow(capture_radius, 2)
return round(sigma, 0)
if __name__ == "__main__":
from doctest import testmod
testmod()
"""
Derivation:
Let: Mt=target mass, Rt=target radius, v=projectile_velocity,
r_0=radius of projectile at instant 0 to CM of target
v_p=v at closest approach,
r_p=radius from projectile to target CM at closest approach,
R_capture= radius of impact for projectile with velocity v
(1)At time=0 the projectile's energy falling from infinity| E=K+U=0.5*m*(v**2)+0
E_initial=0.5*m*(v**2)
(2)at time=0 the angular momentum of the projectile relative to CM target|
L_initial=m*r_0*v*sin(Θ)->m*r_0*v*(R_capture/r_0)->m*v*R_capture
L_i=m*v*R_capture
(3)The energy of the projectile at closest approach will be its kinetic energy
at closest approach plus gravitational potential energy(-(GMm)/R)|
E_p=K_p+U_p->E_p=0.5*m*(v_p**2)-(G*Mt*m)/r_p
E_p=0.0.5*m*(v_p**2)-(G*Mt*m)/r_p
(4)The angular momentum of the projectile relative to the target at closest
approach will be L_p=m*r_p*v_p*sin(Θ), however relative to the target Θ=90°
sin(90°)=1|
L_p=m*r_p*v_p
(5)Using conservation of angular momentum and energy, we can write a quadratic
equation that solves for r_p|
(a)
Ei=Ep-> 0.5*m*(v**2)=0.5*m*(v_p**2)-(G*Mt*m)/r_p-> v**2=v_p**2-(2*G*Mt)/r_p
(b)
Li=Lp-> m*v*R_capture=m*r_p*v_p-> v*R_capture=r_p*v_p-> v_p=(v*R_capture)/r_p
(c) b plugs int a|
v**2=((v*R_capture)/r_p)**2-(2*G*Mt)/r_p->
v**2-(v**2)*(R_c**2)/(r_p**2)+(2*G*Mt)/r_p=0->
(v**2)*(r_p**2)+2*G*Mt*r_p-(v**2)*(R_c**2)=0
(d) Using the quadratic formula, we'll solve for r_p then rearrange to solve to
R_capture
r_p=(-2*G*Mt ± sqrt(4*G^2*Mt^2+ 4(v^4*R_c^2)))/(2*v^2)->
r_p=(-G*Mt ± sqrt(G^2*Mt+v^4*R_c^2))/v^2->
r_p<0 is something we can ignore, as it has no physical meaning for our purposes.->
r_p=(-G*Mt)/v^2 + sqrt(G^2*Mt^2/v^4 + R_c^2)
(e)We are trying to solve for R_c. We are looking for impact, so we want r_p=Rt
Rt + G*Mt/v^2 = sqrt(G^2*Mt^2/v^4 + R_c^2)->
(Rt + G*Mt/v^2)^2 = G^2*Mt^2/v^4 + R_c^2->
Rt^2 + 2*G*Mt*Rt/v^2 + G^2*Mt^2/v^4 = G^2*Mt^2/v^4 + R_c^2->
Rt**2 + 2*G*Mt*Rt/v**2 = R_c**2->
Rt**2 * (1 + 2*G*Mt/Rt *1/v**2) = R_c**2->
escape velocity = sqrt(2GM/R)= v_escape**2=2GM/R->
Rt**2 * (1 + v_esc**2/v**2) = R_c**2->
(6)
R_capture = Rt * sqrt(1 + v_esc**2/v**2)
Source: Problem Set 3 #8 c.Fall_2017|Honors Astronomy|Professor Rachel Bezanson
Source #2: http://www.nssc.ac.cn/wxzygx/weixin/201607/P020160718380095698873.pdf
8.8 Planetary Rendezvous: Pg.368
"""