iframe-proxy | Sunbelt Computer Software

History

177 lines (130 loc) · 5.4 KB

Raw

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

"""

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 sigma, 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 sigma=π*R_capture**2.

This algorithm does not account for an N-body problem.

"""

from math import pow, sqrt # noqa: A004

from scipy.constants import G, c, pi

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")

if target_body_radius < 0:

raise ValueError("Radius cannot be less than 0")

if projectile_velocity > c:

raise ValueError("Cannot go beyond speed of light")

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")

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

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

"""

Sunbelt Computer Software

PL/B Language Development and Support

Navigation Menu

Sunbelt Computer Software

PL/B Language Development and Support

FilesExpand file tree

basic_orbital_capture.py

Latest commit

History

basic_orbital_capture.py

File metadata and controls