1 #include <cstdlib>
2 #include <cstring>
3 #include <map>
4 #include <sstream>
5 using namespace std;
6
7 #include "buildPlanetMap.h"
8 #include "findBodyXYZ.h"
9 #include "keywords.h"
10 #include "PlanetProperties.h"
11 #include "Options.h"
12 #include "readOriginFile.h"
13 #include "Separation.h"
14 #include "xpUtil.h"
15
16 #include "libplanet/Planet.h"
17
18 // if an origin file is being used, set the observer position from it
19 static void
setOriginXYZFromFile(const vector<LBRPoint> & originVector,const vector<LBRPoint>::iterator & iterOriginVector)20 setOriginXYZFromFile(const vector<LBRPoint> &originVector,
21 const vector<LBRPoint>::iterator &iterOriginVector)
22 {
23 Options *options = Options::getInstance();
24
25 // if an origin file is specified, set observer position from it
26 if (!options->OriginFile().empty())
27 {
28 if (options->InterpolateOriginFile())
29 {
30 // interpolate positions in the origin file to the
31 // current time
32 double thisRad, thisLat, thisLon, thisLocalTime = -1;
33 interpolateOriginFile(options->JulianDay(),
34 originVector, thisRad,
35 thisLat, thisLon, thisLocalTime);
36 options->Range(thisRad);
37 options->Latitude(thisLat);
38 options->Longitude(thisLon);
39 options->LocalTime(thisLocalTime);
40 }
41 else
42 {
43 // Use the next time in the origin file
44 options->setTime(iterOriginVector->time);
45
46 // Use the position, if specified, in the origin file
47 if (iterOriginVector->radius > 0)
48 {
49 options->Range(iterOriginVector->radius);
50 options->Latitude(iterOriginVector->latitude);
51 options->Longitude(iterOriginVector->longitude);
52
53 // Use the local time, if specified, in the origin file
54 if (iterOriginVector->localTime > 0)
55 options->LocalTime(iterOriginVector->localTime);
56 }
57 }
58 }
59
60 // if -dynamic_origin was specified, set observer position
61 // from it.
62 if (!options->DynamicOrigin().empty())
63 {
64 LBRPoint originPoint;
65 readDynamicOrigin(options->DynamicOrigin(), originPoint);
66 options->Range(originPoint.radius);
67 options->Latitude(originPoint.latitude);
68 options->Longitude(originPoint.longitude);
69 options->LocalTime(originPoint.localTime);
70 }
71 }
72
73 // set the position of the target. The coordinates are stored in the
74 // Options class.
75 static void
setTargetXYZ(PlanetProperties * planetProperties[])76 setTargetXYZ(PlanetProperties *planetProperties[])
77 {
78 Options *options = Options::getInstance();
79
80 switch (options->TargetMode())
81 {
82 case MAJOR:
83 case RANDOM:
84 {
85 bool no_random_target = true;
86 for (int i = 0; i < RANDOM_BODY; i++)
87 {
88 if (planetProperties[i]->RandomTarget())
89 {
90 no_random_target = false;
91 break;
92 }
93 }
94
95 if (no_random_target)
96 {
97 ostringstream errMsg;
98 errMsg << "random_target is false for all bodies. "
99 << "Check your config file.\n";
100 xpExit(errMsg.str(), __FILE__, __LINE__);
101 }
102
103 bool found_target = false;
104 while (!found_target)
105 {
106 body target = (body) (random() % RANDOM_BODY);
107
108 options->Target(target);
109
110 found_target = planetProperties[target]->RandomTarget();
111
112 if (found_target && options->TargetMode() == MAJOR)
113 {
114 Planet p(options->JulianDay(), options->Target());
115 found_target = (p.Primary() == SUN);
116 }
117 }
118 }
119 // fall through
120 case BODY:
121 {
122 Planet t(options->JulianDay(), options->Target());
123 options->Primary(t.Primary());
124 }
125 break;
126 case XYZ:
127 if (options->Target() == NORAD)
128 {
129 options->Primary(EARTH);
130 }
131 else
132 {
133 options->Primary(SUN);
134 }
135 break;
136 default:
137 xpExit("Unknown target mode?\n", __FILE__, __LINE__);
138 }
139
140 double tX, tY, tZ;
141 if (options->Target() == ALONG_PATH)
142 {
143 findBodyVelocity(options->JulianDay(), options->Origin(),
144 options->OriginID(), options->PathRelativeTo(),
145 options->PathRelativeToID(),
146 tX, tY, tZ);
147
148 tX *= FAR_DISTANCE;
149 tY *= FAR_DISTANCE;
150 tZ *= FAR_DISTANCE;
151 }
152 else
153 {
154 findBodyXYZ(options->JulianDay(), options->Target(),
155 options->TargetID(), tX, tY, tZ);
156 }
157 options->setTarget(tX, tY, tZ);
158 }
159
160 // set the position of the origin. The coordinates are stored in the
161 // Options class.
162 static void
setOriginXYZ(PlanetProperties * planetProperties[])163 setOriginXYZ(PlanetProperties *planetProperties[])
164 {
165 double oX, oY, oZ;
166
167 Options *options = Options::getInstance();
168 switch (options->OriginMode())
169 {
170 case LBR:
171 {
172 if (options->RandomLatLonRot())
173 {
174 double longitude;
175 longitude = random() % 360;
176 longitude *= deg_to_rad;
177 options->Longitude(longitude);
178
179 // Weight random latitudes towards the equator
180 double latitude;
181 latitude = (random() % 2000)/1000.0 - 1;
182 latitude = asin(latitude);
183 options->Latitude(latitude);
184
185 double rotate0;
186 rotate0 = random() % 360;
187 rotate0 *= deg_to_rad;
188 options->Rotate0(rotate0);
189 options->Rotate(0);
190 }
191
192 if (options->Origin() == NAIF || options->Origin() == NORAD)
193 {
194 findBodyXYZ(options->JulianDay(), options->Origin(),
195 options->OriginID(), oX, oY, oZ);
196 }
197 else
198 {
199 // Set the observer's lat & lon. Usually this is in the
200 // frame of the target body, but if the origin file option
201 // is used, the observer is in the origin body's
202 // geographic frame.
203 body referenceBody = options->Target();
204 if (!options->OriginFile().empty()
205 || !options->DynamicOrigin().empty())
206 referenceBody = options->Origin();
207
208 Planet p(options->JulianDay(), referenceBody);
209 p.calcHeliocentricEquatorial();
210
211 if (options->LocalTime() >= 0)
212 {
213 double subSolarLat = 0;
214 double subSolarLon = 0;
215 p.XYZToPlanetographic(0, 0, 0, subSolarLat, subSolarLon);
216
217 double longitude;
218 longitude = (subSolarLon - M_PI
219 + p.Flipped() * options->LocalTime() * M_PI / 12);
220 options->Longitude(longitude);
221 }
222 p.PlanetographicToXYZ(oX, oY, oZ,
223 options->Latitude(), options->Longitude(),
224 options->Range());
225 }
226 }
227 break;
228 case MAJOR:
229 case RANDOM:
230 case SYSTEM:
231 {
232 // check to see that at least one body has random_origin=true
233 bool no_random_origin = true;
234 for (int i = 0; i < RANDOM_BODY; i++)
235 {
236 if (i == options->Target()) continue;
237
238 if (planetProperties[i]->RandomOrigin())
239 {
240 no_random_origin = false;
241 break;
242 }
243 }
244
245 if (no_random_origin)
246 {
247 ostringstream errMsg;
248 errMsg << "Target is "
249 << planetProperties[options->Target()]->Name()
250 << ", random_origin is false for all other bodies. "
251 << "Check your config file.\n";
252 xpExit(errMsg.str(), __FILE__, __LINE__);
253 }
254
255 bool found_origin = false;
256 while (!found_origin)
257 {
258 options->Origin(static_cast<body> (random() % RANDOM_BODY));
259
260 if (options->Verbosity() > 1)
261 {
262 ostringstream msg;
263 msg << "target = " << body_string[options->Target()]
264 << ", origin = " << body_string[options->Origin()] << endl;
265 xpMsg(msg.str(), __FILE__, __LINE__);
266 }
267
268 if (options->Target() == options->Origin()) continue;
269
270 if (options->OriginMode() == RANDOM)
271 {
272 found_origin = true;
273 }
274 else
275 {
276 Planet o(options->JulianDay(), options->Origin());
277
278 if (options->OriginMode() == MAJOR)
279 {
280 found_origin = (o.Primary() == SUN);
281 }
282 else if (options->OriginMode() == SYSTEM)
283 {
284 // SYSTEM means one of three things:
285 // 1) origin is target's primary
286 // 2) target is origin's primary
287 // 3) target and origin have same primary
288 found_origin = ((options->Origin() == options->Primary()
289 || o.Primary() == options->Target()
290 || options->Primary() == o.Primary()));
291 }
292 }
293
294 if (found_origin)
295 found_origin = planetProperties[options->Origin()]->RandomOrigin();
296 } // while (!found_origin)
297 }
298 // fall through
299 case BODY:
300 {
301 findBodyXYZ(options->JulianDay(), options->Origin(),
302 options->OriginID(), oX, oY, oZ);
303
304 if (options->OppositeSide())
305 {
306 double tX, tY, tZ;
307 options->getTarget(tX, tY, tZ);
308 oX = 2*tX - oX;
309 oY = 2*tY - oY;
310 oZ = 2*tZ - oZ;
311 }
312 }
313 break;
314 case ABOVE:
315 case BELOW:
316 {
317 if (options->Target() == SUN)
318 {
319 ostringstream errMsg;
320 errMsg << "-origin above or below not applicable for SUN. "
321 << "Use -latitude 90 or -90 instead\n";
322 xpExit(errMsg.str(), __FILE__, __LINE__);
323 }
324
325 double pX, pY, pZ;
326 findBodyXYZ(options->JulianDay(), options->Primary(), -1, pX, pY, pZ);
327
328 double vX, vY, vZ;
329 findBodyVelocity(options->JulianDay(), options->Target(),
330 options->TargetID(), options->Primary(), -1,
331 vX, vY, vZ);
332
333 // cross product of position and velocity vectors points to the
334 // orbital north pole
335 double tX, tY, tZ;
336 options->getTarget(tX, tY, tZ);
337 double pos[3] = { tX - pX, tY - pY, tZ - pZ };
338 double vel[3] = { vX, vY, vZ };
339
340 double north[3];
341 cross(pos, vel, north);
342
343 double mag = sqrt(dot(north, north));
344 Planet primary(options->JulianDay(), options->Primary());
345 double dist = FAR_DISTANCE * primary.Radius();
346
347 oX = dist * north[0]/mag + pX;
348 oY = dist * north[1]/mag + pY;
349 oZ = dist * north[2]/mag + pZ;
350
351 const double radius = sqrt(pos[0] * pos[0] + pos[1] * pos[1]
352 + pos[2] * pos[2]);
353
354 // will only be zero if -fov or -radius haven't been specified
355 if (options->FieldOfView() < 0)
356 {
357 options->FieldOfView(4*radius/dist); // fit the orbit on the screen
358 options->FOVMode(FOV);
359 }
360
361 if (options->OriginMode() == BELOW)
362 {
363 oX -= 2 * (oX - pX);
364 oY -= 2 * (oY - pX);
365 oZ -= 2 * (oZ - pX);
366 }
367 }
368 break;
369 }
370
371 if (options->RangeSpecified() && options->TargetMode() != XYZ)
372 {
373 Planet p2(options->JulianDay(), options->Target());
374 p2.calcHeliocentricEquatorial();
375
376 double latitude, longitude;
377 p2.XYZToPlanetographic(oX, oY, oZ, latitude, longitude);
378 p2.PlanetographicToXYZ(oX, oY, oZ, latitude, longitude,
379 options->Range());
380 }
381
382 if (options->SeparationTarget() < RANDOM_BODY)
383 {
384 // rectangular coordinates of the separation target
385 double sX, sY, sZ;
386 findBodyXYZ(options->JulianDay(), options->SeparationTarget(), -1,
387 sX, sY, sZ);
388
389 // rectangular coordinates of the primary target
390 double tX, tY, tZ;
391 Planet p(options->JulianDay(), options->Target());
392 p.calcHeliocentricEquatorial();
393 p.getPosition(tX, tY, tZ);
394
395 // Save the range
396 double latitude, longitude, range;
397 p.XYZToPlanetographic(oX, oY, oZ, latitude, longitude, range);
398
399 // place the observer 90 degrees from the separation target in
400 // the primary target's equatorial plane
401 double n[3];
402 p.getBodyNorth(n[0], n[1], n[2]);
403
404 double s[3] = { sX - tX, sY - tY, sZ - tZ };
405 double c[3];
406 cross(s, n, c);
407 c[0] += tX;
408 c[1] += tY;
409 c[2] += tZ;
410
411 p.XYZToPlanetographic(c[0], c[1], c[2], latitude, longitude);
412 p.PlanetographicToXYZ(oX, oY, oZ, latitude, longitude, range);
413
414 Separation sep(oX, oY, oZ, tX, tY, tZ, sX, sY, sZ);
415 sep.setSeparation(options->SeparationDist() * deg_to_rad);
416 sep.getOrigin(oX, oY, oZ);
417 }
418
419 options->setOrigin(oX, oY, oZ);
420 }
421
422 // Set the observer's lat/lon with respect to the target body, if
423 // appropriate
424 void
getObsLatLon(Planet * target,PlanetProperties * planetProperties[])425 getObsLatLon(Planet *target, PlanetProperties *planetProperties[])
426 {
427 Options *options = Options::getInstance();
428
429 // XYZ mode is where the target isn't a planetary body.
430 // (e.g. the Cassini spacecraft)
431 if (options->TargetMode() == XYZ)
432 {
433 if (options->LightTime()) setTargetXYZ(planetProperties);
434 }
435 else
436 {
437 if (target == NULL)
438 xpExit("Can't find target body?\n", __FILE__, __LINE__);
439
440 if (options->LightTime())
441 {
442 double tX, tY, tZ;
443 target->getPosition(tX, tY, tZ);
444 options->setTarget(tX, tY, tZ);
445 }
446
447 // Rectangular coordinates of the observer
448 double oX, oY, oZ;
449 options->getOrigin(oX, oY, oZ);
450
451 // Find the sub-observer lat & lon
452 double obs_lat, obs_lon;
453 target->XYZToPlanetographic(oX, oY, oZ, obs_lat, obs_lon);
454 options->Latitude(obs_lat);
455 options->Longitude(obs_lon);
456
457 // Find the sub-solar lat & lon. This is used for the image
458 // label
459 double sunLat, sunLon;
460 target->XYZToPlanetographic(0, 0, 0, sunLat, sunLon);
461 options->SunLat(sunLat);
462 options->SunLon(sunLon);
463 }
464 }
465
466 // Set the direction of the "Up" vector
467 void
setUpXYZ(const Planet * target,map<double,Planet * > & planetsFromSunMap,double & upX,double & upY,double & upZ)468 setUpXYZ(const Planet *target, map<double, Planet *> &planetsFromSunMap,
469 double &upX, double &upY, double &upZ)
470 {
471 Options *options = Options::getInstance();
472
473 switch (options->North())
474 {
475 default:
476 xpWarn("Unknown value for north, using body\n",
477 __FILE__, __LINE__);
478 case BODY:
479 target->getBodyNorth(upX, upY, upZ);
480 break;
481 case GALACTIC:
482 {
483 const double GN_LAT = 27.4 * deg_to_rad;
484 const double GN_LON = 192.25 * deg_to_rad;
485
486 upX = cos(GN_LAT) * cos(GN_LON);
487 upY = cos(GN_LAT) * sin(GN_LON);
488 upZ = sin(GN_LAT);
489 }
490 break;
491 case ORBIT:
492 {
493 // if it's a moon, pretend its orbital north is the same as
494 // its primary, although in most cases it's the same as its
495 // primary's rotational north
496 if (target->Primary() == SUN)
497 {
498 target->getOrbitalNorth(upX, upY, upZ);
499 }
500 else
501 {
502 Planet *primary = findPlanetinMap(planetsFromSunMap,
503 target->Primary());
504 primary->getOrbitalNorth(upX, upY, upZ);
505 }
506 }
507 break;
508 case PATH:
509 {
510 double vX, vY, vZ;
511 findBodyVelocity(options->JulianDay(),
512 options->Origin(),
513 options->OriginID(),
514 options->PathRelativeTo(),
515 options->PathRelativeToID(),
516 vX, vY, vZ);
517
518 upX = vX;
519 upY = vY;
520 upZ = vZ;
521 }
522 break;
523 case SEPARATION:
524 {
525 if (options->SeparationTarget() < RANDOM_BODY)
526 {
527 double sX, sY, sZ;
528 double tX, tY, tZ;
529 double oX, oY, oZ;
530
531 findBodyXYZ(options->JulianDay(), options->SeparationTarget(),
532 options->TargetID(), sX, sY, sZ);
533 options->getOrigin(oX, oY, oZ);
534 options->getTarget(tX, tY, tZ);
535
536 double t[3] = {tX - oX, tY - oY, tZ - oZ};
537 double s[3] = {sX - oX, sY - oY, sZ - oZ};
538 double c[3];
539 cross(s, t, c);
540
541 if (sin(options->SeparationDist() * deg_to_rad) > 0)
542 {
543 upX = c[0];
544 upY = c[1];
545 upZ = c[2];
546 }
547 else
548 {
549 upX = -c[0];
550 upY = -c[1];
551 upZ = -c[2];
552 }
553 }
554 else
555 {
556 xpWarn("No separation target given, using -north body\n",
557 __FILE__, __LINE__);
558 target->getBodyNorth(upX, upY, upZ);
559 }
560 }
561 break;
562 case TERRESTRIAL:
563 upX = 0;
564 upY = 0;
565 upZ = 1;
566 break;
567 }
568
569 double up[3] = {upX, upY, upZ};
570 normalize(up);
571 upX = up[0];
572 upY = up[1];
573 upZ = up[2];
574 }
575
576 void
setPositions(const vector<LBRPoint> & originVector,const vector<LBRPoint>::iterator & iterOriginVector,Planet * & target,map<double,Planet * > & planetMap,PlanetProperties * planetProperties[])577 setPositions(const vector<LBRPoint> &originVector,
578 const vector<LBRPoint>::iterator &iterOriginVector,
579 Planet *&target, map<double, Planet *> &planetMap,
580 PlanetProperties *planetProperties[])
581 {
582 Options *options = Options::getInstance();
583
584 // set the observer's XYZ coordinates if an origin file has
585 // been specified
586 setOriginXYZFromFile(originVector, iterOriginVector);
587
588 // Set the target's XYZ coordinates.
589 setTargetXYZ(planetProperties);
590
591 // Set the origin's XYZ coordinates.
592 setOriginXYZ(planetProperties);
593
594 // Find the target body in the list
595 body target_body = options->Target();
596 target = findPlanetinMap(planetMap, target_body);
597
598 if (options->LightTime())
599 {
600 double oX, oY, oZ;
601 options->getOrigin(oX, oY, oZ);
602
603 // Destroy the existing planet map and build a new one
604 buildPlanetMap(options->JulianDay(), oX, oY, oZ, true,
605 planetMap);
606
607 target = findPlanetinMap(planetMap, target_body);
608
609 // Set the target's XYZ coordinates.
610 setTargetXYZ(planetProperties);
611
612 // Set the origin's XYZ coordinates.
613 setOriginXYZ(planetProperties);
614 }
615
616 // Now find the observer's lat/lon
617 getObsLatLon(target, planetProperties);
618 }
619