1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.radiation;
18
19 import java.util.ArrayList;
20 import java.util.Collections;
21 import java.util.List;
22
23 import org.hipparchus.CalculusFieldElement;
24 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25 import org.hipparchus.geometry.euclidean.threed.Vector3D;
26 import org.hipparchus.util.FastMath;
27 import org.hipparchus.util.FieldSinCos;
28 import org.hipparchus.util.SinCos;
29 import org.orekit.propagation.FieldSpacecraftState;
30 import org.orekit.propagation.SpacecraftState;
31 import org.orekit.utils.ExtendedPVCoordinatesProvider;
32 import org.orekit.utils.ParameterDriver;
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 public class ECOM2 extends AbstractRadiationForceModel {
74
75
76 public static final String ECOM_COEFFICIENT = "ECOM coefficient";
77
78
79 private static final double MIN_VALUE = Double.NEGATIVE_INFINITY;
80
81
82 private static final double MAX_VALUE = Double.POSITIVE_INFINITY;
83
84
85
86
87
88
89
90 private final double SCALE = FastMath.scalb(1.0, -22);
91
92
93 private final int nD;
94
95
96 private final int nB;
97
98
99
100
101
102
103
104 private final List<ParameterDriver> coefficients;
105
106
107 private final ExtendedPVCoordinatesProvider sun;
108
109
110
111
112
113
114
115
116
117 public ECOM2(final int nD, final int nB, final double value,
118 final ExtendedPVCoordinatesProvider sun, final double equatorialRadius) {
119 super(sun, equatorialRadius);
120 this.nB = nB;
121 this.nD = nD;
122 this.coefficients = new ArrayList<>(2 * (nD + nB) + 3);
123
124
125 coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " B0", value, SCALE, MIN_VALUE, MAX_VALUE));
126 for (int i = 1; i < nB + 1; i++) {
127 coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Bcos" + Integer.toString(i - 1), value, SCALE, MIN_VALUE, MAX_VALUE));
128 }
129 for (int i = nB + 1; i < 2 * nB + 1; i++) {
130 coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Bsin" + Integer.toString(i - (nB + 1)), value, SCALE, MIN_VALUE, MAX_VALUE));
131 }
132
133 coefficients.add(2 * nB + 1, new ParameterDriver(ECOM_COEFFICIENT + " D0", value, SCALE, MIN_VALUE, MAX_VALUE));
134 for (int i = 2 * nB + 2; i < 2 * nB + 2 + nD; i++) {
135 coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Dcos" + Integer.toString(i - (2 * nB + 2)), value, SCALE, MIN_VALUE, MAX_VALUE));
136 }
137 for (int i = 2 * nB + 2 + nD; i < 2 * (nB + nD) + 2; i++) {
138 coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Dsin" + Integer.toString(i - (2 * nB + nD + 2)), value, SCALE, MIN_VALUE, MAX_VALUE));
139 }
140
141 coefficients.add(new ParameterDriver(ECOM_COEFFICIENT + " Y0", value, SCALE, MIN_VALUE, MAX_VALUE));
142
143
144 coefficients.forEach(parameter -> parameter.setSelected(true));
145 this.sun = sun;
146 }
147
148
149 @Override
150 public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {
151
152 final Vector3D Z = s.getPVCoordinates().getMomentum().normalize();
153 final Vector3D sunPos = sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition().normalize();
154 final Vector3D Y = Z.crossProduct(sunPos);
155 final Vector3D X = Y.crossProduct(Z);
156
157 final Vector3D position = s.getPVCoordinates().getPosition().normalize();
158 final Vector3D eD = sunPos.add(-1.0, position);
159 final Vector3D eY = position.crossProduct(eD);
160 final Vector3D eB = eD.crossProduct(eY);
161
162
163 final double delta_u = FastMath.atan2(position.dotProduct(Y), position.dotProduct(X));
164
165
166 double b_u = parameters[0];
167 for (int i = 1; i < nB + 1; i++) {
168 final SinCos sc = FastMath.sinCos(2 * i * delta_u);
169 b_u += parameters[i] * sc.cos() + parameters[i + nB] * sc.sin();
170 }
171
172 double d_u = parameters[2 * nB + 1];
173 for (int i = 1; i < nD + 1; i++) {
174 final SinCos sc = FastMath.sinCos((2 * i - 1) * delta_u);
175 d_u += parameters[2 * nB + 1 + i] * sc.cos() + parameters[2 * nB + 1 + i + nD] * sc.sin();
176 }
177
178 return new Vector3D(d_u, eD, parameters[2 * (nD + nB) + 2], eY, b_u, eB);
179 }
180
181
182 @Override
183 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s, final T[] parameters) {
184
185 final FieldVector3D<T> Z = s.getPVCoordinates().getMomentum().normalize();
186 final FieldVector3D<T> sunPos = sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition().normalize();
187 final FieldVector3D<T> Y = Z.crossProduct(sunPos);
188 final FieldVector3D<T> X = Y.crossProduct(Z);
189
190
191 final FieldVector3D<T> position = s.getPVCoordinates().getPosition().normalize();
192 final FieldVector3D<T> eD = sunPos.add(-1.0, position);
193 final FieldVector3D<T> eY = position.crossProduct(eD);
194 final FieldVector3D<T> eB = eD.crossProduct(eY);
195
196
197 final T delta_u = FastMath.atan2(position.dotProduct(Y), position.dotProduct(X));
198
199
200 T b_u = parameters[0];
201 for (int i = 1; i < nB + 1; i++) {
202 final FieldSinCos<T> sc = FastMath.sinCos(delta_u.multiply(2 * i));
203 b_u = b_u.add(sc.cos().multiply(parameters[i])).add(sc.sin().multiply(parameters[i + nB]));
204 }
205
206 T d_u = parameters[2 * nB + 1];
207
208 for (int i = 1; i < nD + 1; i++) {
209 final FieldSinCos<T> sc = FastMath.sinCos(delta_u.multiply(2 * i - 1));
210 d_u = d_u.add(sc.cos().multiply(parameters[2 * nB + 1 + i])).add(sc.sin().multiply(parameters[2 * nB + 1 + i + nD]));
211 }
212
213 return new FieldVector3D<>(d_u, eD, parameters[2 * (nD + nB) + 2], eY, b_u, eB);
214 }
215
216
217 @Override
218 public List<ParameterDriver> getParametersDrivers() {
219 return Collections.unmodifiableList(coefficients);
220 }
221
222 }
223