1 /* given a Now and an Obj with the object definition portion filled in,
2 * fill in the sky position (s_*) portions.
3 * calculation of positional coordinates reworked by
4 * Michael Sternberg <sternberg@physik.tu-chemnitz.de>
5 * 3/11/98: deflect was using op->s_hlong before being set in cir_pos().
6 * 4/19/98: just edit a comment
7 * 11/22/21: un-swapped arguments "rsn, lsn" in both calls to deflect().
8 */
9
10 #include <stdio.h>
11 #include <math.h>
12 #include <stdlib.h>
13
14 #include "astro.h"
15 #include "preferences.h"
16
17
18 static int obj_planet (Now *np, Obj *op);
19 static int obj_binary (Now *np, Obj *op);
20 static int obj_2binary (Now *np, Obj *op);
21 static int obj_fixed (Now *np, Obj *op);
22 static int obj_elliptical (Now *np, Obj *op);
23 static int obj_hyperbolic (Now *np, Obj *op);
24 static int obj_parabolic (Now *np, Obj *op);
25 static int sun_cir (Now *np, Obj *op);
26 static int moon_cir (Now *np, Obj *op);
27 static double solveKepler (double M, double e);
28 static void binaryStarOrbit (double t, double T, double e, double o, double O,
29 double i, double a, double P, double *thetap, double *rhop);
30 static void cir_sky (Now *np, double lpd, double psi, double rp, double *rho,
31 double lam, double bet, double lsn, double rsn, Obj *op);
32 static void cir_pos (Now *np, double bet, double lam, double *rho, Obj *op);
33 static void elongation (double lam, double bet, double lsn, double *el);
34 static void deflect (double mjd1, double lpd, double psi, double rsn,
35 double lsn, double rho, double *ra, double *dec);
36 static double h_albsize (double H);
37
38 /* given a Now and an Obj, fill in the approprirate s_* fields within Obj.
39 * return 0 if all ok, else -1.
40 */
41 int
obj_cir(Now * np,Obj * op)42 obj_cir (Now *np, Obj *op)
43 {
44 op->o_flags &= ~NOCIRCUM;
45 switch (op->o_type) {
46 case BINARYSTAR: return (obj_binary (np, op));
47 case FIXED: return (obj_fixed (np, op));
48 case ELLIPTICAL: return (obj_elliptical (np, op));
49 case HYPERBOLIC: return (obj_hyperbolic (np, op));
50 case PARABOLIC: return (obj_parabolic (np, op));
51 case EARTHSAT: return (obj_earthsat (np, op));
52 case PLANET: return (obj_planet (np, op));
53 default:
54 printf ("obj_cir() called with type %d %s\n", op->o_type, op->o_name);
55 abort();
56 return (-1); /* just for lint */
57 }
58 }
59
60 static int
obj_planet(Now * np,Obj * op)61 obj_planet (Now *np, Obj *op)
62 {
63 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/
64 double lpd, psi; /* heliocentric ecliptic long and lat */
65 double rp; /* dist from sun */
66 double rho; /* dist from earth */
67 double lam, bet; /* geocentric ecliptic long and lat */
68 double dia, mag; /* angular diameter at 1 AU and magnitude */
69 PLCode p;
70
71 /* validate code and check for a few special cases */
72 p = op->pl_code;
73 if (p == SUN)
74 return (sun_cir (np, op));
75 if (p == MOON)
76 return (moon_cir (np, op));
77 if (op->pl_moon != X_PLANET)
78 return (plmoon_cir (np, op));
79 if (p < 0 || p > MOON) {
80 printf ("unknown planet code: %d\n", p);
81 abort();
82 }
83
84 /* planet itself */
85
86 /* find solar ecliptical longitude and distance to sun from earth */
87 sunpos (mjed, &lsn, &rsn, 0);
88
89 /* find helio long/lat; sun/planet and earth/planet dist; ecliptic
90 * long/lat; diameter and mag.
91 */
92 plans(mjed, p, &lpd, &psi, &rp, &rho, &lam, &bet, &dia, &mag);
93
94 /* fill in all of op->s_* stuff except s_size and s_mag */
95 cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);
96
97 /* set magnitude and angular size */
98 set_smag (op, mag);
99 op->s_size = (float)(dia/rho);
100
101 return (0);
102 }
103
104 static int
obj_binary(Now * np,Obj * op)105 obj_binary (Now *np, Obj *op)
106 {
107 /* always compute circumstances of primary */
108 if (obj_fixed (np, op) < 0)
109 return (0);
110
111 /* compute secondary only if requested, and always reset request flag */
112 if (!op->b_2compute)
113 return (0);
114 op->b_2compute = 0;
115 return (obj_2binary (np, op));
116 }
117
118 /* compute position of secondary component of a BINARYSTAR */
119 static int
obj_2binary(Now * np,Obj * op)120 obj_2binary (Now *np, Obj *op)
121 {
122 if (op->b_nbp > 0) {
123 /* we just have discrete pa/sep, project each from primary */
124 int i;
125 for (i = 0; i < op->b_nbp; i++) {
126 BinPos *bp = &op->b_bp[i];
127 bp->bp_dec = op->s_dec + bp->bp_sep*cos(bp->bp_pa);
128 bp->bp_ra = op->s_ra + bp->bp_sep*sin(bp->bp_pa)/cos(op->s_dec);
129 }
130 } else {
131 BinOrbit *bp = &op->b_bo;
132 double t, theta, rho;
133
134 mjd_year (mjd, &t);
135 binaryStarOrbit (t, bp->bo_T, bp->bo_e, bp->bo_o, bp->bo_O,
136 bp->bo_i, bp->bo_a, bp->bo_P, &theta, &rho);
137 bp->bo_pa = (float)theta;
138 bp->bo_sep = (float)rho;
139 rho = degrad(rho/3600.); /* arc secs to rads */
140 bp->bo_dec = op->s_dec + rho*cos(theta);
141 bp->bo_ra = op->s_ra + rho*sin(theta)/cos(op->s_dec);
142 }
143
144 return (0);
145 }
146
147 /* from W. M. Smart */
148 static void
binaryStarOrbit(double t,double T,double e,double o,double O,double i,double a,double P,double * thetap,double * rhop)149 binaryStarOrbit (
150 double t, /* desired ephemeris epoch, year */
151 double T, /* epoch of periastron, year */
152 double e, /* eccentricity */
153 double o, /* argument of periastron, degrees */
154 double O, /* ascending node, degrees */
155 double i, /* inclination, degrees */
156 double a, /* semi major axis, arcsecs */
157 double P, /* period, years */
158 double *thetap, /* position angle, rads E of N */
159 double *rhop) /* separation, arcsecs */
160 {
161 double M, E, cosE, nu, cosnu, r, rho, theta;
162
163 /* find mean anomaly, insure 0..2*PI */
164 M = 2*PI/P*(t-T);
165 range (&M, 2*PI);
166
167 /* solve for eccentric anomaly */
168 E = solveKepler (M, e);
169 cosE = cos(E);
170
171 /* find true anomaly and separation */
172 cosnu = (cosE - e)/(1.0 - e*cosE);
173 r = a*(1.0 - e*e)/(1.0 + e*cosnu);
174 nu = acos(cosnu);
175 if (E > PI)
176 nu = -nu;
177
178 /* project onto sky */
179 theta = atan(tan(nu+degrad(o))*cos(degrad(i))) + degrad(O);
180 rho = r*cos(nu+degrad(o))/cos(theta-degrad(O));
181 if (rho < 0) {
182 theta += PI;
183 rho = -rho;
184 }
185 range (&theta, 2*PI);
186
187 *thetap = theta;
188 *rhop = rho;
189 }
190
191 /* solve kepler equation using Newton-Raphson search.
192 * Charles and Tatum have shown it always converges starting with PI.
193 */
194 static double
solveKepler(double M,double e)195 solveKepler (double M, double e)
196 {
197 double E, Eprime = PI;
198
199 do {
200 double cosE = cos(Eprime);
201 E = Eprime;
202 Eprime = (M - e*(E*cosE - sin(E)))/(1.0 - e*cosE);
203 } while (fabs(E-Eprime) > 1e-7);
204
205 return (Eprime);
206 }
207
208 static int
obj_fixed(Now * np,Obj * op)209 obj_fixed (Now *np, Obj *op)
210 {
211 double lsn, rsn; /* true geoc lng of sun, dist from sn to earth*/
212 double lam, bet; /* geocentric ecliptic long and lat */
213 double ha; /* local hour angle */
214 double el; /* elongation */
215 double alt, az; /* current alt, az */
216 double ra, dec; /* ra and dec at equinox of date */
217 double rpm, dpm; /* astrometric ra and dec with PM to now */
218 double lst;
219
220 /* on the assumption that the user will stick with their chosen display
221 * epoch for a while, we move the defining values to match and avoid
222 * precession for every call until it is changed again.
223 * N.B. only compare and store jd's to lowest precission (f_epoch).
224 * N.B. maintaining J2k ref (which is arbitrary) helps avoid accum err
225 */
226 if (0 /* disabled in PyEphem */
227 && epoch != EOD && epoch != op->f_epoch) {
228 double pr = op->f_RA, pd = op->f_dec, fe = epoch;
229 /* first bring back to 2k */
230 precess (op->f_epoch, J2000, &pr, &pd);
231 pr += op->f_pmRA*(J2000-op->f_epoch);
232 pd += op->f_pmdec*(J2000-op->f_epoch);
233 /* then to epoch */
234 pr += op->f_pmRA*(fe-J2000);
235 pd += op->f_pmdec*(fe-J2000);
236 precess (J2000, fe, &pr, &pd);
237 op->f_RA = pr;
238 op->f_dec = pd;
239 op->f_epoch = fe;
240 }
241
242 /* apply proper motion .. assume pm epoch reference equals equinox */
243 rpm = op->f_RA + op->f_pmRA*(mjd-op->f_epoch);
244 dpm = op->f_dec + op->f_pmdec*(mjd-op->f_epoch);
245
246 /* set ra/dec to astrometric @ equinox of date */
247 ra = rpm;
248 dec = dpm;
249 if (op->f_epoch != mjed)
250 precess (op->f_epoch, mjed, &ra, &dec);
251
252 /* compute astrometric @ requested equinox */
253 op->s_astrora = rpm;
254 op->s_astrodec = dpm;
255 if (op->f_epoch != epoch)
256 precess (op->f_epoch, epoch, &op->s_astrora, &op->s_astrodec);
257
258 /* convert equatoreal ra/dec to mean geocentric ecliptic lat/long */
259 eq_ecl (mjed, ra, dec, &bet, &lam);
260
261 /* find solar ecliptical long.(mean equinox) and distance from earth */
262 sunpos (mjed, &lsn, &rsn, NULL);
263
264 /* allow for relativistic light bending near the sun */
265 deflect (mjed, lam, bet, rsn, lsn, 1e10, &ra, &dec);
266
267 /* TODO: correction for annual parallax would go here */
268
269 /* correct EOD equatoreal for nutation/aberation to form apparent
270 * geocentric
271 */
272 nut_eq(mjed, &ra, &dec);
273 ab_eq(mjed, lsn, &ra, &dec);
274 op->s_gaera = ra;
275 op->s_gaedec = dec;
276
277 /* set s_ra/dec -- apparent */
278 op->s_ra = ra;
279 op->s_dec = dec;
280
281 /* compute elongation from ecliptic long/lat and sun geocentric long */
282 elongation (lam, bet, lsn, &el);
283 el = raddeg(el);
284 op->s_elong = (float)el;
285
286 /* these are really the same fields ...
287 op->s_mag = op->f_mag;
288 op->s_size = op->f_size;
289 */
290
291 /* alt, az: correct for refraction; use eod ra/dec. */
292 now_lst (np, &lst);
293 ha = hrrad(lst) - ra;
294 hadec_aa (lat, ha, dec, &alt, &az);
295 refract (pressure, temp, alt, &alt);
296 op->s_ha = ha;
297 op->s_alt = alt;
298 op->s_az = az;
299
300 return (0);
301 }
302
303 /* compute sky circumstances of an object in heliocentric elliptic orbit at *np.
304 */
305 static int
obj_elliptical(Now * np,Obj * op)306 obj_elliptical (Now *np, Obj *op)
307 {
308 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/
309 double dt; /* light travel time to object */
310 double lg; /* helio long of earth */
311 double nu; /* true anomaly */
312 double rp=0; /* distance from the sun */
313 double lo, slo, clo; /* angle from ascending node */
314 double inc; /* inclination */
315 double psi=0; /* heliocentric latitude */
316 double spsi=0, cpsi=0; /* trig of heliocentric latitude */
317 double lpd; /* heliocentric longitude */
318 double rho=0; /* distance from the Earth */
319 double om; /* arg of perihelion */
320 double Om; /* long of ascending node. */
321 double lam; /* geocentric ecliptic longitude */
322 double bet; /* geocentric ecliptic latitude */
323 double ll=0, sll, cll; /* helio angle between object and earth */
324 double mag; /* magnitude */
325 double e_n; /* mean daily motion */
326 double tp; /* time from perihelion (days) */
327 double rpd=0;
328 double y;
329 int pass;
330
331 /* find location of earth from sun now */
332 sunpos (mjed, &lsn, &rsn, 0);
333 lg = lsn + PI;
334
335 /* mean daily motion is derived fro mean distance */
336 e_n = 0.9856076686/pow((double)op->e_a, 1.5);
337
338 /* correct for light time by computing position at time mjd, then
339 * again at mjd-dt, where
340 * dt = time it takes light to travel earth-object distance.
341 */
342 dt = 0;
343 for (pass = 0; pass < 2; pass++) {
344
345 reduce_elements (op->e_epoch, mjd-dt, degrad(op->e_inc),
346 degrad (op->e_om), degrad (op->e_Om),
347 &inc, &om, &Om);
348
349 tp = mjed - dt - (op->e_cepoch - op->e_M/e_n);
350 if (vrc (&nu, &rp, tp, op->e_e, op->e_a*(1-op->e_e)) < 0)
351 op->o_flags |= NOCIRCUM;
352 nu = degrad(nu);
353 lo = nu + om;
354 slo = sin(lo);
355 clo = cos(lo);
356 spsi = slo*sin(inc);
357 y = slo*cos(inc);
358 psi = asin(spsi);
359 lpd = atan(y/clo)+Om;
360 if (clo<0) lpd += PI;
361 range (&lpd, 2*PI);
362 cpsi = cos(psi);
363 rpd = rp*cpsi;
364 ll = lpd-lg;
365 rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
366
367 dt = rho*LTAU/3600.0/24.0; /* light travel time, in days / AU */
368 }
369
370 /* compute sin and cos of ll */
371 sll = sin(ll);
372 cll = cos(ll);
373
374 /* find geocentric ecliptic longitude and latitude */
375 if (rpd < rsn)
376 lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
377 else
378 lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
379 range (&lam, 2*PI);
380 bet = atan(rpd*spsi*sin(lam-lpd)/(cpsi*rsn*sll));
381
382 /* fill in all of op->s_* stuff except s_size and s_mag */
383 cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);
384
385 /* compute magnitude and size */
386 if (op->e_mag.whichm == MAG_HG) {
387 /* the H and G parameters from the Astro. Almanac.
388 */
389 hg_mag (op->e_mag.m1, op->e_mag.m2, rp, rho, rsn, &mag);
390 if (op->e_size)
391 op->s_size = (float)(op->e_size / rho);
392 else
393 op->s_size = (float)(h_albsize (op->e_mag.m1)/rho);
394 } else {
395 /* the g/k model of comets */
396 gk_mag (op->e_mag.m1, op->e_mag.m2, rp, rho, &mag);
397 op->s_size = (float)(op->e_size / rho);
398 }
399 set_smag (op, mag);
400
401 return (0);
402 }
403
404 /* compute sky circumstances of an object in heliocentric hyperbolic orbit.
405 */
406 static int
obj_hyperbolic(Now * np,Obj * op)407 obj_hyperbolic (Now *np, Obj *op)
408 {
409 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/
410 double dt; /* light travel time to object */
411 double lg; /* helio long of earth */
412 double nu; /* true anomaly and eccentric anomaly */
413 double rp=0; /* distance from the sun */
414 double lo, slo, clo; /* angle from ascending node */
415 double inc; /* inclination */
416 double psi=0; /* heliocentric latitude */
417 double spsi=0, cpsi=0; /* trig of heliocentric latitude */
418 double lpd; /* heliocentric longitude */
419 double rho=0; /* distance from the Earth */
420 double om; /* arg of perihelion */
421 double Om; /* long of ascending node. */
422 double lam; /* geocentric ecliptic longitude */
423 double bet; /* geocentric ecliptic latitude */
424 double e; /* fast eccentricity */
425 double ll=0, sll, cll; /* helio angle between object and earth */
426 double mag; /* magnitude */
427 double a; /* mean distance */
428 double tp; /* time from perihelion (days) */
429 double rpd=0;
430 double y;
431 int pass;
432
433 /* find solar ecliptical longitude and distance to sun from earth */
434 sunpos (mjed, &lsn, &rsn, 0);
435
436 lg = lsn + PI;
437 e = op->h_e;
438 a = op->h_qp/(e - 1.0);
439
440 /* correct for light time by computing position at time mjd, then
441 * again at mjd-dt, where
442 * dt = time it takes light to travel earth-object distance.
443 */
444 dt = 0;
445 for (pass = 0; pass < 2; pass++) {
446
447 reduce_elements (op->h_epoch, mjd-dt, degrad(op->h_inc),
448 degrad (op->h_om), degrad (op->h_Om),
449 &inc, &om, &Om);
450
451 tp = mjed - dt - op->h_ep;
452 if (vrc (&nu, &rp, tp, op->h_e, op->h_qp) < 0)
453 op->o_flags |= NOCIRCUM;
454 nu = degrad(nu);
455 lo = nu + om;
456 slo = sin(lo);
457 clo = cos(lo);
458 spsi = slo*sin(inc);
459 y = slo*cos(inc);
460 psi = asin(spsi);
461 lpd = atan(y/clo)+Om;
462 if (clo<0) lpd += PI;
463 range (&lpd, 2*PI);
464 cpsi = cos(psi);
465 rpd = rp*cpsi;
466 ll = lpd-lg;
467 rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
468
469 dt = rho*5.775518e-3; /* light travel time, in days */
470 }
471
472 /* compute sin and cos of ll */
473 sll = sin(ll);
474 cll = cos(ll);
475
476 /* find geocentric ecliptic longitude and latitude */
477 if (rpd < rsn)
478 lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
479 else
480 lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
481 range (&lam, 2*PI);
482 bet = atan(rpd*spsi*sin(lam-lpd)/(cpsi*rsn*sll));
483
484 /* fill in all of op->s_* stuff except s_size and s_mag */
485 cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);
486
487 /* compute magnitude and size */
488 gk_mag (op->h_g, op->h_k, rp, rho, &mag);
489 set_smag (op, mag);
490 op->s_size = (float)(op->h_size / rho);
491
492 return (0);
493 }
494
495 /* compute sky circumstances of an object in heliocentric hyperbolic orbit.
496 */
497 static int
obj_parabolic(Now * np,Obj * op)498 obj_parabolic (Now *np, Obj *op)
499 {
500 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/
501 double lam; /* geocentric ecliptic longitude */
502 double bet; /* geocentric ecliptic latitude */
503 double mag; /* magnitude */
504 double inc, om, Om;
505 double lpd, psi, rp, rho;
506 double dt;
507 int pass;
508
509 /* find solar ecliptical longitude and distance to sun from earth */
510 sunpos (mjed, &lsn, &rsn, 0);
511
512 /* two passes to correct lam and bet for light travel time. */
513 dt = 0.0;
514 for (pass = 0; pass < 2; pass++) {
515 reduce_elements (op->p_epoch, mjd-dt, degrad(op->p_inc),
516 degrad(op->p_om), degrad(op->p_Om), &inc, &om, &Om);
517 comet (mjed-dt, op->p_ep, inc, om, op->p_qp, Om,
518 &lpd, &psi, &rp, &rho, &lam, &bet);
519 dt = rho*LTAU/3600.0/24.0; /* light travel time, in days / AU */
520 }
521
522 /* fill in all of op->s_* stuff except s_size and s_mag */
523 cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);
524
525 /* compute magnitude and size */
526 gk_mag (op->p_g, op->p_k, rp, rho, &mag);
527 set_smag (op, mag);
528 op->s_size = (float)(op->p_size / rho);
529
530 return (0);
531 }
532
533 /* find sun's circumstances now.
534 */
535 static int
sun_cir(Now * np,Obj * op)536 sun_cir (Now *np, Obj *op)
537 {
538 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/
539 double bsn; /* true latitude beta of sun */
540 double dhlong;
541
542 sunpos (mjed, &lsn, &rsn, &bsn);/* sun's true coordinates; mean ecl. */
543
544 op->s_sdist = 0.0;
545 op->s_elong = 0.0;
546 op->s_phase = 100.0;
547 set_smag (op, -26.8); /* TODO */
548 dhlong = lsn-PI; /* geo- to helio- centric */
549 range (&dhlong, 2*PI);
550 op->s_hlong = (float)dhlong;
551 op->s_hlat = (float)(-bsn);
552
553 /* fill sun's ra/dec, alt/az in op */
554 cir_pos (np, bsn, lsn, &rsn, op);
555 op->s_edist = (float)rsn;
556 op->s_size = (float)(raddeg(4.65242e-3/rsn)*3600*2);
557
558 return (0);
559 }
560
561 /* find moon's circumstances now.
562 */
563 static int
moon_cir(Now * np,Obj * op)564 moon_cir (Now *np, Obj *op)
565 {
566 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/
567 double lam; /* geocentric ecliptic longitude */
568 double bet; /* geocentric ecliptic latitude */
569 double edistau; /* earth-moon dist, in au */
570 double el; /* elongation, rads east */
571 double ms; /* sun's mean anomaly */
572 double md; /* moon's mean anomaly */
573 double i;
574
575 moon (mjed, &lam, &bet, &edistau, &ms, &md); /* mean ecliptic & EOD*/
576 sunpos (mjed, &lsn, &rsn, NULL); /* mean ecliptic & EOD*/
577
578 op->s_hlong = (float)lam; /* save geo in helio fields */
579 op->s_hlat = (float)bet;
580
581 /* find angular separation from sun */
582 elongation (lam, bet, lsn, &el);
583 op->s_elong = (float)raddeg(el); /* want degrees */
584
585 /* solve triangle of earth, sun, and elongation for moon-sun dist */
586 op->s_sdist = (float) sqrt (edistau*edistau + rsn*rsn
587 - 2.0*edistau*rsn*cos(el));
588
589 /* TODO: improve mag; this is based on a flat moon model. */
590 i = -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el))))
591 + 5*log10(edistau/.0025) /* dist */;
592 set_smag (op, i);
593
594 /* find phase -- allow for projection effects */
595 i = 0.1468*sin(el)*(1 - 0.0549*sin(md))/(1 - 0.0167*sin(ms));
596 op->s_phase = (float)((1+cos(PI-el-degrad(i)))/2*100);
597
598 /* fill moon's ra/dec, alt/az in op and update for topo dist */
599 cir_pos (np, bet, lam, &edistau, op);
600
601 op->s_edist = (float)edistau;
602 op->s_size = (float)(3600*2.0*raddeg(asin(MRAD/MAU/edistau)));
603 /* moon angular dia, seconds */
604
605 return (0);
606 }
607
608 /* fill in all of op->s_* stuff except s_size and s_mag.
609 * this is used for sol system objects (except sun and moon); never FIXED.
610 */
611 static void
cir_sky(Now * np,double lpd,double psi,double rp,double * rho,double lam,double bet,double lsn,double rsn,Obj * op)612 cir_sky (
613 Now *np,
614 double lpd, /* heliocentric ecliptic longitude */
615 double psi, /* heliocentric ecliptic lat */
616 double rp, /* dist from sun */
617 double *rho, /* dist from earth: in as geo, back as geo or topo */
618 double lam, /* true geocentric ecliptic long */
619 double bet, /* true geocentric ecliptic lat */
620 double lsn, /* true geoc lng of sun */
621 double rsn, /* dist from sn to earth*/
622 Obj *op)
623 {
624 double el; /* elongation */
625 double f; /* fractional phase from earth */
626
627 /* compute elongation and phase */
628 elongation (lam, bet, lsn, &el);
629 el = raddeg(el);
630 op->s_elong = (float)el;
631 f = 0.25 * ((rp+ *rho)*(rp+ *rho) - rsn*rsn)/(rp* *rho);
632 op->s_phase = (float)(f*100.0); /* percent */
633
634 /* set heliocentric long/lat; mean ecliptic and EOD */
635 op->s_hlong = (float)lpd;
636 op->s_hlat = (float)psi;
637
638 /* fill solar sys body's ra/dec, alt/az in op */
639 cir_pos (np, bet, lam, rho, op); /* updates rho */
640
641 /* set earth/planet and sun/planet distance */
642 op->s_edist = (float)(*rho);
643 op->s_sdist = (float)rp;
644 }
645
646 /* fill equatoreal and horizontal op-> fields; stern
647 *
648 * input: lam/bet/rho geocentric mean ecliptic and equinox of day
649 *
650 * algorithm at EOD:
651 * ecl_eq --> ra/dec geocentric mean equatoreal EOD (via mean obliq)
652 * deflect --> ra/dec relativistic deflection
653 * nut_eq --> ra/dec geocentric true equatoreal EOD
654 * ab_eq --> ra/dec geocentric apparent equatoreal EOD
655 * if (PREF_GEO) --> output
656 * ta_par --> ra/dec topocentric apparent equatoreal EOD
657 * if (!PREF_GEO) --> output
658 * hadec_aa --> alt/az topocentric horizontal
659 * refract --> alt/az observed --> output
660 *
661 * algorithm at fixed equinox:
662 * ecl_eq --> ra/dec geocentric mean equatoreal EOD (via mean obliq)
663 * deflect --> ra/dec relativistic deflection [for alt/az only]
664 * nut_eq --> ra/dec geocentric true equatoreal EOD [for aa only]
665 * ab_eq --> ra/dec geocentric apparent equatoreal EOD [for aa only]
666 * ta_par --> ra/dec topocentric apparent equatoreal EOD
667 * precess --> ra/dec topocentric equatoreal fixed equinox [eq only]
668 * --> output
669 * hadec_aa --> alt/az topocentric horizontal
670 * refract --> alt/az observed --> output
671 */
672 static void
cir_pos(Now * np,double bet,double lam,double * rho,Obj * op)673 cir_pos (
674 Now *np,
675 double bet, /* geo lat (mean ecliptic of date) */
676 double lam, /* geo long (mean ecliptic of date) */
677 double *rho, /* in: geocentric dist in AU; out: geo- or topocentic dist */
678 Obj *op) /* object to set s_ra/dec as per equinox */
679 {
680 double ra, dec; /* apparent ra/dec, corrected for nut/ab */
681 double tra, tdec; /* astrometric ra/dec, no nut/ab */
682 double lsn, rsn; /* solar geocentric (mean ecliptic of date) */
683 double ha_in, ha_out; /* local hour angle before/after parallax */
684 double dec_out; /* declination after parallax */
685 double dra, ddec; /* parallax correction */
686 double alt, az; /* current alt, az */
687 double lst; /* local sidereal time */
688 double rho_topo; /* topocentric distance in earth radii */
689
690 /* convert to equatoreal [mean equator, with mean obliquity] */
691 ecl_eq (mjed, bet, lam, &ra, &dec);
692 tra = ra; /* keep mean coordinates */
693 tdec = dec;
694
695 /* precess and save astrometric coordinates */
696 if (mjed != epoch)
697 precess (mjed, epoch, &tra, &tdec);
698 op->s_astrora = tra;
699 op->s_astrodec = tdec;
700
701 /* get sun position */
702 sunpos(mjed, &lsn, &rsn, NULL);
703
704 /* allow for relativistic light bending near the sun.
705 * (avoid calling deflect() for the sun itself).
706 */
707 if (!is_planet(op,SUN) && !is_planet(op,MOON))
708 deflect (mjed, op->s_hlong, op->s_hlat, rsn, lsn, *rho, &ra, &dec);
709
710 /* correct ra/dec to form geocentric apparent */
711 nut_eq (mjed, &ra, &dec);
712 if (!is_planet(op,MOON))
713 ab_eq (mjed, lsn, &ra, &dec);
714 op->s_gaera = ra;
715 op->s_gaedec = dec;
716
717 /* find parallax correction for equatoreal coords */
718 now_lst (np, &lst);
719 ha_in = hrrad(lst) - ra;
720 rho_topo = *rho * MAU/ERAD; /* convert to earth radii */
721 ta_par (ha_in, dec, lat, elev, &rho_topo, &ha_out, &dec_out);
722
723 /* transform into alt/az and apply refraction */
724 hadec_aa (lat, ha_out, dec_out, &alt, &az);
725 refract (pressure, temp, alt, &alt);
726 op->s_ha = ha_out;
727 op->s_alt = alt;
728 op->s_az = az;
729
730 /* Get parallax differences and apply to apparent or astrometric place
731 * as needed. For the astrometric place, rotating the CORRECTIONS
732 * back from the nutated equator to the mean equator will be
733 * neglected. This is an effect of about 0.1" at moon distance.
734 * We currently don't have an inverse nutation rotation.
735 */
736 if (pref_get(PREF_EQUATORIAL) == PREF_GEO) {
737 /* no topo corrections to eq. coords */
738 dra = ddec = 0.0;
739 } else {
740 dra = ha_in - ha_out; /* ra sign is opposite of ha */
741 ddec = dec_out - dec;
742 *rho = rho_topo * ERAD/MAU; /* return topocentric distance in AU */
743
744 ra = ra + dra;
745 dec = dec + ddec;
746 }
747 range(&ra, 2*PI);
748 op->s_ra = ra;
749 op->s_dec = dec;
750 }
751
752 /* given geocentric ecliptic longitude and latitude, lam and bet, of some object
753 * and the longitude of the sun, lsn, find the elongation, el. this is the
754 * actual angular separation of the object from the sun, not just the difference
755 * in the longitude. the sign, however, IS set simply as a test on longitude
756 * such that el will be >0 for an evening object <0 for a morning object.
757 * to understand the test for el sign, draw a graph with lam going from 0-2*PI
758 * down the vertical axis, lsn going from 0-2*PI across the hor axis. then
759 * define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and
760 * lam=lsn-PI. the "morning" regions are any values to the lower left of the
761 * first line and bounded within the second pair of lines.
762 * all angles in radians.
763 */
764 static void
elongation(double lam,double bet,double lsn,double * el)765 elongation (double lam, double bet, double lsn, double *el)
766 {
767 *el = acos(cos(bet)*cos(lam-lsn));
768 if (lam>lsn+PI || (lam>lsn-PI && lam<lsn)) *el = - *el;
769 }
770
771 /* apply relativistic light bending correction to ra/dec; stern
772 *
773 * The algorithm is from:
774 * Mean and apparent place computations in the new IAU
775 * system. III - Apparent, topocentric, and astrometric
776 * places of planets and stars
777 * KAPLAN, G. H.; HUGHES, J. A.; SEIDELMANN, P. K.;
778 * SMITH, C. A.; YALLOP, B. D.
779 * Astronomical Journal (ISSN 0004-6256), vol. 97, April 1989, p. 1197-1210.
780 *
781 * This article is a very good collection of formulea for geocentric and
782 * topocentric place calculation in general. The apparent and
783 * astrometric place calculation in this file currently does not follow
784 * the strict algorithm from this paper and hence is not fully correct.
785 * The entire calculation is currently based on the rotating EOD frame and
786 * not the "inertial" J2000 frame.
787 */
788 static void
deflect(double mjd1,double lpd,double psi,double rsn,double lsn,double rho,double * ra,double * dec)789 deflect (
790 double mjd1, /* equinox */
791 double lpd, double psi, /* heliocentric ecliptical long / lat */
792 double rsn, double lsn, /* distance and longitude of sun */
793 double rho, /* geocentric distance */
794 double *ra, double *dec)/* geocentric equatoreal */
795 {
796 double hra, hdec; /* object heliocentric equatoreal */
797 double el; /* HELIOCENTRIC elongation object--earth */
798 double g1, g2; /* relativistic weights */
799 double u[3]; /* object geocentric cartesian */
800 double q[3]; /* object heliocentric cartesian unit vect */
801 double e[3]; /* earth heliocentric cartesian unit vect */
802 double qe, uq, eu; /* scalar products */
803 int i; /* counter */
804
805 #define G 1.32712438e20 /* heliocentric grav const; in m^3*s^-2 */
806 #define c 299792458.0 /* speed of light in m/s */
807
808 elongation(lpd, psi, lsn-PI, &el);
809 el = fabs(el);
810 /* only continue if object is within about 10 deg around the sun,
811 * not obscured by the sun's disc (radius 0.25 deg) and farther away
812 * than the sun.
813 *
814 * precise geocentric deflection is: g1 * tan(el/2)
815 * radially outwards from sun; the vector munching below
816 * just applys this component-wise
817 * Note: el = HELIOCENTRIC elongation.
818 * g1 is always about 0.004 arc seconds
819 * g2 varies from 0 (highest contribution) to 2
820 */
821 if (el<degrad(170) || el>degrad(179.75) || rho<rsn) return;
822
823 /* get cartesian vectors */
824 sphcart(*ra, *dec, rho, u, u+1, u+2);
825
826 ecl_eq(mjd1, psi, lpd, &hra, &hdec);
827 sphcart(hra, hdec, 1.0, q, q+1, q+2);
828
829 ecl_eq(mjd1, 0.0, lsn-PI, &hra, &hdec);
830 sphcart(hra, hdec, 1.0, e, e+1, e+2);
831
832 /* evaluate scalar products */
833 qe = uq = eu = 0.0;
834 for(i=0; i<=2; ++i) {
835 qe += q[i]*e[i];
836 uq += u[i]*q[i];
837 eu += e[i]*u[i];
838 }
839
840 g1 = 2*G/(c*c*MAU)/rsn;
841 g2 = 1 + qe;
842
843 /* now deflect geocentric vector */
844 g1 /= g2;
845 for(i=0; i<=2; ++i)
846 u[i] += g1*(uq*e[i] - eu*q[i]);
847
848 /* back to spherical */
849 cartsph(u[0], u[1], u[2], ra, dec, &rho); /* rho thrown away */
850 }
851
852 /* estimate size in arc seconds @ 1AU from absolute magnitude, H, and assuming
853 * an albedo of 0.1. With this assumption an object with diameter of 1500m
854 * has an absolute mag of 18.
855 */
856 static double
h_albsize(double H)857 h_albsize (double H)
858 {
859 return (3600*raddeg(.707*1500*pow(2.51,(18-H)/2)/MAU));
860 }
861