-
Notifications
You must be signed in to change notification settings - Fork 3
/
IncFunc.m
105 lines (83 loc) · 2.58 KB
/
IncFunc.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
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
%IncFunc Inclination Function using the Allan form.
% See Yan, H., Vadali, S. R., & Alfriend, K. T. (2014). A Simplified Formulation
% of the Satellite Perturbed Relative Motion Problem.
%
% Author: Bharat Mahajan (https://github.com/princemahajan)
function [ F, Fd, Fbysi2, Fdd, Fbysi2d] = IncFunc( n,m,p,i )
% coefficient free of summation
t1 = (-1)^fix((n-m+1)/2)*factorial(n+m)/((2^n)*factorial(n-p)*factorial(p));
si2 = sin(i/2);
ci2 = cos(i/2);
% j range
ji = max([0, n-m-2*p]);
jf = min([n-m, 2*n-2*p]);
% j summation
jsum = 0;
jsumbysi2 = 0;
jsumd = 0;
jsumd2 = 0;
jsumbysi2d = 0;
for jctr = ji:1:jf
% binomial coeff-1
try
n1 = 2*n-2*p;
k1 = jctr;
bc1 = nchoosek(n1,k1);
catch
bc1 = 0;
end
% binomial coeff-2
try
n2 = 2*p;
k2 = n-m-jctr;
bc2 = nchoosek(n2,k2);
catch
bc2 = 0;
end
% sum for F
jsum = jsum + ((-1)^jctr)*bc1*bc2*(si2.^(m-n+2*p+2*jctr)).*ci2.^(3*n-m-2*jctr-2*p);
% sum for F/s(i/2)
jsumbysi2 = jsumbysi2 + ((-1)^jctr)*bc1*bc2*(si2.^(m-n+2*p+2*jctr-1)).*ci2.^(3*n-m-2*jctr-2*p);
% First Derivative
if (m-n+2*p+2*jctr) == 0
sterm = 0;
else
sterm = (m-n+2*p+2*jctr)*(si2.^(m-n+2*p+2*jctr-1)).*ci2/2.*(ci2.^(3*n-m-2*jctr-2*p));
end
if (3*n-m-2*jctr-2*p) == 0
cterm = 0;
else
cterm = -(si2.^(m-n+2*p+2*jctr))*(3*n-m-2*jctr-2*p).*(ci2.^(3*n-m-2*jctr-2*p-1)).*si2/2;
end
jsumd = jsumd + ((-1)^jctr)*bc1*bc2*(sterm + cterm);
% derivative of F/si2
if (m-n+2*p+2*jctr-1) == 0
sterm = 0;
else
sterm = (m-n+2*p+2*jctr-1)*(si2.^(m-n+2*p+2*jctr-2)).*ci2/2.*(ci2.^(3*n-m-2*jctr-2*p));
end
if (3*n-m-2*jctr-2*p) == 0
cterm = 0;
else
cterm = -(si2.^(m-n+2*p+2*jctr-1))*(3*n-m-2*jctr-2*p).*(ci2.^(3*n-m-2*jctr-2*p-1)).*si2/2;
end
jsumbysi2d = jsumbysi2d + ((-1)^jctr)*bc1*bc2*(sterm + cterm);
% Second Derivative
s = m - n + 2*p + 2*jctr;
csterm1 = 0;
csterm2 = 0;
if s ~= 0 && s ~= 1
csterm2 = s*(s-1)*ci2.^(2*n-s+2)*si2.^(s-2);
end
if (2*n-s) ~= 0 && (2*n-s) ~= 1
csterm1 = (2*n-s)*(2*n-s-1)*si2.^(s+2).*ci2.^(2*n-s-2);
end
jsumd2 = jsumd2 + ((-1)^jctr)*bc1*bc2*1/4*(ci2.^(2*n-s).*si2.^s.*(-(2*n-s)*(s+1) - (2*n-s+1)*s) + csterm1 + csterm2);
end
Fbysi2 = t1*jsumbysi2;
F = t1*jsum;
% Derivatives of Inclination Function
Fd = t1*jsumd;
Fdd = t1*jsumd2;
Fbysi2d = t1*jsumbysi2d;
end