1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.models.earth.weather;
18
19 import java.util.List;
20 import java.util.function.ToDoubleFunction;
21
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.util.FastMath;
24 import org.hipparchus.util.FieldSinCos;
25 import org.hipparchus.util.MathUtils;
26 import org.hipparchus.util.SinCos;
27 import org.orekit.errors.OrekitException;
28 import org.orekit.errors.OrekitMessages;
29 import org.orekit.utils.Constants;
30 import org.orekit.utils.units.Unit;
31
32
33
34
35
36
37 class Grid {
38
39
40 private final Indexer latitudeIndexer;
41
42
43 private final Indexer longitudeIndexer;
44
45
46 private final GridEntry[][] entries;
47
48
49
50
51
52 Grid(final List<GridEntry> loadedEntries, final String name) {
53
54
55 latitudeIndexer = new Indexer(loadedEntries, GridEntry::getLatitude, name);
56 longitudeIndexer = new Indexer(loadedEntries, GridEntry::getLongitude, name);
57
58
59 entries = new GridEntry[latitudeIndexer.n][longitudeIndexer.n + 1];
60 for (final GridEntry entry : loadedEntries) {
61 final int ia = latitudeIndexer.closeIndex(entry.getLatitude());
62 final int io = longitudeIndexer.closeIndex(entry.getLongitude());
63 entries[ia][io] = entry;
64 }
65
66
67 for (int ia = 0; ia < latitudeIndexer.n; ia++) {
68 if (entries[ia][0] != null) {
69 entries[ia][longitudeIndexer.n] = entries[ia][0].buildWrappedEntry();
70 }
71 }
72
73
74 for (final GridEntry[] row : entries) {
75 for (final GridEntry entry : row) {
76 if (entry == null) {
77 throw new OrekitException(OrekitMessages.IRREGULAR_OR_INCOMPLETE_GRID, name);
78 }
79 }
80 }
81
82 }
83
84
85
86
87
88 private int getSouthIndex(final double latitude) {
89
90 return FastMath.min(latitudeIndexer.lowIndex(latitude), latitudeIndexer.n - 2);
91 }
92
93
94
95
96
97 private int getWestIndex(final double longitude) {
98
99 return longitudeIndexer.lowIndex(longitude);
100 }
101
102
103
104
105
106
107
108
109 CellInterpolator getInterpolator(final double latitude, final double longitude,
110 final double altitude, final double deltaRef) {
111
112
113 final double normalizedLongitude =
114 MathUtils.normalizeAngle(longitude,
115 entries[0][0].getLongitude() + FastMath.PI);
116
117
118 final int southIndex = getSouthIndex(latitude);
119 final int westIndex = getWestIndex(normalizedLongitude);
120
121 final double coef = (deltaRef / Constants.JULIAN_YEAR) * 2 * FastMath.PI;
122 final SinCos sc1 = FastMath.sinCos(coef);
123 final SinCos sc2 = FastMath.sinCos(2.0 * coef);
124
125
126 return new CellInterpolator(latitude, normalizedLongitude,
127 entries[southIndex ][westIndex ].evaluate(sc1, sc2, altitude),
128 entries[southIndex ][westIndex + 1].evaluate(sc1, sc2, altitude),
129 entries[southIndex + 1][westIndex ].evaluate(sc1, sc2, altitude),
130 entries[southIndex + 1][westIndex + 1].evaluate(sc1, sc2, altitude));
131
132 }
133
134
135
136
137
138
139
140
141
142 <T extends CalculusFieldElement<T>> FieldCellInterpolator<T> getInterpolator(final T latitude, final T longitude,
143 final T altitude, final T deltaRef) {
144
145
146 final T normalizedLongitude =
147 MathUtils.normalizeAngle(longitude,
148 longitude.newInstance(entries[0][0].getLongitude() + FastMath.PI));
149
150
151 final int southIndex = getSouthIndex(latitude.getReal());
152 final int westIndex = getWestIndex(normalizedLongitude.getReal());
153
154 final T coef = deltaRef.multiply(2 * FastMath.PI / Constants.JULIAN_YEAR);
155 final FieldSinCos<T> sc1 = FastMath.sinCos(coef);
156 final FieldSinCos<T> sc2 = FastMath.sinCos(coef.multiply(2));
157
158
159 return new FieldCellInterpolator<>(latitude, normalizedLongitude,
160 entries[southIndex ][westIndex ].evaluate(sc1, sc2, altitude),
161 entries[southIndex ][westIndex + 1].evaluate(sc1, sc2, altitude),
162 entries[southIndex + 1][westIndex ].evaluate(sc1, sc2, altitude),
163 entries[southIndex + 1][westIndex + 1].evaluate(sc1, sc2, altitude));
164
165 }
166
167
168
169
170
171 boolean hasModels(final SeasonalModelType... types) {
172 boolean hasAll = true;
173 for (final SeasonalModelType type : types) {
174 hasAll &= entries[0][0].hasModel(type);
175 }
176 return hasAll;
177 }
178
179
180
181
182 private static class Indexer {
183
184
185 private final double min;
186
187
188 private final double step;
189
190
191 private final int n;
192
193
194
195
196
197
198 Indexer(final List<GridEntry> entries, final ToDoubleFunction<GridEntry> extractor, final String name) {
199
200 final double tolerance = Unit.parse("mas").toSI(1.0);
201
202
203 double inf = Double.POSITIVE_INFINITY;
204 double sup = Double.NEGATIVE_INFINITY;
205 for (final GridEntry entry : entries) {
206 final double coordinate = extractor.applyAsDouble(entry);
207 inf = FastMath.min(inf, coordinate);
208 sup = FastMath.max(sup, coordinate);
209 }
210
211
212 double firstStep = Double.POSITIVE_INFINITY;
213 for (final GridEntry entry : entries) {
214 final double delta = extractor.applyAsDouble(entry) - inf;
215 if (delta > tolerance) {
216
217 firstStep = FastMath.min(firstStep, delta);
218 }
219 }
220
221
222 this.min = inf;
223 this.step = firstStep;
224 this.n = 1 + (int) FastMath.rint((sup - inf) / firstStep);
225
226
227 for (final GridEntry entry : entries) {
228 final double coordinate = extractor.applyAsDouble(entry);
229 final double rebuilt = min + closeIndex(coordinate) * step;
230 if (FastMath.abs(coordinate - rebuilt) > tolerance) {
231 throw new OrekitException(OrekitMessages.IRREGULAR_OR_INCOMPLETE_GRID, name);
232 }
233 }
234
235 }
236
237
238
239
240
241 public int lowIndex(final double coordinate) {
242 return (int) FastMath.floor((coordinate - min) / step);
243 }
244
245
246
247
248
249 public int closeIndex(final double coordinate) {
250 return (int) FastMath.rint((coordinate - min) / step);
251 }
252
253 }
254
255 }