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

Last change on this file since 1783 was 1719, checked in by cmv, 24 years ago

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