1 #include <clocale>
2 #include <cstring>
3 #include <fstream>
4 #include <map>
5 #include <sstream>
6 #include <SpiceUsr.h>
7 using namespace std;
8
9 #include "findFile.h"
10 #include "keywords.h"
11 #include "Options.h"
12 #include "parse.h"
13 #include "sphericalToPixel.h"
14 #include "View.h"
15 #include "xpUtil.h"
16
17 #include "libannotate/libannotate.h"
18 #include "libplanet/Planet.h"
19 #include "libprojection/ProjectionBase.h"
20
21 bool
calculateSpicePosition(double jd,const int naifInt,Planet * relative,const int relativeInt,double & X,double & Y,double & Z)22 calculateSpicePosition(double jd, const int naifInt,
23 Planet *relative, const int relativeInt,
24 double &X, double &Y, double &Z)
25 {
26 SpiceInt naifId = static_cast<SpiceInt> (naifInt);
27 SpiceInt relativeTo = static_cast<SpiceInt> (relativeInt);
28
29 // seconds past J2000
30 const SpiceDouble et = ((jd - 2451545.0) * 86400);
31 SpiceDouble pos[3];
32 SpiceDouble lt;
33
34 spkgps_c(naifId, et, "J2000", relativeTo, pos, <);
35 if (return_c()) // true if spkgps_c fails
36 {
37 reset_c(); // reset the SPICE error flags
38 return(false);
39 }
40
41 // convert from km to AU
42 for (int i = 0; i < 3; i++) pos[i] /= AU_to_km;
43
44 relative->getPosition(X, Y, Z);
45 X += pos[0];
46 Y += pos[1];
47 Z += pos[2];
48
49 Options *options = Options::getInstance();
50 if (options->LightTime())
51 {
52 // Rectangular coordinates of the observer
53 double oX, oY, oZ;
54 options->getOrigin(oX, oY, oZ);
55
56 // Now get the position relative to the origin
57 double dX = X - oX;
58 double dY = Y - oY;
59 double dZ = Z - oZ;
60 double dist = sqrt(dX*dX + dY*dY + dZ*dZ);
61
62 double light_time = dist * AU_to_km / 299792.458;
63 spkgps_c(naifId, et - light_time, "J2000", relativeTo, pos, <);
64 if (return_c())
65 {
66 reset_c();
67 return(false);
68 }
69
70 for (int i = 0; i < 3; i++) pos[i] /= AU_to_km;
71
72 relative->getPosition(X, Y, Z);
73 X += pos[0];
74 Y += pos[1];
75 Z += pos[2];
76 }
77 return(true);
78 }
79
80 static void
readSpiceFile(const char * line,map<double,Planet * > & planetsFromSunMap,View * view,ProjectionBase * projection,multimap<double,Annotation * > & annotationMap)81 readSpiceFile(const char *line,
82 map<double, Planet *> &planetsFromSunMap,
83 View *view, ProjectionBase *projection,
84 multimap<double, Annotation *> &annotationMap)
85 {
86 int i = 0;
87 while (isDelimiter(line[i]))
88 {
89 i++;
90 if (static_cast<unsigned int> (i) > strlen(line)) return;
91 }
92 if (isEndOfLine(line[i])) return;
93
94 Options *options = Options::getInstance();
95
96 int align = AUTO;
97 string name("");
98
99 unsigned char color[3] = { 255, 0, 0 };
100 string font("");
101 int fontSize = -1;
102 string image("");
103 int symbolSize = 2;
104 bool syntaxError = false;
105
106 int thickness = 1;
107
108 double trailStart = 0;
109 double trailEnd = 0;
110 double trailInterval = 1;
111
112 bool transparency = false;
113 unsigned char transparent_pixel[3];
114
115 bool haveId = false;
116 int naifInt = 0; // Solar system barycenter
117 int relativeInt = 10; // SUN
118 SpiceChar relativeName[128];
119
120 SpiceBoolean found;
121 bodc2n_c(static_cast<SpiceInt> (relativeInt),
122 128, relativeName, &found);
123
124 while (static_cast<unsigned int> (i) < strlen(line))
125 {
126 char *returnString = NULL;
127 int val = parse(i, line, returnString);
128
129 switch (val)
130 {
131 case ALIGN:
132 if (returnString == NULL) break;
133 switch (returnString[0])
134 {
135 case 'r':
136 case 'R':
137 align = RIGHT;
138 break;
139 case 'l':
140 case 'L':
141 align = LEFT;
142 break;
143 case 'a':
144 case 'A':
145 align = ABOVE;
146 break;
147 case 'b':
148 case 'B':
149 align = BELOW;
150 break;
151 case 'c':
152 case 'C':
153 align = CENTER;
154 break;
155 default:
156 xpWarn("Unrecognized option for align in spice file\n",
157 __FILE__, __LINE__);
158 syntaxError = true;
159 break;
160 }
161 break;
162 case COLOR:
163 {
164 int r, g, b;
165 if (sscanf(returnString, "%d,%d,%d", &r, &g, &b) == 3)
166 {
167 color[0] = r & 0xff;
168 color[1] = g & 0xff;
169 color[2] = b & 0xff;
170 }
171 else
172 {
173 xpWarn("need three values for color\n", __FILE__, __LINE__);
174 syntaxError = true;
175 }
176 }
177 break;
178 case FONT:
179 font.assign(returnString);
180 break;
181 case FONTSIZE:
182 sscanf(returnString, "%d", &fontSize);
183 if (fontSize <= 0)
184 {
185 xpWarn("fontSize must be positive.\n", __FILE__, __LINE__);
186 syntaxError = true;
187 }
188 break;
189 case IMAGE:
190 image.assign(returnString);
191 break;
192 case LATLON:
193 {
194 sscanf(returnString, "%d", &naifInt);
195
196 SpiceChar bodyName[128];
197 bodc2n_c(static_cast<SpiceInt> (naifInt), 128, bodyName, &found);
198
199 if (found != SPICETRUE)
200 {
201 ostringstream errStr;
202 errStr << "Invalid NAIF id code: " << naifInt << ".\n";
203 xpWarn(errStr.str(), __FILE__, __LINE__);
204 syntaxError = true;
205 }
206 else
207 {
208 if (name.empty()) name.assign(bodyName);
209 haveId = true;
210 }
211 }
212 break;
213 case NAME:
214 name.assign(returnString);
215 break;
216 case ORIGIN: // relative_to keyword
217 {
218 sscanf(returnString, "%d", &relativeInt);
219 bodc2n_c(static_cast<SpiceInt> (relativeInt), 128,
220 relativeName, &found);
221 if (found != SPICETRUE)
222 {
223 ostringstream errStr;
224 errStr << "Invalid NAIF id code: " << relativeInt << ".\n";
225 xpWarn(errStr.str(), __FILE__, __LINE__);
226 syntaxError = true;
227 }
228 }
229 break;
230 case SYMBOLSIZE:
231 sscanf(returnString, "%d", &symbolSize);
232 if (symbolSize < 0) symbolSize = 2;
233 break;
234 case THICKNESS:
235 sscanf(returnString, "%d", &thickness);
236 if (thickness < 1)
237 {
238 xpWarn("thickness must be positive.\n",
239 __FILE__, __LINE__);
240 syntaxError = true;
241 }
242 break;
243 case TRAIL:
244 {
245 checkLocale(LC_NUMERIC, "C");
246 if (!sscanf(returnString, "%lf,%lf,%lf", &trailStart, &trailEnd,
247 &trailInterval) == 3)
248 {
249 xpWarn("Need three values for trail{}!\n",
250 __FILE__, __LINE__);
251 syntaxError = true;
252 }
253 else
254 {
255 if (trailInterval < 1e-4) trailInterval = 1e-4;
256 }
257 checkLocale(LC_NUMERIC, "");
258 }
259 break;
260 case TRANSPARENT:
261 {
262 int r, g, b;
263 if (sscanf(returnString, "%d,%d,%d", &r, &g, &b) == 3)
264 {
265 transparent_pixel[0] = r & 0xff;
266 transparent_pixel[1] = g & 0xff;
267 transparent_pixel[2] = b & 0xff;
268 }
269 else
270 {
271 xpWarn("Need three values for transparency pixel!\n",
272 __FILE__, __LINE__);
273 syntaxError = true;
274 }
275
276 transparency = true;
277 }
278 break;
279 case UNKNOWN:
280 syntaxError = true;
281 default:
282 case DELIMITER:
283 break;
284 }
285
286 if (val != DELIMITER && options->Verbosity() > 3)
287 {
288 ostringstream msg;
289 msg << "value is " << keyWordString[val - '?'];
290 if (returnString != NULL)
291 msg << ", returnString is " << returnString;
292 msg << endl;
293 xpMsg(msg.str(), __FILE__, __LINE__);
294 }
295
296 delete [] returnString;
297
298 if (syntaxError)
299 {
300 ostringstream errStr;
301 errStr << "Syntax error in spice file\n";
302 errStr << "line is \"" << line << "\"" << endl;
303 xpWarn(errStr.str(), __FILE__, __LINE__);
304 return;
305 }
306
307 if (val == ENDOFLINE) break;
308 }
309
310 if (!haveId)
311 {
312 ostringstream errStr;
313 errStr << "Can't compute position of " << name
314 << " relative to " << relativeName
315 << ", try a major planet or satellite.\n";
316 xpWarn(errStr.str(), __FILE__, __LINE__);
317 return;
318 }
319
320 Planet *relative = NULL;
321 for (map<double, Planet *>::iterator it0 = planetsFromSunMap.begin();
322 it0 != planetsFromSunMap.end(); it0++)
323 {
324 relative = it0->second;
325 if (relativeInt == naif_id[relative->Index()]) break;
326 relative = NULL;
327 }
328
329 if (relative == NULL) return;
330
331 const double jd = options->JulianDay();
332 double X, Y, Z;
333 if (calculateSpicePosition(jd, naifInt, relative, relativeInt, X, Y, Z))
334 {
335 if (options->Verbosity() > 1)
336 {
337 ostringstream xpStr;
338 xpStr << "Calculating position of " << name
339 << " relative to " << relativeName << endl;
340 xpMsg(xpStr.str(), __FILE__, __LINE__);
341 }
342 }
343 else
344 {
345 ostringstream errStr;
346 errStr << "Can't compute position of " << name
347 << " relative to " << relativeName << endl;
348 xpWarn(errStr.str(), __FILE__, __LINE__);
349 return;
350 }
351
352 bool plotThis = false;
353 if (view != NULL)
354 {
355 double pX, pY, pZ;
356 view->XYZToPixel(X, Y, Z, pX, pY, pZ);
357 pX += options->CenterX();
358 pY += options->CenterY();
359
360 plotThis = (pZ > 0);
361
362 // Rectangular coordinates of the observer
363 double oX, oY, oZ;
364 options->getOrigin(oX, oY, oZ);
365
366 // Now get the position relative to the origin
367 double dX = X - oX;
368 double dY = Y - oY;
369 double dZ = Z - oZ;
370 double dist = sqrt(dX*dX + dY*dY + dZ*dZ);
371
372 // don't plot this point if it's too close to a planet
373 for (map<double, Planet *>::iterator it0 = planetsFromSunMap.begin();
374 it0 != planetsFromSunMap.end(); it0++)
375 {
376 Planet *planet = it0->second;
377 double rX, rY, rZ;
378 planet->getPosition(rX, rY, rZ);
379 view->XYZToPixel(rX, rY, rZ, rX, rY, rZ);
380 rX += options->CenterX();
381 rY += options->CenterY();
382
383 double pixelDist = sqrt((rX - pX)*(rX - pX)
384 + (rY - pY)*(rY - pY));
385 if (pixelDist < 1)
386 {
387 double planetRadius = planet->Radius() / dist;
388 plotThis = (planetRadius / options->FieldOfView() > 0.01);
389 break;
390 }
391 }
392
393 if (options->Verbosity() > 0 && plotThis)
394 {
395 ostringstream msg;
396 char buffer[256];
397 char shortName[10];
398 memcpy(shortName, name.c_str(), 9);
399 shortName[9] = '\0';
400 snprintf(buffer, 256, "%10s%10.4f%8.1f%8.1f\n",
401 shortName, dist, pX, pY);
402 msg << buffer;
403 xpMsg(msg.str(), __FILE__, __LINE__);
404 }
405
406 X = pX;
407 Y = pY;
408 Z = pZ;
409 }
410 else
411 {
412 double lat, lon;
413 relative->XYZToPlanetographic(X, Y, Z, lat, lon);
414 plotThis = projection->sphericalToPixel(lon * relative->Flipped(),
415 lat, X, Y);
416 Z = 0;
417 }
418
419 if (plotThis)
420 {
421 const int ix = static_cast<int> (floor(X + 0.5));
422 const int iy = static_cast<int> (floor(Y + 0.5));
423
424 int xOffset = 0;
425 int yOffset = 0;
426 if (image.empty())
427 {
428 Symbol *sym = new Symbol(color, ix, iy, symbolSize);
429 annotationMap.insert(pair<const double, Annotation*>(Z, sym));
430 xOffset = symbolSize;
431 yOffset = symbolSize;
432 }
433 else if (image.compare("none") != 0)
434 {
435 unsigned char *transparent = (transparency ? transparent_pixel : NULL);
436 Icon *icon = new Icon(ix, iy, image, transparent);
437 annotationMap.insert(pair<const double, Annotation*>(Z, icon));
438 xOffset = icon->Width() / 2;
439 yOffset = icon->Height() / 2;
440 }
441
442 if (!name.empty())
443 {
444 Text *t = new Text(color, ix, iy, xOffset, yOffset,
445 align, name);
446 annotationMap.insert(pair<const double, Annotation*>(Z, t));
447 }
448 }
449
450 if (trailInterval > 0)
451 {
452 double X0, Y0, Z0;
453 double X1, Y1, Z1;
454 calculateSpicePosition(jd + trailStart, naifInt,
455 relative, relativeInt,
456 X0, Y0, Z0);
457 for (double et = trailStart + trailInterval;
458 et <= trailEnd; et += trailInterval)
459 {
460 calculateSpicePosition(jd + et, naifInt,
461 relative, relativeInt,
462 X1, Y1, Z1);
463 double newX0 = X1;
464 double newY0 = Y1;
465 double newZ0 = Z1;
466 if (view != NULL)
467 {
468 view->XYZToPixel(X0, Y0, Z0, X0, Y0, Z0);
469 X0 += options->CenterX();
470 Y0 += options->CenterY();
471
472 view->XYZToPixel(X1, Y1, Z1, X1, Y1, Z1);
473 X1 += options->CenterX();
474 Y1 += options->CenterY();
475 }
476 else
477 {
478 double lat, lon;
479 relative->XYZToPlanetographic(X0, Y0, Z0, lat, lon);
480 if (projection->sphericalToPixel(lon * relative->Flipped(),
481 lat, X0, Y0))
482 Z0 = 0;
483 else
484 Z0 = -1;
485
486 relative->XYZToPlanetographic(X1, Y1, Z1, lat, lon);
487 if (projection->sphericalToPixel(lon * relative->Flipped(),
488 lat, X1, Y1))
489 Z1 = 0;
490 else
491 Z1 = -1;
492 }
493
494 if (Z0 >= 0 && Z1 >= 0)
495 {
496 double Z = 0.5 * (Z0 + Z1);
497 LineSegment *ls = new LineSegment(color, thickness,
498 X1, Y1, X0, Y0);
499 annotationMap.insert(pair<const double, Annotation*>(Z, ls));
500 }
501 X0 = newX0;
502 Y0 = newY0;
503 Z0 = newZ0;
504 }
505 }
506 }
507
508 void
processSpiceKernels(const bool load)509 processSpiceKernels(const bool load)
510 {
511 // Set the SPICELIB error response action to "RETURN":
512 erract_c ( "SET", 200, "RETURN" );
513
514 // Output ALL CSPICE error messages on error:
515 errprt_c ( "SET", 200, "NONE, ALL" );
516
517 Options *options = Options::getInstance();
518 vector<string> spiceFiles = options->SpiceFiles();
519 for (unsigned int i = 0; i < spiceFiles.size(); i++)
520 {
521 string kernelFile(spiceFiles[i]);
522 kernelFile += ".krn";
523 if (findFile(kernelFile, "spice"))
524 {
525 ifstream inFile(kernelFile.c_str());
526 char *line = new char[MAX_LINE_LENGTH];
527 while (inFile.getline(line, MAX_LINE_LENGTH, '\n'))
528 {
529 int ii = 0;
530 while (isDelimiter(line[ii]))
531 {
532 ii++;
533 if (static_cast<unsigned int> (ii) > strlen(line))
534 continue;
535 }
536 if (isEndOfLine(line[ii])) continue;
537
538 char *ptr = &line[ii];
539 while (!(isDelimiter(*ptr) || isEndOfLine(*ptr)))
540 ptr++;
541 *ptr = '\0';
542
543 string spiceFile(&line[ii]);
544 if (findFile(spiceFile, "spice"))
545 {
546 if (load)
547 furnsh_c(spiceFile.c_str());
548 else
549 unload_c(spiceFile.c_str());
550 }
551 }
552
553 inFile.close();
554 delete [] line;
555 }
556 else
557 {
558 ostringstream errStr;
559 errStr << "Can't find spice kernel file " << kernelFile << endl;
560 xpWarn(errStr.str(), __FILE__, __LINE__);
561 }
562 }
563 }
564
565 void
addSpiceObjects(map<double,Planet * > & planetsFromSunMap,View * view,ProjectionBase * projection,multimap<double,Annotation * > & annotationMap)566 addSpiceObjects(map<double, Planet *> &planetsFromSunMap,
567 View *view, ProjectionBase *projection,
568 multimap<double, Annotation *> &annotationMap)
569 {
570 Options *options = Options::getInstance();
571 vector<string> spiceFiles = options->SpiceFiles();
572 for (unsigned int i = 0; i < spiceFiles.size(); i++)
573 {
574 string spiceFile(spiceFiles[i]);
575 if (findFile(spiceFile, "spice"))
576 {
577 ifstream inFile(spiceFile.c_str());
578 char *line = new char[MAX_LINE_LENGTH];
579 while (inFile.getline(line, MAX_LINE_LENGTH, '\n'))
580 readSpiceFile(line, planetsFromSunMap, view, projection,
581 annotationMap);
582 inFile.close();
583 delete [] line;
584 }
585 else
586 {
587 ostringstream errStr;
588 errStr << "Can't load spice file " << spiceFile << endl;
589 xpWarn(errStr.str(), __FILE__, __LINE__);
590 }
591 }
592 }
593