source: Sophya/trunk/SophyaExt/XephemAstroLib/circum.c@ 4024

Last change on this file since 4024 was 4017, checked in by cmv, 14 years ago

fichiers de Xephem 3.7.5 update, cmv 21/09/2011
File size: 26.7 KB
Line 
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 */
8
9#include <stdio.h>
10#include <math.h>
11#include <stdlib.h>
12
13#include "astro.h"
14#include "preferences.h"
15
16
17static int obj_planet (Now *np, Obj *op);
18static int obj_binary (Now *np, Obj *op);
19static int obj_2binary (Now *np, Obj *op);
20static int obj_fixed (Now *np, Obj *op);
21static int obj_elliptical (Now *np, Obj *op);
22static int obj_hyperbolic (Now *np, Obj *op);
23static int obj_parabolic (Now *np, Obj *op);
24static int sun_cir (Now *np, Obj *op);
25static int moon_cir (Now *np, Obj *op);
26static double solveKepler (double M, double e);
27static void binaryStarOrbit (double t, double T, double e, double o, double O,
28 double i, double a, double P, double *thetap, double *rhop);
29static void cir_sky (Now *np, double lpd, double psi, double rp, double *rho,
30 double lam, double bet, double lsn, double rsn, Obj *op);
31static void cir_pos (Now *np, double bet, double lam, double *rho, Obj *op);
32static void elongation (double lam, double bet, double lsn, double *el);
33static void deflect (double mjd1, double lpd, double psi, double rsn,
34 double lsn, double rho, double *ra, double *dec);
35static double h_albsize (double H);
36
37/* given a Now and an Obj, fill in the approprirate s_* fields within Obj.
38 * return 0 if all ok, else -1.
39 */
40int
41obj_cir (Now *np, Obj *op)
42{
43 op->o_flags &= ~NOCIRCUM;
44 switch (op->o_type) {
45 case BINARYSTAR: return (obj_binary (np, op));
46 case FIXED: return (obj_fixed (np, op));
47 case ELLIPTICAL: return (obj_elliptical (np, op));
48 case HYPERBOLIC: return (obj_hyperbolic (np, op));
49 case PARABOLIC: return (obj_parabolic (np, op));
50 case EARTHSAT: return (obj_earthsat (np, op));
51 case PLANET: return (obj_planet (np, op));
52 default:
53 printf ("obj_cir() called with type %d %s\n", op->o_type, op->o_name);
54 abort();
55 return (-1); /* just for lint */
56 }
57}
58
59static int
60obj_planet (Now *np, Obj *op)
61{
62 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/
63 double lpd, psi; /* heliocentric ecliptic long and lat */
64 double rp; /* dist from sun */
65 double rho; /* dist from earth */
66 double lam, bet; /* geocentric ecliptic long and lat */
67 double dia, mag; /* angular diameter at 1 AU and magnitude */
68 PLCode p;
69
70 /* validate code and check for a few special cases */
71 p = op->pl_code;
72 if (p == SUN)
73 return (sun_cir (np, op));
74 if (p == MOON)
75 return (moon_cir (np, op));
76 if (op->pl_moon != X_PLANET)
77 return (plmoon_cir (np, op));
78 if (p < 0 || p > MOON) {
79 printf ("unknown planet code: %d\n", p);
80 abort();
81 }
82
83 /* planet itself */
84
85 /* find solar ecliptical longitude and distance to sun from earth */
86 sunpos (mjed, &lsn, &rsn, 0);
87
88 /* find helio long/lat; sun/planet and earth/planet dist; ecliptic
89 * long/lat; diameter and mag.
90 */
91 plans(mjed, p, &lpd, &psi, &rp, &rho, &lam, &bet, &dia, &mag);
92
93 /* fill in all of op->s_* stuff except s_size and s_mag */
94 cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);
95
96 /* set magnitude and angular size */
97 set_smag (op, mag);
98 op->s_size = (float)(dia/rho);
99
100 return (0);
101}
102
103static int
104obj_binary (Now *np, Obj *op)
105{
106 /* always compute circumstances of primary */
107 if (obj_fixed (np, op) < 0)
108 return (0);
109
110 /* compute secondary only if requested, and always reset request flag */
111 if (!op->b_2compute)
112 return (0);
113 op->b_2compute = 0;
114 return (obj_2binary (np, op));
115}
116
117/* compute position of secondary component of a BINARYSTAR */
118static int
119obj_2binary (Now *np, Obj *op)
120{
121 if (op->b_nbp > 0) {
122 /* we just have discrete pa/sep, project each from primary */
123 int i;
124 for (i = 0; i < op->b_nbp; i++) {
125 BinPos *bp = &op->b_bp[i];
126 bp->bp_dec = op->s_dec + bp->bp_sep*cos(bp->bp_pa);
127 bp->bp_ra = op->s_ra + bp->bp_sep*sin(bp->bp_pa)/cos(op->s_dec);
128 }
129 } else {
130 BinOrbit *bp = &op->b_bo;
131 double t, theta, rho;
132
133 mjd_year (mjd, &t);
134 binaryStarOrbit (t, bp->bo_T, bp->bo_e, bp->bo_o, bp->bo_O,
135 bp->bo_i, bp->bo_a, bp->bo_P, &theta, &rho);
136 bp->bo_pa = (float)theta;
137 bp->bo_sep = (float)rho;
138 rho = degrad(rho/3600.); /* arc secs to rads */
139 bp->bo_dec = op->s_dec + rho*cos(theta);
140 bp->bo_ra = op->s_ra + rho*sin(theta)/cos(op->s_dec);
141 }
142
143 return (0);
144}
145
146/* from W. M. Smart */
147static void
148binaryStarOrbit (
149double t, /* desired ephemeris epoch, year */
150double T, /* epoch of periastron, year */
151double e, /* eccentricity */
152double o, /* argument of periastron, degrees */
153double O, /* ascending node, degrees */
154double i, /* inclination, degrees */
155double a, /* semi major axis, arcsecs */
156double P, /* period, years */
157double *thetap, /* position angle, rads E of N */
158double *rhop) /* separation, arcsecs */
159{
160 double M, E, cosE, nu, cosnu, r, rho, theta;
161
162 /* find mean anomaly, insure 0..2*PI */
163 M = 2*PI/P*(t-T);
164 range (&M, 2*PI);
165
166 /* solve for eccentric anomaly */
167 E = solveKepler (M, e);
168 cosE = cos(E);
169
170 /* find true anomaly and separation */
171 cosnu = (cosE - e)/(1.0 - e*cosE);
172 r = a*(1.0 - e*e)/(1.0 + e*cosnu);
173 nu = acos(cosnu);
174 if (E > PI)
175 nu = -nu;
176
177 /* project onto sky */
178 theta = atan(tan(nu+degrad(o))*cos(degrad(i))) + degrad(O);
179 rho = r*cos(nu+degrad(o))/cos(theta-degrad(O));
180 if (rho < 0) {
181 theta += PI;
182 rho = -rho;
183 }
184 range (&theta, 2*PI);
185
186 *thetap = theta;
187 *rhop = rho;
188}
189
190/* solve kepler equation using Newton-Raphson search.
191 * Charles and Tatum have shown it always converges starting with PI.
192 */
193static double
194solveKepler (double M, double e)
195{
196 double E, Eprime = PI;
197
198 do {
199 double cosE = cos(Eprime);
200 E = Eprime;
201 Eprime = (M - e*(E*cosE - sin(E)))/(1.0 - e*cosE);
202 } while (fabs(E-Eprime) > 1e-7);
203
204 return (Eprime);
205}
206
207static int
208obj_fixed (Now *np, Obj *op)
209{
210 double lsn, rsn; /* true geoc lng of sun, dist from sn to earth*/
211 double lam, bet; /* geocentric ecliptic long and lat */
212 double ha; /* local hour angle */
213 double el; /* elongation */
214 double alt, az; /* current alt, az */
215 double ra, dec; /* ra and dec at equinox of date */
216 double rpm, dpm; /* astrometric ra and dec with PM to now */
217 double lst;
218
219 /* on the assumption that the user will stick with their chosen display
220 * epoch for a while, we move the defining values to match and avoid
221 * precession for every call until it is changed again.
222 * N.B. only compare and store jd's to lowest precission (f_epoch).
223 * N.B. maintaining J2k ref (which is arbitrary) helps avoid accum err
224 */
225 if (epoch != EOD && (float)epoch != (float)op->f_epoch) {
226 double pr = op->f_RA, pd = op->f_dec, fe = (float)epoch;
227 /* first bring back to 2k */
228 precess (op->f_epoch, J2000, &pr, &pd);
229 pr += op->f_pmRA*(J2000-op->f_epoch);
230 pd += op->f_pmdec*(J2000-op->f_epoch);
231 /* then to epoch */
232 pr += op->f_pmRA*(fe-J2000);
233 pd += op->f_pmdec*(fe-J2000);
234 precess (J2000, fe, &pr, &pd);
235 op->f_RA = (float)pr;
236 op->f_dec = (float)pd;
237 op->f_epoch = (float)fe;
238 }
239
240 /* apply proper motion .. assume pm epoch reference equals equinox */
241 rpm = op->f_RA + op->f_pmRA*(mjd-op->f_epoch);
242 dpm = op->f_dec + op->f_pmdec*(mjd-op->f_epoch);
243
244 /* set ra/dec to astrometric @ equinox of date */
245 ra = rpm;
246 dec = dpm;
247 precess (op->f_epoch, mjed, &ra, &dec);
248
249 /* convert equatoreal ra/dec to mean geocentric ecliptic lat/long */
250 eq_ecl (mjed, ra, dec, &bet, &lam);
251
252 /* find solar ecliptical long.(mean equinox) and distance from earth */
253 sunpos (mjed, &lsn, &rsn, NULL);
254
255 /* allow for relativistic light bending near the sun */
256 deflect (mjed, lam, bet, lsn, rsn, 1e10, &ra, &dec);
257
258 /* TODO: correction for annual parallax would go here */
259
260 /* correct EOD equatoreal for nutation/aberation to form apparent
261 * geocentric
262 */
263 nut_eq(mjed, &ra, &dec);
264 ab_eq(mjed, lsn, &ra, &dec);
265 op->s_gaera = (float)ra;
266 op->s_gaedec = (float)dec;
267
268 /* set s_ra/dec -- apparent if EOD else astrometric */
269 if (epoch == EOD) {
270 op->s_ra = (float)ra;
271 op->s_dec = (float)dec;
272 } else {
273 /* annual parallax at time mjd is to be added here, too, but
274 * technically in the frame of equinox (usually different from mjd)
275 */
276 op->s_ra = rpm;
277 op->s_dec = dpm;
278 }
279
280 /* compute elongation from ecliptic long/lat and sun geocentric long */
281 elongation (lam, bet, lsn, &el);
282 el = raddeg(el);
283 op->s_elong = (float)el;
284
285 /* these are really the same fields ...
286 op->s_mag = op->f_mag;
287 op->s_size = op->f_size;
288 */
289
290 /* alt, az: correct for refraction; use eod ra/dec. */
291 now_lst (np, &lst);
292 ha = hrrad(lst) - ra;
293 hadec_aa (lat, ha, dec, &alt, &az);
294 refract (pressure, temp, alt, &alt);
295 op->s_alt = alt;
296 op->s_az = az;
297
298 return (0);
299}
300
301/* compute sky circumstances of an object in heliocentric elliptic orbit at *np.
302 */
303static int
304obj_elliptical (Now *np, Obj *op)
305{
306 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/
307 double dt; /* light travel time to object */
308 double lg; /* helio long of earth */
309 double nu; /* true anomaly */
310 double rp=0; /* distance from the sun */
311 double lo, slo, clo; /* angle from ascending node */
312 double inc; /* inclination */
313 double psi=0; /* heliocentric latitude */
314 double spsi=0, cpsi=0; /* trig of heliocentric latitude */
315 double lpd; /* heliocentric longitude */
316 double rho=0; /* distance from the Earth */
317 double om; /* arg of perihelion */
318 double Om; /* long of ascending node. */
319 double lam; /* geocentric ecliptic longitude */
320 double bet; /* geocentric ecliptic latitude */
321 double ll=0, sll, cll; /* helio angle between object and earth */
322 double mag; /* magnitude */
323 double e_n; /* mean daily motion */
324 double tp; /* time from perihelion (days) */
325 double rpd=0;
326 double y;
327 int pass;
328
329 /* find location of earth from sun now */
330 sunpos (mjed, &lsn, &rsn, 0);
331 lg = lsn + PI;
332
333 /* mean daily motion is derived fro mean distance */
334 e_n = 0.9856076686/pow((double)op->e_a, 1.5);
335
336 /* correct for light time by computing position at time mjd, then
337 * again at mjd-dt, where
338 * dt = time it takes light to travel earth-object distance.
339 */
340 dt = 0;
341 for (pass = 0; pass < 2; pass++) {
342
343 reduce_elements (op->e_epoch, mjd-dt, degrad(op->e_inc),
344 degrad (op->e_om), degrad (op->e_Om),
345 &inc, &om, &Om);
346
347 tp = mjed - dt - (op->e_cepoch - op->e_M/e_n);
348 if (vrc (&nu, &rp, tp, op->e_e, op->e_a*(1-op->e_e)) < 0)
349 op->o_flags |= NOCIRCUM;
350 nu = degrad(nu);
351 lo = nu + om;
352 slo = sin(lo);
353 clo = cos(lo);
354 spsi = slo*sin(inc);
355 y = slo*cos(inc);
356 psi = asin(spsi);
357 lpd = atan(y/clo)+Om;
358 if (clo<0) lpd += PI;
359 range (&lpd, 2*PI);
360 cpsi = cos(psi);
361 rpd = rp*cpsi;
362 ll = lpd-lg;
363 rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
364
365 dt = rho*LTAU/3600.0/24.0; /* light travel time, in days / AU */
366 }
367
368 /* compute sin and cos of ll */
369 sll = sin(ll);
370 cll = cos(ll);
371
372 /* find geocentric ecliptic longitude and latitude */
373 if (rpd < rsn)
374 lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
375 else
376 lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
377 range (&lam, 2*PI);
378 bet = atan(rpd*spsi*sin(lam-lpd)/(cpsi*rsn*sll));
379
380 /* fill in all of op->s_* stuff except s_size and s_mag */
381 cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op);
382
383 /* compute magnitude and size */
384 if (op->e_mag.whichm == MAG_HG) {
385 /* the H and G parameters from the Astro. Almanac.
386 */
387 if (op->e_size)
388 op->s_size = (float)(op->e_size / rho);
389 else {
390 hg_mag (op->e_mag.m1, op->e_mag.m2, rp, rho, rsn, &mag);
391 op->s_size = (float)(h_albsize (op->e_mag.m1)/rho);
392
393 }
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 */
406static int
407obj_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 */
497static int
498obj_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 */
535static int
536sun_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 */
563static int
564moon_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 */
611static void
612cir_sky (
613Now *np,
614double lpd, /* heliocentric ecliptic longitude */
615double psi, /* heliocentric ecliptic lat */
616double rp, /* dist from sun */
617double *rho, /* dist from earth: in as geo, back as geo or topo */
618double lam, /* true geocentric ecliptic long */
619double bet, /* true geocentric ecliptic lat */
620double lsn, /* true geoc lng of sun */
621double rsn, /* dist from sn to earth*/
622Obj *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 */
672static void
673cir_pos (
674Now *np,
675double bet, /* geo lat (mean ecliptic of date) */
676double lam, /* geo long (mean ecliptic of date) */
677double *rho, /* in: geocentric dist in AU; out: geo- or topocentic dist */
678Obj *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 /* get sun position */
696 sunpos(mjed, &lsn, &rsn, NULL);
697
698 /* allow for relativistic light bending near the sun.
699 * (avoid calling deflect() for the sun itself).
700 */
701 if (!is_planet(op,SUN) && !is_planet(op,MOON))
702 deflect (mjed, op->s_hlong, op->s_hlat, lsn, rsn, *rho, &ra, &dec);
703
704 /* correct ra/dec to form geocentric apparent */
705 nut_eq (mjed, &ra, &dec);
706 if (!is_planet(op,MOON))
707 ab_eq (mjed, lsn, &ra, &dec);
708 op->s_gaera = (float)ra;
709 op->s_gaedec = (float)dec;
710
711 /* find parallax correction for equatoreal coords */
712 now_lst (np, &lst);
713 ha_in = hrrad(lst) - ra;
714 rho_topo = *rho * MAU/ERAD; /* convert to earth radii */
715 ta_par (ha_in, dec, lat, elev, &rho_topo, &ha_out, &dec_out);
716
717 /* transform into alt/az and apply refraction */
718 hadec_aa (lat, ha_out, dec_out, &alt, &az);
719 refract (pressure, temp, alt, &alt);
720 op->s_alt = alt;
721 op->s_az = az;
722
723 /* Get parallax differences and apply to apparent or astrometric place
724 * as needed. For the astrometric place, rotating the CORRECTIONS
725 * back from the nutated equator to the mean equator will be
726 * neglected. This is an effect of about 0.1" at moon distance.
727 * We currently don't have an inverse nutation rotation.
728 */
729 if (pref_get(PREF_EQUATORIAL) == PREF_GEO) {
730 /* no topo corrections to eq. coords */
731 dra = ddec = 0.0;
732 } else {
733 dra = ha_in - ha_out; /* ra sign is opposite of ha */
734 ddec = dec_out - dec;
735 *rho = rho_topo * ERAD/MAU; /* return topocentric distance in AU */
736 }
737
738 /* fill in ra/dec fields */
739 if (epoch == EOD) { /* apparent geo/topocentric */
740 ra = ra + dra;
741 dec = dec + ddec;
742 } else { /* astrometric geo/topocent */
743 ra = tra + dra;
744 dec = tdec + ddec;
745 precess (mjed, epoch, &ra, &dec);
746 }
747 range(&ra, 2*PI);
748 op->s_ra = (float)ra;
749 op->s_dec = (float)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 */
764static void
765elongation (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 */
788static void
789deflect (
790double mjd1, /* equinox */
791double lpd, double psi, /* heliocentric ecliptical long / lat */
792double rsn, double lsn, /* distance and longitude of sun */
793double rho, /* geocentric distance */
794double *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 */
856static double
857h_albsize (double H)
858{
859 return (3600*raddeg(.707*1500*pow(2.51,(18-H)/2)/MAU));
860}
861
862/* For RCS Only -- Do Not Edit */
863static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: circum.c,v $ $Date: 2011-09-21 16:17:48 $ $Revision: 1.10 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.