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, &lt);
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, &lt);
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