-
Notifications
You must be signed in to change notification settings - Fork 1
/
coast.m
48 lines (33 loc) · 905 Bytes
/
coast.m
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
function dZ = coast(~,Z)
global CRAFT PLANET ATMOSPHERE
T = CRAFT(1);
II = CRAFT(2);
IF = CRAFT(3);
G = PLANET(1);
R = PLANET(2);
RS = PLANET(3);
S = PLANET(4);
D = ATMOSPHERE(1);
H = ATMOSPHERE(3);
x = Z(1);
y = Z(2);
vx = Z(3);
vy = Z(4);
m = Z(5);
p = [x,y];
pang = atan2(y,x)-pi/2;
vr = [vx-cos(pang)*RS, vy-sin(pang)*RS];
d = norm(p);
h = d-R;
ap = exp(-h/H);
grav = -S/d^3*p;
drag = -D*ap*norm(vr)*vr;
drag_ang = atan2(drag(2),drag(1));
throttle = m*norm(drag)/T; % throttle thrust to counteract drag exactly
thrust = throttle*T/m*[cos(drag_ang+180), sin(drag_ang+180)];
a = grav+drag+thrust;
dm = throttle*T/(((IF-II)*ap-IF)*G);
ax = a(1);
ay = a(2);
dZ = [vx; vy; ax; ay; dm];
end