forked from poliastro/validation
-
Notifications
You must be signed in to change notification settings - Fork 0
/
sidereal.m
87 lines (79 loc) · 3.31 KB
/
sidereal.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
% ----------------------------------------------------------------------------
%
% function sidereal
%
% this function calulates the transformation matrix that accounts for the
% effects of sidereal time. Notice that deltaspi should not be mod'ed to a
% positive number because it is multiplied rather than used in a
% trigonometric argument.
%
% author : david vallado 719-573-2600 25 jun 2002
%
% revisions
% vallado - fix units on kinematic terms 5 sep 2002
% vallado - add terms 30 sep 2002
% vallado - consolidate with iau 2000 14 feb 2005
%
% inputs description range / units
% jdut1 - julian centuries of ut1 days
% deltapsi - nutation angle rad
% meaneps - mean obliquity of the ecliptic rad
% omega - long of asc node of moon rad
% lod - length of day sec
% eqeterms - terms for ast calculation 0,2
%
% outputs :
% st - transformation matrix for pef - tod
% stdot - transformation matrix for pef - tod rate
%
% locals :
% gmst - mean greenwich sidereal time 0 to 2pi rad
% ast - apparent gmst 0 to 2pi rad
% hr - hour hr
% min - minutes min
% sec - seconds sec
% temp - temporary vector
% tempval - temporary variable
%
% coupling :
%
% references :
% vallado 2013, 223-224
%
% [st, stdot] = sidereal (jdut1, deltapsi, meaneps, omega, lod, eqeterms );
% ----------------------------------------------------------------------------
function [st,stdot] = sidereal (jdut1, deltapsi, meaneps, omega, lod, eqeterms )
% ------------------------ find gmst --------------------------
gmst= gstime( jdut1 );
% ------------------------ find mean ast ----------------------
% after 1997, kinematic terms apply as well as gemoetric in eqe
if (jdut1 > 2450449.5 ) && (eqeterms > 0)
ast= gmst + deltapsi* cos(meaneps) ...
+ 0.00264*pi /(3600*180)*sin(omega) ...
+ 0.000063*pi /(3600*180)*sin(2.0 *omega);
else
ast= gmst + deltapsi* cos(meaneps);
end
ast = rem (ast, 2.0*pi);
thetasa = 7.29211514670698e-05 * (1.0 - lod/86400.0 );
omegaearth = thetasa;
%fprintf(1,'st gmst %11.8f ast %11.8f ome %11.8f \n', gmst*180/pi, ast*180/pi, omegaearth*180/pi );
st(1,1) = cos(ast);
st(1,2) = -sin(ast);
st(1,3) = 0.0;
st(2,1) = sin(ast);
st(2,2) = cos(ast);
st(2,3) = 0.0;
st(3,1) = 0.0;
st(3,2) = 0.0;
st(3,3) = 1.0;
% compute sidereal time rate matrix
stdot(1,1) = -omegaearth * sin(ast);
stdot(1,2) = -omegaearth * cos(ast);
stdot(1,3) = 0.0;
stdot(2,1) = omegaearth * cos(ast);
stdot(2,2) = -omegaearth * sin(ast);
stdot(2,3) = 0.0;
stdot(3,1) = 0.0;
stdot(3,2) = 0.0;
stdot(3,3) = 0.0;