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

Last change on this file since 1518 was 1457, checked in by cmv, 25 years ago

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