#include "SolSystem.h"#include <iostream>#include <cmath>#include <string>#include "vsop87_data.h"#include "chap95_data.h"Go to the source code of this file.
Compounds | |
| struct | plantbl |
Defines | |
| #define | A2000 2451545.00000000 |
| #define | Conv_Rad 0.017453292 |
| #define | Conv_Grad 57.29578122 |
| #define | Conv_Ora_Rad 0.261799387 |
| #define | SPD (24.0*3600.0) |
| #define | MJD0 2415020.0 |
| #define | J2000 (2451545.0 - MJD0) |
| #define | MAU (1.4959787e11) |
| #define | ERAD (6.37816e6) |
| #define | SIDRATE .9972695677 |
| #define | TMACC (10./3600./24.0) |
| #define | G 1.32712438e20 |
| #define | c 299792458.0 |
| #define | MAXPASSES 20 |
| #define | FIRSTSTEP (1.0/60.0/24.0) |
| #define | MAXLOOPS 10 |
| #define | MAXERR (0.25/60.) |
| #define | EPS (1e-9) |
| #define | VSOP_ASCALE 1e8 |
| #define | VSOP_SPHERICAL 1 |
| #define | VSOP_GETRATE 0 |
| #define | VSOP_A1000 365250.0 |
| #define | VSOP_MAXALPHA 5 |
| #define | DCOS(x) cos(degrad(x)) |
| #define | DSIN(x) sin(degrad(x)) |
| #define | DASIN(x) raddeg(asin(x)) |
| #define | DATAN2(y, x) raddeg(atan2((y),(x))) |
| #define | ERRLMT 0.0001 |
| #define | TWOPI (2*PI) |
| #define | STOPERR (1e-8) |
| #define | CHAP_MAXTPOW 2 |
| #define | CHAR short |
| #define | NARGS 18 |
| #define | MOSHIER_J2000 (2451545.0) |
| #define | MOSHIER_BEGIN (1221000.5 - MJD0) |
| #define | MOSHIER_END (2798525.5 - MJD0) |
| #define | DTR 1.7453292519943295769e-2 |
| #define | RTD 5.7295779513082320877e1 |
| #define | RTS 2.0626480624709635516e5 |
| #define | STR 4.8481368110953599359e-6 |
| #define | AUKM 1.4959787e8 |
| #define | EarthRadius 6378.16 |
| #define | NUT_SCALE 1e4 |
| #define | NUT_SERIES 106 |
| #define | NUT_MAXMUL 4 |
| #define | SECPERCIRC (3600.*360.) |
| #define | LTLIM 14.5 |
| #define | GELIM 15.5 |
| #define | MAXRERR degrad(0.1/3600.) |
| #define | ABERR_CONST (20.49552/3600./180.*PI) |
| #define | AB_ECL_EOD 0 |
| #define | AB_EQ_EOD 1 |
| #define | EQtoECL 1 |
| #define | ECLtoEQ (-1) |
| #define | TABSTART 1620.0 |
| #define | TABEND 2004.0 |
| #define | TABSIZ 385 |
Functions | |
| void | precess (double mjd1, double mjd2, double *ra, double *dec) |
| void | range (double *v, double r) |
| int | Conv_GST_UT (double jd, double gst, double *tu) |
| int | Sun_Pos (double jd, double *ra, double *decs, double *lon, double *r) |
| int | Conv_lb_ad (double jd, double l, double b, double *ra, double *dec) |
| void | riset (double ra, double dec, double lat, double dis, double *lstr, double *lsts, double *azr, double *azs, int *status) |
| int | vsop87 (double mjd, int obj, double prec, double *ret) |
| void | sunpos (double mjd, double *lsn, double *rsn, double *bsn) |
| void | cal_mjd (int mn, double dy, int yr, double *mjd) |
| double | mjd_day (double jd) |
| void | rnd_second (double *t) |
| void | year_mjd (double y, double *mjd) |
| void | mjd_year (double mjd, double *yr) |
| void | mjd_dpm (double mjd, int *ndays) |
| int | mjd_dow (double mjd, int *dow) |
| void | mjd_cal (double mjd, int *mn, double dy, int *yr) |
| void | cartsph (double x, double y, double z, double *l, double *b, double *r) |
| void | sphcart (double l, double b, double r, double *x, double *y, double *z) |
| void | zero_mem (void *loc, unsigned len) |
| void | reduce_elements (double mjd0, double mjd, double inc0, double ap0, double om0, double *inc, double *ap, double *om) |
| void | obliquity (double mjd, double *eps) |
| void | anomaly (double ma, double s, double *nu, double *ea) |
| void | comet (double mjd, double ep, double inc, double ap, double qp, double om, double *lpd, double *psi, double *rp, double *rho, double *lam, double *bet) |
| int | chap95 (double mjd, int obj, double prec, double *ret) |
| void | moon (double mjd, double *lam, double *bet, double *rho, double *msp, double *mdp) |
| void | nut_eq (double mjd, double *ra, double *dec) |
| void | nutation (double mjd, double *deps, double *dpsi) |
| void | ta_par (double tha, double tdec, double phi, double ht, double *rho, double *aha, double *adec) |
| void | refract (double pr, double tr, double ta, double *aa) |
| void | ab_eq (double mjd, double lsn, double *ra, double *dec) |
| void | ab_ecl (double mjd, double lsn, double *lam, double *bet) |
| void | unrefract (double pr, double tr, double aa, double *ta) |
| void | eq_ecl (double mjd, double ra, double dec, double *lat, double *lng) |
| void | ecl_eq (double mjd, double lat, double lng, double *ra, double *dec) |
| void | aa_hadec (double lat, double alt, double az, double *ha, double *dec) |
| void | hadec_aa (double lat, double ha, double dec, double *alt, double *az) |
| void | solve_sphere (double A, double b, double cc, double sc, double *cap, double *Bp) |
| double | deltat (double mjd) |
| double | mm_mjed (double mjd) |
| double | Tempo_Siderale (double J_D, double Ora_Un_Dec) |
| void | now_lst (double mdj, double lng, double *lstp) |
| void | plans (double mjd, int p, double *lpd0, double *psi0, double *rp0, double *rho0, double *lam, double *bet, double *dia, double *mag) |
| void | mjd_cal (double mjd, int *mn, double *dy, int *yr) |
| double | mjd_hr (double jd) |
| double | geoc_lat (double phi) |
Variables | |
| double | ss [NARGS][30] |
| double | cc [NARGS][30] |
| plantbl | moonlr |
| plantbl | moonlat |
|
|
Definition at line 13 of file SolSystem.cxx. |
|
|
Definition at line 5945 of file SolSystem.cxx. |
|
|
Definition at line 5946 of file SolSystem.cxx. |
|
|
Definition at line 5947 of file SolSystem.cxx. |
|
|
Definition at line 2074 of file SolSystem.cxx. |
|
|
Definition at line 1871 of file SolSystem.cxx. |
|
|
Definition at line 2024 of file SolSystem.cxx. |
|
|
Definition at line 16 of file SolSystem.cxx. |
|
|
Definition at line 17 of file SolSystem.cxx. |
|
|
Definition at line 15 of file SolSystem.cxx. |
|
|
Definition at line 1109 of file SolSystem.cxx. |
|
|
Definition at line 1110 of file SolSystem.cxx. |
|
|
Definition at line 1107 of file SolSystem.cxx. |
|
|
Definition at line 1108 of file SolSystem.cxx. |
|
|
Definition at line 2070 of file SolSystem.cxx. |
|
|
Definition at line 6064 of file SolSystem.cxx. |
|
|
|
|
|
Definition at line 6063 of file SolSystem.cxx. |
|
|
Definition at line 23 of file SolSystem.cxx. |
|
|
|
|
|
Definition at line 5365 of file SolSystem.cxx. |
|
|
|
|
|
|
|
|
|
|
|
Definition at line 21 of file SolSystem.cxx. |
|
|
|
|
|
Definition at line 22 of file SolSystem.cxx. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Definition at line 20 of file SolSystem.cxx. |
|
|
Definition at line 2051 of file SolSystem.cxx. |
|
|
Definition at line 2052 of file SolSystem.cxx. |
|
|
Definition at line 2049 of file SolSystem.cxx. |
|
|
Definition at line 2026 of file SolSystem.cxx. |
|
|
Definition at line 5410 of file SolSystem.cxx. |
|
|
Definition at line 5408 of file SolSystem.cxx. |
|
|
Definition at line 5409 of file SolSystem.cxx. |
|
|
Definition at line 2071 of file SolSystem.cxx. |
|
|
Definition at line 2072 of file SolSystem.cxx. |
|
|
Definition at line 5411 of file SolSystem.cxx. |
|
|
Definition at line 24 of file SolSystem.cxx. |
|
|
Definition at line 19 of file SolSystem.cxx. |
|
|
Definition at line 1798 of file SolSystem.cxx. |
|
|
Definition at line 2073 of file SolSystem.cxx. |
|
|
Definition at line 6236 of file SolSystem.cxx. |
|
|
Definition at line 6237 of file SolSystem.cxx. |
|
|
Definition at line 6235 of file SolSystem.cxx. |
|
|
Definition at line 79 of file SolSystem.cxx. |
|
|
Definition at line 1797 of file SolSystem.cxx. |
|
|
Definition at line 912 of file SolSystem.cxx. |
|
|
Definition at line 906 of file SolSystem.cxx. |
|
|
Definition at line 910 of file SolSystem.cxx. |
|
|
Definition at line 913 of file SolSystem.cxx. |
|
|
Definition at line 909 of file SolSystem.cxx. |
|
|
|
|
||||||||||||||||
|
Definition at line 779 of file SolSystem.cxx. Referenced by SIDRATE().
00781 {
00782 double t,t0;
00783 t=(jd-A2000)/36525.;
00784 t0=fmod(6.697374558 + 2400.051336 * t + 0.000025862 * t * t,24.);
00785 if(t0<0.)t0+=24.;
00786 *tu=fmod((gst-t0),24.);
00787 if(*tu<0.)*tu+=24.;
00788 (*tu)=(*tu)*0.9972695663;
00789 return 0;
00790 }
|
|
||||||||||||||||||||||||
|
Definition at line 808 of file SolSystem.cxx. Referenced by SIDRATE(), and Sun_Pos().
00811 {
00812 double e,ras,a,c1;
00813
00814 e=0.409087723-0.069813e-9*(jd-A2000);/*obliquita eclittica*/
00815 a=cos(e)*sin(l*Conv_Rad)-tan(b*Conv_Rad)*sin(e);
00816 c1=cos(l*Conv_Rad);
00817 ras=atan(a/c1);
00818 if(l>90.)ras+=PI;
00819 if(l>270.)ras+=PI;
00820 *ra=(ras/Conv_Ora_Rad);
00821 a=sin(e)*sin(l*Conv_Rad)*cos(b*Conv_Rad);
00822 c1=sin(b*Conv_Rad)*cos(e) ;
00823 *dec=asin( a + c1 )/Conv_Rad;
00824 return 0;
00825 }
|
|
||||||||||||||||||||||||
|
Definition at line 792 of file SolSystem.cxx. Referenced by SIDRATE().
00793 {
00794 double TW=2.*PI,t,m,tl,tl0,m1;
00795 t=jd-A2000;
00796 tl=4.8949504+0.017202792*t; /*Long media sole*/
00797 while(tl<0.)tl+=TW;
00798 tl0= 6.240040768+.01720197*t; /*Anomalia media sole*/
00799 while(tl0<0.)tl0+=TW;
00800 m=tl+0.033423055*sin(tl0)+0.000349065*sin(2.*tl0); /*Oss.Long sole*/
00801 m1=fmod(m,TW);
00802 *r=1.00014-0.01671*cos(tl0)-0.00014*cos(2.*tl0);
00803 *lon=(m1 -9.93e-5)/Conv_Rad;
00804 Conv_lb_ad(jd,*lon,0.,ra,decs);
00805 return 0;
00806 }
|
|
||||||||||||
|
Definition at line 723 of file SolSystem.cxx. Referenced by SIDRATE().
00724 {
00725 double M,T,T1,Tempo_Siderale_0,Tempo_Siderale_Ora,Tempo_Siderale_Loc;
00726
00727 T = ((J_D) - A2000) / 36525.;
00728 T1 = (24110.54841 + 8640184.812866 * T + 0.0093103 * T * T)/86400.0;
00729 Tempo_Siderale_0 = modf(T1,&M) * 24.;
00730 Tempo_Siderale_Ora = Tempo_Siderale_0 + Ora_Un_Dec * 1.00273790935;
00731 if (Tempo_Siderale_Ora < 0.) Tempo_Siderale_Ora = Tempo_Siderale_Ora + 24.;
00732 if (Tempo_Siderale_Ora >= 24.) Tempo_Siderale_Ora = Tempo_Siderale_Ora - 24.;
00733 Tempo_Siderale_Loc = Tempo_Siderale_Ora;
00734 if (Tempo_Siderale_Loc < 0.) Tempo_Siderale_Loc = Tempo_Siderale_Loc + 24.;
00735 if (Tempo_Siderale_Loc >= 24.) Tempo_Siderale_Loc = Tempo_Siderale_Loc - 24.;
00736 return Tempo_Siderale_Loc;
00737 } /* Calcolo_Tempo_Siderale */
|
|
||||||||||||||||||||||||
|
Definition at line 6127 of file SolSystem.cxx. Referenced by SIDRATE().
06128 {
06129 aaha_aux (lat, az, alt, ha, dec);
06130 if (*ha > PI)
06131 *ha -= 2*PI;
06132 }
|
|
||||||||||||||||||||
|
Definition at line 5956 of file SolSystem.cxx. Referenced by SIDRATE().
05957 {
05958 ab_aux(mjd, lam, bet, lsn, AB_ECL_EOD);
05959 }
|
|
||||||||||||||||||||
|
Definition at line 5965 of file SolSystem.cxx. Referenced by SIDRATE().
05966 {
05967 ab_aux(mjd, ra, dec, lsn, AB_EQ_EOD);
05968 }
|
|
||||||||||||||||||||
|
Definition at line 1804 of file SolSystem.cxx. Referenced by SIDRATE().
01805 {
01806 double m, fea, corr;
01807
01808 if (s < 1.0) {
01809 /* elliptical */
01810 double dla;
01811
01812 m = ma-TWOPI*(long)(ma/TWOPI);
01813 if (m > PI) m -= TWOPI;
01814 if (m < -PI) m += TWOPI;
01815 fea = m;
01816
01817 for (;;) {
01818 dla = fea-(s*sin(fea))-m;
01819 if (fabs(dla)<STOPERR)
01820 break;
01821 /* avoid runnaway corrections for e>.97 and M near 0*/
01822 corr = 1-(s*cos(fea));
01823 if (corr < .1) corr = .1;
01824 dla /= corr;
01825 fea -= dla;
01826 }
01827 *nu = 2*atan(sqrt((1+s)/(1-s))*tan(fea/2));
01828 } else {
01829 /* hyperbolic */
01830 double fea1;
01831
01832 m = fabs(ma);
01833 fea = m / (s-1.);
01834 fea1 = pow(6*m/(s*s),1./3.);
01835 /* whichever is smaller is the better initial guess */
01836 if (fea1 < fea) fea = fea1;
01837
01838 corr = 1;
01839 while (fabs(corr) > STOPERR) {
01840 corr = (m - s * sinh(fea) + fea) / (s*cosh(fea) - 1);
01841 fea += corr;
01842 }
01843 if (ma < 0.) fea = -fea;
01844 *nu = 2*atan(sqrt((s+1)/(s-1))*tanh(fea/2));
01845 }
01846 *ea = fea;
01847 }
|
|
||||||||||||||||||||
|
Definition at line 1227 of file SolSystem.cxx. Referenced by SIDRATE(), mjd_year(), and year_mjd().
01228 {
01229 static double last_mjd, last_dy;
01230 static int last_mn, last_yr;
01231 int b, d, m, y;
01232 long c;
01233
01234 if (mn == last_mn && yr == last_yr && dy == last_dy) {
01235 *mjd = last_mjd;
01236 return;
01237 }
01238
01239 m = mn;
01240 y = (yr < 0) ? yr + 1 : yr;
01241 if (mn < 3) {
01242 m += 12;
01243 y -= 1;
01244 }
01245
01246 if (yr < 1582 || (yr == 1582 && (mn < 10 || (mn == 10 && dy < 15))))
01247 b = 0;
01248 else {
01249 int a;
01250 a = y/100;
01251 b = 2 - a + a/4;
01252 }
01253
01254 if (y < 0)
01255 c = (long)((365.25*y) - 0.75) - 694025L;
01256 else
01257 c = (long)(365.25*y) - 694025L;
01258
01259 d = (int)30.6001*(m+1);
01260
01261 *mjd = b + c + d + dy - 0.5;
01262
01263 last_mn = mn;
01264 last_dy = dy;
01265 last_yr = yr;
01266 last_mjd = *mjd;
01267 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 1436 of file SolSystem.cxx. Referenced by SIDRATE(), nut_eq(), plans(), and ta_par().
01437 : rectangular coordinates
01438 double *l, *b, *r; result: spherical coordinates */
01439 {
01440 double rho = x*x + y*y;
01441
01442 if (rho > 1e-35) { /* standard case: off axis */
01443 *l = atan2(y, x);
01444 range (l, 2*PI);
01445 *b = atan2(z, sqrt(rho));
01446 *r = sqrt(rho + z*z);
01447 } else { /* degenerate case; avoid math error */
01448 *l = 0.0;
01449 if (z == 0.0)
01450 *b = 0.0;
01451 else
01452 *b = (z > 0.0) ? PI/2. : -PI/2.;
01453 *r = fabs(z);
01454 }
01455 }
|
|
||||||||||||||||||||
|
Definition at line 1895 of file SolSystem.cxx. Referenced by SIDRATE().
01896 {
01897 static double a0[] = { /* semimajor axes for precision ctrl */
01898 0.39, 0.72, 1.5, 5.2, 9.6, 19.2, 30.1, 39.5, 1.0
01899 };
01900 double sum[CHAP_MAXTPOW+1][6]; /* [T^0, ..][X,Y,Z,X',Y',Z'] */
01901 double T, t; /* time in centuries and years */
01902 double ca, sa, Nu; /* aux vars for terms */
01903 double precT[CHAP_MAXTPOW+1]; /* T-augmented precision threshold */
01904 chap95_rec *rec; /* term coeffs */
01905 int cooidx;
01906
01907 /* check parameters */
01908 if (mjd < CHAP_BEGIN || mjd > CHAP_END)
01909 return (1);
01910
01911 if (obj < JUPITER || obj > PLUTO)
01912 return (2);
01913
01914 if (prec < 0.0 || prec > 1e-3)
01915 return (3);
01916
01917 /* init the sums */
01918 zero_mem ((void *)sum, sizeof(sum));
01919
01920 T = (mjd - J2000)/36525.0; /* centuries since J2000.0 */
01921
01922 /* modify precision treshold for
01923 * a) term storing scale
01924 * b) convert radians to au
01925 * c) account for skipped terms (more terms needed for better prec)
01926 * threshold empirically established similar to VSOP; stern
01927 * d) augment for secular terms
01928 */
01929 precT[0] = prec * CHAP_SCALE /* a) */
01930 * a0[obj] /* b) */
01931 / (10. * (-log10(prec + 1e-35) - 2)); /* c) */
01932 t = 1./(fabs(T) + 1e-35); /* d) */
01933 precT[1] = precT[0]*t;
01934 precT[2] = precT[1]*t;
01935
01936 t = T * 100.0; /* YEARS since J2000.0 */
01937
01938 ca = sa = Nu = 0.; /* shut up compiler warning 'uninitialised' */
01939
01940 switch (obj) { /* set initial term record pointer */
01941 case JUPITER: rec = chap95_jupiter; break;
01942 case SATURN: rec = chap95_saturn; break;
01943 case URANUS: rec = chap95_uranus; break;
01944 case NEPTUNE: rec = chap95_neptune; break;
01945 case PLUTO: rec = chap95_pluto; break;
01946 default:
01947 return (2); /* wrong object: severe internal trouble */
01948 }
01949
01950 /* do the term summation into sum[T^n] slots */
01951 for (; rec->n >= 0; ++rec) {
01952 double *amp;
01953
01954 /* NOTE: The formula
01955 * X = SUM[i=1,Records] T**n_i*(CX_i*cos(Nu_k*t)+SX_i*sin(Nu_k*t))
01956 * could be rewritten as SUM( ... A sin (B + C*t) )
01957 * "saving" trigonometric calls. However, e.g. for Pluto,
01958 * there are only 65 distinct angles NU_k (130 trig calls).
01959 * With that manipulation, EVERY arg_i would be different for X,
01960 * Y and Z, which is 3*96 terms. Hence, the formulation as
01961 * given is good (optimal?).
01962 */
01963
01964 for (cooidx = 0, amp = rec->amp; cooidx < 3; ++cooidx) {
01965 double C, S, term, termdot;
01966 short n; /* fast access */
01967
01968 C = *amp++;
01969 S = *amp++;
01970 n = rec->n;
01971
01972 /* drop term if too small
01973 * this is quite expensive: 17% of loop time
01974 */
01975 if (fabs(C) + fabs(S) < precT[n])
01976 continue;
01977
01978 if (n == 0 && cooidx == 0) { /* new Nu only here */
01979 double arg;
01980
01981 Nu = rec->Nu;
01982 arg = Nu * t;
01983 arg -= floor(arg/(2.*PI))*(2.*PI);
01984 ca = cos(arg); /* blast it - even for Nu = 0.0 */
01985 sa = sin(arg);
01986 }
01987
01988 term = C * ca + S * sa;
01989 sum[n][cooidx] += term;
01990 #if CHAP_GETRATE
01991 termdot = (-C * sa + S * ca) * Nu;
01992 sum[n][cooidx+3] += termdot;
01993 if (n > 0) sum[n - 1][cooidx+3] += n/100.0 * term;
01994 #endif
01995 } /* cooidx */
01996 } /* records */
01997
01998 /* apply powers of time and sum up */
01999 for (cooidx = 0; cooidx < 6; ++cooidx) {
02000 ret[cooidx] = (sum[0][cooidx] +
02001 T * (sum[1][cooidx] +
02002 T * (sum[2][cooidx] )) )/CHAP_SCALE;
02003 }
02004
02005 /* TEST: if the MAIN terms are dropped, get angular residue
02006 ret[0] = sqrt(ret[0]*ret[0] + ret[1]*ret[1] + ret[2]*ret[2])/a0[obj];
02007 */
02008
02009 #if CHAP_GETRATE
02010 for (cooidx = 3; cooidx < 6; ++cooidx) {
02011 ret[cooidx] /= 365.25; /* yearly to daily rate */
02012 }
02013 #endif
02014
02015 return (0);
02016 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 1551 of file SolSystem.cxx. Referenced by SIDRATE().
01552 {
01553 double w, s, s2;
01554 double l, sl, cl, y;
01555 double spsi, cpsi;
01556 double rd, lsn, rsn;
01557 double lg, re, ll;
01558 double cll, sll;
01559 double nu;
01560
01561 #define ERRLMT 0.0001
01562 w = ((mjd-ep)*3.649116e-02)/(qp*sqrt(qp));
01563 s = w/3;
01564 for (;;) {
01565 double d;
01566 s2 = s*s;
01567 d = (s2+3)*s-w;
01568 if (fabs(d) <= ERRLMT)
01569 break;
01570 s = ((2*s*s2)+w)/(3*(s2+1));
01571 }
01572
01573 nu = 2*atan(s);
01574 *rp = qp*(1+s2);
01575 l = nu+ap;
01576 sl = sin(l);
01577 cl = cos(l);
01578 spsi = sl*sin(inc);
01579 *psi = asin(spsi);
01580 y = sl*cos(inc);
01581 *lpd = atan(y/cl)+om;
01582 cpsi = cos(*psi);
01583 if (cl<0) *lpd += PI;
01584 range (lpd, 2*PI);
01585 rd = *rp * cpsi;
01586 sunpos (mjd, &lsn, &rsn, 0);
01587 lg = lsn+PI;
01588 re = rsn;
01589 ll = *lpd - lg;
01590 cll = cos(ll);
01591 sll = sin(ll);
01592 *rho = sqrt((re * re)+(*rp * *rp)-(2*re*rd*cll));
01593 if (rd<re)
01594 *lam = atan((-1*rd*sll)/(re-(rd*cll)))+lg+PI;
01595 else
01596 *lam = atan((re*sll)/(rd-(re*cll)))+*lpd;
01597 range (lam, 2*PI);
01598 *bet = atan((rd*spsi*sin(*lam-*lpd))/(cpsi*re*sll));
01599 }
|
|
|
Definition at line 6309 of file SolSystem.cxx. Referenced by SIDRATE(), and mm_mjed().
06310 {
06311 double Y;
06312 double p, B;
06313 int d[6];
06314 int i, iy, k;
06315 static double ans;
06316 static double lastmjd = -10000;
06317
06318 if (mjd == lastmjd) {
06319 return(ans);
06320 }
06321 lastmjd = mjd;
06322
06323 Y = 2000.0 + (mjd - J2000)/365.25;
06324
06325 if( Y > TABEND ) {
06326 /* linear interpolation from table end; stern */
06327 B = Y - TABEND;
06328 ans = dt[TABSIZ-1] + B * (dt[TABSIZ-1] - dt[TABSIZ-2]);
06329 ans *= 0.01;
06330 return(ans);
06331 }
06332
06333 if( Y < TABSTART ) {
06334 if( Y >= 948.0 ) {
06335 /* Stephenson and Morrison, stated domain is 948 to 1600:
06336 * 25.5(centuries from 1800)^2 - 1.9159(centuries from 1955)^2
06337 */
06338 B = 0.01*(Y - 2000.0);
06339 ans = (23.58 * B + 100.3)*B + 101.6;
06340 } else {
06341 /* Borkowski */
06342 B = 0.01*(Y - 2000.0) + 3.75;
06343 ans = 35.0 * B * B + 40.;
06344 }
06345 return(ans);
06346 }
06347
06348 /* Besselian interpolation from tabulated values.
06349 * See AA page K11.
06350 */
06351
06352 /* Index into the table.
06353 */
06354 double a=Y;
06355 p = floor(a);
06356 iy =(int)( p - TABSTART);
06357 /* Zeroth order estimate is value at start of year
06358 */
06359 ans = dt[iy];
06360 k = iy + 1;
06361 if( k >= TABSIZ )
06362 goto done; /* No data, can't go on. */
06363
06364 /* The fraction of tabulation interval
06365 */
06366 p = Y - p;
06367
06368 /* First order interpolated value
06369 */
06370 ans += p*(dt[k] - dt[iy]);
06371 if( (iy-1 < 0) || (iy+2 >= TABSIZ) )
06372 goto done; /* can't do second differences */
06373
06374 /* Make table of first differences
06375 */
06376 k = iy - 2;
06377 for( i=0; i<5; i++ ) {
06378 if( (k < 0) || (k+1 >= TABSIZ) )
06379 d[i] = 0;
06380 else d[i] = dt[k+1] - dt[k];
06381 k += 1;
06382 }
06383
06384 /* Compute second differences
06385 */
06386 for( i=0; i<4; i++ )
06387 d[i] = d[i+1] - d[i];
06388 B = 0.25*p*(p-1.0);
06389 ans += B*(d[1] + d[2]);
06390 if( iy+2 >= TABSIZ )
06391 goto done;
06392
06393 /* Compute third differences
06394 */
06395 for( i=0; i<3; i++ )
06396 d[i] = d[i+1] - d[i];
06397 B = 2.0*B/3.0;
06398 ans += (p-0.5)*B*d[1];
06399 if( (iy-2 < 0) || (iy+3 > TABSIZ) )
06400 goto done;
06401
06402 /* Compute fourth differences
06403 */
06404 for( i=0; i<2; i++ )
06405 d[i] = d[i+1] - d[i];
06406 B = 0.125*B*(p+1.0)*(p-2.0);
06407 ans += B*(d[0] + d[1]);
06408
06409 done:
06410 /* Astronomical Almanac table is corrected by adding the expression
06411 * -0.000091 (ndot + 26)(year-1955)^2 seconds
06412 * to entries prior to 1955 (AA page K8), where ndot is the secular
06413 * tidal term in the mean motion of the Moon.
06414 *
06415 * Entries after 1955 are referred to atomic time standards and
06416 * are not affected by errors in Lunar or planetary theory.
06417 */
06418 ans *= 0.01;
06419 if( Y < 1955.0 ) {
06420 B = (Y - 1955.0);
06421 ans += -0.000091 * (-25.8 + 26.0) * B * B;
06422 }
06423 return( ans );
06424 }
|
|
||||||||||||||||||||||||
|
Definition at line 6084 of file SolSystem.cxx. Referenced by SIDRATE().
06085 {
06086 ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
06087 }
|
|
||||||||||||||||||||||||
|
Definition at line 6073 of file SolSystem.cxx. Referenced by SIDRATE().
06074 {
06075 ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
06076 }
|
|
|
|
|
||||||||||||||||||||||||
|
Definition at line 6138 of file SolSystem.cxx. Referenced by SIDRATE().
06139 {
06140 aaha_aux (lat, ha, dec, az, alt);
06141 }
|
|
||||||||||||||||||||
|
Definition at line 1274 of file SolSystem.cxx. Referenced by SIDRATE(), mjd_dpm(), and mjd_year().
01275 {
01276 static double last_mjd, last_dy;
01277 static int last_mn, last_yr;
01278 double d, f;
01279 double i, a, b, ce, g;
01280
01281 /* we get called with 0 quite a bit from unused epoch fields.
01282 * 0 is noon the last day of 1899.
01283 */
01284 if (mjd == 0.0) {
01285 *mn = 12;
01286 *dy = 31.5;
01287 *yr = 1899;
01288 return;
01289 }
01290
01291 if (mjd == last_mjd) {
01292 *mn = last_mn;
01293 *yr = last_yr;
01294 *dy = last_dy;
01295 return;
01296 }
01297
01298 d = mjd + 0.5;
01299 i = floor(d);
01300 f = d-i;
01301 if (f == 1) {
01302 f = 0;
01303 i += 1;
01304 }
01305
01306 if (i > -115860.0) {
01307 a = floor((i/36524.25)+.99835726)+14;
01308 i += 1 + a - floor(a/4.0);
01309 }
01310
01311 b = floor((i/365.25)+.802601);
01312 ce = i - floor((365.25*b)+.750001)+416;
01313 g = floor(ce/30.6001);
01314 *mn = (int)(g - 1);
01315 *dy = ce - floor(30.6001*g)+f;
01316 *yr = (int)(b + 1899);
01317
01318 if (g > 13.5)
01319 *mn = (int)(g - 13);
01320 if (*mn < 2.5)
01321 *yr = (int)(b + 1900);
01322 if (*yr < 1)
01323 *yr -= 1;
01324
01325 last_mn = *mn;
01326 last_dy = *dy;
01327 last_yr = *yr;
01328 last_mjd = mjd;
01329 }
|
|
||||||||||||||||||||
|
|
|
|
Definition at line 1414 of file SolSystem.cxx. Referenced by SIDRATE(), SolSystem::ephi_riset(), and mjd_hr().
01416 {
01417 return (floor(jd-0.5)+0.5);
01418 }
|
|
||||||||||||
|
Definition at line 1336 of file SolSystem.cxx. Referenced by SIDRATE().
01338 {
01339 /* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
01340 * (Pope Gregory XIII dropped 10 days, Oct 5..14, and improved the leap-
01341 * year algorithm). however, Great Britian and the colonies did not
01342 * adopt it until Sept 14, 1752 (they dropped 11 days, Sept 3-13,
01343 * due to additional accumulated error). leap years before 1752 thus
01344 * can not easily be accounted for from the cal_mjd() number...
01345 */
01346 if (mjd < -53798.5) {
01347 /* pre sept 14, 1752 too hard to correct |:-S */
01348 return (-1);
01349 }
01350 *dow = ((long)floor(mjd-.5) + 1) % 7;/* 1/1/1900 (mjd 0.5) is a Monday*/
01351 if (*dow < 0)
01352 *dow += 7;
01353 return (0);
01354 }
|
|
||||||||||||
|
Definition at line 1358 of file SolSystem.cxx. Referenced by SIDRATE().
01359 {
01360 static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
01361 int m, y;
01362 double d;
01363
01364 mjd_cal (mjd, &m, &d, &y);
01365 *ndays = (m==2 && ((y%4==0 && y%100!=0)||y%400==0)) ? 29 : dpm[m-1];
01366 }
|
|
|
Definition at line 1421 of file SolSystem.cxx. 01422 {
01423 return ((jd-mjd_day(jd))*24.0);
01424 }
|
|
||||||||||||
|
Definition at line 1370 of file SolSystem.cxx. Referenced by SIDRATE().
01371 {
01372 static double last_mjd, last_yr;
01373 int m, y;
01374 double d;
01375 double e0, e1; /* mjd of start of this year, start of next year */
01376
01377 if (mjd == last_mjd) {
01378 *yr = last_yr;
01379 return;
01380 }
01381
01382 mjd_cal (mjd, &m, &d, &y);
01383 if (y == -1) y = -2;
01384 cal_mjd (1, 1.0, y, &e0);
01385 cal_mjd (1, 1.0, y+1, &e1);
01386 *yr = y + (mjd - e0)/(e1 - e0);
01387
01388 last_mjd = mjd;
01389 last_yr = *yr;
01390 }
|
|
|
Definition at line 773 of file SolSystem.cxx. Referenced by SIDRATE(), SolSystem::ephi_moon(), SolSystem::ephi_planet(), and SolSystem::ephi_sun().
00774 {
00775 return (mjd + deltat(mjd)/86400.0);
00776 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 5383 of file SolSystem.cxx. Referenced by SIDRATE(), and SolSystem::ephi_moon().
05384 {
05385 double pobj[3], dt;
05386 double hp;
05387
05388 if (mjd >= MOSHIER_BEGIN && mjd <= MOSHIER_END) {
05389 /* retard for light time */
05390 moon_fast (mjd, lam, bet, &hp, msp, mdp);
05391 *rho = EarthRadius/AUKM/sin(hp);
05392 dt = *rho * 5.7755183e-3; /* speed of light in a.u/day */
05393 gecmoon(mjd + MJD0 - dt, &moonlr, &moonlat, pobj);
05394
05395 *lam = pobj[0];
05396 range (lam, 2*PI);
05397 *bet = pobj[1];
05398 *rho = pobj[2];
05399 *msp = STR * Args[11]; /* don't need range correction here */
05400 *mdp = STR * Args[12];
05401 } else {
05402 moon_fast (mjd, lam, bet, &hp, msp, mdp);
05403 *rho = EarthRadius/AUKM/sin(hp);
05404
05405 }
05406 }
|
|
||||||||||||||||
|
Definition at line 745 of file SolSystem.cxx. Referenced by SIDRATE(), SolSystem::find_0alt(), and SolSystem::find_transit().
00749 {
00750 static double last_mjd = -23243, last_lng = 121212, last_lst;
00751 double eps, lst, deps, dpsi;
00752
00753 if (last_mjd == mjd && last_lng == lng) {
00754 *lstp = last_lst;
00755 return;
00756 }
00757
00758 //utc_gst (mjd_day(mjd), mjd_hr(mjd), &lst);
00759 lst=Tempo_Siderale(mjd_day(mjd)+MJD0,mjd_hr(mjd));
00760 lst += radhr(lng);
00761
00762 obliquity(mjd, &eps);
00763 nutation(mjd, &deps, &dpsi);
00764 lst += radhr(dpsi*cos(eps+deps));
00765
00766 range (&lst, 24.0);
00767
00768 last_mjd = mjd;
00769 last_lng = lng;
00770 *lstp = last_lst = lst;
00771 }
|
|
||||||||||||||||
|
Definition at line 5769 of file SolSystem.cxx. Referenced by SIDRATE().
05770 {
05771 static double lastmjd = -10000;
05772 static double a[3][3]; /* rotation matrix */
05773 double xold, yold, zold, x, y, z;
05774
05775 if (mjd != lastmjd) {
05776 double epsilon, dpsi, deps;
05777 double se, ce, sp, cp, sede, cede;
05778
05779 obliquity(mjd, &epsilon);
05780 nutation(mjd, &deps, &dpsi);
05781
05782 /* the rotation matrix a applies the nutation correction to
05783 * a vector of equatoreal coordinates Xeq to Xeq' by 3 subsequent
05784 * rotations: R1 - from equatoreal to ecliptic system by
05785 * rotation of angle epsilon about x, R2 - rotate ecliptic
05786 * system by -dpsi about its z, R3 - from ecliptic to equatoreal
05787 * by rotation of angle -(epsilon + deps)
05788 *
05789 * Xeq' = A * Xeq = R3 * R2 * R1 * Xeq
05790 *
05791 * [ 1 0 0 ]
05792 * R1 = [ 0 cos(eps) sin(eps) ]
05793 * [ 0 - sin(eps) cos(eps) ]
05794 *
05795 * [ cos(dpsi) - sin(dpsi) 0 ]
05796 * R2 = [ sin(dpsi) cos(dpsi) 0 ]
05797 * [ 0 0 1 ]
05798 *
05799 * [ 1 0 0 ]
05800 * R3 = [ 0 cos(eps + deps) - sin(eps + deps) ]
05801 * [ 0 sin(eps + deps) cos(eps + deps) ]
05802 *
05803 * for efficiency, here is a explicitely:
05804 */
05805
05806 se = sin(epsilon);
05807 ce = cos(epsilon);
05808 sp = sin(dpsi);
05809 cp = cos(dpsi);
05810 sede = sin(epsilon + deps);
05811 cede = cos(epsilon + deps);
05812
05813 a[0][0] = cp;
05814 a[0][1] = -sp*ce;
05815 a[0][2] = -sp*se;
05816
05817 a[1][0] = cede*sp;
05818 a[1][1] = cede*cp*ce+sede*se;
05819 a[1][2] = cede*cp*se-sede*ce;
05820
05821 a[2][0] = sede*sp;
05822 a[2][1] = sede*cp*ce-cede*se;
05823 a[2][2] = sede*cp*se+cede*ce;
05824
05825 lastmjd = mjd;
05826 }
05827
05828 sphcart(*ra, *dec, 1.0, &xold, &yold, &zold);
05829 x = a[0][0] * xold + a[0][1] * yold + a[0][2] * zold;
05830 y = a[1][0] * xold + a[1][1] * yold + a[1][2] * zold;
05831 z = a[2][0] * xold + a[2][1] * yold + a[2][2] * zold;
05832 cartsph(x, y, z, ra, dec, &zold); /* radius should be 1.0 */
05833 if (*ra < 0.) *ra += 2.*PI; /* make positive for display */
05834 }
|
|
||||||||||||||||
|
Definition at line 5671 of file SolSystem.cxx. Referenced by SIDRATE(), and nut_eq().
05673 : precision parameter in arc seconds
05674 double *dpsi;*/
05675 {
05676 static double lastmjd = -10000, lastdeps, lastdpsi;
05677 double T, T2, T3, T10; /* jul cent since J2000 */
05678 double prec; /* series precis in arc sec */
05679 int i, isecul; /* index in term table */
05680 static double delcache[5][2*NUT_MAXMUL+1];
05681 /* cache for multiples of delaunay args
05682 * [M',M,F,D,Om][-min*x, .. , 0, .., max*x]
05683 * make static to have unfilled fields cleared on init
05684 */
05685
05686 if (mjd == lastmjd) {
05687 *deps = lastdeps;
05688 *dpsi = lastdpsi;
05689 return;
05690 }
05691
05692 prec = 0.0;
05693
05694 #if 0 /* this is if deps should contain a precision value */
05695 prec =* deps;
05696 if (prec < 0.0 || prec > 1.0) /* accept only sane value */
05697 prec = 1.0;
05698 #endif
05699
05700 /* augment for abundance of small terms */
05701 prec *= NUT_SCALE/10;
05702
05703 T = (mjd - J2000)/36525.;
05704 T2 = T * T;
05705 T3 = T2 * T;
05706 T10 = T/10.;
05707
05708 /* calculate delaunay args and place in cache */
05709 for (i = 0; i < 5; ++i) {
05710 double x;
05711 short j;
05712
05713 x = delaunay[i][0] +
05714 delaunay[i][1] * T +
05715 delaunay[i][2] * T2 +
05716 delaunay[i][3] * T3;
05717
05718 /* convert to radians */
05719 x /= SECPERCIRC;
05720 x -= floor(x);
05721 x *= 2.*PI;
05722
05723 /* fill cache table */
05724 for (j = 0; j <= 2*NUT_MAXMUL; ++j)
05725 delcache[i][j] = (j - NUT_MAXMUL) * x;
05726 }
05727
05728 /* find dpsi and deps */
05729 lastdpsi = lastdeps = 0.;
05730 for (i = isecul = 0; i < NUT_SERIES ; ++i) {
05731 double arg = 0., ampsin, ampcos;
05732 short j;
05733
05734 if (ampconst[i][0] || ampconst[i][1]) {
05735 /* take non-secular terms from simple array */
05736 ampsin = ampconst[i][0];
05737 ampcos = ampconst[i][1];
05738 } else {
05739 /* secular terms from different array */
05740 ampsin = ampsecul[isecul][1] + ampsecul[isecul][2] * T10;
05741 ampcos = ampsecul[isecul][3] + ampsecul[isecul][4] * T10;
05742 ++isecul;
05743 }
05744
05745 for (j = 0; j < 5; ++j)
05746 arg += delcache[j][NUT_MAXMUL + multarg[i][j]];
05747
05748 if (fabs(ampsin) >= prec)
05749 lastdpsi += ampsin * sin(arg);
05750
05751 if (fabs(ampcos) >= prec)
05752 lastdeps += ampcos * cos(arg);
05753
05754 }
05755
05756 /* convert to radians.
05757 */
05758 lastdpsi = degrad(lastdpsi/3600./NUT_SCALE);
05759 lastdeps = degrad(lastdeps/3600./NUT_SCALE);
05760
05761 lastmjd = mjd;
05762 *deps = lastdeps;
05763 *dpsi = lastdpsi;
05764 }
|
|
||||||||||||
|
Definition at line 1601 of file SolSystem.cxx. Referenced by SIDRATE(), and nut_eq().
01602 {
01603 static double lastmjd = -16347, lasteps;
01604
01605 if (mjd != lastmjd) {
01606 double t = (mjd - J2000)/36525.; /* centuries from J2000 */
01607 lasteps = degrad(23.4392911 + /* 23^ 26' 21".448 */
01608 t * (-46.8150 +
01609 t * ( -0.00059 +
01610 t * ( 0.001813 )))/3600.0);
01611 lastmjd = mjd;
01612 }
01613 *eps = lasteps;
01614 }
|
|
||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 1735 of file SolSystem.cxx. Referenced by SIDRATE(), and SolSystem::ephi_planet().
01736 {
01737 static double lastmjd = -10000;
01738 static double lsn, bsn, rsn; /* geometric geocentric coords of sun */
01739 static double xsn, ysn, zsn;
01740 double lp, bp, rp; /* heliocentric coords of planet */
01741 double xp, yp, zp, rho; /* rect. coords and geocentric dist. */
01742 double dt; /* light time */
01743 int pass;
01744
01745 /* get sun cartesian; needed only once at mjd */
01746 if (mjd != lastmjd) {
01747 sunpos (mjd, &lsn, &rsn, &bsn);
01748 sphcart (lsn, bsn, rsn, &xsn, &ysn, &zsn);
01749 lastmjd = mjd;
01750 }
01751
01752 /* first find the true position of the planet at mjd.
01753 * then repeat a second time for a slightly different time based
01754 * on the position found in the first pass to account for light-travel
01755 * time.
01756 */
01757 dt = 0.0;
01758 for (pass = 0; pass < 2; pass++) {
01759 double ret[6];
01760
01761 /* get spherical coordinates of planet from precision routines,
01762 * retarded for light time in second pass;
01763 * alternative option: vsop allows calculating rates.
01764 */
01765 planpos(mjd - dt, p, 0.0, ret);
01766
01767 lp = ret[0];
01768 bp = ret[1];
01769 rp = ret[2];
01770
01771 sphcart (lp, bp, rp, &xp, &yp, &zp);
01772 cartsph (xp + xsn, yp + ysn, zp + zsn, lam, bet, &rho);
01773
01774 if (pass == 0) {
01775 /* save heliocentric coordinates at first pass since, being
01776 * true, they are NOT to be corrected for light-travel time.
01777 */
01778 *lpd0 = lp;
01779 range (lpd0, 2.*PI);
01780 *psi0 = bp;
01781 *rp0 = rp;
01782 *rho0 = rho;
01783 }
01784
01785 /* when we view a planet we see it in the position it occupied
01786 * dt days ago, where rho is the distance between it and earth,
01787 * in AU. use this as the new time for the next pass.
01788 */
01789 dt = rho * 5.7755183e-3;
01790 }
01791
01792 *dia = vis_elements[p][0];
01793 *mag = vis_elements[p][1];
01794 }
|
|
||||||||||||||||||||
|
Definition at line 1119 of file SolSystem.cxx. Referenced by SIDRATE().
01122 {
01123 precess_hiprec (mjd1, mjd2, ra, dec);
01124 }
|
|
||||||||||||
|
Definition at line 900 of file SolSystem.cxx. Referenced by SIDRATE(), car2sph(), cartsph(), comet(), SolSystem::ephi_sun(), astro::EarthOrbit::initialize(), main(), moon(), plans(), astro::EarthOrbit::position(), riset(), solve_sphere(), sunpos(), and ta_par().
00901 {
00902 *v -= r*floor(*v/r);
00903 }
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 1459 of file SolSystem.cxx. Referenced by SIDRATE().
01468 {
01469 double t0, t1;
01470 double tt, tt2, t02, tt3;
01471 double eta, th, th0;
01472 double a, b;
01473 double dap;
01474 double cinc, sinc;
01475 double ot, sot, cot, ot1;
01476 double seta, ceta;
01477
01478 if (fabs(mjd - mjd0) < 1e-5) {
01479 /* sin(eta) blows for inc < 10 degrees -- anyway, no need */
01480 *inc = inc0;
01481 *ap = ap0;
01482 *om = om0;
01483 return;
01484 }
01485
01486 t0 = mjd0/365250.0;
01487 t1 = mjd/365250.0;
01488
01489 tt = t1-t0;
01490 tt2 = tt*tt;
01491 t02 = t0*t0;
01492 tt3 = tt*tt2;
01493 eta = (471.07-6.75*t0+.57*t02)*tt+(.57*t0-3.37)*tt2+.05*tt3;
01494 th0 = 32869.0*t0+56*t02-(8694+55*t0)*tt+3*tt2;
01495 eta = degrad(eta/3600.0);
01496 th0 = degrad((th0/3600.0)+173.950833);
01497 th = (50256.41+222.29*t0+.26*t02)*tt+(111.15+.26*t0)*tt2+.1*tt3;
01498 th = th0+degrad(th/3600.0);
01499 cinc = cos(inc0);
01500 sinc = sin(inc0);
01501 ot = om0-th0;
01502 sot = sin(ot);
01503 cot = cos(ot);
01504 seta = sin(eta);
01505 ceta = cos(eta);
01506 a = sinc*sot;
01507 b = ceta*sinc*cot-seta*cinc;
01508 ot1 = atan(a/b);
01509 if (b<0) ot1 += PI;
01510 b = sinc*ceta-cinc*seta*cot;
01511 a = -1*seta*sot;
01512 dap = atan(a/b);
01513 if (b<0) dap += PI;
01514
01515 *ap = ap0+dap;
01516 range (ap, 2*PI);
01517 *om = ot1+th;
01518 range (om, 2*PI);
01519
01520 if (inc0<.175)
01521 *inc = asin(a/sin(dap));
01522 else
01523 *inc = 1.570796327-asin((cinc*ceta)+(sinc*seta*cot));
01524 }
|
|
||||||||||||||||||||
|
Definition at line 5916 of file SolSystem.cxx. Referenced by SIDRATE().
05917 {
05918 #define MAXRERR degrad(0.1/3600.) /* desired accuracy, rads */
05919
05920 double d, t, t0, a;
05921
05922 /* first guess of error is to go backwards.
05923 * make use that we know delta-apparent is always < delta-true.
05924 */
05925 unrefract (pr, tr, ta, &t);
05926 d = 0.8*(ta - t);
05927 t0 = t;
05928 a = ta;
05929
05930 /* use secant method to discover a value that unrefracts to ta.
05931 * max=7 ave=2.4 loops in hundreds of test cases.
05932 */
05933 do {
05934 a += d;
05935 unrefract (pr, tr, a, &t);
05936 d *= -(ta - t)/(t0 - t);
05937 t0 = t;
05938 } while (fabs(ta-t) > MAXRERR);
05939
05940 *aa = a;
05941
05942 #undef MAXRERR
05943 }
|
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 827 of file SolSystem.cxx. Referenced by SIDRATE(), and SolSystem::ephi_riset().
00829 {
00830 #define EPS (1e-9) /* math rounding fudge - always the way, eh? */
00831 double h; /* hour angle */
00832 double cos_h; /* cos h */
00833 double z; /* zenith angle */
00834 double zmin, zmax; /* Minimum and maximum zenith angles */
00835 double xaz, yaz; /* components of az */
00836 int shemi; /* flag for southern hemisphere reflection */
00837
00838 /* reflect lat and dec if in southern hemisphere, then az back later */
00839 if ((shemi= (lat < 0.)) != 0) {
00840 lat = -lat;
00841 dec = -dec;
00842 }
00843
00844 /* establish zenith angle, and its extrema */
00845 z = (PI/2.) + dis;
00846 zmin = fabs (dec - lat);
00847 zmax = PI - fabs(dec + lat);
00848
00849 /* first consider special cases.
00850 * these also avoid any boundary problems in subsequent computations.
00851 */
00852 if (zmax <= z + EPS) {
00853 *status = -1; /* never sets */
00854 return;
00855 }
00856 if (zmin >= z - EPS) {
00857 *status = 1; /* never rises */
00858 return;
00859 }
00860
00861 /* compute rising hour angle -- beware found off */
00862 cos_h = (cos(z)-sin(lat)*sin(dec))/(cos(lat)*cos(dec));
00863 if (cos_h >= 1.)
00864 h = 0.;
00865 else if (cos_h <= -1.)
00866 h = PI;
00867 else
00868 h = acos (cos_h);
00869
00870 /* compute setting azimuth -- beware found off */
00871 xaz = sin(dec)*cos(lat)-cos(dec)*cos(h)*sin(lat);
00872 yaz = -1.*cos(dec)*sin(h);
00873 if (xaz == 0.) {
00874 if (yaz > 0)
00875 *azs = PI/2;
00876 else
00877 *azs = -PI/2;
00878 } else
00879 *azs = atan2 (yaz, xaz);
00880
00881 /* reflect az back if southern */
00882 if (shemi)
00883 *azs = PI - *azs;
00884 range(azs, 2.*PI);
00885
00886 /* rising is just the opposite side */
00887 *azr = 2.*PI - *azs;
00888 range(azr, 2.*PI);
00889
00890 /* rise and set are just ha either side of ra */
00891 *lstr = radhr(ra-h);
00892 range(lstr,24.0);
00893 *lsts = radhr(ra+h);
00894 range(lsts,24.0);
00895
00896 /* OK */
00897 *status = 0;
00898 }
|
|
|
Definition at line 1407 of file SolSystem.cxx. Referenced by SIDRATE().
|
|
||||||||||||||||||||||||||||
|
Definition at line 6192 of file SolSystem.cxx. Referenced by SIDRATE().
06193 {
06194 double cb = cos(b), sb = sin(b);
06195 double cA = cos(A);
06196 double ca;
06197 double B;
06198
06199 ca = cb*cc + sb*sc*cA;
06200 if (ca > 1.0) ca = 1.0;
06201 if (ca < -1.0) ca = -1.0;
06202 if (cap)
06203 *cap = ca;
06204
06205 if (!Bp)
06206 return;
06207
06208 if (cc > .99999) {
06209 /* as c approaches 0, B approaches pi - A */
06210 B = PI - A;
06211 } else if (cc < -.99999) {
06212 /* as c approaches PI, B approaches A */
06213 B = A;
06214 } else {
06215 /* compute cB and sB and remove common factor of sa from quotient.
06216 * be careful where B causes atan to blow.
06217 */
06218 double sA = sin(A);
06219 double x, y;
06220
06221 y = sA*sb*sc;
06222 x = cb - ca*cc;
06223
06224 if (fabs(x) < 1e-5)
06225 B = y < 0 ? 3*PI/2 : PI/2;
06226 else
06227 B = atan2 (y, x);
06228 }
06229
06230 *Bp = B;
06231 range (Bp, 2*PI);
06232 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 1426 of file SolSystem.cxx. Referenced by SIDRATE(), nut_eq(), plans(), and ta_par().
01427 {
01428 double rcb = r * cos(b);
01429
01430 *x = rcb * cos(l);
01431 *y = rcb * sin(l);
01432 *z = r * sin(b);
01433 }
|
|
||||||||||||||||||||
|
Definition at line 1079 of file SolSystem.cxx. Referenced by SIDRATE(), comet(), SolSystem::ephi_moon(), SolSystem::ephi_planet(), SolSystem::ephi_sun(), and plans().
01080 {
01081 static double last_mjd = -3691, last_lsn, last_rsn, last_bsn;
01082 double ret[6];
01083
01084 if (mjd == last_mjd) {
01085 *lsn = last_lsn;
01086 *rsn = last_rsn;
01087 if (bsn) *bsn = last_bsn;
01088 return;
01089 }
01090
01091 vsop87(mjd, SUN, 0.0, ret); /* full precision earth pos */
01092
01093 *lsn = ret[0] - PI; /* revert to sun pos */
01094 range (lsn, 2*PI); /* normalise */
01095
01096 last_lsn = *lsn; /* memorise */
01097 last_rsn = *rsn = ret[2];
01098 last_bsn = -ret[1];
01099 last_mjd = mjd;
01100
01101 if (bsn) *bsn = last_bsn; /* assign only if non-NULL pointer */
01102
01103 }
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 5836 of file SolSystem.cxx. Referenced by SIDRATE().
05837 {
05838 static double last_phi = 1000.0, last_ht = -1000.0, xobs, zobs;
05839 double x, y, z; /* obj cartesian coord, in Earth radii */
05840
05841 /* avoid calcs involving the same phi and ht */
05842 if (phi != last_phi || ht != last_ht) {
05843 double cphi, sphi, robs, e2 = (2 - 1/298.257)/298.257;
05844 cphi = cos(phi);
05845 sphi = sin(phi);
05846 robs = 1/sqrt(1 - e2 * sphi * sphi);
05847
05848 /* observer coordinates: x to meridian, y east, z north */
05849 xobs = (robs + ht) * cphi;
05850 zobs = (robs*(1-e2) + ht) * sphi;
05851 last_phi = phi;
05852 last_ht = ht;
05853 }
05854
05855 sphcart(-tha, tdec, *rho, &x, &y, &z);
05856 cartsph(x - xobs, y, z - zobs, aha, adec, rho);
05857 *aha *= -1;
05858 range (aha, 2*PI);
05859 }
|
|
||||||||||||||||||||
|
Definition at line 5865 of file SolSystem.cxx. Referenced by SIDRATE(), and refract().
05866 {
05867 #define LTLIM 14.5
05868 #define GELIM 15.5
05869
05870 double aadeg = raddeg(aa);
05871
05872 if (aadeg < LTLIM)
05873 unrefractLT15 (pr, tr, aa, ta);
05874 else if (aadeg >= GELIM)
05875 unrefractGE15 (pr, tr, aa, ta);
05876 else {
05877 /* smooth blend -- important for inverse */
05878 double taLT, taGE, p;
05879
05880 unrefractLT15 (pr, tr, aa, &taLT);
05881 unrefractGE15 (pr, tr, aa, &taGE);
05882 p = (aadeg - LTLIM)/(GELIM - LTLIM);
05883 *ta = taLT + (taGE - taLT)*p;
05884 }
05885
05886 }
|
|
||||||||||||||||||||
|
Definition at line 971 of file SolSystem.cxx. Referenced by SIDRATE(), and sunpos().
00972 {
00973 static double (*vx_map[])[3] = { /* data tables */
00974 vx_mercury, vx_venus, vx_mars, vx_jupiter,
00975 vx_saturn, vx_uranus, vx_neptune, 0, vx_earth,
00976 };
00977 static int (*vn_map[])[3] = { /* indexes */
00978 vn_mercury, vn_venus, vn_mars, vn_jupiter,
00979 vn_saturn, vn_uranus, vn_neptune, 0, vn_earth,
00980 };
00981 static double a0[] = { /* semimajor axes; for precision ctrl only */
00982 0.39, 0.72, 1.5, 5.2, 9.6, 19.2, 30.1, 39.5, 1.0,
00983 };
00984 double (*vx_obj)[3] = vx_map[obj]; /* VSOP87 data and indexes */
00985 int (*vn_obj)[3] = vn_map[obj];
00986
00987 double t[VSOP_MAXALPHA+1]; /* powers of time */
00988 double t_abs[VSOP_MAXALPHA+1]; /* powers of abs(time) */
00989 double q; /* aux for precision control */
00990 int i, cooidx, alpha; /* misc indexes */
00991
00992 if (obj == PLUTO || obj > SUN)
00993 return (2);
00994
00995 if (prec < 0.0 || prec > 1e-3)
00996 return(3);
00997
00998 /* zero result array */
00999 for (i = 0; i < 6; ++i) ret[i] = 0.0;
01000
01001 /* time and its powers */
01002 t[0] = 1.0;
01003 t[1] = (mjd - J2000)/VSOP_A1000;
01004 for (i = 2; i <= VSOP_MAXALPHA; ++i) t[i] = t[i-1] * t[1];
01005 t_abs[0] = 1.0;
01006 for (i = 1; i <= VSOP_MAXALPHA; ++i) t_abs[i] = fabs(t[i]);
01007
01008 /* precision control */
01009 q = -log10(prec + 1e-35) - 2; /* decades below 1e-2 */
01010 q = VSOP_ASCALE * prec / 10.0 / q; /* reduce threshold progressively
01011 * for higher precision */
01012
01013 /* do the term summation; first the spatial dimensions */
01014 for (cooidx = 0; cooidx < 3; ++cooidx) {
01015
01016 /* then the powers of time */
01017 for (alpha = 0; vn_obj[alpha+1][cooidx] ; ++alpha) {
01018 double p, term, termdot;
01019
01020 /* precision threshold */
01021 p = q/(t_abs[alpha] + alpha * t_abs[alpha-1] * 1e-4 + 1e-35);
01022 #if VSOP_SPHERICAL
01023 if (cooidx == 2) /* scale by semimajor axis for radius */
01024 #endif
01025 p *= a0[obj];
01026
01027 term = termdot = 0.0;
01028 for (i = vn_obj[alpha][cooidx]; i < vn_obj[alpha+1][cooidx]; ++i) {
01029 double a, b, c, arg;
01030
01031 a = vx_obj[i][0];
01032 if (a < p) continue; /* ignore small terms */
01033
01034 b = vx_obj[i][1];
01035 c = vx_obj[i][2];
01036
01037 arg = b + c * t[1];
01038 term += a * cos(arg);
01039 #if VSOP_GETRATE
01040 termdot += -c * a * sin(arg);
01041 #endif
01042 }
01043
01044 ret[cooidx] += t[alpha] * term;
01045 #if VSOP_GETRATE
01046 ret[cooidx + 3] += t[alpha] * termdot +
01047 ((alpha > 0) ? alpha * t[alpha - 1] * term : 0.0);
01048 #endif
01049 } /* alpha */
01050 } /* cooidx */
01051
01052 for (i = 0; i < 6; ++i) ret[i] /= VSOP_ASCALE;
01053
01054 #if VSOP_SPHERICAL
01055 /* reduce longitude to 0..2pi */
01056 ret[0] -= floor(ret[0]/(2.*PI)) * (2.*PI);
01057 #endif
01058
01059 #if VSOP_GETRATE
01060 /* convert millenium rate to day rate */
01061 for (i = 3; i < 6; ++i) ret[i] /= VSOP_A1000;
01062 #endif
01063
01064 #if VSOP_SPHERICAL
01065 /* reduction from dynamical equinox of VSOP87 to FK5;
01066 */
01067 if (prec < 5e-7) { /* 5e-7 rad = 0.1 arc seconds */
01068 double L1, c1, s1;
01069 L1 = ret[0] - degrad(13.97 * t[1] - 0.031 * t[2]);
01070 c1 = cos(L1); s1 = sin(L1);
01071 ret[0] += degrad(-0.09033 + 0.03916 * (c1 + s1) * tan(ret[1]))/3600.0;
01072 ret[1] += degrad(0.03916 * (c1 - s1))/3600.0;
01073 }
01074 #endif
01075
01076 return (0);
01077 }
|
|
||||||||||||
|
Definition at line 1394 of file SolSystem.cxx. Referenced by SIDRATE().
|
|
||||||||||||
|
Definition at line 2018 of file SolSystem.cxx. Referenced by SIDRATE(), and chap95().
02019 {
02020 (void) memset (loc, 0, len);
02021 }
|
|
|
Definition at line 2064 of file SolSystem.cxx. |
|
|
Initial value: {
{ 0, 26, 29, 8, 3, 5, 0, 0, 0, 6, 5, 3, 5, 1, 0, 0, 0, 0,},
3,
bargs,
btabl,
btabb,
btabr,
0.0000000000000000e+00,
3.6525000000000000e+06,
1.0000000000000000e-04,
}Definition at line 4770 of file SolSystem.cxx. |
|
|
Initial value: {
{ 3, 26, 29, 23, 5, 10, 0, 0, 0, 8, 4, 4, 6, 2, 0, 0, 0, 0,},
3,
lrargs,
lrtabl,
lrtabb,
lrtabr,
2.5735686895300000e-03,
3.6525000000000000e+06,
1.0000000000000000e-04,
}Definition at line 4758 of file SolSystem.cxx. |
|
|
Definition at line 2063 of file SolSystem.cxx. |
1.2.11.1 written by Dimitri van Heesch,
© 1997-2001