![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
Changeset 2551 in Sophya for trunk/SophyaExt/XephemAstroLib/circum.c
- Timestamp:
- Jun 15, 2004, 6:54:12 PM (21 years ago)
- File:
-
- 1 edited
-
trunk/SophyaExt/XephemAstroLib/circum.c (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaExt/XephemAstroLib/circum.c
r1719 r2551 9 9 #include <stdio.h> 10 10 #include <math.h> 11 #if defined(__STDC__)12 11 #include <stdlib.h> 13 #endif 14 15 #include "P_.h" 12 16 13 #include "astro.h" 17 #include "circum.h"18 14 #include "preferences.h" 19 15 20 16 21 static int obj_planet P_((Now *np, Obj *op)); 22 static int obj_fixed P_((Now *np, Obj *op)); 23 static int obj_elliptical P_((Now *np, Obj *op)); 24 static int obj_hyperbolic P_((Now *np, Obj *op)); 25 static int obj_parabolic P_((Now *np, Obj *op)); 26 static int sun_cir P_((Now *np, Obj *op)); 27 static int moon_cir P_((Now *np, Obj *op)); 28 static 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)); 30 static void cir_pos P_((Now *np, double bet, double lam, double *rho, Obj *op)); 31 static void elongation P_((double lam, double bet, double lsn, double *el)); 32 static void deflect P_((double mjd1, double lpd, double psi, double rsn, 33 double lsn, double rho, double *ra, double *dec)); 34 static double h_albsize P_((double H)); 17 static int obj_planet (Now *np, Obj *op); 18 static int obj_binary (Now *np, Obj *op); 19 static int obj_2binary (Now *np, Obj *op); 20 static int obj_fixed (Now *np, Obj *op); 21 static int obj_elliptical (Now *np, Obj *op); 22 static int obj_hyperbolic (Now *np, Obj *op); 23 static int obj_parabolic (Now *np, Obj *op); 24 static int sun_cir (Now *np, Obj *op); 25 static int moon_cir (Now *np, Obj *op); 26 static double solveKepler (double M, double e); 27 static void binaryStarOrbit (double t, double T, double e, double o, double O, 28 double i, double a, double P, double *thetap, double *rhop); 29 static void cir_sky (Now *np, double lpd, double psi, double rp, double *rho, 30 double lam, double bet, double lsn, double rsn, Obj *op); 31 static void cir_pos (Now *np, double bet, double lam, double *rho, Obj *op); 32 static void elongation (double lam, double bet, double lsn, double *el); 33 static void deflect (double mjd1, double lpd, double psi, double rsn, 34 double lsn, double rho, double *ra, double *dec); 35 static double h_albsize (double H); 35 36 36 37 /* given a Now and an Obj, fill in the approprirate s_* fields within Obj. … … 38 39 */ 39 40 int 40 obj_cir (np, op) 41 Now *np; 42 Obj *op; 43 { 41 obj_cir (Now *np, Obj *op) 42 { 43 op->o_flags &= ~NOCIRCUM; 44 44 switch (op->o_type) { 45 case BINARYSTAR: return (obj_binary (np, op)); 45 46 case FIXED: return (obj_fixed (np, op)); 46 47 case ELLIPTICAL: return (obj_elliptical (np, op)); … … 50 51 case PLANET: return (obj_planet (np, op)); 51 52 default: 52 printf ("obj_cir() called with type %d \n", op->o_type);53 exit(1);53 printf ("obj_cir() called with type %d %s\n", op->o_type, op->o_name); 54 abort(); 54 55 return (-1); /* just for lint */ 55 56 } … … 57 58 58 59 static int 59 obj_planet (np, op) 60 Now *np; 61 Obj *op; 60 obj_planet (Now *np, Obj *op) 62 61 { 63 62 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ … … 67 66 double lam, bet; /* geocentric ecliptic long and lat */ 68 67 double dia, mag; /* angular diameter at 1 AU and magnitude */ 69 intp;68 PLCode p; 70 69 71 70 /* validate code and check for a few special cases */ 72 p = op->pl.pl_code; 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)); 73 78 if (p < 0 || p > MOON) { 74 79 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)); 80 abort(); 81 } 82 83 /* planet itself */ 81 84 82 85 /* find solar ecliptical longitude and distance to sun from earth */ … … 99 102 100 103 static int 101 obj_fixed (np, op) 102 Now *np; 103 Obj *op; 104 obj_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 */ 118 static int 119 obj_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 */ 147 static void 148 binaryStarOrbit ( 149 double t, /* desired ephemeris epoch, year */ 150 double T, /* epoch of periastron, year */ 151 double e, /* eccentricity */ 152 double o, /* argument of periastron, degrees */ 153 double O, /* ascending node, degrees */ 154 double i, /* inclination, degrees */ 155 double a, /* semi major axis, arcsecs */ 156 double P, /* period, years */ 157 double *thetap, /* position angle, rads E of N */ 158 double *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 */ 193 static double 194 solveKepler (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 207 static int 208 obj_fixed (Now *np, Obj *op) 104 209 { 105 210 double lsn, rsn; /* true geoc lng of sun, dist from sn to earth*/ … … 108 213 double el; /* elongation */ 109 214 double alt, az; /* current alt, az */ 110 double ra, dec; /* ra and dec at epoch of date */ 215 double ra, dec; /* ra and dec at equinox of date */ 216 double rpm, dpm; /* astrometric ra and dec with PM to now */ 111 217 double lst; 112 218 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; 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; 129 247 precess (op->f_epoch, mjed, &ra, &dec); 130 248 … … 154 272 } else { 155 273 /* annual parallax at time mjd is to be added here, too, but 156 * technically in the frame of e poch(usually different from mjd)274 * technically in the frame of equinox (usually different from mjd) 157 275 */ 158 op->s_ra = op->f_RA; /* already precessed */159 op->s_dec = op->f_dec;276 op->s_ra = rpm; 277 op->s_dec = dpm; 160 278 } 161 279 … … 184 302 */ 185 303 static int 186 obj_elliptical (np, op) 187 Now *np; 188 Obj *op; 304 obj_elliptical (Now *np, Obj *op) 189 305 { 190 306 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ … … 230 346 231 347 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)); 348 if (vrc (&nu, &rp, tp, op->e_e, op->e_a*(1-op->e_e)) < 0) 349 op->o_flags |= NOCIRCUM; 233 350 nu = degrad(nu); 234 351 lo = nu + om; … … 288 405 */ 289 406 static int 290 obj_hyperbolic (np, op) 291 Now *np; 292 Obj *op; 407 obj_hyperbolic (Now *np, Obj *op) 293 408 { 294 409 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ … … 309 424 double e; /* fast eccentricity */ 310 425 double ll=0, sll, cll; /* helio angle between object and earth */ 311 double n; /* mean daily motion */312 426 double mag; /* magnitude */ 313 427 double a; /* mean distance */ … … 323 437 e = op->h_e; 324 438 a = op->h_qp/(e - 1.0); 325 n = .98563/sqrt(a*a*a);326 439 327 440 /* correct for light time by computing position at time mjd, then … … 337 450 338 451 tp = mjed - dt - op->h_ep; 339 vrc (&nu, &rp, tp, op->h_e, op->h_qp); 452 if (vrc (&nu, &rp, tp, op->h_e, op->h_qp) < 0) 453 op->o_flags |= NOCIRCUM; 340 454 nu = degrad(nu); 341 455 lo = nu + om; … … 382 496 */ 383 497 static int 384 obj_parabolic (np, op) 385 Now *np; 386 Obj *op; 498 obj_parabolic (Now *np, Obj *op) 387 499 { 388 500 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ … … 422 534 */ 423 535 static int 424 sun_cir (np, op) 425 Now *np; 426 Obj *op; 536 sun_cir (Now *np, Obj *op) 427 537 { 428 538 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ … … 452 562 */ 453 563 static int 454 moon_cir (np, op) 455 Now *np; 456 Obj *op; 564 moon_cir (Now *np, Obj *op) 457 565 { 458 566 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ … … 480 588 481 589 /* 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))))); 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); 483 593 484 594 /* find phase -- allow for projection effects */ … … 500 610 */ 501 611 static void 502 cir_sky (np, lpd, psi, rp, rho, lam, bet, lsn, rsn, op) 503 Now *np; 504 double lpd, psi; /* heliocentric ecliptic long and lat */ 505 double rp; /* dist from sun */ 506 double *rho; /* dist from earth: in as geo, back as geo or topo */ 507 double lam, bet; /* true geocentric ecliptic long and lat */ 508 double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ 509 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) 510 623 { 511 624 double el; /* elongation */ … … 546 659 * refract --> alt/az observed --> output 547 660 * 548 * algorithm at fixed e poch:661 * algorithm at fixed equinox: 549 662 * ecl_eq --> ra/dec geocentric mean equatoreal EOD (via mean obliq) 550 663 * deflect --> ra/dec relativistic deflection [for alt/az only] … … 558 671 */ 559 672 static void 560 cir_pos (np, bet, lam, rho, op) 561 Now *np; 562 double bet, lam;/* geo lat/long (mean ecliptic of date) */ 563 double *rho; /* in: geocentric dist in AU; out: geo- or topocentic dist */ 564 Obj *op; /* object to set s_ra/dec as per epoch */ 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 */ 565 679 { 566 680 double ra, dec; /* apparent ra/dec, corrected for nut/ab */ … … 649 763 */ 650 764 static void 651 elongation (lam, bet, lsn, el) 652 double lam, bet, lsn; 653 double *el; 765 elongation (double lam, double bet, double lsn, double *el) 654 766 { 655 767 *el = acos(cos(bet)*cos(lam-lsn)); … … 675 787 */ 676 788 static void 677 deflect ( mjd1, lpd, psi, lsn, rsn, rho, ra, dec)678 double mjd1 ; /* epoch*/679 double lpd, psi;/* heliocentric ecliptical long / lat */680 double rsn, lsn;/* distance and longitude of sun */681 double rho ;/* geocentric distance */682 double *ra, *dec;/* geocentric equatoreal */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 */ 683 795 { 684 796 double hra, hdec; /* object heliocentric equatoreal */ … … 742 854 */ 743 855 static double 744 h_albsize (H) 745 double H; 856 h_albsize (double H) 746 857 { 747 858 return (3600*raddeg(.707*1500*pow(2.51,(18-H)/2)/MAU)); … … 749 860 750 861 /* For RCS Only -- Do Not Edit */ 751 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: circum.c,v $ $Date: 200 1-10-22 12:08:26 $ $Revision: 1.2$ $Name: not supported by cvs2svn $"};862 static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: circum.c,v $ $Date: 2004-06-15 16:52:37 $ $Revision: 1.3 $ $Name: not supported by cvs2svn $"};
Note:
See TracChangeset
for help on using the changeset viewer.
