source: ntrip/trunk/BNC/bnctides.cpp@ 2581

Last change on this file since 2581 was 2581, checked in by mervart, 14 years ago
File size: 7.0 KB
Line 
1
2#include <cmath>
3#include <iostream>
4#include <iomanip>
5
6#include "bnctides.h"
7#include "bncutils.h"
8
9using namespace std;
10
11// Auxiliary Functions
12///////////////////////////////////////////////////////////////////////////
13namespace {
14
15 static const double RHO_DEG = 180.0 / M_PI;
16 static const double RHO_SEC = 3600.0 * RHO_DEG;
17 static const double MJD_J2000 = 51544.5;
18
19 double Frac (double x) { return x-floor(x); };
20 double Modulo (double x, double y) { return y*Frac(x/y); }
21
22 Matrix rotX(double Angle) {
23 const double C = cos(Angle);
24 const double S = sin(Angle);
25 Matrix UU(3,3);
26 UU[0][0] = 1.0; UU[0][1] = 0.0; UU[0][2] = 0.0;
27 UU[1][0] = 0.0; UU[1][1] = +C; UU[1][2] = +S;
28 UU[2][0] = 0.0; UU[2][1] = -S; UU[2][2] = +C;
29 return UU;
30 }
31
32 Matrix rotY(double Angle) {
33 const double C = cos(Angle);
34 const double S = sin(Angle);
35 Matrix UU(3,3);
36 UU[0][0] = +C; UU[0][1] = 0.0; UU[0][2] = -S;
37 UU[1][0] = 0.0; UU[1][1] = 1.0; UU[1][2] = 0.0;
38 UU[2][0] = +S; UU[2][1] = 0.0; UU[2][2] = +C;
39 return UU;
40 }
41
42 Matrix rotZ(double Angle) {
43 const double C = cos(Angle);
44 const double S = sin(Angle);
45 Matrix UU(3,3);
46 UU[0][0] = +C; UU[0][1] = +S; UU[0][2] = 0.0;
47 UU[1][0] = -S; UU[1][1] = +C; UU[1][2] = 0.0;
48 UU[2][0] = 0.0; UU[2][1] = 0.0; UU[2][2] = 1.0;
49 return UU;
50 }
51}
52
53// Greenwich Mean Sidereal Time
54///////////////////////////////////////////////////////////////////////////
55double GMST(double Mjd_UT1) {
56
57 const double Secs = 86400.0;
58
59 double Mjd_0 = floor(Mjd_UT1);
60 double UT1 = Secs*(Mjd_UT1-Mjd_0);
61 double T_0 = (Mjd_0 -MJD_J2000)/36525.0;
62 double T = (Mjd_UT1-MJD_J2000)/36525.0;
63
64 double gmst = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1
65 + (0.093104-6.2e-6*T)*T*T;
66
67 return 2.0*M_PI*Frac(gmst/Secs);
68}
69
70// Nutation Matrix
71///////////////////////////////////////////////////////////////////////////
72Matrix NutMatrix(double Mjd_TT) {
73
74 const double T = (Mjd_TT-MJD_J2000)/36525.0;
75
76 double ls = 2.0*M_PI*Frac(0.993133+ 99.997306*T);
77 double D = 2.0*M_PI*Frac(0.827362+1236.853087*T);
78 double F = 2.0*M_PI*Frac(0.259089+1342.227826*T);
79 double N = 2.0*M_PI*Frac(0.347346- 5.372447*T);
80
81 double dpsi = ( -17.200*sin(N) - 1.319*sin(2*(F-D+N)) - 0.227*sin(2*(F+N))
82 + 0.206*sin(2*N) + 0.143*sin(ls) ) / RHO_SEC;
83 double deps = ( + 9.203*cos(N) + 0.574*cos(2*(F-D+N)) + 0.098*cos(2*(F+N))
84 - 0.090*cos(2*N) ) / RHO_SEC;
85
86 double eps = 0.4090928-2.2696E-4*T;
87
88 return rotX(-eps-deps)*rotZ(-dpsi)*rotX(+eps);
89}
90
91// Precession Matrix
92///////////////////////////////////////////////////////////////////////////
93Matrix PrecMatrix (double Mjd_1, double Mjd_2) {
94
95 const double T = (Mjd_1-MJD_J2000)/36525.0;
96 const double dT = (Mjd_2-Mjd_1)/36525.0;
97
98 double zeta = ( (2306.2181+(1.39656-0.000139*T)*T)+
99 ((0.30188-0.000344*T)+0.017998*dT)*dT )*dT/RHO_SEC;
100 double z = zeta + ( (0.79280+0.000411*T)+0.000205*dT)*dT*dT/RHO_SEC;
101 double theta = ( (2004.3109-(0.85330+0.000217*T)*T)-
102 ((0.42665+0.000217*T)+0.041833*dT)*dT )*dT/RHO_SEC;
103
104 return rotZ(-z) * rotY(theta) * rotZ(-zeta);
105}
106
107// Sun's position
108///////////////////////////////////////////////////////////////////////////
109ColumnVector Sun(double Mjd_TT) {
110
111 const double eps = 23.43929111/RHO_DEG;
112 const double T = (Mjd_TT-MJD_J2000)/36525.0;
113
114 double M = 2.0*M_PI * Frac ( 0.9931267 + 99.9973583*T);
115 double L = 2.0*M_PI * Frac ( 0.7859444 + M/2.0*M_PI +
116 (6892.0*sin(M)+72.0*sin(2.0*M)) / 1296.0e3);
117 double r = 149.619e9 - 2.499e9*cos(M) - 0.021e9*cos(2*M);
118
119 ColumnVector r_Sun(3);
120 r_Sun << r*cos(L) << r*sin(L) << 0.0; r_Sun = rotX(-eps) * r_Sun;
121
122 return rotZ(GMST(Mjd_TT))
123 * NutMatrix(Mjd_TT)
124 * PrecMatrix(MJD_J2000, Mjd_TT)
125 * r_Sun;
126}
127
128// Moon's position
129///////////////////////////////////////////////////////////////////////////
130ColumnVector Moon(double Mjd_TT) {
131
132 const double eps = 23.43929111/RHO_DEG;
133 const double T = (Mjd_TT-MJD_J2000)/36525.0;
134
135 double L_0 = Frac ( 0.606433 + 1336.851344*T );
136 double l = 2.0*M_PI*Frac ( 0.374897 + 1325.552410*T );
137 double lp = 2.0*M_PI*Frac ( 0.993133 + 99.997361*T );
138 double D = 2.0*M_PI*Frac ( 0.827361 + 1236.853086*T );
139 double F = 2.0*M_PI*Frac ( 0.259086 + 1342.227825*T );
140
141 double dL = +22640*sin(l) - 4586*sin(l-2*D) + 2370*sin(2*D) + 769*sin(2*l)
142 -668*sin(lp) - 412*sin(2*F) - 212*sin(2*l-2*D)- 206*sin(l+lp-2*D)
143 +192*sin(l+2*D) - 165*sin(lp-2*D) - 125*sin(D) - 110*sin(l+lp)
144 +148*sin(l-lp) - 55*sin(2*F-2*D);
145
146 double L = 2.0*M_PI * Frac( L_0 + dL/1296.0e3 );
147
148 double S = F + (dL+412*sin(2*F)+541*sin(lp)) / RHO_SEC;
149 double h = F-2*D;
150 double N = -526*sin(h) + 44*sin(l+h) - 31*sin(-l+h) - 23*sin(lp+h)
151 +11*sin(-lp+h) - 25*sin(-2*l+F) + 21*sin(-l+F);
152
153 double B = ( 18520.0*sin(S) + N ) / RHO_SEC;
154
155 double cosB = cos(B);
156
157 double R = 385000e3 - 20905e3*cos(l) - 3699e3*cos(2*D-l) - 2956e3*cos(2*D)
158 -570e3*cos(2*l) + 246e3*cos(2*l-2*D) - 205e3*cos(lp-2*D)
159 -171e3*cos(l+2*D) - 152e3*cos(l+lp-2*D);
160
161 ColumnVector r_Moon(3);
162 r_Moon << R*cos(L)*cosB << R*sin(L)*cosB << R*sin(B);
163 r_Moon = rotX(-eps) * r_Moon;
164
165 return rotZ(GMST(Mjd_TT))
166 * NutMatrix(Mjd_TT)
167 * PrecMatrix(MJD_J2000, Mjd_TT)
168 * r_Moon;
169}
170
171// Tidal Correction
172////////////////////////////////////////////////////////////////////////////
173void tides(const bncTime& time, ColumnVector& xyz) {
174
175 static double lastMjd = 0.0;
176 static ColumnVector xSun;
177 static ColumnVector xMoon;
178 static double rSun;
179 static double rMoon;
180
181 double Mjd = time.mjd() + time.daysec() / 86400.0;
182
183 if (Mjd != lastMjd) {
184 lastMjd = Mjd;
185 xSun = Sun(Mjd);
186 rSun = sqrt(DotProduct(xSun,xSun));
187 xSun /= rSun;
188 xMoon = Moon(Mjd);
189 rMoon = sqrt(DotProduct(xMoon,xMoon));
190 xMoon /= rMoon;
191 }
192
193 double rRec = sqrt(DotProduct(xyz, xyz));
194 ColumnVector xyzUnit = xyz / rRec;
195
196 // Love's Numbers
197 // --------------
198 const double H2 = 0.6090;
199 const double L2 = 0.0852;
200
201 // Tidal Displacement
202 // ------------------
203 double scSun = DotProduct(xyzUnit, xSun);
204 double scMoon = DotProduct(xyzUnit, xMoon);
205
206 double p2Sun = 3.0 * (H2/2.0-L2) * scSun * scSun - H2/2.0;
207 double p2Moon = 3.0 * (H2/2.0-L2) * scMoon * scMoon - H2/2.0;
208
209 double x2Sun = 3.0 * L2 * scSun;
210 double x2Moon = 3.0 * L2 * scMoon;
211
212 const double gmWGS = 398.6005e12;
213 const double gms = 1.3271250e20;
214 const double gmm = 4.9027890e12;
215
216 double facSun = gms / gmWGS *
217 (rRec * rRec * rRec * rRec) / (rSun * rSun * rSun);
218
219 double facMoon = gmm / gmWGS *
220 (rRec * rRec * rRec * rRec) / (rMoon * rMoon * rMoon);
221
222 ColumnVector dX = facSun * (x2Sun * xSun + p2Sun * xyzUnit) +
223 facMoon * (x2Moon * xMoon + p2Moon * xyzUnit);
224
225 xyz += dX;
226}
Note: See TracBrowser for help on using the repository browser.