1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.models.earth.atmosphere;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.util.FastMath;
21 import org.hipparchus.util.FieldSinCos;
22 import org.hipparchus.util.MathUtils;
23 import org.orekit.annotation.DefaultDataContext;
24 import org.orekit.bodies.BodyShape;
25 import org.orekit.data.DataContext;
26 import org.orekit.time.AbsoluteDate;
27 import org.orekit.time.TimeScale;
28 import org.orekit.utils.ExtendedPositionProvider;
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 public class JB2006 extends AbstractJacchiaBowmanModel {
66
67
68 private static final double[] FZM = { 0.111613e+00, -0.159000e-02, 0.126190e-01, -0.100064e-01, -0.237509e-04, 0.260759e-04};
69
70
71 private static final double[] GTM = {-0.833646e+00, -0.265450e+00, 0.467603e+00, -0.299906e+00, -0.105451e+00,
72 -0.165537e-01, -0.380037e-01, -0.150991e-01, -0.541280e-01, 0.119554e-01,
73 0.437544e-02, -0.369016e-02, 0.206763e-02, -0.142888e-02, -0.867124e-05,
74 0.189032e-04, 0.156988e-03, 0.491286e-03, -0.391484e-04, -0.126854e-04,
75 0.134078e-04, -0.614176e-05, 0.343423e-05};
76
77
78 private final JB2006InputParameters inputParams;
79
80
81
82
83
84
85
86
87 @DefaultDataContext
88 public JB2006(final JB2006InputParameters parameters, final ExtendedPositionProvider sun, final BodyShape earth) {
89 this(parameters, sun, earth, DataContext.getDefault().getTimeScales().getUTC());
90 }
91
92
93
94
95
96
97
98
99
100 public JB2006(final JB2006InputParameters parameters, final ExtendedPositionProvider sun,
101 final BodyShape earth, final TimeScale utc) {
102 super(sun, utc, earth, parameters.getMinDate(), parameters.getMaxDate());
103 this.inputParams = parameters;
104 }
105
106
107 @Override
108 protected double computeTInf(final AbsoluteDate date, final double tsubl, final double dtclst) {
109 return tsubl + getDtg(date) + dtclst;
110 }
111
112
113 @Override
114 protected <T extends CalculusFieldElement<T>> T computeTInf(final AbsoluteDate date, final T tsubl, final T dtclst) {
115 return tsubl.add(getDtg(date)).add(dtclst);
116 }
117
118
119 @Override
120 protected double computeTc(final AbsoluteDate date) {
121 final double f10 = inputParams.getF10(date);
122 final double f10B = inputParams.getF10B(date);
123 final double s10 = inputParams.getS10(date);
124 final double s10B = inputParams.getS10B(date);
125 final double xm10 = inputParams.getXM10(date);
126 final double xm10B = inputParams.getXM10B(date);
127 return 379 + 3.353 * f10B + 0.358 * (f10 - f10B) + 2.094 * (s10 - s10B) + 0.343 * (xm10 - xm10B);
128 }
129
130
131 @Override
132 protected double getF10(final AbsoluteDate date) {
133 return inputParams.getF10(date);
134 }
135
136
137 @Override
138 protected double getF10B(final AbsoluteDate date) {
139 return inputParams.getF10B(date);
140 }
141
142
143 @Override
144 protected double semian(final AbsoluteDate date, final double day, final double height) {
145
146 final double f10Bar = inputParams.getF10B(date);
147 final double f10Bar2 = f10Bar * f10Bar;
148 final double htz = height / 1000.0;
149
150
151 final double fzz = FZM[0] + FZM[1] * f10Bar + FZM[2] * f10Bar * htz + FZM[3] * f10Bar * htz * htz + FZM[4] * f10Bar * f10Bar * htz + FZM[5] * f10Bar * f10Bar * htz * htz;
152
153
154 final double tau = MathUtils.TWO_PI * (day - 1.0) / 365;
155 final double sin1P = FastMath.sin(tau);
156 final double cos1P = FastMath.cos(tau);
157 final double sin2P = FastMath.sin(2.0 * tau);
158 final double cos2P = FastMath.cos(2.0 * tau);
159 final double sin3P = FastMath.sin(3.0 * tau);
160 final double cos3P = FastMath.cos(3.0 * tau);
161 final double sin4P = FastMath.sin(4.0 * tau);
162 final double cos4P = FastMath.cos(4.0 * tau);
163 final double gtz = GTM[0] +
164 GTM[1] * sin1P +
165 GTM[2] * cos1P +
166 GTM[3] * sin2P +
167 GTM[4] * cos2P +
168 GTM[5] * sin3P +
169 GTM[6] * cos3P +
170 GTM[7] * sin4P +
171 GTM[8] * cos4P +
172 GTM[9] * f10Bar +
173 GTM[10] * f10Bar * sin1P +
174 GTM[11] * f10Bar * cos1P +
175 GTM[12] * f10Bar * sin2P +
176 GTM[13] * f10Bar * cos2P +
177 GTM[14] * f10Bar * sin3P +
178 GTM[15] * f10Bar * cos3P +
179 GTM[16] * f10Bar * sin4P +
180 GTM[17] * f10Bar * cos4P +
181 GTM[18] * f10Bar2 +
182 GTM[19] * f10Bar2 * sin1P +
183 GTM[20] * f10Bar2 * cos1P +
184 GTM[21] * f10Bar2 * sin2P +
185 GTM[22] * f10Bar2 * cos2P;
186
187 return FastMath.max(1.0e-6, fzz) * gtz;
188 }
189
190 @Override
191 protected <T extends CalculusFieldElement<T>> T semian(final AbsoluteDate date, final T doy, final T height) {
192
193 final double f10Bar = getF10B(date);
194 final double f10Bar2 = f10Bar * f10Bar;
195 final T htz = height.divide(1000.0);
196
197
198 final T fzz = htz.multiply(FZM[2] * f10Bar).add(htz.square().multiply(FZM[3] * f10Bar)).add(htz.multiply(FZM[4] * f10Bar * f10Bar)).add(htz.square().multiply(FZM[5] * f10Bar * f10Bar)).add(FZM[0] + FZM[1] * f10Bar);
199
200
201 final T tau = doy.subtract(1).divide(365).multiply(MathUtils.TWO_PI);
202 final FieldSinCos<T> sc1P = FastMath.sinCos(tau);
203 final FieldSinCos<T> sc2P = FastMath.sinCos(tau.multiply(2.0));
204 final FieldSinCos<T> sc3P = FastMath.sinCos(tau.multiply(3.0));
205 final FieldSinCos<T> sc4P = FastMath.sinCos(tau.multiply(4.0));
206 final T gtz = sc1P.sin().multiply(GTM[1]).add(
207 sc1P.cos().multiply(GTM[2])).add(
208 sc2P.sin().multiply(GTM[3])).add(
209 sc2P.cos().multiply(GTM[4])).add(
210 sc3P.sin().multiply(GTM[5])).add(
211 sc3P.cos().multiply(GTM[6])).add(
212 sc4P.sin().multiply(GTM[7])).add(
213 sc4P.cos().multiply(GTM[8])).add(
214 GTM[9] * f10Bar).add(
215 sc1P.sin().multiply(f10Bar).multiply(GTM[10])).add(
216 sc1P.cos().multiply(f10Bar).multiply(GTM[11])).add(
217 sc2P.sin().multiply(f10Bar).multiply(GTM[12])).add(
218 sc2P.cos().multiply(f10Bar).multiply(GTM[13])).add(
219 sc3P.sin().multiply(f10Bar).multiply(GTM[14])).add(
220 sc3P.cos().multiply(f10Bar).multiply(GTM[15])).add(
221 sc4P.sin().multiply(f10Bar).multiply(GTM[16])).add(
222 sc4P.cos().multiply(f10Bar).multiply(GTM[17])).add(
223 GTM[18] * f10Bar2).add(
224 sc1P.sin().multiply(f10Bar2).multiply(GTM[19])).add(
225 sc1P.cos().multiply(f10Bar2).multiply(GTM[20])).add(
226 sc2P.sin().multiply(f10Bar2).multiply(GTM[21])).add(
227 sc2P.cos().multiply(f10Bar2).multiply(GTM[22])).add(GTM[0]);
228
229 return fzz.getReal() > 1.0e-6 ? gtz.multiply(fzz) : gtz.multiply(1.0e-6);
230 }
231
232
233
234
235
236 private double getDtg(final AbsoluteDate date) {
237
238 final double ap = inputParams.getAp(date);
239 final double expAp = FastMath.exp(-0.08 * ap);
240 return ap + 100. * (1. - expAp);
241 }
242 }