/* given a Now and an Obj with the object definition portion filled in, * fill in the sky position (s_*) portions. * calculation of positional coordinates reworked by * Michael Sternberg * 3/11/98: deflect was using op->s_hlong before being set in cir_pos(). * 4/19/98: just edit a comment */ #include #include #if defined(__STDC__) #include #endif #include "P_.h" #include "astro.h" #include "circum.h" #include "preferences.h" static int obj_planet P_((Now *np, Obj *op)); static int obj_fixed P_((Now *np, Obj *op)); static int obj_elliptical P_((Now *np, Obj *op)); static int obj_hyperbolic P_((Now *np, Obj *op)); static int obj_parabolic P_((Now *np, Obj *op)); static int sun_cir P_((Now *np, Obj *op)); static int moon_cir P_((Now *np, Obj *op)); static void cir_sky P_((Now *np, double lpd, double psi, double rp, double *rho, double lam, double bet, double lsn, double rsn, Obj *op)); static void cir_pos P_((Now *np, double bet, double lam, double *rho, Obj *op)); static void elongation P_((double lam, double bet, double lsn, double *el)); static void deflect P_((double mjd1, double lpd, double psi, double rsn, double lsn, double rho, double *ra, double *dec)); static double h_albsize P_((double H)); /* given a Now and an Obj, fill in the approprirate s_* fields within Obj. * return 0 if all ok, else -1. */ int obj_cir (np, op) Now *np; Obj *op; { switch (op->o_type) { case FIXED: return (obj_fixed (np, op)); case ELLIPTICAL: return (obj_elliptical (np, op)); case HYPERBOLIC: return (obj_hyperbolic (np, op)); case PARABOLIC: return (obj_parabolic (np, op)); case EARTHSAT: return (obj_earthsat (np, op)); case PLANET: return (obj_planet (np, op)); default: printf ("obj_cir() called with type %d\n", op->o_type); exit(1); return (-1); /* just for lint */ } } static int obj_planet (np, op) Now *np; Obj *op; { double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ double lpd, psi; /* heliocentric ecliptic long and lat */ double rp; /* dist from sun */ double rho; /* dist from earth */ double lam, bet; /* geocentric ecliptic long and lat */ double dia, mag; /* angular diameter at 1 AU and magnitude */ double f; /* fractional phase from earth */ int p; /* validate code and check for a few special cases */ p = op->pl.pl_code; if (p < 0 || p > MOON) { printf ("unknown planet code: %d\n", p); exit(1); } else if (p == SUN) return (sun_cir (np, op)); else if (p == MOON) return (moon_cir (np, op)); /* find solar ecliptical longitude and distance to sun from earth */ sunpos (mjed, &lsn, &rsn, 0); /* find helio long/lat; sun/planet and earth/plant dist; ecliptic * long/lat; diameter and mag. */ plans(mjed, p, &lpd, &psi, &rp, &rho, &lam, &bet, &dia, &mag); /* fill in all of op->s_* stuff except s_size and s_mag */ cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op); /* compute magnitude and angular size */ f = op->s_phase ? 5*log10(rp*rho) - 5*log10(op->s_phase/100) : 100; set_smag (op, mag+f); op->s_size = (float)(dia/rho); return (0); } static int obj_fixed (np, op) Now *np; Obj *op; { double lsn, rsn; /* true geoc lng of sun, dist from sn to earth*/ double lam, bet; /* geocentric ecliptic long and lat */ double ha; /* local hour angle */ double el; /* elongation */ double alt, az; /* current alt, az */ double ra, dec; /* ra and dec at epoch of date */ double lst; if (epoch != EOD && (float)epoch != op->f_epoch) { /* want a certain epoch -- if it's not what the database is at * we change the original to save time next time assuming the * user is likely to stick with this for a while. */ double tra = op->f_RA, tdec = op->f_dec; float tepoch = (float)epoch; /* compare w/float precision */ precess (op->f_epoch, tepoch, &tra, &tdec); op->f_epoch = tepoch; op->f_RA = (float)tra; op->f_dec = (float)tdec; } /* set ra/dec to astrometric @ epoch of date */ ra = op->f_RA; dec = op->f_dec; precess (op->f_epoch, mjd, &ra, &dec); /* convert equatoreal ra/dec to mean geocentric ecliptic lat/long */ eq_ecl (mjd, ra, dec, &bet, &lam); /* find solar ecliptical long.(mean equinox) and distance from earth */ sunpos (mjed, &lsn, &rsn, NULL); /* allow for relativistic light bending near the sun */ deflect (mjd, lam, bet, lsn, rsn, 1e10, &ra, &dec); /* TODO: correction for annual parallax would go here */ /* correct EOD equatoreal for nutation/aberation to form apparent * geocentric */ nut_eq(mjd, &ra, &dec); ab_eq(mjd, lsn, &ra, &dec); op->s_gaera = (float)ra; op->s_gaedec = (float)dec; /* set s_ra/dec -- apparent if EOD else astrometric */ if (epoch == EOD) { op->s_ra = (float)ra; op->s_dec = (float)dec; } else { /* annual parallax at time mjd is to be added here, too, but * technically in the frame of epoch (usually different from mjd) */ op->s_ra = op->f_RA; /* already precessed */ op->s_dec = op->f_dec; } /* compute elongation from ecliptic long/lat and sun geocentric long */ elongation (lam, bet, lsn, &el); el = raddeg(el); op->s_elong = (float)el; /* these are really the same fields ... op->s_mag = op->f_mag; op->s_size = op->f_size; */ /* alt, az: correct for refraction; use eod ra/dec. */ now_lst (np, &lst); ha = hrrad(lst) - ra; hadec_aa (lat, ha, dec, &alt, &az); refract (pressure, temp, alt, &alt); op->s_alt = alt; op->s_az = az; return (0); } /* compute sky circumstances of an object in heliocentric elliptic orbit at *np. */ static int obj_elliptical (np, op) Now *np; Obj *op; { double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ double dt; /* light travel time to object */ double lg; /* helio long of earth */ double nu, ea; /* true anomaly and eccentric anomaly */ double ma; /* mean anomaly */ double rp=0; /* distance from the sun */ double lo, slo, clo; /* angle from ascending node */ double inc; /* inclination */ double psi=0; /* heliocentric latitude */ double spsi=0, cpsi=0; /* trig of heliocentric latitude */ double lpd; /* heliocentric longitude */ double rho=0; /* distance from the Earth */ double om; /* arg of perihelion */ double Om; /* long of ascending node. */ double lam; /* geocentric ecliptic longitude */ double bet; /* geocentric ecliptic latitude */ double e; /* fast eccentricity */ double ll=0, sll, cll; /* helio angle between object and earth */ double mag; /* magnitude */ double e_n; /* mean daily motion */ double rpd=0; double y; int pass; /* find location of earth from sun now */ sunpos (mjed, &lsn, &rsn, 0); lg = lsn + PI; /* faster access to eccentricty */ e = op->e_e; /* mean daily motion is derived fro mean distance */ e_n = 0.9856076686/pow((double)op->e_a, 1.5); /* correct for light time by computing position at time mjd, then * again at mjd-dt, where * dt = time it takes light to travel earth-object distance. */ dt = 0; for (pass = 0; pass < 2; pass++) { reduce_elements (op->e_epoch, mjd-dt, degrad(op->e_inc), degrad (op->e_om), degrad (op->e_Om), &inc, &om, &Om); ma = degrad (op->e_M + (mjed - op->e_cepoch - dt) * e_n); anomaly (ma, e, &nu, &ea); rp = op->e_a * (1-e*e) / (1+e*cos(nu)); lo = nu + om; slo = sin(lo); clo = cos(lo); spsi = slo*sin(inc); y = slo*cos(inc); psi = asin(spsi); lpd = atan(y/clo)+Om; if (clo<0) lpd += PI; range (&lpd, 2*PI); cpsi = cos(psi); rpd = rp*cpsi; ll = lpd-lg; rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll)); dt = rho*LTAU/3600.0/24.0; /* light travel time, in days / AU */ } /* compute sin and cos of ll */ sll = sin(ll); cll = cos(ll); /* find geocentric ecliptic longitude and latitude */ if (rpd < rsn) lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI; else lam = atan(rsn*sll/(rpd-rsn*cll))+lpd; range (&lam, 2*PI); bet = atan(rpd*spsi*sin(lam-lpd)/(cpsi*rsn*sll)); /* fill in all of op->s_* stuff except s_size and s_mag */ cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op); /* compute magnitude and size */ if (op->e_mag.whichm == MAG_HG) { /* the H and G parameters from the Astro. Almanac. */ if (op->e_size) op->s_size = (float)(op->e_size / rho); else { hg_mag (op->e_mag.m1, op->e_mag.m2, rp, rho, rsn, &mag); op->s_size = (float)(h_albsize (op->e_mag.m1)/rho); } } else { /* the g/k model of comets */ gk_mag (op->e_mag.m1, op->e_mag.m2, rp, rho, &mag); op->s_size = (float)(op->e_size / rho); } set_smag (op, mag); return (0); } /* compute sky circumstances of an object in heliocentric hyperbolic orbit. */ static int obj_hyperbolic (np, op) Now *np; Obj *op; { double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ double dt; /* light travel time to object */ double lg; /* helio long of earth */ double nu, ea; /* true anomaly and eccentric anomaly */ double ma; /* mean anomaly */ double rp=0; /* distance from the sun */ double lo, slo, clo; /* angle from ascending node */ double inc; /* inclination */ double psi=0; /* heliocentric latitude */ double spsi=0, cpsi=0; /* trig of heliocentric latitude */ double lpd; /* heliocentric longitude */ double rho=0; /* distance from the Earth */ double om; /* arg of perihelion */ double Om; /* long of ascending node. */ double lam; /* geocentric ecliptic longitude */ double bet; /* geocentric ecliptic latitude */ double e; /* fast eccentricity */ double ll=0, sll, cll; /* helio angle between object and earth */ double n; /* mean daily motion */ double mag; /* magnitude */ double a; /* mean distance */ double rpd=0; double y; int pass; /* find solar ecliptical longitude and distance to sun from earth */ sunpos (mjed, &lsn, &rsn, 0); lg = lsn + PI; e = op->h_e; a = op->h_qp/(e - 1.0); n = .98563/sqrt(a*a*a); /* correct for light time by computing position at time mjd, then * again at mjd-dt, where * dt = time it takes light to travel earth-object distance. */ dt = 0; for (pass = 0; pass < 2; pass++) { reduce_elements (op->h_epoch, mjd-dt, degrad(op->h_inc), degrad (op->h_om), degrad (op->h_Om), &inc, &om, &Om); ma = degrad ((mjed - op->h_ep - dt) * n); anomaly (ma, e, &nu, &ea); rp = a * (e*e-1.0) / (1.0+e*cos(nu)); lo = nu + om; slo = sin(lo); clo = cos(lo); spsi = slo*sin(inc); y = slo*cos(inc); psi = asin(spsi); lpd = atan(y/clo)+Om; if (clo<0) lpd += PI; range (&lpd, 2*PI); cpsi = cos(psi); rpd = rp*cpsi; ll = lpd-lg; rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll)); dt = rho*5.775518e-3; /* light travel time, in days */ } /* compute sin and cos of ll */ sll = sin(ll); cll = cos(ll); /* find geocentric ecliptic longitude and latitude */ if (rpd < rsn) lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI; else lam = atan(rsn*sll/(rpd-rsn*cll))+lpd; range (&lam, 2*PI); bet = atan(rpd*spsi*sin(lam-lpd)/(cpsi*rsn*sll)); /* fill in all of op->s_* stuff except s_size and s_mag */ cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op); /* compute magnitude and size */ gk_mag (op->h_g, op->h_k, rp, rho, &mag); set_smag (op, mag); op->s_size = (float)(op->h_size / rho); return (0); } /* compute sky circumstances of an object in heliocentric hyperbolic orbit. */ static int obj_parabolic (np, op) Now *np; Obj *op; { double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ double lam; /* geocentric ecliptic longitude */ double bet; /* geocentric ecliptic latitude */ double mag; /* magnitude */ double inc, om, Om; double lpd, psi, rp, rho; double dt; int pass; /* find solar ecliptical longitude and distance to sun from earth */ sunpos (mjed, &lsn, &rsn, 0); /* two passes to correct lam and bet for light travel time. */ dt = 0.0; for (pass = 0; pass < 2; pass++) { reduce_elements (op->p_epoch, mjd-dt, degrad(op->p_inc), degrad(op->p_om), degrad(op->p_Om), &inc, &om, &Om); comet (mjed-dt, op->p_ep, inc, om, op->p_qp, Om, &lpd, &psi, &rp, &rho, &lam, &bet); dt = rho*LTAU/3600.0/24.0; /* light travel time, in days / AU */ } /* fill in all of op->s_* stuff except s_size and s_mag */ cir_sky (np, lpd, psi, rp, &rho, lam, bet, lsn, rsn, op); /* compute magnitude and size */ gk_mag (op->p_g, op->p_k, rp, rho, &mag); set_smag (op, mag); op->s_size = (float)(op->p_size / rho); return (0); } /* find sun's circumstances now. */ static int sun_cir (np, op) Now *np; Obj *op; { double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ double bsn; /* true latitude beta of sun */ double dhlong; sunpos (mjed, &lsn, &rsn, &bsn);/* sun's true coordinates; mean ecl. */ op->s_sdist = 0.0; op->s_elong = 0.0; op->s_phase = 100.0; set_smag (op, -26.8); /* TODO */ dhlong = lsn-PI; /* geo- to helio- centric */ range (&dhlong, 2*PI); op->s_hlong = (float)dhlong; op->s_hlat = (float)(-bsn); /* fill sun's ra/dec, alt/az in op */ cir_pos (np, bsn, lsn, &rsn, op); op->s_edist = (float)rsn; op->s_size = (float)(raddeg(4.65242e-3/rsn)*3600*2); return (0); } /* find moon's circumstances now. */ static int moon_cir (np, op) Now *np; Obj *op; { double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ double lam; /* geocentric ecliptic longitude */ double bet; /* geocentric ecliptic latitude */ double edistau; /* earth-moon dist, in au */ double el; /* elongation, rads east */ double ms; /* sun's mean anomaly */ double md; /* moon's mean anomaly */ double i; moon (mjed, &lam, &bet, &edistau, &ms, &md); /* mean ecliptic & EOD*/ sunpos (mjed, &lsn, &rsn, NULL); /* mean ecliptic & EOD*/ op->s_hlong = (float)lam; /* save geo in helio fields */ op->s_hlat = (float)bet; /* find angular separation from sun */ elongation (lam, bet, lsn, &el); op->s_elong = (float)raddeg(el); /* want degrees */ /* solve triangle of earth, sun, and elongation for moon-sun dist */ op->s_sdist = (float) sqrt (edistau*edistau + rsn*rsn - 2.0*edistau*rsn*cos(el)); /* TODO: improve mag; this is based on a flat moon model. */ set_smag (op, -12.7 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el))))); /* find phase -- allow for projection effects */ i = 0.1468*sin(el)*(1 - 0.0549*sin(md))/(1 - 0.0167*sin(ms)); op->s_phase = (float)((1+cos(PI-el-degrad(i)))/2*100); /* fill moon's ra/dec, alt/az in op and update for topo dist */ cir_pos (np, bet, lam, &edistau, op); op->s_edist = (float)edistau; op->s_size = (float)(3600*2.0*raddeg(asin(MRAD/MAU/edistau))); /* moon angular dia, seconds */ return (0); } /* fill in all of op->s_* stuff except s_size and s_mag. * this is used for sol system objects (except sun and moon); never FIXED. */ static void cir_sky (np, lpd, psi, rp, rho, lam, bet, lsn, rsn, op) Now *np; double lpd, psi; /* heliocentric ecliptic long and lat */ double rp; /* dist from sun */ double *rho; /* dist from earth: in as geo, back as geo or topo */ double lam, bet; /* true geocentric ecliptic long and lat */ double lsn, rsn; /* true geoc lng of sun; dist from sn to earth*/ Obj *op; { double el; /* elongation */ double f; /* fractional phase from earth */ /* compute elongation and phase */ elongation (lam, bet, lsn, &el); el = raddeg(el); op->s_elong = (float)el; f = 0.25 * ((rp+ *rho)*(rp+ *rho) - rsn*rsn)/(rp* *rho); op->s_phase = (float)(f*100.0); /* percent */ /* set heliocentric long/lat; mean ecliptic and EOD */ op->s_hlong = (float)lpd; op->s_hlat = (float)psi; /* fill solar sys body's ra/dec, alt/az in op */ cir_pos (np, bet, lam, rho, op); /* updates rho */ /* set earth/planet and sun/planet distance */ op->s_edist = (float)(*rho); op->s_sdist = (float)rp; } /* fill equatoreal and horizontal op-> fields; stern * * input: lam/bet/rho geocentric mean ecliptic and equinox of day * * algorithm at EOD: * ecl_eq --> ra/dec geocentric mean equatoreal EOD (via mean obliq) * deflect --> ra/dec relativistic deflection * nut_eq --> ra/dec geocentric true equatoreal EOD * ab_eq --> ra/dec geocentric apparent equatoreal EOD * if (PREF_GEO) --> output * ta_par --> ra/dec topocentric apparent equatoreal EOD * if (!PREF_GEO) --> output * hadec_aa --> alt/az topocentric horizontal * refract --> alt/az observed --> output * * algorithm at fixed epoch: * ecl_eq --> ra/dec geocentric mean equatoreal EOD (via mean obliq) * deflect --> ra/dec relativistic deflection [for alt/az only] * nut_eq --> ra/dec geocentric true equatoreal EOD [for aa only] * ab_eq --> ra/dec geocentric apparent equatoreal EOD [for aa only] * ta_par --> ra/dec topocentric apparent equatoreal EOD * precess --> ra/dec topocentric equatoreal fixed equinox [eq only] * --> output * hadec_aa --> alt/az topocentric horizontal * refract --> alt/az observed --> output */ static void cir_pos (np, bet, lam, rho, op) Now *np; double bet, lam;/* geo lat/long (mean ecliptic of date) */ double *rho; /* in: geocentric dist in AU; out: geo- or topocentic dist */ Obj *op; /* object to set s_ra/dec as per epoch */ { double ra, dec; /* apparent ra/dec, corrected for nut/ab */ double tra, tdec; /* astrometric ra/dec, no nut/ab */ double lsn, rsn; /* solar geocentric (mean ecliptic of date) */ double ha_in, ha_out; /* local hour angle before/after parallax */ double dec_out; /* declination after parallax */ double dra, ddec; /* parallax correction */ double alt, az; /* current alt, az */ double lst; /* local sidereal time */ double rho_topo; /* topocentric distance in earth radii */ /* convert to equatoreal [mean equator, with mean obliquity] */ ecl_eq (mjd, bet, lam, &ra, &dec); tra = ra; /* keep mean coordinates */ tdec = dec; /* get sun position */ sunpos(mjed, &lsn, &rsn, NULL); /* allow for relativistic light bending near the sun. * (avoid calling deflect() for the sun itself). */ if (!is_planet(op,SUN) && !is_planet(op,MOON)) deflect (mjd, op->s_hlong, op->s_hlat, lsn, rsn, *rho, &ra, &dec); /* correct ra/dec to form geocentric apparent */ nut_eq (mjd, &ra, &dec); if (!is_planet(op,MOON)) ab_eq (mjd, lsn, &ra, &dec); op->s_gaera = (float)ra; op->s_gaedec = (float)dec; /* find parallax correction for equatoreal coords */ now_lst (np, &lst); ha_in = hrrad(lst) - ra; rho_topo = *rho * MAU/ERAD; /* convert to earth radii */ ta_par (ha_in, dec, lat, elev, &rho_topo, &ha_out, &dec_out); /* transform into alt/az and apply refraction */ hadec_aa (lat, ha_out, dec_out, &alt, &az); refract (pressure, temp, alt, &alt); op->s_alt = alt; op->s_az = az; /* Get parallax differences and apply to apparent or astrometric place * as needed. For the astrometric place, rotating the CORRECTIONS * back from the nutated equator to the mean equator will be * neglected. This is an effect of about 0.1" at moon distance. * We currently don't have an inverse nutation rotation. */ if (pref_get(PREF_EQUATORIAL) == PREF_GEO) { /* no topo corrections to eq. coords */ dra = ddec = 0.0; } else { dra = ha_in - ha_out; /* ra sign is opposite of ha */ ddec = dec_out - dec; *rho = rho_topo * ERAD/MAU; /* return topocentric distance in AU */ } /* fill in ra/dec fields */ if (epoch == EOD) { /* apparent geo/topocentric */ ra = ra + dra; dec = dec + ddec; } else { /* astrometric geo/topocent */ ra = tra + dra; dec = tdec + ddec; precess (mjd, epoch, &ra, &dec); } range(&ra, 2*PI); op->s_ra = (float)ra; op->s_dec = (float)dec; } /* given geocentric ecliptic longitude and latitude, lam and bet, of some object * and the longitude of the sun, lsn, find the elongation, el. this is the * actual angular separation of the object from the sun, not just the difference * in the longitude. the sign, however, IS set simply as a test on longitude * such that el will be >0 for an evening object <0 for a morning object. * to understand the test for el sign, draw a graph with lam going from 0-2*PI * down the vertical axis, lsn going from 0-2*PI across the hor axis. then * define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and * lam=lsn-PI. the "morning" regions are any values to the lower left of the * first line and bounded within the second pair of lines. * all angles in radians. */ static void elongation (lam, bet, lsn, el) double lam, bet, lsn; double *el; { *el = acos(cos(bet)*cos(lam-lsn)); if (lam>lsn+PI || (lam>lsn-PI && lamdegrad(179.75)) return; /* get cartesian vectors */ sphcart(*ra, *dec, rho, u, u+1, u+2); ecl_eq(mjd1, psi, lpd, &hra, &hdec); sphcart(hra, hdec, 1.0, q, q+1, q+2); ecl_eq(mjd1, 0.0, lsn-PI, &hra, &hdec); sphcart(hra, hdec, 1.0, e, e+1, e+2); /* evaluate scalar products */ qe = uq = eu = 0.0; for(i=0; i<=2; ++i) { qe += q[i]*e[i]; uq += u[i]*q[i]; eu += e[i]*u[i]; } g1 = 2*G/(c*c*MAU)/rsn; g2 = 1 + qe; /* now deflect geocentric vector */ g1 /= g2; for(i=0; i<=2; ++i) u[i] += g1*(uq*e[i] - eu*q[i]); /* back to spherical */ cartsph(u[0], u[1], u[2], ra, dec, &rho); /* rho thrown away */ } /* estimate size in arc seconds @ 1AU from absolute magnitude, H, and assuming * an albedo of 0.1. With this assumption an object with diameter of 1500m * has an absolute mag of 18. */ static double h_albsize (H) double H; { return (3600*raddeg(.707*1500*pow(2.51,(18-H)/2)/MAU)); } /* For RCS Only -- Do Not Edit */ static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: circum.c,v $ $Date: 2000/10/07 05:12:17 $ $Revision: 1.4 $ $Name: $"};