-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathvplanet.c
2012 lines (1724 loc) · 68.4 KB
/
vplanet.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*
Solar System Live -- Main server program
by John Walker
http://www.fourmilab.ch/
*/
#include "vplanet.h"
/* #define TESTMODE */ /* Direct subsequent requests to Solar.t script */
/* If CACHE_WARNING is defined, generated HTML will contain a warning
to users not to link to the ephemeral cache files. */
/* #define CACHE_WARNING */ /**/
#define TimeLimit 25 /* CPU time limit in seconds */
#ifdef TimeLimit
#define __USE_XOPEN_EXTENDED /* Required to define sigset() with GCC library */
#include <signal.h>
#include <sys/time.h>
#include <sys/resource.h>
#endif
double siteLat = 47.0, siteLon = -7.0;
char orbitalElements[maxOrbitalElements][132]; /* Elements of body being tracked */
char savedOrbitalElements[maxOrbitalElements][132]; /* Elements submitted with dynamic request */
char *progpath; /* Program invocation path (argv[0]) */
static int directCall = FALSE; /* Were we invoked directly, without uncgi first ? */
static int embeddedImage = FALSE; /* Is this an embedded image specified with a "?di=" block ? */
char *elementSpecification = NULL; /* Contents of orbit elements textbox in the request */
static int dto = 0; /* Date/time option */
static BITMAPFILEHEADER bfh;
static struct bmih {
BITMAPINFOHEADER bmiHeader;
RGBQUAD bmiColors[256];
} bmi, *tmbmap;
static int bmRowLen;
static unsigned char *bmBits = NULL, *tmbits;
static unsigned char ct[4];
static double jt, sunra, sundec, sunlong, gt;
static struct tm ctim;
static char viewdesc[132], dtutc[80];
static int transparent = -1, cBlue, cGreen, cBorder = -1, cWhite,
cGrey, cLtGrey, cDkRed, cRed, cBlack;
static int dumpelem = FALSE; /* Dump orbital elements ? */
static int stereoView = FALSE;
static int skyColourScheme = SCS_COLOUR; /* Colour scheme */
static double stereoOcular = 6; /* Ocular separation - cross, + wall */
#define stereoSepPixels 10
static int useimage = FALSE; /* Use images for planets ? */
#ifdef SAVESET
static int bookmark = FALSE; /* Save parameters for bookmark ? */
static int gotten = FALSE; /* Was "get" method used ? */
#endif
#ifdef DBLOG
static FILE *db;
#define DB(x) fprintf(db, "%s\n", x)
#else
#define DB(x)
#endif
extern char **environ; /* Environment variables */
#ifdef TimeLimit
/* DINGALING -- Handle SIGALRM when time runs out. */
static void dingaling(int sig/*, int code, struct sigcontext *scp, char *addr*/)
{
struct rusage u;
if (sig == SIGALRM) {
getrusage(RUSAGE_SELF, &u);
if ((u.ru_utime.tv_sec + u.ru_stime.tv_sec) > TimeLimit) {
printf("<h1><font color=\"#FF0000\">Time limit (%d seconds) exceeded.</font></h1>\n", TimeLimit);
exit(0);
}
alarm(TimeLimit); /* Wind the cat */
}
}
#endif
/* LOAD_BITMAP -- Load a Microsoft Windows .BMP file into memory. */
#define rB(x) (x) = getc(fp)
#define rS(x) dstat = fread(ct, 2, 1, fp); (x) = ((ct[1] << 8) | ct[0])
#define rL(x) dstat = fread(ct, 4, 1, fp); (x) = (((long) ct[3] << 24) | ((long) ct[2] << 16) | (ct[1] << 8) | ct[0])
static void load_bitmap(const char *bmpfile)
{
FILE *fp = fopen(bmpfile, "r");
int dstat;
if (bmBits != NULL) {
free(bmBits);
bmBits = NULL;
}
if (fp != NULL) {
rS(bfh.bfType);
rL(bfh.bfSize);
rS(bfh.bfReserved1);
rS(bfh.bfReserved2);
rL(bfh.bfOffBits);
rL(bmi.bmiHeader.biSize);
rL(bmi.bmiHeader.biWidth);
rL(bmi.bmiHeader.biHeight);
rS(bmi.bmiHeader.biPlanes);
rS(bmi.bmiHeader.biBitCount);
rL(bmi.bmiHeader.biCompression);
rL(bmi.bmiHeader.biSizeImage);
rL(bmi.bmiHeader.biXPelsPerMeter);
rL(bmi.bmiHeader.biYPelsPerMeter);
rL(bmi.bmiHeader.biClrUsed);
rL(bmi.bmiHeader.biClrImportant);
assert(bmi.bmiHeader.biBitCount == 8);
assert(bmi.bmiHeader.biCompression == BI_RGB);
if (bmi.bmiHeader.biClrUsed == 0) {
bmi.bmiHeader.biClrUsed = 1 << bmi.bmiHeader.biBitCount;
}
/* Read the colour palette table. */
dstat = fread(&bmi.bmiColors[0], sizeof(RGBQUAD), bmi.bmiHeader.biClrUsed, fp);
/* Compute the row length and read the bitmap. */
bmRowLen = (bmi.bmiHeader.biWidth + 3) & ~3;
fseek(fp, bfh.bfOffBits, SEEK_SET);
bmBits = (unsigned char *) malloc(bmRowLen * bmi.bmiHeader.biHeight);
dstat = fread(bmBits, bmRowLen * bmi.bmiHeader.biHeight, 1, fp);
fclose(fp);
dstat++; /* Purely to get rid of "set but not used" warning */
}
}
static void dump_ppm(const char *fname, const struct bmih *bi, const unsigned char *bits)
{
FILE *fp;
int i, j;
if (fname[0] == '-') {
fp = stdout;
} else {
fp = fopen(fname, "w");
}
assert(fp != NULL);
fprintf(fp, "P6 %ld %ld %d\n",
(long)( bi->bmiHeader.biWidth), (long) (bi->bmiHeader.biHeight), 255);
for (i = 0; i < bi->bmiHeader.biHeight; i++) {
unsigned const char *px = bits +
(((bi->bmiHeader.biHeight - 1) - i) *
((bi->bmiHeader.biWidth + 3) & ~3));
for (j = 0; j < bi->bmiHeader.biWidth; j++) {
int p = *px++;
putc(bi->bmiColors[p].rgbRed, fp);
putc(bi->bmiColors[p].rgbGreen, fp);
putc(bi->bmiColors[p].rgbBlue, fp);
}
}
if (fp != stdout) {
fclose(fp);
}
}
/* PARSEDEG -- Parse latitude and longitude degree specifications. */
static double parseDeg(const char *dec)
{
if ((strchr(dec, 'd') != NULL) ||
(strchr(dec, 'D') != NULL) ||
(strchr(dec, ' ') != NULL) ||
(strchr(dec, '°') != NULL)) {
double dd = 0, mm = 0, ss = 0;
char c1, c2;
sscanf(dec, "%lf%c%lf%c%lf", &dd, &c1, &mm, &c2, &ss);
return sgn(dd) * (abs(dd) + (mm / 60) + (ss / 3600));
} else {
return atof(dec);
}
}
/* LOCFILE -- Expand local file name to full directory name from
which we're executing this program. */
static char *locfile(const char *basename)
{
static char pn[PATH_MAX];
strcpy(pn, progpath);
if (strrchr(pn, '/') != NULL) {
*(strrchr(pn, '/') + 1) = EOS;
} else {
strcpy(pn, "./");
}
if (directCall || embeddedImage) {
strcat(pn, "../cgi-executables/");
}
strcat(pn, basename);
return pn;
}
/* EDDEG -- Edit degrees and minutes. */
static char *eddeg(const double ds)
{
static char buf[2][20];
static int n = 0;
long a = (long) ((abs(ds) + 0.00001) * 3600);
int d, m, s;
n = (n + 1) % 2;
d = a / 3600;
m = (a / 60) % 60;
s = (a % 60);
sprintf(buf[n], "%s%d\260", ds < 0 ? "-" : "", d);
if (m > 0 || s > 0) {
sprintf(buf[n] + strlen(buf[n]), "%d'", m);
if (s > 0) {
sprintf(buf[n] + strlen(buf[n]), "%d"", s);
}
}
return buf[n];
}
/* EDLAT -- Edit latitude. */
static char *edlat(const double lat)
{
static char slat[30];
double ulat = fixangle(rtd(lat));
if (ulat >= 180) {
ulat = -(360 - ulat);
}
sprintf(slat, "%s%s", eddeg(abs(ulat)), ulat < 0 ? "S" : "N");
return slat;
}
/* EDLON -- Edit longitude. */
static char *edlon(const double lon)
{
static char slon[30];
double ulon = fixangle(rtd(lon));
if (ulon >= 180) {
ulon = -(360 - ulon);
}
sprintf(slon, "%s%s", eddeg(abs(ulon)), ulon < 0 ? "E" : "W");
return slon;
}
/* EDVPOS -- Edit viewing position and return string. */
static char *edvpos(const double tlat, const double tlon)
{
static char vpos[80];
sprintf(vpos, "%s %s",
edlat(tlat), edlon(tlon));
return vpos;
}
/* PSRVECTOR -- Draw a vector into the current image array. */
#define PixA(x, y) ((tpixel + ((x) + (((DWORD) (tmhm1 - (y))) * (tmrowlen)))))
#define Opix(x, y) ((opixel + ((x) + (((DWORD) ((adiheight - 1) - (y))) * (rowlen)))))
static int rowlen, tmrowlen, adiheight;
static unsigned char *opixel;
static unsigned char *tpixel;
void psrvector(const unsigned int fx, const unsigned int fy,
const unsigned int tx, const unsigned int ty, const unsigned int colour)
{
unsigned int cfx, cfy, ctx, cty;
int x, y, m, f, xinc, loopcnt, dx, dy;
#ifdef FlipY
fy = (iarray[1] - 1) - fy;
ty = (iarray[1] - 1) - ty;
#endif
if (fy < ty) {
cfx = fx;
cfy = fy;
ctx = tx;
cty = ty;
} else {
cfx = tx;
cfy = ty;
ctx = fx;
cty = fy;
}
dy = cty - cfy;
dx = ctx - cfx;
if (dx < 0) {
xinc = -1;
dx = -dx;
} else {
xinc = 1;
}
x = cfx;
y = cfy;
if (dx > dy) {
f = dy - dx / 2;
m = dx - dy;
loopcnt = dx + 1;
while (loopcnt--) {
*Opix(x, y) = colour;
x += xinc;
if (f > 0) {
f -= m;
y++;
} else {
f += dy;
}
}
} else {
f = -(dy / 2 - dx);
m = dy - dx;
loopcnt = dy + 1;
while (loopcnt--) {
*Opix(x, y) = colour;
y++;
if (f > 0) {
x += xinc;
f -= m;
} else {
f += dx;
}
}
}
}
static double orrLat = 90; /* Orrery latitude */
static double orrLong = 0; /* Orrery longitude */
static int orrInner = TRUE; /* Orrery inner solar system only ? */
static int orrScale = 2; /* Distance scale */
static double earthlong, earthrv; /* Heliocentric coordinates of the Earth */
static double minperh = 0.28; /* Minimum perihelion (usually Mercury, but possible ast./comet) */
/* The fastest-changing aspect of Solar system orbits is the perihelion shift in
the orbit of Mercury, and that amounts to only about 1.5 degrees per century (
0.559974 +/= 0.000041 arc seconds per century, to be precise). So, we only
bother to recalculate the orbits when they're a century out of date. */
/* Sidereal rotation period of planets in days */
static double rotationPeriod[] = {
87.969,
224.701,
365.256,
686.980,
4332.589,
10759.22,
30685.4,
60189.0,
90465.0,
0.0 /* Asteroid period calculated from orbital elements */
};
/* SETPROJECT -- Define projection matrix for a given viewpoint. */
static void setproject(const double orrlat, const double orrlong, const double stereo)
{
matident(cT);
if (stereo != 0.0) {
rot(dtr(stereo), Y); /* Stereo ocular offset */
}
rot(dtr(180.0 - (orrlat - 90.0)), X); /* Tilt to latitide */
rot(dtr(90.0 - orrlong), Z); /* Spin to longitude */
}
/* PROJECT -- Project a longitude, latitude, and radius vector into the
display space. */
static void project(const double along, const double alat, const double arv, const int cr,
const int plno, const int nplan, const double diau,
int *u, int *v, int *z)
{
double lrv, acoslat, vx, vy, vz;
vector p, uvw;
#define MercOffset 16
if (orrScale == 0) {
/* Log projection */
lrv = log10(1.0 + 9.0 * ((arv - minperh) / (diau)));
if (lrv < 0) {
lrv = 0;
} else {
lrv = (int) (MercOffset + (cr - MercOffset) * lrv);
}
} else if (orrScale == 1) {
/* Real projection */
lrv = (int) (MercOffset + (cr - MercOffset) * (arv / diau));
} else {/* orrScale == 2 */
/* Equal projection */
lrv = (int) (MercOffset + (cr - MercOffset) * ((plno == 10 ? 4.5 : ((double) plno)) / ((double) nplan)));
}
/* Get rectangular co-ordinates of point. */
acoslat = abs(cos(dtr(alat)));
vecget(p, lrv * sin(dtr(along)) * acoslat,
lrv * cos(dtr(along)) * acoslat,
-lrv * sin(dtr(alat)));
/* Transform by viewing matrix */
vecxmat(uvw, p, cT);
/* Store projected co-ordinates. */
vecput(&vx, &vy, &vz, uvw);
/* Return screen co-ordinates (relative to the centre of
the screen. The Z co-ordinate represents the distance
of the point behind (if negative) or above (if positive)
the screen. This is used to sort the order in which the
planet icons are drawn to that any obscurations are done
in the correct order. */
*u = (int) vx;
*v = (int) vy;
*z = (int) (vz * 2000);
}
/* ORRERY -- Generate Orrery view. */
static void orrery(const int isize, const int tmwidth, const int tmheight,
const int diwidth, const int diheight, const double faketime,
const int stereo)
{
int i, x, y, ix, iy, o = -1, t, tmhm1, cr, cpx = 0, cpy = 0,
nplan = aTracked ? 10 : (orrInner ? 4 : 9),
siwidth = diwidth, sbias, pass,
planX[11], planY[11], planZ[11], planord[11];
double diau = orrInner ? (aTracked ?
max(1.7, min(5.5, ((ast_info.SemiMajorAU * (1.0 + ast_info.Eccentricity)) * 1.1)))
: 1.7) :
(aTracked ?
max(50, min(150, ((ast_info.SemiMajorAU * (1.0 + ast_info.Eccentricity)) * 1.1)))
: 50),
eplan = aTracked ? (orrInner ? 5 : 9) : (orrInner ? 4 : 9);
double dialim = diau * 2;
BITMAPINFOHEADER *bh, *obmap;
WORD ncol, offbits;
double orrlat = orrLat, orrlong = orrLong;
#ifdef SaveOrbits
FILE *fp = fopen("/tmp/orbits.c", "w");
#endif
/*
if (aTracked) {
fprintf(stderr, "Diau = %.2f, acomp = %.2f\n", diau,
(ast_info.SemiMajorAU * (1.0 + ast_info.Eccentricity)) * 1.1);
}
*/
if (stereo) {
siwidth += diwidth + stereoSepPixels;
}
adiheight = diheight;
bh = &bmi.bmiHeader;
ncol = (WORD) (bh->biClrUsed == 0 ? (1L << bh->biBitCount) : bh->biClrUsed);
offbits = (WORD) (bh->biSize + ncol * sizeof(RGBQUAD));
tpixel = bmBits;
tmrowlen = ((tmwidth + (sizeof(LONG) - 1)) / sizeof(LONG)) * sizeof(LONG);
tmhm1 = tmheight - 1;
if (tmbmap != NULL) {
free(tmbmap);
tmbmap = NULL;
}
rowlen = ((siwidth + (sizeof(LONG) - 1)) / sizeof(LONG)) * sizeof(LONG);
obmap = (BITMAPINFOHEADER *) malloc(offbits + diheight * rowlen);
memset((char *) obmap, 0, offbits + diheight * rowlen);
memcpy((char *) obmap, (char *) bh, offbits);
obmap->biWidth = siwidth;
obmap->biHeight = diheight;
obmap->biSizeImage = ((DWORD) diheight) * rowlen;
opixel = ((unsigned char *) (obmap)) + offbits;
tmbmap = (struct bmih *) obmap;
tmbits = (unsigned char *) opixel;
if ((skyColourScheme == SCS_BLACK_ON_WHITE) ||
(skyColourScheme == SCS_COLOUR_ON_WHITE)) {
for (y = 0; y < diheight; y++) {
memset(Opix(0, y), cWhite, stereo ? siwidth : diwidth);
}
}
/* If we're generating a stereo pair, draw the border between
the two images. This provides a reference which makes them
easier to fuse. */
if (stereo) {
for (y = 0; y < diheight; y++) {
for (x = diwidth; x < diwidth + stereoSepPixels; x++) {
*Opix(x, y) = cBorder;
}
}
}
for (pass = 0; pass < (stereo ? 2 : 1); pass++) {
sbias = (pass == 0) ? 0 : (diwidth + stereoSepPixels);
setproject(orrlat, orrlong,
stereo ? (stereoOcular * (pass == 0 ? 1 : -1)) : 0.0);
{
double plong, plat, prv, sjd, timestep;
int cx, cy, vx, vy, vz, fvx = 0, fvy = 0, ftf = TRUE, nsteps, colour = 0;
char spi[sizeof(planet_info)];
#define inLim(x) ((x) >= 0 && (x) < diheight)
#define MoveTo(x, y) cpx = (x), cpy = (y)
#define LineTo(x, y) if (inLim(cpx) && inLim(cpy) && inLim(x) && inLim(y)) { \
psrvector(sbias + cpx, cpy, \
sbias + (x), (y), colour); \
} cpx = (x), cpy = (y)
#define nSteps 64
/* Constrain use of canned orbits to decent interval around J2000.0 */
if (abs((faketime - J2000) / JulianCentury) <= 5) {
o = 0;
}
cx = diwidth / 2;
cy = diheight / 2;
cr = (min(diwidth, diheight) / 2) - (MercOffset / 2);
for (i = 1; i <= min(nplan, 9); i++) {
ftf = TRUE;
#ifdef SaveOrbits
fprintf(fp, "\n /* Planet %d */ \n\n", i);
#endif
/* Test for zero radius vector which indicates planetary theory not
known for this date (Pluto only). If that's the case, skip
drawing the planet. */
if ((orrInner && (i > 4)) ||
((i != 3) && (planet_info[i].hrv == 0.0))) {
continue;
}
plutoPrecise = FALSE; /* Allow suspect values to close Pluto orbit */
memcpy(spi, (char *) planet_info, sizeof spi);
nsteps = nSteps;
/* Draw the orbital path for the planet by stepping forward in
time one revolution period. (N.b. for the planets whose
orbits are almost perfectly circular, it would be *enormously*
faster to just draw circles of the mean radius. But for the
moment we do it the hard way for all of them. This may be
slow, but it's *right*.) */
timestep = rotationPeriod[i - 1] / nsteps;
for (sjd = faketime; sjd <= faketime + rotationPeriod[i - 1]; sjd += timestep) {
if (o >= 0 && i <= 9) {
plong = (((double) currentOrbits[o].ulong) * 360.0) / 65535.0;
plat = (((double) currentOrbits[o].ulat) * 180.0) / 32767.0;
prv = (((double) currentOrbits[o++].urv) * 50.0) / 65535.0;
} else {
if (i == 3) {
double sunra, sundec;
sunpos(sjd, TRUE, &sunra, &sundec, &prv, &plong);
plong += 180;
plat = 0;
} else {
planets(sjd, 1 << i);
plong = planet_info[i].hlong;
plat = planet_info[i].hlat;
prv = planet_info[i].hrv;
}
}
#ifdef SaveOrbits
if (fp != NULL && i <= 9) {
unsigned int ulong, urv;
int ulat;
ulong = (unsigned int) ((fixangle(plong) * 65535.0) / 360.0);
ulat = (int) ((plat * 32767.0) / 180.0);
urv = (unsigned int) ((prv * 65535.0) / 50.0);
fprintf(fp, " { %6u, %6d, %6u },\n", ulong, ulat, urv);
}
#endif
if (prv > 0) {
project(plong, plat, prv, cr, i, eplan, diau,
&vx, &vy, &vz);
vx += cx;
vy = cy - vy;
if (ftf) {
MoveTo(fvx = vx, fvy = vy);
ftf = FALSE;
} else {
if (plat < 0) {
colour = Cscheme(cGreen, cGreen, cLtGrey, cGrey, cDkRed);
} else {
colour = Cscheme(cBlue, cBlue, cBlack, cWhite, cRed);
}
LineTo(vx, vy);
}
}
}
if (!ftf) {
/* Draw line to guarantee orbit is closed */
LineTo(fvx, fvy);
}
plutoPrecise = TRUE;
memcpy((char *) planet_info, spi, sizeof spi);
}
/* If an asteroid or comet is being tracked, we have to handle it
specially. First of all, we can never use the canned planetary
orbits to avoid calculation. Second, unlike planetary orbits,
asteroidal and cometary orbits can have large eccentricity,
which means that we can't use a naive approximation to the
number of steps needed to draw the orbit smoothly, but must
instead take into account the orbital velocity at various
locations on the orbit and adjust the time step accordingly.
Finally, the orbit may not be closed either because it is
genuinely parabolic or hyperbolic, or merely such an eccentric
ellipse that the curve is truncated at the edge of the display. */
if (aTracked) {
double eSmA, lprv, aprv = 0;
int arm, avx = 0, avy = 0, lvx = 0, lvy = 0;
/* int aclosed = FALSE; */
i = 10;
memcpy(spi, (char *) planet_info, sizeof spi);
/* Because the orbit may not be closed, we start at the
perihelion date and plot into the future until we
either arrive at the aphelion or we run into the
cutoff based on the scale of the display. Then we go
back to the perihelion and plot backward until we
encounter the equivalent event on the other branch. */
eSmA = ast_info.Eccentricity < 1.0 ? ast_info.SemiMajorAU :
40;
for (arm = -1; arm <= 1; arm += 2) {
ftf = TRUE; /* No vector on first point */
sjd = ast_info.PeriDate;
lprv = 0;
while (TRUE) {
double ov = 0.5;
planets(sjd, 1 << 10);
plong = planet_info[i].hlong;
plat = planet_info[i].hlat;
prv = planet_info[i].hrv;
if (prv > 0) {
lvx = vx;
lvy = vy;
project(plong, plat, prv, cr, 10, eplan, diau,
&vx, &vy, &vz);
vx += cx;
vy = cy - vy;
if (ftf) {
MoveTo(vx, vy);
ftf = FALSE;
} else {
if (plat < 0) {
colour = cGreen;
} else {
colour = cBlue;
}
LineTo(vx, vy);
}
}
#ifdef AST_DEBUG
else {
fprintf(stderr, "Prv = %.4f sjd = %.4f\n", prv, sjd);
break;
}
if (!finite(prv)) {
fprintf(stderr, "Not finite Prv = %.4f sjd = %.4f\n", prv, sjd);
break;
}
#endif
if ((prv == 0) || (prv > dialim) || (prv < lprv)) {
if (prv > 0) {
if (arm > 0) {
/* If orbit closed, connect ends of orbit */
if (aprv <= lprv && prv <= lprv) {
LineTo(avx, avy);
LineTo(avx, avy); /* Set last pixel */
}
} else {
/* aclosed = prv < lprv; */
aprv = prv;
avx = lvx;
avy = lvy;
}
}
break;
}
lprv = prv;
/* Orbital velocity. Note that since we plot both branches of
the orbit outward from the perihelion, we always overestimate
the average velocity for the segment, guaranteeing that the
curve will be smooth at the expense of a few extra segments.
If the orbit of the body extends into the outer solar system,
the calculation of its orbital velocity may result in taking
the square root of a negative number due to the hokey way we
we estimate it. If this situation arises, we simply use the
last valid orbital velocity computed (which will result in
more segments than needed) for the remaining part of the orbit. */
#ifdef AST_DEBUG
if (prv == 0) {
fprintf(stderr, "Doooh! OV calc Prv = %.4f sjd = %.4f\n", prv, sjd);
break;
}
#endif
if ((1 / (2 * eSmA)) < (1 / prv)) {
/* fprintf(stderr, "Limiting OV, was %.4f at prv %.4f\n", ov, prv); */
ov = 42.1219 * sqrt((1 / prv) - (1 / (2 * eSmA)));
}
/* Time step to move 0.1 AU in inner solar system, correspondingly
more in full system. */
sjd += arm * (orrInner ? 1 : (prv < 1.0 ? 1 : max(1, (eSmA / 5)))) *
(((AstronomicalUnit / ov) / (24.0 * 60 * 60)) / 10);
#ifdef AST_DEBUG
if (!finite(sjd)) {
fprintf(stderr, "Doooh! SJD increment arm = %d prv = %.4f ov = %.4f\n", arm, prv, ov);
break;
}
#endif
}
}
memcpy((char *) planet_info, spi, sizeof spi);
}
}
/* Plug in well-known values for the heliocentric
position of the Sun. We only do this for the
first image of a stereo pair, since the existing
value will be re-used for the second. */
if (pass == 0) {
planord[0] = planX[0] = planY[0] = planZ[0] = 0;
}
/* Calculate transformed plot positions for the planet
icons. */
for (i = 1; i <= 10; i++) {
int o;
/* On the first image of a stereo pair, fill the table
of planets in linear order. For the second image,
where we want to re-use the sort order of the
original image (see comment below for why), find the
planet in the order table and use that slot. Note
that we have to include a bail-out in case this
planet wasn't rendered in the selected view. */
if (pass == 0) {
o = i;
} else {
for (o = 0; o <= 10; o++) {
if (planord[o] == i) {
break;
}
}
if (o > 10) {
continue;
}
}
planord[o] = i;
if (i <= nplan) {
double plong, plat, prv;
if (i == 3) {
plong = earthlong;
plat = 0;
prv = earthrv;
} else {
plong = planet_info[i].hlong;
plat = planet_info[i].hlat;
prv = planet_info[i].hrv;
}
if ((orrInner && (i > 4 && i <= 9)) || (prv == 0.0)) {
planord[o] = planX[o] = planY[o] = planZ[o] = -1;
continue; /* Planetary theory for Pluto invalid at this time */
}
project(plong, plat, prv, cr, i, eplan, diau,
&planX[o], &planY[o], &planZ[o]);
} else {
planord[o] = planX[o] = planY[o] = planZ[o] = -1;
}
}
/* Sort planets by projected distance. We do this only
for the first image of a stereo pair. The second stereo
image is always drawn with the same order of planets
as the first to avoid peekaboo planets when the
slight difference in viewpoint changes the Z-depth. */
if (pass == 0) {
for (i = 0; i <= 10; i++) {
int j;
for (j = i + 1; j <= 10; j++) {
if (planZ[j] < planZ[i]) {
int t;
#define Swop(p) t = p[i]; p[i] = p[j]; p[j] = t
Swop(planX);
Swop(planY);
Swop(planZ);
Swop(planord);
#undef Swop
}
}
}
}
/* Plot the planet icons in order of projected distance
from the view: most distant to nearest. This performs
a "painter's algorithm" rendering which ensures all
obscurations are done in the right order. */
for (i = 0; i <= 10; i++) {
if (planord[i] >= 0) {
ix = diwidth / 2 + planX[i];
iy = diheight / 2 - planY[i];
t = planord[i] * 32; /* Select proper icon */
if (planord[i] == 10 && ast_info.cometary) {
t += 32;
}
for (x = 0; x < 32; x++) {
int ax = (ix - 16) + x;
if (ax >= 0 && ax < diwidth) {
for (y = 0; y < 32; y++) {
int ay = (iy - 16) + y;
if (ay >= 0 && ay < diheight) {
unsigned char p = *PixA(x + t, y);
if (p != transparent) {
*Opix(sbias + ax, ay) = p;
}
}
}
}
}
}
}
}
#ifdef SaveOrbits
fclose(fp);
#endif
}
/* WRITE_CACHE_IMAGE_WARNING -- Include warning in HTML source to deter
users from linking to images in the cache. */
#ifdef CACHE_WARNING
static void write_cache_image_warning(FILE *fp)
{
#define W(x) fprintf(fp, "%s\n", x)
W("<!--");
W(" WARNING!!!");
W("");
W("You may be looking at the source for this Solar System Live result");
W("document in order to figure out where the image file comes from");
W("so you can link to it from your site. DON'T DO THAT!");
W("");
W("The image files generated by Solar System Live are stored in");
W("a temporary \"cache\" directory, and are kept only long enough");
W("to guarantee that users will be able to download them to their");
W("browsers. After 10 minutes or so, images are automatically deleted.");
W("So, if you link to one of these images, it may work if you try it");
W("right away, but the link is sure to break before long, leaving an");
W("ugly broken image on your page.");
W("");
W("Rather than linking to the image, load it in your browser and then");
W("save a local copy on your site and reference it from there.");
W("If you want to link to an image which will be updated every time");
W("users request it, create a custom request with:");
W(" http://www.fourmilab.ch/solar/custom.html");
W("then copy the URL it generates into a link on your page.");
W("");
W("Thank you for your understanding and cooperation.");
W("-->");
#undef W
}
#endif
/* WRITEPOST -- Write balance of HTML file, including controls
which permit regenerating the view. */
static void writepost(FILE *fp, const double vlat, const double vlon, const double alt,
const char *bmpname, const int innersystem, const int orbittype, const int tmsize, const double jt)
{
int i, j;
double lat, lon;
static const char checked[] = " checked=\"checked\"";
lon = fixangle(rtd(vlon));
if (lon >= 180) {
lon = -(360 - lon);
}
lat = fixangle(rtd(vlat));
if (lat >= 180) {
lat = -(360 - lat);
}
#define W(x) fprintf(fp, "%s\n", x)
#define L(x, y) fprintf(fp, "<a href=\"/solar/help/%s.html\"><b>%s:</b></a>", y, x);
#ifdef TESTMODE
/* We always use the "get" method in TESTMODE so that the arguments are visible
and can be edited for debugging. */
W("<form method=\"get\" name=\"request\" action=\"/cgi-bin/Solar.t\">");
#else
#ifdef SAVESET
fprintf(fp, "<form method=\"%s\" name=\"request\" action=\"/cgi-bin/Solar\">\n", bookmark ? "get" : "post");
if (bookmark) {
W("<input name=\"got\" type=\"hidden\" value=\"-xh\" />\n");
}
#else
fprintf(fp, "<form method=\"post\" name=\"request\" action=\"/cgi-bin/Solar\">\n");
#endif
#endif
W("<center>");
W("<p />");
L("Time", "timedate");
fprintf(fp, "<input type=\"radio\" name=\"date\" onclick=\"0\" value=\"0\"%s /> <a href=\"/solar/help/timedate.html#Now\">Now</a>\n",
dto == 0 ? checked : "");
fprintf(fp, "<input type=\"radio\" name=\"date\" onclick=\"0\" value=\"1\"%s /> <a href=\"/solar/help/timedate.html#UTC\">UTC:</a> ",
dto == 1 ? checked : "");
fprintf(fp, "<input type=\"text\" name=\"utc\" value=\"%s\" size=\"20\" onchange=\"document.request.date[1].checked=true;\" />\n",
dtutc);
fprintf(fp, "<input type=\"radio\" name=\"date\" onclick=\"0\" value=\"2\"%s /> <a href=\"/solar/help/timedate.html#Julian\">Julian:</a> ",
dto == 2 ? checked : "");
fprintf(fp, "<input type=\"text\" name=\"jd\" value=\"%.5f\" size=\"15\" onchange=\"document.request.date[2].checked=true;\" />\n",
jt);
W("<br />");
#ifdef SAVESET
fprintf(fp, " <input type=\"checkbox\" name=\"bmark\" value=\"-xs\"%s /> ",
(bookmark && !gotten) ? "checked" : "");
fprintf(fp, "<a href=\"/solar/help/saveset.html\"><b>Save settings</b></a>");
#endif
W( "<input type=\"submit\" value=\"Update\" />");
W("<br />");
L("Show", "icons");
fprintf(fp, "<label><input type=\"radio\" name=\"img\" value=\"-k0\"%s /> Icons</label>\n",
!useimage ? checked : "");
fprintf(fp, "<label><input type=\"radio\" name=\"img\" value=\"-k1\"%s /> Images</label>\n",
useimage ? checked : "");
W("<br />");
L("Display", "in-out");
fprintf(fp, "<label><input type=\"radio\" name=\"sys\" value=\"-Sf\"%s /> Full system</label>\n",
!innersystem ? checked : "");
fprintf(fp, "<label><input type=\"radio\" name=\"sys\" value=\"-Si\"%s /> Inner system</label>\n",
innersystem ? checked : "");