00001
00002
00004
00005 #include <iostream>
00006 #include <cmath>
00007 #include <string>
00008 #include "vsop87_data.h"
00009 #include "chap95_data.h"
00010
00011
00012
00013 #define A2000 2451545.00000000
00014
00015 #define Conv_Rad 0.017453292
00016 #define Conv_Grad 57.29578122
00017 #define Conv_Ora_Rad 0.261799387
00018
00019 #define SPD (24.0*3600.0)
00020 #define MJD0 2415020.0
00021 #define J2000 (2451545.0 - MJD0)
00022 #define MAU (1.4959787e11)
00023 #define ERAD (6.37816e6)
00024 #define SIDRATE .9972695677
00025
00026
00027
00028 extern "C" {void precess (double mjd1, double mjd2, double *ra, double *dec);
00029 void range (double *v, double r);
00030 int Conv_GST_UT(double jd,double gst,double *tu);
00031 int Sun_Pos(double jd,double *ra,double *decs,double *lon,double *r);
00032 int Conv_lb_ad(double jd,double l,double b,double *ra,double *dec);
00033 void riset (double ra, double dec, double lat, double dis, double *lstr,double * lsts,
00034 double *azr, double *azs, int *status);
00035 int vsop87 (double mjd, int obj, double prec, double *ret);
00036 void sunpos (double mjd, double *lsn, double *rsn, double *bsn);
00037
00038 void cal_mjd (int mn, double dy, int yr, double *mjd);
00039 double mjd_day(double jd);
00040 void rnd_second (double *t);
00041 void year_mjd (double y, double *mjd);
00042 void mjd_year (double mjd, double *yr);
00043 void mjd_dpm (double mjd, int *ndays);
00044 int mjd_dow (double mjd,int *dow);
00045 void mjd_cal (double mjd, int *mn, double dy, int *yr);
00046
00047 void cartsph (double x,double y,double z,double *l, double *b, double *r);
00048 void sphcart (double l, double b, double r, double *x, double *y, double *z);
00049 void zero_mem (void *loc, unsigned len);
00050
00051 void reduce_elements (double mjd0, double mjd, double inc0, double ap0, double om0, double *inc, double *ap, double *om);
00052 void obliquity (double mjd, double *eps);
00053 void anomaly (double ma, double s, double *nu, double *ea);
00054 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);
00055 int chap95 (double mjd, int obj, double prec, double *ret);
00056 void moon (double mjd, double *lam, double *bet, double *rho,
00057 double *msp, double *mdp);
00058 void nut_eq (double mjd, double *ra, double *dec);
00059 void nutation (double mjd, double *deps, double *dpsi);
00060 void ta_par (double tha, double tdec, double phi, double ht, double *rho, double *aha, double *adec);
00061 void refract (double pr,double tr, double ta,double *aa);
00062 void ab_eq (double mjd, double lsn, double *ra, double *dec);
00063 void ab_ecl (double mjd, double lsn, double *lam, double *bet);
00064 void unrefract (double pr, double tr, double aa, double *ta);
00065
00066 void eq_ecl (double mjd, double ra, double dec, double *lat, double *lng);
00067 void ecl_eq (double mjd, double lat, double lng, double *ra, double *dec);
00068
00069 void aa_hadec (double lat, double alt,double az, double *ha, double *dec);
00070 void hadec_aa (double lat, double ha, double dec, double *alt, double *az);
00071
00072 void solve_sphere (double A, double b, double cc, double sc, double *cap, double *Bp);
00073 double deltat(double mjd);
00074 double mm_mjed (double mjd);
00075 double Tempo_Siderale(double J_D,double Ora_Un_Dec);
00076 void now_lst(double mdj,double lng, double *lstp);
00077 void plans (double mjd, int p, double *lpd0, double *psi0, double *rp0, double *rho0, double *lam, double *bet, double *dia, double *mag);
00078 }
00079 #define TMACC (10./3600./24.0)
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00095
00097
00098 SolSystem::SolSystem()
00099 {
00100 LOCSET=-1;
00101 dis=9.89e-3;
00102 lat=degrad(43.);
00103 lng=0.0;
00104 elev=0.0;
00105 GG=A2000;
00106 Obj=8;
00107 }
00108
00109 SolSystem::~SolSystem()
00110 {
00111 LOCSET=-1;
00112 }
00113
00114 void SolSystem::RiseSet()
00115 {
00116 ephi_riset ();
00117
00118 }
00119
00120 void SolSystem::RiseSet(double diso)
00121 {
00122 dis=degrad(diso);
00123 ephi_riset ();
00124
00125 }
00126
00127
00128 void SolSystem::CalculatePos(double jd, int obo)
00129 {
00130 if(LOCSET){
00131 GG=jd;
00132 if(obo>=0 && obo<=9)
00133 {
00134 Obj=obo;
00135 switch(Obj){
00136 case 8: ephi_sun();
00137 break;
00138 case 9:ephi_moon();
00139 break;
00140 default:
00141 ephi_planet();
00142 break;
00143 }
00144 }
00145 }
00146 }
00147
00148 void SolSystem::CalculatePos(double JD)
00149 {
00150 if(LOCSET){
00151
00152 SetJD(JD);
00153 switch(Obj){
00154 case 8: ephi_sun();
00155 break;
00156 case 9:ephi_moon();
00157 break;
00158 default:
00159 ephi_planet();
00160 break;
00161 }
00162 }
00163 }
00164
00165 void SolSystem::SetLocation(double lati, double loni, double eleva)
00166 {
00167 lat=degrad(lati);
00168 lng=degrad(loni);
00169 elev=eleva;
00170
00171 LOCSET=1;
00172 }
00173 void SolSystem::SetJD(double jd)
00174 {
00175 GG=jd;
00176 }
00177
00178 void SolSystem::SetObj(int nob)
00179 {
00180 if(nob>=0 && nob<=9)
00181 Obj=nob;
00182 }
00183
00184 void SolSystem::elongation (double lam, double bet, double lsn, double *el)
00185 {
00186 *el = acos(cos(bet)*cos(lam-lsn));
00187 if (lam>lsn+PI || (lam>lsn-PI && lam<lsn)) *el = - *el;
00188 }
00189
00190 void SolSystem::deflect (double mjd1, double lpd, double psi, double lsn, double rsn, double rho, double *ra,double *dec)
00191
00192
00193
00194
00195
00196 {
00197 double hra, hdec;
00198 double el;
00199 double g1, g2;
00200 double u[3];
00201 double q[3];
00202 double e[3];
00203 double qe, uq, eu;
00204 int i;
00205
00206 #define G 1.32712438e20
00207 #define c 299792458.0
00208
00209 elongation(lpd, psi, lsn-PI, &el);
00210 el = fabs(el);
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221 if (el<degrad(170) || el>degrad(179.75)) return;
00222
00223
00224 sphcart(*ra, *dec, rho, u, u+1, u+2);
00225
00226 ecl_eq(mjd1, psi, lpd, &hra, &hdec);
00227 sphcart(hra, hdec, 1.0, q, q+1, q+2);
00228
00229 ecl_eq(mjd1, 0.0, lsn-PI, &hra, &hdec);
00230 sphcart(hra, hdec, 1.0, e, e+1, e+2);
00231
00232
00233 qe = uq = eu = 0.0;
00234 for(i=0; i<=2; ++i) {
00235 qe += q[i]*e[i];
00236 uq += u[i]*q[i];
00237 eu += e[i]*u[i];
00238 }
00239
00240 g1 = 2*G/(c*c*MAU)/rsn;
00241 g2 = 1 + qe;
00242
00243
00244 g1 /= g2;
00245 for(i=0; i<=2; ++i)
00246 u[i] += g1*(uq*e[i] - eu*q[i]);
00247
00248
00249 cartsph(u[0], u[1], u[2], ra, dec, &rho);
00250 #undef G
00251 #undef c
00252 }
00253
00254
00255 int SolSystem::ephi_sun()
00256 {
00257 double lsn, rsn;
00258 double bsn;
00259 double dhlong;
00260 double mjd, mjed;
00261 double ra,dec;
00262
00263 mjd=GG-MJD0;
00264 mjed=mm_mjed(mjd);
00265
00266
00267 sunpos (mjed, &lsn, &rsn, &bsn);
00268
00269 Sdist = 0.0;
00270 Elong = 0.0;
00271 Phase = 100.0;
00272 Mag= -26.8;
00273 dhlong = lsn-PI;
00274 range (&dhlong, 2*PI);
00275 Hlong = dhlong;
00276 Hlat = -bsn;
00277 ephi_pos(mjd,SUN,bsn, lsn, &rsn, &ra,&dec);
00278
00279
00280 Edist = rsn;
00281 Size = raddeg(4.65242e-3/rsn)*3600*2;
00282 HRa=radhr(ra);
00283 DDec=raddeg(dec);
00284 Ra=ra;
00285 Dec=dec;
00286 return (0);
00287 }
00288
00289
00290
00291
00292 int SolSystem::ephi_moon()
00293
00294 {
00295 double lsn, rsn;
00296 double lam;
00297 double bet;
00298 double edistau;
00299 double el;
00300 double ms;
00301 double md;
00302 double i;
00303 double mjd, mjed;
00304 double ra,dec;
00305
00306 mjd=GG-MJD0;
00307 mjed=mm_mjed(mjd);
00308
00309 moon (mjed, &lam, &bet, &edistau, &ms, &md);
00310 sunpos (mjed, &lsn, &rsn, NULL);
00311
00312
00313
00314 Hlong=lam;
00315 Hlat=bet;
00316
00317 elongation (lam, bet, lsn, &el);
00318 Elong = raddeg(el);
00319
00320
00321 Sdist= sqrt (edistau*edistau + rsn*rsn - 2.0*edistau*rsn*cos(el));
00322
00323
00324
00325
00326
00327
00328
00329 i = 0.1468*sin(el)*(1 - 0.0549*sin(md))/(1 - 0.0167*sin(ms));
00330
00331
00332 Phase = (1+cos(PI-el-degrad(i)))/2*100;
00333
00334
00335
00336 ephi_pos(mjd,MOON,bet, lam, &edistau, &ra,&dec);
00337
00338 HRa=radhr(ra);
00339 DDec=raddeg(dec);
00340 Ra=ra;
00341 Dec=dec;
00342 Edist=edistau;
00343
00344
00345
00346 return (0);
00347 }
00348
00349 void SolSystem::ephi_planet()
00350 {
00351 double el;
00352 double f;
00353 double mjd, mjed;
00354
00355
00356 double lsn, rsn;
00357 double lpd, psi;
00358 double rp;
00359 double rho;
00360 double lam, bet;
00361 double dia, mag;
00362 double ra, dec;
00363 mjd=GG-MJD0;
00364 mjed=mm_mjed(mjd);
00365
00366 sunpos (mjed, &lsn, &rsn, 0);
00367
00368
00369
00370
00371 plans(mjed, Obj, &lpd, &psi, &rp, &rho, &lam, &bet, &dia, &mag);
00372 Size=dia;
00373 Mag=mag;
00374
00375 elongation (lam, bet, lsn, &el);
00376 Elong=el = raddeg(el);
00377
00378 f = 0.25 * ((rp+rho)*(rp+ rho) - rsn*rsn)/(rp*rho);
00379 Phase= f*100.0;
00380
00381
00382 Hlong = lpd;
00383 Hlat = psi;
00384
00385
00386
00387 ephi_pos(mjd,Obj,bet, lam, &rho, &ra,&dec);
00388 HRa=radhr(ra);
00389 DDec=raddeg(dec);
00390 Ra=ra;
00391 Dec=dec;
00392
00393 Edist = rho;
00394
00395 }
00396
00397
00398
00399 void SolSystem::ephi_pos(double mjd,int obj, double bet, double lam, double *rho, double *ra1, double *dec1)
00400
00401
00402
00403
00404
00405 {
00406 double ra, dec;
00407 double tra, tdec;
00408 double lsn, rsn;
00409 double ha_in, ha_out;
00410 double dec_out;
00411 double dra, ddec;
00412 double alt, az;
00413 double lst;
00414 double rho_topo;
00415 double temp_ra,temp_dec;
00416 double mjed;
00417
00418
00419 ecl_eq (mjd, bet, lam, &ra, &dec);
00420 tra = ra;
00421 tdec = dec;
00422
00423
00424 mjed=mm_mjed(mjd);
00425 sunpos(mjed, &lsn, &rsn, NULL);
00426
00427
00428
00429
00430
00431 if((obj!=SUN)&&(obj!=MOON))
00432 deflect (mjd, Hlong, Hlat, lsn, rsn, *rho, &ra, &dec);
00433
00434
00435 nut_eq (mjd, &ra, &dec);
00436
00437 if(obj!=MOON)
00438 ab_eq (mjd, lsn, &ra, &dec);
00439 GRa = raddeg(ra);
00440 GDec = raddeg(dec);
00441 temp_ra = ra;
00442 temp_dec = dec;
00443
00444
00445
00446
00447
00448
00449
00450 now_lst(mjd,lng,&lst);
00451 ha_in = hrrad(lst) - ra;
00452 rho_topo = *rho * MAU/ERAD;
00453 double dd=(elev/1000.)/ERAD;
00454 ta_par (ha_in, dec, lat, dd, &rho_topo, &ha_out, &dec_out);
00455
00456
00457 ha_in = hrrad(lst) - ra;
00458 hadec_aa (lat, ha_out, dec_out, &alt, &az);
00459
00460
00461 Alt = alt;
00462 Az = az;
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474 dra = ha_in - ha_out;
00475 ddec = dec_out - dec;
00476 *rho = rho_topo * ERAD/MAU;
00477
00478
00479
00480
00481 ra = ra + dra;
00482 dec = dec + ddec;
00483
00484
00485
00486
00487
00488 range(&ra, 2*PI);
00489
00490
00491 *ra1=ra;
00492 *dec1=dec;
00493
00494 }
00495
00496
00497
00498
00499
00500
00501 void SolSystem::ephi_riset ()
00502 {
00503 double mjdn;
00504 double lstn;
00505 double lr, ls;
00506 double ar, as;
00507 double ran;
00508 int rss;
00509
00510 double mjd=GG-MJD0;
00511
00512 Flags = 0;
00513
00514
00515 mjdn = mjd_day(mjd)+0.5;
00516 now_lst (mjdn,lng, &lstn);
00517
00518 SetJD(mjdn+MJD0);
00519
00520 switch(Obj){
00521 case 8: ephi_sun();
00522 break;
00523 case 9: ephi_moon();
00524 break;
00525 default:
00526 ephi_planet();
00527 break;
00528 }
00529 ran=Ra;
00530
00531
00532
00533
00534
00535
00536
00537 riset (Ra,Dec, lat, dis+0.01, &lr, &ls, &ar, &as, &rss);
00538 switch (rss) {
00539 case 0: break;
00540 case 1: Flags = RS_NEVERUP; RiseTu=0.0; RiseAz=0.0;return;
00541 case -1: Flags = RS_CIRCUMPOLAR; goto dotransit;
00542 default: Flags = RS_ERROR; RiseTu=0.0; RiseAz=0.0;return;
00543 }
00544
00545
00546
00547 double mjdt;
00548 switch (find_0alt (mjdn,(lr - lstn)/SIDRATE, dis, &ar, &mjdt)) {
00549 case 0:
00550 RiseTu = mjd_hr(mjdt);
00551 RiseAz = raddeg(ar);
00552 break;
00553 case -1:
00554 Flags |= RS_RISERR;
00555 RiseTu=0.0; RiseAz=0.0;
00556 break;
00557 case -2:
00558 case -3:
00559 Flags |= RS_NORISE;
00560 RiseTu=0.0; RiseAz=0.0;
00561 break;
00562 }
00563
00564
00565 switch (find_0alt (mjdn, (ls - lstn)/SIDRATE, dis, &ar, &mjdt)) {
00566 case 0:
00567 SetTu = mjd_hr(mjdt);
00568 SetAz = raddeg(ar);
00569 break;
00570 case -1:
00571 Flags |= RS_SETERR;
00572 SetTu=-1.0; SetAz=0.0;
00573 break;
00574 case -2:
00575 case -3:
00576 Flags |= RS_NOSET;
00577 SetTu=-1.0; SetAz=0.0;
00578 break;
00579 }
00580
00581
00582 dotransit:
00583 switch (find_transit (mjdn,(radhr(ran) - lstn)/SIDRATE, dis, &as, &mjdt)) {
00584 case 0:
00585 TranTu = mjd_hr(mjdt);
00586 TranAlt = raddeg(as);
00587 break;
00588 case -1:
00589 Flags |= RS_TRANSERR;
00590 TranTu=-1.0;TranAlt=0.0;
00591 break;
00592 case -2:
00593 Flags |= RS_NOTRANS;
00594 TranTu=-1.0;TranAlt=0.0;
00595 break;
00596 }
00597 }
00598
00599 int SolSystem::find_0alt (double mjd, double dt, double dis, double *az, double *mjdt)
00600 {
00601 #define MAXPASSES 20
00602 #define FIRSTSTEP (1.0/60.0/24.0)
00603
00604 double a0 = 0, alt=0;
00605 double mjdn = mjd;
00606 int npasses;
00607
00608 if (dt < -12.0)
00609 dt += 24.0;
00610 if (dt > 12.0)
00611 dt -= 24.0;
00612
00613
00614 dt /= 24.0;
00615
00616
00617
00618 npasses = 0;
00619 do {
00620 double a1,lst,ha;
00621
00622 mjd += dt;
00623 *mjdt=mjd;
00624
00625 SetJD(mjd+MJD0);
00626 switch(Obj){
00627 case 8: ephi_sun();
00628 break;
00629 case 9:ephi_moon();
00630 break;
00631 default:
00632 ephi_planet();
00633 break;
00634 }
00635 now_lst (mjd,lng, &lst);
00636
00637 ha= hrrad(lst) - Ra;
00638
00639 hadec_aa (lat, ha, Dec, &alt, az);
00640
00641 a1 = alt;
00642
00643 dt = (npasses == 0) ? FIRSTSTEP : (dis+a1)*dt/(a0-a1);
00644 a0 = a1;
00645
00646 } while (++npasses < MAXPASSES && fabs(dt) > TMACC);
00647
00648
00649 if (npasses == MAXPASSES)
00650 return (-3);
00651
00652 return (fabs(mjdn-mjd) < .5 ? 0 : -2);
00653
00654 #undef MAXPASSES
00655 #undef FIRSTSTEP
00656 }
00657
00658
00659
00660
00661
00662
00663
00664 int SolSystem::find_transit (double mjd, double dt, double dis, double *alt, double *mjdt)
00665 {
00666 #define MAXLOOPS 10
00667 #define MAXERR (0.25/60.)
00668 double mjdn = mjd;
00669 double lst,az;
00670 int i;
00671
00672
00673 if (dt < -12.0)
00674 dt += 24.0;
00675 if (dt > 12.0)
00676 dt -= 24.0;
00677
00678 i = 0;
00679 double ha;
00680 do {
00681 mjd += dt/24.0;
00682 SetJD(mjd+MJD0);
00683 switch(Obj){
00684 case 8: ephi_sun();
00685 break;
00686 case 9:ephi_moon();
00687 break;
00688 default:
00689 ephi_planet();
00690 break;
00691 }
00692 now_lst (mjd,lng, &lst);
00693 ha= hrrad(lst) - Ra;
00694 hadec_aa (lat, ha, Dec, alt, &az);
00695 dt = (radhr(Ra) - lst);
00696 if (dt < -12.0)
00697 dt += 24.0;
00698 if (dt > 12.0)
00699 dt -= 24.0;
00700 } while (++i < MAXLOOPS && fabs(dt) > MAXERR);
00701
00702
00703 if (i == MAXLOOPS)
00704 return (-1);
00705 *mjdt=mjd;
00706 return (fabs(mjd - mjdn) < 0.5 ? 0 : -2);
00707
00708 #undef MAXLOOPS
00709 #undef MAXERR
00710 }
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723 double Tempo_Siderale(double J_D,double Ora_Un_Dec)
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 }
00738
00739
00740
00741
00742
00743
00744
00745 void now_lst(double mjd,double lng, double *lstp)
00746
00747
00748
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
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 }
00773 double mm_mjed (double mjd)
00774 {
00775 return (mjd + deltat(mjd)/86400.0);
00776 }
00777
00778
00779 int Conv_GST_UT(double jd,double gst,double *tu)
00780
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 }
00791
00792 int Sun_Pos(double jd,double *ra,double *decs,double *lon,double *r)
00793 {
00794 double TW=2.*PI,t,m,tl,tl0,m1;
00795 t=jd-A2000;
00796 tl=4.8949504+0.017202792*t;
00797 while(tl<0.)tl+=TW;
00798 tl0= 6.240040768+.01720197*t;
00799 while(tl0<0.)tl0+=TW;
00800 m=tl+0.033423055*sin(tl0)+0.000349065*sin(2.*tl0);
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 }
00807
00808 int Conv_lb_ad(double jd,double l,double b,double *ra,double *dec)
00809
00810
00811 {
00812 double e,ras,a,c1;
00813
00814 e=0.409087723-0.069813e-9*(jd-A2000);
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 }
00826
00827 void riset (double ra, double dec, double lat, double dis, double *lstr,double * lsts,
00828 double *azr, double *azs, int *status)
00829 {
00830 #define EPS (1e-9)
00831 double h;
00832 double cos_h;
00833 double z;
00834 double zmin, zmax;
00835 double xaz, yaz;
00836 int shemi;
00837
00838
00839 if ((shemi= (lat < 0.)) != 0) {
00840 lat = -lat;
00841 dec = -dec;
00842 }
00843
00844
00845 z = (PI/2.) + dis;
00846 zmin = fabs (dec - lat);
00847 zmax = PI - fabs(dec + lat);
00848
00849
00850
00851
00852 if (zmax <= z + EPS) {
00853 *status = -1;
00854 return;
00855 }
00856 if (zmin >= z - EPS) {
00857 *status = 1;
00858 return;
00859 }
00860
00861
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
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
00882 if (shemi)
00883 *azs = PI - *azs;
00884 range(azs, 2.*PI);
00885
00886
00887 *azr = 2.*PI - *azs;
00888 range(azr, 2.*PI);
00889
00890
00891 *lstr = radhr(ra-h);
00892 range(lstr,24.0);
00893 *lsts = radhr(ra+h);
00894 range(lsts,24.0);
00895
00896
00897 *status = 0;
00898 }
00899
00900 void range (double *v, double r)
00901 {
00902 *v -= r*floor(*v/r);
00903 }
00904
00905
00906 #define VSOP_ASCALE 1e8
00907
00908
00909 #define VSOP_SPHERICAL 1
00910 #define VSOP_GETRATE 0
00911
00912 #define VSOP_A1000 365250.0
00913 #define VSOP_MAXALPHA 5
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971 int vsop87 (double mjd, int obj, double prec, double *ret)
00972 {
00973 static double (*vx_map[])[3] = {
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] = {
00978 vn_mercury, vn_venus, vn_mars, vn_jupiter,
00979 vn_saturn, vn_uranus, vn_neptune, 0, vn_earth,
00980 };
00981 static double a0[] = {
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];
00985 int (*vn_obj)[3] = vn_map[obj];
00986
00987 double t[VSOP_MAXALPHA+1];
00988 double t_abs[VSOP_MAXALPHA+1];
00989 double q;
00990 int i, cooidx, alpha;
00991
00992 if (obj == PLUTO || obj > SUN)
00993 return (2);
00994
00995 if (prec < 0.0 || prec > 1e-3)
00996 return(3);
00997
00998
00999 for (i = 0; i < 6; ++i) ret[i] = 0.0;
01000
01001
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
01009 q = -log10(prec + 1e-35) - 2;
01010 q = VSOP_ASCALE * prec / 10.0 / q;
01011
01012
01013
01014 for (cooidx = 0; cooidx < 3; ++cooidx) {
01015
01016
01017 for (alpha = 0; vn_obj[alpha+1][cooidx] ; ++alpha) {
01018 double p, term, termdot;
01019
01020
01021 p = q/(t_abs[alpha] + alpha * t_abs[alpha-1] * 1e-4 + 1e-35);
01022 #if VSOP_SPHERICAL
01023 if (cooidx == 2)
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;
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 }
01050 }
01051
01052 for (i = 0; i < 6; ++i) ret[i] /= VSOP_ASCALE;
01053
01054 #if VSOP_SPHERICAL
01055
01056 ret[0] -= floor(ret[0]/(2.*PI)) * (2.*PI);
01057 #endif
01058
01059 #if VSOP_GETRATE
01060
01061 for (i = 3; i < 6; ++i) ret[i] /= VSOP_A1000;
01062 #endif
01063
01064 #if VSOP_SPHERICAL
01065
01066
01067 if (prec < 5e-7) {
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 }
01078
01079 void sunpos (double mjd, double *lsn, double *rsn, double *bsn)
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);
01092
01093 *lsn = ret[0] - PI;
01094 range (lsn, 2*PI);
01095
01096 last_lsn = *lsn;
01097 last_rsn = *rsn = ret[2];
01098 last_bsn = -ret[1];
01099 last_mjd = mjd;
01100
01101 if (bsn) *bsn = last_bsn;
01102
01103 }
01104
01105
01106
01107 #define DCOS(x) cos(degrad(x))
01108 #define DSIN(x) sin(degrad(x))
01109 #define DASIN(x) raddeg(asin(x))
01110 #define DATAN2(y,x) raddeg(atan2((y),(x)))
01111
01112 static void precess_hiprec (double mjd1, double mjd2, double *ra,
01113 double *dec);
01114
01115
01116
01117
01118
01119 void precess (double mjd1, double mjd2, double *ra, double *dec)
01120
01121
01122 {
01123 precess_hiprec (mjd1, mjd2, ra, dec);
01124 }
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141 static void precess_hiprec (double mjd1, double mjd2, double *ra, double *dec)
01142
01143
01144 {
01145 static double last_mjd1 = -213.432, last_from;
01146 static double last_mjd2 = -213.432, last_to;
01147 double zeta_A, z_A, theta_A;
01148 double T;
01149 double A, B, C;
01150 double alpha, delta;
01151 double alpha_in, delta_in;
01152 double from_equinox, to_equinox;
01153 double alpha2000, delta2000;
01154
01155
01156
01157
01158 if (last_mjd1 == mjd1)
01159 from_equinox = last_from;
01160 else {
01161 mjd_year (mjd1, &from_equinox);
01162 last_mjd1 = mjd1;
01163 last_from = from_equinox;
01164 }
01165 if (last_mjd2 == mjd2)
01166 to_equinox = last_to;
01167 else {
01168 mjd_year (mjd2, &to_equinox);
01169 last_mjd2 = mjd2;
01170 last_to = to_equinox;
01171 }
01172
01173
01174 alpha_in = raddeg(*ra);
01175 delta_in = raddeg(*dec);
01176
01177
01178
01179 if (fabs (from_equinox-2000.0) > .02) {
01180 T = (from_equinox - 2000.0)/100.0;
01181 zeta_A = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T;
01182 z_A = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T;
01183 theta_A = 0.5567530* T - 0.0001185* T*T - 0.0000116* T*T*T;
01184
01185 A = DSIN(alpha_in - z_A) * DCOS(delta_in);
01186 B = DCOS(alpha_in - z_A) * DCOS(theta_A) * DCOS(delta_in)
01187 + DSIN(theta_A) * DSIN(delta_in);
01188 C = -DCOS(alpha_in - z_A) * DSIN(theta_A) * DCOS(delta_in)
01189 + DCOS(theta_A) * DSIN(delta_in);
01190
01191 alpha2000 = DATAN2(A,B) - zeta_A;
01192 range (&alpha2000, 360.0);
01193 delta2000 = DASIN(C);
01194 } else {
01195
01196 alpha2000 = alpha_in;
01197 delta2000 = delta_in;
01198 };
01199
01200
01201
01202 if (fabs (to_equinox - 2000.0) > .02) {
01203 T = (to_equinox - 2000.0)/100.0;
01204 zeta_A = 0.6406161* T + 0.0000839* T*T + 0.0000050* T*T*T;
01205 z_A = 0.6406161* T + 0.0003041* T*T + 0.0000051* T*T*T;
01206 theta_A = 0.5567530* T - 0.0001185* T*T - 0.0000116* T*T*T;
01207
01208 A = DSIN(alpha2000 + zeta_A) * DCOS(delta2000);
01209 B = DCOS(alpha2000 + zeta_A) * DCOS(theta_A) * DCOS(delta2000)
01210 - DSIN(theta_A) * DSIN(delta2000);
01211 C = DCOS(alpha2000 + zeta_A) * DSIN(theta_A) * DCOS(delta2000)
01212 + DCOS(theta_A) * DSIN(delta2000);
01213
01214 alpha = DATAN2(A,B) + z_A;
01215 range(&alpha, 360.0);
01216 delta = DASIN(C);
01217 } else {
01218
01219 alpha = alpha2000;
01220 delta = delta2000;
01221 };
01222
01223 *ra = degrad(alpha);
01224 *dec = degrad(delta);
01225 }
01226
01227 void cal_mjd (int mn, double dy, int yr, double *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 }
01268
01269
01270
01271
01272
01273
01274 void mjd_cal (double mjd, int *mn, double *dy, int *yr)
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
01282
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 }
01330
01331
01332
01333
01334
01335
01336 int mjd_dow (double mjd,int *dow)
01337
01338 {
01339
01340
01341
01342
01343
01344
01345
01346 if (mjd < -53798.5) {
01347
01348 return (-1);
01349 }
01350 *dow = ((long)floor(mjd-.5) + 1) % 7;
01351 if (*dow < 0)
01352 *dow += 7;
01353 return (0);
01354 }
01355
01356
01357
01358 void mjd_dpm (double mjd, int *ndays)
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 }
01367
01368
01369
01370 void mjd_year (double mjd, double *yr)
01371 {
01372 static double last_mjd, last_yr;
01373 int m, y;
01374 double d;
01375 double e0, e1;
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 }
01391
01392
01393
01394 void year_mjd (double y, double *mjd)
01395 {
01396 double e0, e1;
01397 int yf =(int) floor (y);
01398 if (yf == -1) yf = -2;
01399
01400 cal_mjd (1, 1.0, yf, &e0);
01401 cal_mjd (1, 1.0, yf+1, &e1);
01402 *mjd = e0 + (y - yf)*(e1-e0);
01403 }
01404
01405
01406
01407 void rnd_second (double *t)
01408 {
01409 *t = floor(*t*SPD+0.5)/SPD;
01410 }
01411
01412
01413
01414 double mjd_day(double jd)
01415
01416 {
01417 return (floor(jd-0.5)+0.5);
01418 }
01419
01420
01421 double mjd_hr(double jd)
01422 {
01423 return ((jd-mjd_day(jd))*24.0);
01424 }
01425
01426 void sphcart (double l, double b, double r, double *x, double *y, double *z)
01427 {
01428 double rcb = r * cos(b);
01429
01430 *x = rcb * cos(l);
01431 *y = rcb * sin(l);
01432 *z = r * sin(b);
01433 }
01434
01435
01436 void cartsph (double x,double y,double z,double *l, double *b, double *r)
01437
01438
01439 {
01440 double rho = x*x + y*y;
01441
01442 if (rho > 1e-35) {
01443 *l = atan2(y, x);
01444 range (l, 2*PI);
01445 *b = atan2(z, sqrt(rho));
01446 *r = sqrt(rho + z*z);
01447 } else {
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 }
01456
01457
01458
01459 void reduce_elements (double mjd0, double mjd, double inc0, double ap0, double om0, double *inc, double *ap, double *om)
01460
01461
01462
01463
01464
01465
01466
01467
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
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 }
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551 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)
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 }
01600
01601 void obliquity (double mjd, double *eps)
01602 {
01603 static double lastmjd = -16347, lasteps;
01604
01605 if (mjd != lastmjd) {
01606 double t = (mjd - J2000)/36525.;
01607 lasteps = degrad(23.4392911 +
01608 t * (-46.8150 +
01609 t * ( -0.00059 +
01610 t * ( 0.001813 )))/3600.0);
01611 lastmjd = mjd;
01612 }
01613 *eps = lasteps;
01614 }
01615
01616 static void pluto_ell (double mjd, double *ret);
01617 static void chap_trans (double mjd, double *ret);
01618 static void planpos (double mjd, int obj, double prec, double *ret);
01619
01620
01621 static void chap_trans (double mjd, double *ret)
01622
01623
01624 {
01625 double ra, dec, r, eps;
01626 double sr, cr, sd, cd, se, ce;
01627
01628 cartsph(ret[0], ret[1], ret[2], &ra, &dec, &r);
01629 precess(J2000, mjd, &ra, &dec);
01630 obliquity(mjd, &eps);
01631 sr = sin(ra); cr = cos(ra);
01632 sd = sin(dec); cd = cos(dec);
01633 se = sin(eps); ce = cos(eps);
01634 ret[0] = atan2( sr * ce + sd/cd * se, cr);
01635 ret[1] = asin( sd * ce - cd * se * sr);
01636 ret[2] = r;
01637 }
01638
01639
01640
01641
01642 static void pluto_ell (double mjd, double *ret)
01643
01644
01645 {
01646
01647
01648
01649 double a = 39.543,
01650 e = 0.2490,
01651 inc0 = 17.140,
01652 Om0 = 110.307,
01653 omeg0 = 113.768,
01654 mjdp = 2448045.539 - MJD0,
01655 mjdeq = J2000,
01656 n = 144.9600/36525.;
01657
01658 double inc, Om, omeg;
01659 double ma, ea, nu;
01660 double lo, slo, clo;
01661
01662 reduce_elements(mjdeq, mjd, degrad(inc0), degrad(omeg0), degrad(Om0),
01663 &inc, &omeg, &Om);
01664 ma = degrad((mjd - mjdp) * n);
01665 anomaly(ma, e, &nu, &ea);
01666 ret[2] = a * (1.0 - e*cos(ea));
01667 lo = omeg + nu;
01668 slo = sin(lo);
01669 clo = cos(lo);
01670 ret[1] = asin(slo * sin(inc));
01671 ret[0] = atan2(slo * cos(inc), clo) + Om;
01672 }
01673
01674
01675
01676
01677
01678
01679 static void planpos (double mjd, int obj, double prec, double *ret)
01680 {
01681 if (mjd >= CHAP_BEGIN && mjd <= CHAP_END) {
01682 if (obj >= JUPITER) {
01683 chap95(mjd, obj, prec, ret);
01684 chap_trans (mjd, ret);
01685 } else {
01686 vsop87(mjd, obj, prec, ret);
01687 }
01688 } else {
01689 if (obj != PLUTO) {
01690 vsop87(mjd, obj, prec, ret);
01691 } else {
01692 pluto_ell(mjd, ret);
01693 }
01694 }
01695 }
01696
01697
01698
01699
01700
01701
01702
01703 static double vis_elements[8][2] = {
01704 { 6.74, -0.42, },
01705 { 16.92, -4.34, },
01706 { 9.36, -1.20, },
01707 { 196.74, -9.4, },
01708 { 165.6, -8.88, },
01709 { 65.8, -7.19, },
01710 { 62.2, -6.87, },
01711 { 8.2, -1.0, }
01712 };
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735 void plans (double mjd, int p, double *lpd0, double *psi0, double *rp0, double *rho0, double *lam, double *bet, double *dia, double *mag)
01736 {
01737 static double lastmjd = -10000;
01738 static double lsn, bsn, rsn;
01739 static double xsn, ysn, zsn;
01740 double lp, bp, rp;
01741 double xp, yp, zp, rho;
01742 double dt;
01743 int pass;
01744
01745
01746 if (mjd != lastmjd) {
01747 sunpos (mjd, &lsn, &rsn, &bsn);
01748 sphcart (lsn, bsn, rsn, &xsn, &ysn, &zsn);
01749 lastmjd = mjd;
01750 }
01751
01752
01753
01754
01755
01756
01757 dt = 0.0;
01758 for (pass = 0; pass < 2; pass++) {
01759 double ret[6];
01760
01761
01762
01763
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
01776
01777
01778 *lpd0 = lp;
01779 range (lpd0, 2.*PI);
01780 *psi0 = bp;
01781 *rp0 = rp;
01782 *rho0 = rho;
01783 }
01784
01785
01786
01787
01788
01789 dt = rho * 5.7755183e-3;
01790 }
01791
01792 *dia = vis_elements[p][0];
01793 *mag = vis_elements[p][1];
01794 }
01795
01796
01797 #define TWOPI (2*PI)
01798 #define STOPERR (1e-8)
01799
01800
01801
01802
01803
01804 void anomaly (double ma, double s, double *nu, double *ea)
01805 {
01806 double m, fea, corr;
01807
01808 if (s < 1.0) {
01809
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
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
01830 double fea1;
01831
01832 m = fabs(ma);
01833 fea = m / (s-1.);
01834 fea1 = pow(6*m/(s*s),1./3.);
01835
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 }
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871 #define CHAP_MAXTPOW 2
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895 int chap95 (double mjd, int obj, double prec, double *ret)
01896 {
01897 static double a0[] = {
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];
01901 double T, t;
01902 double ca, sa, Nu;
01903 double precT[CHAP_MAXTPOW+1];
01904 chap95_rec *rec;
01905 int cooidx;
01906
01907
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
01918 zero_mem ((void *)sum, sizeof(sum));
01919
01920 T = (mjd - J2000)/36525.0;
01921
01922
01923
01924
01925
01926
01927
01928
01929 precT[0] = prec * CHAP_SCALE
01930 * a0[obj]
01931 / (10. * (-log10(prec + 1e-35) - 2));
01932 t = 1./(fabs(T) + 1e-35);
01933 precT[1] = precT[0]*t;
01934 precT[2] = precT[1]*t;
01935
01936 t = T * 100.0;
01937
01938 ca = sa = Nu = 0.;
01939
01940 switch (obj) {
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);
01948 }
01949
01950
01951 for (; rec->n >= 0; ++rec) {
01952 double *amp;
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964 for (cooidx = 0, amp = rec->amp; cooidx < 3; ++cooidx) {
01965 double C, S, term, termdot;
01966 short n;
01967
01968 C = *amp++;
01969 S = *amp++;
01970 n = rec->n;
01971
01972
01973
01974
01975 if (fabs(C) + fabs(S) < precT[n])
01976 continue;
01977
01978 if (n == 0 && cooidx == 0) {
01979 double arg;
01980
01981 Nu = rec->Nu;
01982 arg = Nu * t;
01983 arg -= floor(arg/(2.*PI))*(2.*PI);
01984 ca = cos(arg);
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 }
01996 }
01997
01998
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
02006
02007
02008
02009 #if CHAP_GETRATE
02010 for (cooidx = 3; cooidx < 6; ++cooidx) {
02011 ret[cooidx] /= 365.25;
02012 }
02013 #endif
02014
02015 return (0);
02016 }
02017
02018 void zero_mem (void *loc, unsigned len)
02019 {
02020 (void) memset (loc, 0, len);
02021 }
02022
02023
02024 #define CHAR short
02025
02026 #define NARGS 18
02027
02028 struct plantbl {
02029 char max_harmonic[NARGS];
02030 char max_power_of_t;
02031 CHAR *arg_tbl;
02032 long *lon_tbl;
02033 long *lat_tbl;
02034 long *rad_tbl;
02035 double distance;
02036 double timescale;
02037 double trunclvl;
02038 };
02039
02040 static double mods3600 (double x);
02041 static void mean_elements (double JED);
02042 static int sscc (int k, double arg, int n);
02043 static int g2plan (double J, struct plantbl *plan, double *pobj, int flag);
02044 static double g1plan (double J, struct plantbl *plan);
02045 static int gecmoon (double J, struct plantbl *lrtab,
02046 struct plantbl *lattab, double *pobj);
02047
02048
02049 #define MOSHIER_J2000 (2451545.0)
02050
02051 #define MOSHIER_BEGIN (1221000.5 - MJD0)
02052 #define MOSHIER_END (2798525.5 - MJD0)
02053
02054
02055 static double Args[NARGS];
02056 static double LP_equinox;
02057 static double NF_arcsec;
02058 static double Ea_arcsec;
02059 static double pA_precession;
02060
02061
02062
02063 double ss[NARGS][30];
02064 double cc[NARGS][30];
02065
02066
02067 static double T;
02068
02069
02070 #define DTR 1.7453292519943295769e-2
02071 #define RTD 5.7295779513082320877e1
02072 #define RTS 2.0626480624709635516e5
02073 #define STR 4.8481368110953599359e-6
02074 #define AUKM 1.4959787e8
02075
02076
02077 static long lrtabl[] = {
02078 175667, 66453, 5249, -42,
02079 20057, 403, -2360, 6148,
02080 -7644, 24646, -1273, 9127,
02081 -1395, 1958,
02082 232, -289,
02083 -97, 553, 69, 130,
02084 -80, 6,
02085 129, -868, 26, -89,
02086 1042, 1172, 194, -112,
02087 -47433, -241666, 224626, -103752,
02088 63419, 127606,
02089 2294, -691, -1827, -1254,
02090 -1, -119,
02091 1057, 324,
02092 505, -195, 254, -641,
02093 -36, 1008, -1082, -3,
02094 -87, 122,
02095 161, 11,
02096 2, -106,
02097 29, -123,
02098 -32, 41,
02099 -524, -35,
02100 133, -595,
02101 225, 837, -108, -191,
02102 -2294, 841, -340, -394,
02103 -351, -1039, 238, -108,
02104 -66, 21,
02105 1405, 869, 520, 2776,
02106 -174, 71,
02107 425, 652, -1260, -80,
02108 249, 77,
02109 -192, -17,
02110 -97, 134,
02111 -7, -54,
02112 -802, -7436, -2824, 70869,
02113 -35, 2481,
02114 1865, 1749, -2166, 2415,
02115 33, -183,
02116 -835, 283,
02117 27, -45,
02118 56, 235,
02119 2, 718,
02120 -1206, 275, -87, -158,
02121 -7, -2534, 0, 10774,
02122 1, -324,
02123 -208, 821,
02124 281, 1340, -797, 440,
02125 224, 72,
02126 -65, -5,
02127 -7, -44,
02128 -48, 66,
02129 -151, -40,
02130 -41, -45,
02131 76, -108,
02132 -18, 1202, 0, -2501,
02133 1438, -595, 900, 3040,
02134 -3435, -5,
02135 -100, -26,
02136 0, -13714,
02137 -183, 68,
02138 453, -83,
02139 -228, 325,
02140 97, 13,
02141 2, 105,
02142 -61, 257,
02143 0, 57,
02144 88, -11,
02145 -1, -8220,
02146 0, 275,
02147 -43, -10,
02148 -199, 105,
02149 1, -5849, 2, 24887,
02150 -128, 48,
02151 712, 970, -1407, 845,
02152 -266, 378,
02153 311, 1526, -1751, 27,
02154 0, -185858,
02155 133, 6383,
02156 -108, 25,
02157 -7, 1944,
02158 5, 390,
02159 -11, 31,
02160 277, -384, 158, 72,
02161 -81, -41, -13, -111,
02162 -2332, -65804, -698, 505812,
02163 34, 1676716, -72, -6664384,
02164 154, -57, 52, 95,
02165 -4, -5,
02166 -7, 37,
02167 -63, -32,
02168 4, 3349, 1, -14370,
02169 16, -83,
02170 0, -401,
02171 13, 3013,
02172 48, -20,
02173 0, 250,
02174 51, -79,
02175 -7, -146,
02176 148, 9,
02177 0, -64,
02178 -17, -59,
02179 -67, -492,
02180 -2, 2116601,
02181 -12, -1848,
02182 8, -436,
02183 -6, 324, 0, -1363,
02184 -163, 9,
02185 0, -74,
02186 63, 8167, -29, 37587,
02187 -22, -74501,
02188 -71, 497,
02189 -1, 551747,
02190 -87, -22,
02191 0, -51,
02192 -1, -463,
02193 0, -444,
02194 3, 89,
02195 15, -84,
02196 -36, -6829, -5, -21663,
02197 0, 86058,
02198 0, -298,
02199 -2, 751, -2, -1015,
02200 0, 69,
02201 1, -4989, 0, 21458,
02202 0, -330,
02203 0, -7,
02204 0, -226,
02205 0, -1407, 0, 2942,
02206 0, 66,
02207 0, 667,
02208 0, -155,
02209 0, 105,
02210 0, -107,
02211 0, -74,
02212 0, -52,
02213 0, 91,
02214 0, 59,
02215 0, 235,
02216 -1, -1819, 0, 2470,
02217 71, 13,
02218 0, 1026,
02219 14, -54,
02220 0, -174,
02221 -121, -19,
02222 0, -200,
02223 0, 3008,
02224 -16, -8043, -10, -37136,
02225 -3, 73724,
02226 -157, -5,
02227 0, -854,
02228 8, 147,
02229 -13, -893,
02230 0, 11869,
02231 -23, -172,
02232 89, 14,
02233 -1, 872, 0, -3744,
02234 11, 1606,
02235 0, -559,
02236 -1, -2530,
02237 0, 454,
02238 0, -193,
02239 -60, -10,
02240 -82, -13,
02241 -75, 6,
02242 36, 81,
02243 354, -162836, 148, -516569,
02244 4, 2054441,
02245 4, -94,
02246 39, 38,
02247 61, -30,
02248 2, 121,
02249 -11, 590,
02250 62, 2108,
02251 0, -12242,
02252 -476, -42,
02253 -84, 113,
02254 -394, 236,
02255 0, 276,
02256 -49, 31,
02257 0, 86,
02258 1, -1313,
02259 1, 69,
02260 -60, 88,
02261 -46, 18,
02262 0, -63818,
02263 14, -93,
02264 113, 547, -618, 17,
02265 -7, 12290, -1, -25679,
02266 0, 92,
02267 -115, 50,
02268 -48, 233,
02269 4, 1311, 1, -5567,
02270 3, 1251,
02271 29, 548,
02272 -244, 257,
02273 -2, 1825,
02274 42, 637,
02275 -46, 68,
02276 -62, 8,
02277 3, 110,
02278 445, -100, -316, -202,
02279 2925, -621, 763, 1495,
02280 -169, -184, 20, -76,
02281 -475, -138, 8, -141,
02282 -197, 1351, -1284, 422,
02283 -129, 1879, -102, 8382,
02284 -9, 45864958,
02285 -215, 1350, -1285, 422,
02286 -481, -136, 8, -140,
02287 40, -53,
02288 2622, -543, 700, 1406,
02289 402, -95, -318, -194,
02290 122, 13,
02291 -30, 147,
02292 -121, -902,
02293 61, -23,
02294 -63, 7,
02295 69, 479,
02296 -224, 228,
02297 -7, 500,
02298 0, -429,
02299 -42, 193,
02300 -92, 37,
02301 67, 5,
02302 -350, -31,
02303 0, 67,
02304 -55, -5,
02305 0, 47,
02306 -36, 53,
02307 5, 561,
02308 0, -126,
02309 0, 871,
02310 -52, 4,
02311 -201, 116922, -22, 371352,
02312 -12, -1473285,
02313 0, 87,
02314 -164, 84,
02315 -3, 422,
02316 30, 1434,
02317 -26, 38,
02318 2, -1249943,
02319 -404, -34,
02320 -57, 79,
02321 5, 509,
02322 1, 131,
02323 -344, 168,
02324 112, 22540, 30, 71218,
02325 18, -283983,
02326 0, -851,
02327 0, -1538,
02328 0, 1360,
02329 -12, 51,
02330 -48, 68,
02331 88, -20,
02332 1, 63,
02333 0, -568,
02334 303, 25,
02335 0, -122,
02336 87, 586, -606, -14,
02337 0, -100,
02338 -85, 8,
02339 -165, 54,
02340 -45, 140,
02341 0, -54,
02342 4, -831, 1, 3495,
02343 31, 116,
02344 -46, -11,
02345 -371, 190,
02346 -507, 399,
02347 -2, 57,
02348 -60, 36,
02349 -198, -1174, -613, 4988,
02350 -87, -4,
02351 141, 560, -276, 187,
02352 1876, 1379, 778, 4386,
02353 24, -15,
02354 167, -774,
02355 -71, -9,
02356 -62, 90,
02357 98, 580, -663, -7,
02358 34, -112,
02359 57, 15,
02360 -355, -214,
02361 -3240, -13605, 12229, -5723,
02362 3496, 7063,
02363 33, -51,
02364 1908, 1160, -226, 715,
02365 964, 1170, -1264, 623,
02366 14071, 5280, 5614, 3026,
02367 488, 1576, -2, 226395859,
02368 824, 1106, -1287, 617,
02369 1917, 1156, -214, 718,
02370 90, -97,
02371 12078, -2366, 3282, 6668,
02372 -219, 9179, 593, 2015,
02373 -282, -186,
02374 57, 25,
02375 31, -102,
02376 -77, -4,
02377 -268, -341, -7, -45,
02378 -3, 74,
02379 15, -615,
02380 -88, -7,
02381 234, -353,
02382 1, -119,
02383 -163, -1159, -601, 4969,
02384 22, -58,
02385 -17, -11434,
02386 17, 54,
02387 348, 348, -460, 434,
02388 -371, 175,
02389 -11, -204,
02390 4, -6440,
02391 -5, -53,
02392 -4, -14388, -37, -45231,
02393 -7, 179562,
02394 -44, 136,
02395 -160, 49,
02396 -101, 81,
02397 -1, -188,
02398 0, 2,
02399 -4, 12124, -11, -25217,
02400 71, 543, -557, -14,
02401 -75, 526,
02402 0, 395274,
02403 -233, -16,
02404 93, -20,
02405 -43, 61,
02406 0, -1275,
02407 0, -824,
02408 1, -415, 0, 1762,
02409 -261, 131,
02410 -45, 64,
02411 -297, -25,
02412 0, -17533,
02413 -6, -56,
02414 21, 1100,
02415 1, 327,
02416 1, 66,
02417 23, -217,
02418 -83, -7,
02419 83, 86847, 49, 275754,
02420 -4, -1093857,
02421 -46, 2,
02422 0, -24,
02423 0, -419,
02424 0, -5833,
02425 1, 506,
02426 0, -827,
02427 -1, -377,
02428 -11, -78,
02429 0, 131945,
02430 -2, -334,
02431 1, -75,
02432 0, -72,
02433 0, -213,
02434 -6, 5564, -2, -11618,
02435 0, 1790,
02436 0, -131,
02437 0, 6,
02438 0, -76,
02439 0, -130,
02440 0, -1115, 0, 4783,
02441 0, -195,
02442 0, -627,
02443 0, -55,
02444 0, -83,
02445 0, 163,
02446 0, -54,
02447 0, 82,
02448 0, 149,
02449 0, -754, 0, 1578,
02450 0, 138,
02451 0, 68,
02452 2, -2506, 0, 3399,
02453 0, -125,
02454 86, 16,
02455 0, -6350, 0, 27316,
02456 18, -63,
02457 0, -169,
02458 -1, 46,
02459 -136, -21,
02460 0, -239,
02461 -30, -8788, -15, -40549,
02462 -4, 80514,
02463 -46, -8,
02464 -168, -6,
02465 -1, 536, 0, -2314,
02466 9, 148,
02467 -13, -842,
02468 -1, 307713,
02469 -23, -175,
02470 95, 15,
02471 0, -297,
02472 11, 1341,
02473 0, -106,
02474 0, 5,
02475 -4, 68,
02476 -114, 10,
02477 32, 75,
02478 159, -130487, 98, -413967,
02479 2, 1647339,
02480 -4, -85,
02481 100, -46,
02482 2, 95,
02483 -11, 461,
02484 51, 1647,
02485 0, -32090,
02486 -375, -33,
02487 -65, 86,
02488 -300, 180,
02489 0, 836, 0, -3576,
02490 0, -222,
02491 0, -993,
02492 -41, 60,
02493 0, -4537,
02494 -431, -34,
02495 2, 927, 0, -1931,
02496 -79, 33,
02497 -31, 144,
02498 -1, 284, 0, -1207,
02499 0, 88,
02500 -11, 315,
02501 -178, 177,
02502 -1, 144,
02503 -58, 986,
02504 11, 86,
02505 -228, -110,
02506 2636, -494, 718, 1474,
02507 28, -35,
02508 -24, 782, -797, 277,
02509 2142, -1231, 856, 1853,
02510 74, 10797, 0, 23699298,
02511 -21, 786, -796, 277,
02512 27, -34,
02513 2615, -494, 712, 1461,
02514 -226, -109,
02515 -11, 663,
02516 0, -123,
02517 -169, 157,
02518 -54, 266,
02519 0, -76,
02520 1, -634, 0, 2738,
02521 -25, 106,
02522 -63, 24,
02523 0, -372,
02524 -221, -24,
02525 0, -5356,
02526 0, -219,
02527 0, 91,
02528 -28, 7684, -6, 24391,
02529 -1, -96795,
02530 -77, 43,
02531 2, 95,
02532 -47, -3,
02533 0, -84530,
02534 2, 310,
02535 1, 88,
02536 111, 19331, 32, 61306,
02537 4, -243595,
02538 0, 770,
02539 0, -103,
02540 0, 160,
02541 0, 356,
02542 0, 236,
02543 -41, 354,
02544 39, 303,
02545 12, -56,
02546 873, -143, 238, 482,
02547 -28, 35,
02548 -93, 31,
02549 -3, 7690211,
02550 -91, 33,
02551 -34, 43,
02552 824, -130, 226, 450,
02553 -39, 341,
02554 -1, -687,
02555 0, -303,
02556 11, -2935, 1, 12618,
02557 121, 924, 9, -1836,
02558 -268, -1144, -678, 3685,
02559 -69, -261,
02560 0, -4115951,
02561 -69, -261,
02562 5, -151,
02563 0, -88,
02564 0, 91,
02565 0, 187,
02566 0, -1281,
02567 1, 77,
02568 1, 6059, 3, 19238,
02569 0, -76305,
02570 0, -90,
02571 0, -238,
02572 0, -962, 0, 4133,
02573 0, 96,
02574 0, 9483,
02575 0, 85,
02576 0, -688,
02577 0, -5607,
02578 0, 55,
02579 0, -752,
02580 0, 71,
02581 0, 303,
02582 0, -288,
02583 0, 57,
02584 0, 45,
02585 0, 189,
02586 0, 401,
02587 0, -1474, 0, 3087,
02588 0, -71,
02589 0, 2925,
02590 0, -75,
02591 0, 359,
02592 0, 55,
02593 1, -10155, 0, 43735,
02594 0, -572,
02595 0, -49,
02596 0, -660,
02597 0, -3591, 0, 7516,
02598 0, 668,
02599 -1, -53,
02600 -2, 384259,
02601 0, -163,
02602 0, -93,
02603 1, 112,
02604 -95, -11528, -22, -36505,
02605 -1, 145308,
02606 5, 145,
02607 0, 4047,
02608 1, 1483, 0, -6352,
02609 0, 991, 0, -4262,
02610 0, -93,
02611 0, -334,
02612 0, -160,
02613 0, -153,
02614 -10, 127,
02615 51, 185,
02616 -77, 18,
02617 56, 1217, 6, 1919574,
02618 -74, 17,
02619 50, 180,
02620 -5, 93,
02621 0, -104,
02622 0, -58,
02623 -3, -353, -1, 1499,
02624 0, -229,
02625 -15, 86,
02626 0, -93657,
02627 0, 1561, 0, -6693,
02628 0, -5839,
02629 1, 6791, 0, -29143,
02630 1, -701, 0, 3015,
02631 0, 2543,
02632 0, 693,
02633 -1, 361233,
02634 0, -50,
02635 0, 946,
02636 -1, -140,
02637 -70, 407,
02638 0, -450995,
02639 0, -368,
02640 0, 54,
02641 0, -802,
02642 0, -96,
02643 0, 1274, 0, -5459,
02644 0, -614, 0, 2633,
02645 0, 685,
02646 0, -915,
02647 0, -85,
02648 0, 88,
02649 0, 106,
02650 0, 928,
02651 0, -726, 0, 1523,
02652 0, 5715,
02653 0, -4338, 0, 18706,
02654 0, -135,
02655 0, -132,
02656 0, -158,
02657 0, -98,
02658 0, 680,
02659 -1, 138968,
02660 0, -192,
02661 0, -1698,
02662 0, -2734, 0, 11769,
02663 0, 4,
02664 0, 673, 0, -2891,
02665 0, 889, 0, -3821,
02666 0, 121,
02667 -1, 143783,
02668 0, 231,
02669 -9, 51,
02670 0, -57413,
02671 0, -483,
02672 0, -407,
02673 0, 676, 0, -2902,
02674 0, 531,
02675 0, 445,
02676 0, 672,
02677 0, 19336,
02678 0, 70,
02679 0, -39976,
02680 0, -68,
02681 0, 4203,
02682 0, -406,
02683 0, 446,
02684 0, -108,
02685 0, 79,
02686 0, 84,
02687 0, 734,
02688 0, 255,
02689 0, 3944,
02690 0, -655, 0, 2825,
02691 0, -109,
02692 0, -234,
02693 0, 57,
02694 0, 19773,
02695 0, -2013,
02696 0, 958,
02697 0, -521,
02698 0, -757,
02699 0, 10594,
02700 0, -9901,
02701 0, 199,
02702 0, -275,
02703 0, 64,
02704 0, 54,
02705 0, 165,
02706 0, 1110,
02707 0, -3286,
02708 0, 909,
02709 0, 54,
02710 0, 87,
02711 0, 258,
02712 0, 1261,
02713 0, -51,
02714 0, 336,
02715 0, -114,
02716 0, 2185,
02717 0, -850,
02718 0, 75,
02719 0, -69,
02720 0, -103,
02721 0, 776,
02722 0, -1238,
02723 0, 137,
02724 0, 67,
02725 0, -260,
02726 0, 130,
02727 0, 49,
02728 0, 228,
02729 0, 215,
02730 0, -178,
02731 0, 57,
02732 0, -133,
02733 };
02734 static long lrtabb[] = {-1};
02735 static long lrtabr[] = {
02736 -5422, -2120, 1077, 772,
02737 39, 75, 3, 10,
02738 -468, -326, -113, -78,
02739 -4, -2,
02740 1, 3,
02741 29, 24, 4, 2,
02742 1, 0,
02743 -9, 7, -2, 0,
02744 -32, -13, -3, -3,
02745 233, 126, 89, 77,
02746 -33, 16,
02747 3, -3, 0, -1,
02748 2, 0,
02749 0, 1,
02750 4, 9, 1, 1,
02751 16, -1, 0, 18,
02752 3, 2,
02753 0, 0,
02754 0, 0,
02755 0, 0,
02756 0, 0,
02757 0, -1,
02758 -22, -5,
02759 10, 3, 1, 1,
02760 -15, 7, -2, 1,
02761 -8, -11, -1, -2,
02762 -1, 1,
02763 46, -58, 126, -23,
02764 4, 8,
02765 35, 8, 10, -17,
02766 0, 0,
02767 0, 0,
02768 -10, -7,
02769 0, 0,
02770 -23, 3, 151, 10,
02771 -327, 0,
02772 4, -5, 6, 5,
02773 1, 0,
02774 -1, -3,
02775 0, 0,
02776 0, 1,
02777 -185, 0,
02778 -3, -24, -5, -2,
02779 -1062, 3, 4560, 0,
02780 -3, 0,
02781 4, 1,
02782 8, -1, 2, 4,
02783 0, 1,
02784 0, -1,
02785 0, 0,
02786 -1, 0,
02787 0, 1,
02788 0, 0,
02789 -1, -1,
02790 277, 3, -583, 1,
02791 -1, 4, -32, 7,
02792 0, -34,
02793 1, -1,
02794 -23685, 0,
02795 -1, -2,
02796 -1, -7,
02797 -5, -4,
02798 0, 2,
02799 -2, 0,
02800 -5, -1,
02801 35, 0,
02802 0, 2,
02803 202, 0,
02804 180, 0,
02805 0, -1,
02806 -3, -6,
02807 -193, 0, 770, -1,
02808 -2, -4,
02809 -32, 23, -28, -46,
02810 -13, -9,
02811 -54, 10, -1, -61,
02812 -44895, 0,
02813 -230, 5,
02814 -1, -4,
02815 -71, 0,
02816 -15, 0,
02817 1, 0,
02818 15, 11, -3, 6,
02819 2, -3, 4, -1,
02820 2576, -138, -19881, -47,
02821 -65906, -1, 261925, -4,
02822 -2, -7, 4, -2,
02823 0, 0,
02824 -1, 0,
02825 1, -3,
02826 172, -2, -727, 0,
02827 4, 1,
02828 324, 0,
02829 -139, 1,
02830 1, 3,
02831 -276, 0,
02832 5, 3,
02833 9, 0,
02834 -1, 10,
02835 -37, 0,
02836 5, -1,
02837 76, -10,
02838 1318810, 1,
02839 12, -1,
02840 -38, 1,
02841 -141, 0, 611, 0,
02842 0, -11,
02843 4, 0,
02844 -627, 2, -2882, -3,
02845 5711, -2,
02846 -48, -7,
02847 55294, 0,
02848 2, -7,
02849 31, 0,
02850 34, 0,
02851 -259, 0,
02852 -55, 2,
02853 6, 3,
02854 -4273, 20, -13554, 3,
02855 53878, 0,
02856 -46, 0,
02857 -85, 0, 114, 0,
02858 -45, 0,
02859 -818, 0, 3520, 0,
02860 34, 0,
02861 -157, 0,
02862 29, 0,
02863 -878, 0, 1838, 0,
02864 -428, 0,
02865 161, 0,
02866 24, 0,
02867 65, 0,
02868 19, 0,
02869 15, 0,
02870 12, 0,
02871 -26, 0,
02872 -14, 0,
02873 -149, 0,
02874 584, 0, -793, 0,
02875 4, -23,
02876 -238, 0,
02877 -18, -5,
02878 45, 0,
02879 -7, 42,
02880 79, 0,
02881 -1723, 0,
02882 2895, -6, 13362, -4,
02883 -26525, -1,
02884 -2, 57,
02885 291, 0,
02886 52, -3,
02887 -327, 5,
02888 -2755, 0,
02889 -63, 9,
02890 5, -33,
02891 -261, -1, 1122, 0,
02892 621, -4,
02893 -227, 0,
02894 1077, 0,
02895 -167, 0,
02896 85, 0,
02897 -4, 23,
02898 -5, 32,
02899 3, 30,
02900 -32, 14,
02901 64607, 141, 204958, 59,
02902 -815115, 2,
02903 -37, -1,
02904 15, -15,
02905 12, 24,
02906 48, -1,
02907 235, 4,
02908 843, -25,
02909 4621, 0,
02910 -17, 191,
02911 45, 34,
02912 95, 159,
02913 -132, 0,
02914 13, 20,
02915 32, 0,
02916 -540, 0,
02917 29, 0,
02918 37, 25,
02919 8, 19,
02920 22127, 0,
02921 -35, -5,
02922 232, -48, 7, 262,
02923 5428, 3, -11342, 1,
02924 -45, 0,
02925 -21, -49,
02926 -100, -21,
02927 -626, 1, 2665, 0,
02928 532, -2,
02929 235, -12,
02930 -111, -105,
02931 774, 1,
02932 -283, 17,
02933 29, 20,
02934 3, 27,
02935 47, -2,
02936 -43, -192, -87, 136,
02937 -269, -1264, 646, -330,
02938 -79, 73, -33, -9,
02939 60, -205, 61, 4,
02940 -584, -85, -182, -555,
02941 -780, -57, -3488, -45,
02942 -19818328, -4,
02943 583, 93, 182, 555,
02944 -59, 208, -60, -4,
02945 23, 17,
02946 235, 1133, -608, 302,
02947 41, 174, 84, -137,
02948 6, -53,
02949 63, 13,
02950 -392, 52,
02951 -10, -27,
02952 -3, -27,
02953 199, -31,
02954 99, 97,
02955 -218, -3,
02956 209, 0,
02957 84, 18,
02958 16, 40,
02959 2, -30,
02960 14, -154,
02961 30, 0,
02962 -2, 24,
02963 -108, 0,
02964 -24, -16,
02965 262, -2,
02966 55, 0,
02967 -304, 0,
02968 2, 25,
02969 55112, 95, 175036, 11,
02970 -694477, 5,
02971 41, 0,
02972 -38, -76,
02973 199, 1,
02974 679, -14,
02975 -17, -12,
02976 582619, 1,
02977 -16, 191,
02978 38, 27,
02979 -234, 2,
02980 -60, 0,
02981 80, 163,
02982 -10296, 48, -32526, 13,
02983 129703, 8,
02984 -1366, 0,
02985 -741, 0,
02986 -646, 0,
02987 25, 6,
02988 33, 23,
02989 10, 43,
02990 -31, 0,
02991 -6, 0,
02992 -12, 147,
02993 59, 0,
02994 287, -42, -7, 297,
02995 -59, 0,
02996 -4, -42,
02997 -27, -81,
02998 -69, -22,
02999 27, 0,
03000 -423, -2, 1779, -1,
03001 -57, 15,
03002 5, -23,
03003 94, 182,
03004 -197, -250,
03005 24, 1,
03006 -18, -30,
03007 581, -98, -2473, -303,
03008 -2, 43,
03009 -277, 70, -92, -136,
03010 -681, 925, -2165, 384,
03011 -8, -12,
03012 382, 82,
03013 -4, 35,
03014 -45, -31,
03015 -286, 48, 3, -328,
03016 -55, -17,
03017 8, -28,
03018 -106, 175,
03019 -6735, 1601, -2832, -6052,
03020 3495, -1730,
03021 -25, -17,
03022 -574, 944, -354, -112,
03023 -579, 476, -308, -625,
03024 -2411, 7074, -1529, 2828,
03025 -1335, 247,-112000844, -1,
03026 545, -409, 305, 637,
03027 572, -950, 356, 106,
03028 48, 44,
03029 1170, 5974, -3298, 1624,
03030 -4538, -106, -996, 294,
03031 92, -139,
03032 -12, 28,
03033 50, 16,
03034 2, -38,
03035 169, -133, 22, -3,
03036 38, 1,
03037 305, 7,
03038 4, -44,
03039 175, 116,
03040 59, 1,
03041 -573, 81, 2453, 297,
03042 29, 11,
03043 5674, -8,
03044 -27, 9,
03045 173, -173, 215, 228,
03046 -87, -184,
03047 102, -5,
03048 3206, 2,
03049 -53, 2,
03050 7159, -7, 22505, -19,
03051 -89344, -3,
03052 67, 22,
03053 24, 79,
03054 -40, -50,
03055 94, 0,
03056 186, 0,
03057 -6063, 0, 12612, -5,
03058 -271, 35, 7, -278,
03059 -479, -74,
03060 426754, 0,
03061 8, -116,
03062 -10, -47,
03063 -31, -22,
03064 645, 0,
03065 426, 0,
03066 -213, 0, 903, 0,
03067 -67, -133,
03068 -33, -23,
03069 13, -152,
03070 -9316, 0,
03071 29, -3,
03072 -564, 11,
03073 -167, 0,
03074 -34, 0,
03075 114, 12,
03076 4, -44,
03077 -44561, 42, -141493, 25,
03078 561256, -2,
03079 -1, -24,
03080 -261, 0,
03081 211, 0,
03082 -4263, 0,
03083 -262, 1,
03084 1842, 0,
03085 202, 0,
03086 41, -6,
03087 77165, 0,
03088 176, -1,
03089 39, 1,
03090 -24, 0,
03091 118, 0,
03092 -2991, -4, 6245, -1,
03093 46886, 0,
03094 -75, 0,
03095 -100, 0,
03096 40, 0,
03097 75, 0,
03098 -618, 0, 2652, 0,
03099 112, 0,
03100 1780, 0,
03101 30, 0,
03102 49, 0,
03103 86, 0,
03104 33, 0,
03105 -30, 0,
03106 -95, 0,
03107 277, 0, -580, 0,
03108 -35, 0,
03109 -319, 0,
03110 1622, 1, -2201, 0,
03111 79, 0,
03112 10, -57,
03113 2363, 0, -10162, 0,
03114 -41, -12,
03115 62, 0,
03116 30, 1,
03117 -14, 89,
03118 -2721, 0,
03119 5780, -19, 26674, -10,
03120 -52964, -2,
03121 -5, 30,
03122 -4, 111,
03123 -317, -1, 1369, 0,
03124 93, -6,
03125 -564, 9,
03126 -115913, 0,
03127 -113, 15,
03128 10, -62,
03129 99, 0,
03130 891, -7,
03131 36, 0,
03132 108, 0,
03133 -42, -2,
03134 7, 75,
03135 -50, 21,
03136 86822, 104, 275441, 65,
03137 -1096109, 1,
03138 -56, 3,
03139 31, 66,
03140 63, -1,
03141 307, 7,
03142 1097, -34,
03143 17453, 0,
03144 -22, 250,
03145 57, 43,
03146 120, 200,
03147 -297, 0, 1269, 0,
03148 166, 0,
03149 -662, 0,
03150 40, 28,
03151 1521, 0,
03152 -23, 288,
03153 351, -2, -729, 0,
03154 -22, -52,
03155 -96, -21,
03156 -139, -1, 589, 0,
03157 35, 0,
03158 210, 7,
03159 -118, -119,
03160 62, 0,
03161 -583, -26,
03162 -42, 5,
03163 -73, 152,
03164 -330, -1759, 983, -479,
03165 -23, -19,
03166 -522, -15, -185, -533,
03167 739, 1559, -1300, 614,
03168 -7332, 52, -15836758, 0,
03169 524, 16, 185, 532,
03170 23, 18,
03171 330, 1751, -978, 476,
03172 73, -151,
03173 519, 18,
03174 38, 0,
03175 105, 113,
03176 -178, -37,
03177 26, 0,
03178 262, 1, -1139, 0,
03179 71, 17,
03180 16, 42,
03181 151, 0,
03182 16, -148,
03183 4147, 0,
03184 149, 0,
03185 -30, 0,
03186 2980, 9, 9454, 2,
03187 -37519, 0,
03188 -28, -49,
03189 37, -1,
03190 2, -31,
03191 33870, 0,
03192 -208, 1,
03193 -59, 1,
03194 -13105, 68, -41564, 21,
03195 165148, 3,
03196 -1022, 0,
03197 -40, 0,
03198 -132, 0,
03199 -228, 0,
03200 95, 0,
03201 -138, -16,
03202 -126, 16,
03203 24, 5,
03204 -57, -346, 191, -94,
03205 -14, -11,
03206 -12, -37,
03207 -3053364, -1,
03208 13, 36,
03209 17, 13,
03210 51, 327, -179, 90,
03211 138, 16,
03212 233, 0,
03213 62, 0,
03214 1164, 0, -5000, 0,
03215 -407, 117, 770, 9,
03216 -4, 1, 21, 2,
03217 1, 0,
03218 -16869, 0,
03219 -1, 0,
03220 1, 0,
03221 35, 0,
03222 -78, 0,
03223 78, 0,
03224 -533, 0,
03225 -31, 1,
03226 -2448, -6, -7768, -1,
03227 30812, 0,
03228 37, 0,
03229 -227, 0,
03230 197, 0, -846, 0,
03231 -77, 0,
03232 4171, 0,
03233 -67, 0,
03234 287, 0,
03235 2532, 0,
03236 -19, 0,
03237 -40, 0,
03238 -56, 0,
03239 128, 0,
03240 83, 0,
03241 -45, 0,
03242 -36, 0,
03243 -92, 0,
03244 -134, 0,
03245 714, 0, -1495, 0,
03246 32, 0,
03247 -981, 0,
03248 15, 0,
03249 -166, 0,
03250 -59, 0,
03251 4923, 0, -21203, 0,
03252 246, 0,
03253 15, 0,
03254 104, 0,
03255 1683, 0, -3523, 0,
03256 -865, 0,
03257 -25, 1,
03258 -186329, -1,
03259 10, 0,
03260 50, 0,
03261 53, 0,
03262 5455, -45, 17271, -10,
03263 -68747, 0,
03264 69, -2,
03265 -7604, 0,
03266 -724, 1, 3101, 0,
03267 -46, 0, 200, 0,
03268 -44, 0,
03269 97, 0,
03270 -53, 0,
03271 62, 0,
03272 -54, -4,
03273 88, -24,
03274 -9, -36,
03275 -581, 27, -914711, 3,
03276 8, 35,
03277 -86, 24,
03278 51, 3,
03279 48, 0,
03280 26, 0,
03281 133, 1, -577, 0,
03282 105, 0,
03283 -3, -1,
03284 3194, 0,
03285 528, 0, -2263, 0,
03286 2028, 0,
03287 -3266, 1, 14016, 0,
03288 10, 0, -41, 0,
03289 -100, 0,
03290 -32, 0,
03291 -124348, 0,
03292 16, 0,
03293 -325, 0,
03294 50, -1,
03295 1, 0,
03296 -553, 0,
03297 0, 0,
03298 0, 0,
03299 2, 0,
03300 -34, 0,
03301 -444, 0, 1902, 0,
03302 9, 0, -37, 0,
03303 254, 0,
03304 156, 0,
03305 -2, 0,
03306 -35, 0,
03307 -48, 0,
03308 -368, 0,
03309 327, 0, -686, 0,
03310 -2263, 0,
03311 1952, 0, -8418, 0,
03312 -13, 0,
03313 52, 0,
03314 9, 0,
03315 21, 0,
03316 -261, 0,
03317 -62404, 0,
03318 0, 0,
03319 79, 0,
03320 1056, 0, -4547, 0,
03321 -351, 0,
03322 -305, 0, 1310, 0,
03323 -1, 0, 6, 0,
03324 0, 0,
03325 -55953, 0,
03326 -80, 0,
03327 0, 0,
03328 168, 0,
03329 -147, 0,
03330 127, 0,
03331 -265, 0, 1138, 0,
03332 -1, 0,
03333 -9, 0,
03334 -8, 0,
03335 -5984, 0,
03336 -22, 0,
03337 -5, 0,
03338 0, 0,
03339 0, 0,
03340 127, 0,
03341 -2, 0,
03342 10, 0,
03343 -31, 0,
03344 -29, 0,
03345 -286, 0,
03346 -98, 0,
03347 -1535, 0,
03348 252, 0, -1087, 0,
03349 43, 0,
03350 4, 0,
03351 -19, 0,
03352 -7620, 0,
03353 29, 0,
03354 -322, 0,
03355 203, 0,
03356 0, 0,
03357 -3587, 0,
03358 10, 0,
03359 0, 0,
03360 94, 0,
03361 0, 0,
03362 -1, 0,
03363 -1, 0,
03364 -315, 0,
03365 1, 0,
03366 0, 0,
03367 0, 0,
03368 -30, 0,
03369 -94, 0,
03370 -460, 0,
03371 1, 0,
03372 -114, 0,
03373 0, 0,
03374 -746, 0,
03375 4, 0,
03376 -23, 0,
03377 24, 0,
03378 0, 0,
03379 -237, 0,
03380 1, 0,
03381 0, 0,
03382 -18, 0,
03383 0, 0,
03384 0, 0,
03385 -16, 0,
03386 -76, 0,
03387 -67, 0,
03388 0, 0,
03389 -16, 0,
03390 0, 0,
03391 };
03392 static CHAR lrargs[] = {
03393 0, 3,
03394 3, 4, 3, -8, 4, 3, 5, 1,
03395 2, 2, 5, -5, 6, 2,
03396 5, -1, 10, 2, 13, -1, 11, 3, 3, -7, 4, 0,
03397 3, 1, 13, -1, 11, 2, 5, 1,
03398 2, 4, 5,-10, 6, 0,
03399 4, 2, 10, -2, 13, 14, 3,-23, 4, 1,
03400 3, 3, 2, -7, 3, 4, 4, 1,
03401 3, -1, 13, 18, 2,-16, 3, 2,
03402 2, 8, 2,-13, 3, 1,
03403 5, 2, 10, -2, 13, 2, 3, -3, 5, 1, 6, 0,
03404 3, -1, 13, 26, 2,-29, 3, 0,
03405 3, 1, 10, -1, 11, 2, 4, 1,
03406 4, 1, 10, -1, 13, 3, 2, -4, 3, 1,
03407 4, 1, 10, -1, 13, 3, 3, -4, 4, 0,
03408 3, -1, 10, 15, 2,-12, 3, 0,
03409 4, 2, 10, -3, 13, 24, 2,-24, 3, 0,
03410 3, -1, 10, 23, 2,-25, 3, 0,
03411 4, 1, 10, -1, 11, 1, 3, 1, 6, 0,
03412 4, 2, 10, -2, 11, 5, 2, -6, 3, 0,
03413 4, 2, 10, -2, 13, 6, 2, -8, 3, 0,
03414 4, -2, 10, 1, 13, 12, 2, -8, 3, 1,
03415 5, -1, 10, 1, 13, -1, 11, 20, 2,-20, 3, 1,
03416 4, -2, 10, 1, 13, 3, 1, -1, 3, 1,
03417 5, 2, 10, -2, 13, 2, 3, -5, 5, 5, 6, 0,
03418 4, 2, 10, -2, 13, 2, 3, -3, 5, 1,
03419 4, 2, 10, -2, 13, 6, 3, -8, 4, 0,
03420 4, -2, 10, 1, 13, 20, 2,-21, 3, 1,
03421 4, 1, 10, -1, 11, 1, 3, 1, 5, 0,
03422 1, 1, 6, 0,
03423 4, 2, 10, -2, 13, 5, 3, -6, 4, 0,
03424 3, 3, 2, -5, 3, 2, 5, 0,
03425 2, -1, 11, 1, 14, 1,
03426 4, 2, 10, -2, 13, 2, 3, -2, 5, 0,
03427 2, 1, 3, -2, 4, 1,
03428 4, 1, 10, -1, 11, 5, 2, -7, 3, 0,
03429 1, 1, 5, 0,
03430 2, 7, 3,-13, 4, 0,
03431 4, -2, 10, 1, 13, 15, 2,-13, 3, 0,
03432 4, 2, 10, -2, 13, 3, 2, -3, 3, 0,
03433 2, -2, 11, 2, 14, 1,
03434 3, 1, 10, 1, 12, -1, 13, 1,
03435 3, -1, 13, 21, 2,-21, 3, 0,
03436 2, 3, 2, -5, 3, 0,
03437 2, 2, 3, -4, 4, 1,
03438 2, 5, 2, -8, 3, 0,
03439 3, -1, 13, 23, 2,-24, 3, 0,
03440 2, 6, 3,-11, 4, 0,
03441 1, 2, 5, 0,
03442 2, 3, 3, -6, 4, 0,
03443 2, 5, 3, -9, 4, 0,
03444 4, 1, 10, -1, 11, 1, 3, -2, 5, 0,
03445 3, 2, 10, 2, 12, -2, 13, 1,
03446 2, 2, 2, -3, 3, 2,
03447 2, 4, 3, -7, 4, 0,
03448 2, 2, 13, -2, 11, 0,
03449 2, 3, 3, -5, 4, 0,
03450 2, 1, 2, -2, 3, 0,
03451 2, 2, 3, -3, 4, 0,
03452 4, 1, 10, -1, 11, 4, 2, -5, 3, 0,
03453 2, 1, 3, -1, 4, 0,
03454 2, 4, 2, -6, 3, 0,
03455 4, 2, 10, -2, 13, 2, 2, -2, 3, 0,
03456 3, 1, 10, -1, 11, 1, 2, 0,
03457 2, 1, 2, -1, 3, 0,
03458 3, 1, 12, 2, 13, -2, 11, 0,
03459 2, 5, 3, -8, 4, 0,
03460 2, 1, 3, -3, 5, 0,
03461 3, 2, 10, 1, 12, -2, 13, 1,
03462 2, 4, 3, -6, 4, 0,
03463 2, 1, 3, -2, 5, 1,
03464 2, 3, 3, -4, 4, 0,
03465 2, 3, 2, -4, 3, 1,
03466 2, 1, 10, -1, 13, 0,
03467 2, 1, 3, -1, 5, 0,
03468 2, 1, 3, -2, 6, 0,
03469 2, 2, 3, -2, 4, 0,
03470 2, 1, 3, -1, 6, 0,
03471 2, 8, 2,-14, 3, 0,
03472 3, 1, 3, 2, 5, -5, 6, 1,
03473 3, 5, 3, -8, 4, 3, 5, 1,
03474 1, 1, 12, 3,
03475 3, 3, 3, -8, 4, 3, 5, 1,
03476 3, 1, 3, -2, 5, 5, 6, 0,
03477 2, 8, 2,-12, 3, 0,
03478 2, 1, 3, 1, 5, 0,
03479 3, 2, 10, 1, 12, -2, 11, 1,
03480 2, 5, 2, -7, 3, 0,
03481 3, 1, 10, 1, 13, -2, 11, 0,
03482 2, 2, 2, -2, 3, 0,
03483 2, 5, 3, -7, 4, 0,
03484 3, 1, 12, -2, 13, 2, 11, 0,
03485 2, 4, 3, -5, 4, 0,
03486 2, 3, 3, -3, 4, 0,
03487 1, 1, 2, 0,
03488 3, 3, 10, 1, 12, -3, 13, 0,
03489 2, 2, 3, -4, 5, 0,
03490 2, 2, 3, -3, 5, 0,
03491 2, 2, 10, -2, 13, 0,
03492 2, 2, 3, -2, 5, 0,
03493 2, 3, 2, -3, 3, 0,
03494 3, 1, 10, -1, 12, -1, 13, 1,
03495 2, 2, 3, -1, 5, 0,
03496 2, 2, 3, -2, 6, 0,
03497 1, 2, 12, 2,
03498 3, -2, 10, 1, 11, 1, 14, 0,
03499 2, 2, 10, -2, 11, 0,
03500 2, 2, 2, -1, 3, 0,
03501 4, -2, 10, 2, 13, 1, 2, -1, 3, 0,
03502 2, 4, 2, -4, 3, 0,
03503 2, 3, 10, -3, 13, 0,
03504 4, -2, 10, 2, 13, 1, 3, -1, 5, 0,
03505 2, 3, 3, -3, 5, 0,
03506 3, 2, 10, -1, 12, -2, 13, 2,
03507 3, 3, 10, -1, 13, -2, 11, 0,
03508 1, 3, 12, 1,
03509 4, -2, 10, 2, 13, 2, 2, -2, 3, 0,
03510 3, 2, 10, -1, 12, -2, 11, 1,
03511 2, 5, 2, -5, 3, 0,
03512 2, 4, 10, -4, 13, 0,
03513 2, 6, 2, -6, 3, 0,
03514 3, 2, 10, -2, 12, -2, 13, 1,
03515 3, 4, 10, -2, 13, -2, 11, 0,
03516 3, 2, 10, -2, 12, -2, 11, 0,
03517 2, 7, 2, -7, 3, 0,
03518 3, 2, 10, -3, 12, -2, 13, 0,
03519 2, 8, 2, -8, 3, 0,
03520 2, 9, 2, -9, 3, 0,
03521 2, 10, 2,-10, 3, 0,
03522 3, 2, 10, -4, 12, -1, 13, 0,
03523 3, 4, 10, -2, 12, -3, 13, 0,
03524 4, 4, 10, -1, 12, -1, 13, -2, 11, 0,
03525 3, 2, 10, -3, 12, -1, 13, 1,
03526 4, -2, 10, 1, 13, 3, 3, -2, 5, 0,
03527 3, 4, 10, -1, 12, -3, 13, 0,
03528 4, -2, 10, 1, 13, 3, 3, -3, 5, 0,
03529 4, 2, 10, -2, 12, 1, 13, -2, 11, 0,
03530 4, -2, 10, 1, 13, 2, 2, -1, 3, 0,
03531 3, 3, 10, -1, 12, -2, 11, 0,
03532 3, 4, 10, -1, 13, -2, 11, 0,
03533 3, 2, 10, -2, 12, -1, 13, 2,
03534 4, -2, 10, 1, 13, 2, 3, -1, 5, 0,
03535 3, 3, 10, -1, 12, -2, 13, 0,
03536 4, -2, 10, 1, 13, 3, 2, -3, 3, 0,
03537 4, -2, 10, 1, 13, 2, 3, -2, 5, 0,
03538 2, 4, 10, -3, 13, 0,
03539 4, -2, 10, 1, 13, 2, 3, -3, 5, 0,
03540 3, -2, 10, 1, 13, 1, 2, 0,
03541 4, 2, 10, -1, 12, 1, 13, -2, 11, 1,
03542 4, -2, 10, 1, 13, 2, 2, -2, 3, 0,
03543 2, 3, 12, -1, 13, 0,
03544 2, 3, 10, -2, 11, 0,
03545 2, 1, 10, -2, 12, 0,
03546 4, 4, 10, 1, 12, -1, 13, -2, 11, 0,
03547 3, -1, 13, 3, 2, -2, 3, 0,
03548 3, -1, 13, 3, 3, -2, 5, 0,
03549 3, -2, 10, 18, 2,-15, 3, 0,
03550 5, 2, 10, -1, 13, 3, 3, -8, 4, 3, 5, 0,
03551 3, 2, 10, -1, 12, -1, 13, 2,
03552 5, -2, 10, 1, 13, 5, 3, -8, 4, 3, 5, 0,
03553 5, -2, 10, 1, 13, 1, 3, 2, 5, -5, 6, 0,
03554 4, 2, 10, -2, 13, 18, 2,-17, 3, 0,
03555 4, -2, 10, 1, 13, 1, 3, -1, 6, 0,
03556 4, -2, 10, 1, 13, 2, 3, -2, 4, 0,
03557 4, -2, 10, 1, 13, 1, 3, -1, 5, 0,
03558 2, 3, 10, -2, 13, 0,
03559 4, -2, 10, 1, 13, 3, 2, -4, 3, 0,
03560 4, -2, 10, 1, 13, 3, 3, -4, 4, 0,
03561 4, -2, 10, 1, 13, 1, 3, -2, 5, 0,
03562 3, 4, 10, 1, 12, -3, 13, 0,
03563 4, -2, 10, 1, 13, 1, 3, -3, 5, 0,
03564 3, -1, 13, 4, 2, -4, 3, 0,
03565 4, -2, 10, 1, 13, 1, 2, -1, 3, 0,
03566 4, -2, 10, 1, 13, 1, 3, -1, 4, 0,
03567 4, -2, 10, 1, 13, 2, 3, -3, 4, 0,
03568 4, -2, 10, 1, 13, 3, 3, -5, 4, 0,
03569 3, 2, 10, 1, 13, -2, 11, 0,
03570 4, -2, 10, -1, 13, 1, 11, 1, 14, 0,
03571 4, -2, 10, 1, 13, 2, 2, -3, 3, 1,
03572 2, 2, 12, -1, 13, 1,
03573 3, 3, 10, 1, 12, -2, 11, 0,
03574 4, 2, 10, -1, 13, 2, 3, -4, 4, 0,
03575 4, 2, 10, -1, 13, 3, 2, -5, 3, 0,
03576 2, 1, 10, -1, 12, 1,
03577 3, -1, 13, 3, 2, -3, 3, 0,
03578 3, -2, 10, 1, 13, 1, 5, 0,
03579 4, 2, 10, -1, 13, 1, 3, -2, 4, 0,
03580 3, -1, 13, 2, 3, -2, 5, 0,
03581 4, 2, 10, -1, 13, -1, 11, 1, 14, 0,
03582 3, -1, 13, 5, 3, -6, 4, 0,
03583 3, -2, 10, 1, 13, 1, 6, 0,
03584 3, -1, 10, 1, 3, -1, 5, 0,
03585 4, -2, 10, 1, 13, 8, 2,-13, 3, 1,
03586 3, -2, 10, 18, 2,-16, 3, 1,
03587 5, -2, 10, 1, 13, 3, 2, -7, 3, 4, 4, 1,
03588 4, 2, 10, -1, 13, 2, 5, -5, 6, 1,
03589 5, 2, 10, -1, 13, 4, 3, -8, 4, 3, 5, 1,
03590 2, 2, 10, -1, 13, 2,
03591 5, -2, 10, 1, 13, 4, 3, -8, 4, 3, 5, 1,
03592 4, -2, 10, 1, 13, 2, 5, -5, 6, 1,
03593 5, 2, 10, -1, 13, 3, 2, -7, 3, 4, 4, 0,
03594 4, 2, 10, -2, 13, 18, 2,-16, 3, 1,
03595 4, 2, 10, -1, 13, 8, 2,-13, 3, 1,
03596 3, -1, 10, 3, 2, -4, 3, 0,
03597 3, -1, 13, 6, 2, -8, 3, 0,
03598 3, -1, 13, 2, 3, -3, 5, 0,
03599 3, -1, 13, 6, 3, -8, 4, 0,
03600 3, 2, 10, -1, 13, 1, 6, 0,
03601 4, -2, 10, 1, 13, -1, 11, 1, 14, 0,
03602 4, -2, 10, 1, 13, 1, 3, -2, 4, 0,
03603 3, 2, 10, -1, 13, 1, 5, 0,
03604 3, 3, 10, 1, 12, -2, 13, 0,
03605 4, -2, 10, 1, 13, 3, 2, -5, 3, 0,
03606 4, -2, 10, 1, 13, 2, 3, -4, 4, 0,
03607 2, -1, 13, 1, 2, 0,
03608 4, 2, 10, -1, 13, 2, 2, -3, 3, 0,
03609 3, -1, 10, 1, 2, -1, 3, 0,
03610 3, -1, 13, 4, 2, -5, 3, 0,
03611 3, 2, 10, -3, 13, 2, 11, 0,
03612 4, 2, 10, -1, 13, 2, 3, -3, 4, 0,
03613 3, -1, 13, 2, 2, -2, 3, 0,
03614 4, 2, 10, -1, 13, 1, 2, -1, 3, 0,
03615 4, 2, 10, 1, 12, 1, 13, -2, 11, 0,
03616 3, -2, 13, 18, 2,-15, 3, 0,
03617 2, 1, 12, -1, 13, 2,
03618 3, -1, 13, 1, 3, -1, 6, 0,
03619 4, 2, 10, -1, 13, 1, 3, -2, 5, 0,
03620 3, -1, 13, 2, 3, -2, 4, 0,
03621 3, -1, 13, 1, 3, -1, 5, 0,
03622 4, 2, 10, -1, 13, 3, 3, -4, 4, 0,
03623 1, 1, 10, 0,
03624 3, -1, 13, 3, 2, -4, 3, 0,
03625 3, -1, 13, 3, 3, -4, 4, 0,
03626 4, 2, 10, -1, 13, 1, 3, -1, 5, 0,
03627 4, 2, 10, -1, 13, 2, 3, -2, 4, 0,
03628 3, -1, 13, 1, 3, -2, 5, 0,
03629 3, 2, 10, 1, 12, -1, 13, 2,
03630 3, 1, 12, 1, 13, -2, 11, 0,
03631 3, -1, 13, 1, 2, -1, 3, 0,
03632 4, 2, 10, -1, 13, 2, 2, -2, 3, 0,
03633 3, -1, 13, 4, 2, -6, 3, 0,
03634 3, -1, 13, 2, 3, -3, 4, 0,
03635 3, 1, 13, 1, 2, -2, 3, 0,
03636 4, 2, 10, -1, 13, 3, 3, -3, 4, 0,
03637 2, 3, 13, -2, 11, 0,
03638 4, 2, 10, -1, 13, 4, 2, -5, 3, 0,
03639 3, 1, 10, 1, 2, -1, 3, 0,
03640 3, -1, 13, 2, 2, -3, 3, 1,
03641 3, 2, 10, 2, 12, -3, 13, 0,
03642 3, 2, 10, -1, 13, 1, 2, 0,
03643 3, 1, 13, 2, 3, -4, 4, 0,
03644 3, 1, 13, 3, 2, -5, 3, 0,
03645 2, 21, 2,-21, 3, 0,
03646 3, 1, 10, 1, 12, -2, 13, 1,
03647 4, 2, 10, -1, 13, 2, 3, -4, 5, 0,
03648 4, 2, 10, -1, 13, 7, 3,-10, 4, 0,
03649 2, -1, 13, 1, 5, 0,
03650 3, 1, 13, 1, 3, -2, 4, 0,
03651 4, 2, 10, -3, 13, 2, 3, -2, 5, 0,
03652 3, 1, 10, 1, 3, -2, 5, 0,
03653 3, 1, 13, -1, 11, 1, 14, 1,
03654 2, -1, 13, 1, 6, 0,
03655 4, 2, 10, -1, 13, 6, 3, -8, 4, 1,
03656 4, 2, 10, -1, 13, 2, 3, -3, 5, 1,
03657 3, -1, 13, 8, 3,-15, 4, 0,
03658 4, 2, 10, -1, 13, 6, 2, -8, 3, 0,
03659 5, 2, 10, -1, 13, -2, 11, 5, 2, -6, 3, 0,
03660 3, 1, 10, 3, 3, -4, 4, 0,
03661 3, 1, 10, 3, 2, -4, 3, 1,
03662 4, 1, 10, -1, 13, -1, 11, 2, 4, 0,
03663 3, -2, 13, 26, 2,-29, 3, 0,
03664 3, -1, 13, 8, 2,-13, 3, 0,
03665 3, -2, 13, 18, 2,-16, 3, 2,
03666 4, -1, 13, 3, 2, -7, 3, 4, 4, 0,
03667 3, 1, 13, 2, 5, -5, 6, 1,
03668 4, 1, 13, 4, 3, -8, 4, 3, 5, 1,
03669 1, 1, 13, 3,
03670 4, -1, 13, 4, 3, -8, 4, 3, 5, 1,
03671 3, -1, 13, 2, 5, -5, 6, 1,
03672 4, 1, 13, 3, 2, -7, 3, 4, 4, 0,
03673 2, 18, 2,-16, 3, 1,
03674 3, 1, 13, 8, 2,-13, 3, 2,
03675 2, 26, 2,-29, 3, 0,
03676 4, 1, 10, 1, 13, -1, 11, 2, 4, 0,
03677 5, 2, 10, 1, 13, -2, 11, 5, 2, -6, 3, 0,
03678 3, 1, 13, 8, 3,-15, 4, 1,
03679 4, 2, 10, -3, 13, 2, 3, -3, 5, 0,
03680 3, 1, 10, 1, 3, -1, 5, 0,
03681 2, 1, 13, 1, 6, 0,
03682 4, 2, 10, -1, 13, 5, 3, -6, 4, 0,
03683 3, 1, 10, 2, 3, -2, 4, 0,
03684 3, -1, 13, -1, 11, 1, 14, 1,
03685 4, 2, 10, -1, 13, 2, 3, -5, 6, 0,
03686 4, 2, 10, -1, 13, 2, 3, -2, 5, 0,
03687 5, 2, 10, -1, 13, 2, 3, -4, 5, 5, 6, 0,
03688 3, -1, 13, 1, 3, -2, 4, 1,
03689 2, 1, 13, 1, 5, 0,
03690 4, 2, 10, -1, 13, 4, 3, -4, 4, 0,
03691 4, 2, 10, -1, 13, 3, 2, -3, 3, 0,
03692 4, 2, 10, 2, 12, -1, 13, -2, 11, 0,
03693 2, 1, 10, 1, 12, 2,
03694 3, -1, 13, 3, 2, -5, 3, 0,
03695 3, -1, 13, 2, 3, -4, 4, 0,
03696 4, 2, 10, -1, 13, 2, 3, -1, 5, 0,
03697 4, 2, 10, -1, 13, 2, 3, -2, 6, 0,
03698 3, 1, 10, 1, 12, -2, 11, 0,
03699 3, 2, 10, 2, 12, -1, 13, 1,
03700 3, 1, 13, 2, 2, -3, 3, 1,
03701 3, -1, 13, 1, 11, 1, 14, 0,
03702 2, 1, 13, -2, 11, 0,
03703 4, 2, 10, -1, 13, 5, 2, -6, 3, 0,
03704 3, -1, 13, 1, 2, -2, 3, 0,
03705 3, 1, 13, 2, 3, -3, 4, 0,
03706 3, 1, 13, 1, 2, -1, 3, 0,
03707 4, 2, 10, -1, 13, 4, 2, -4, 3, 0,
03708 3, 2, 10, 1, 12, -3, 13, 1,
03709 3, 1, 13, 1, 3, -2, 5, 0,
03710 3, 1, 13, 3, 3, -4, 4, 0,
03711 3, 1, 13, 3, 2, -4, 3, 0,
03712 2, 1, 10, -2, 13, 0,
03713 4, 2, 10, -1, 13, 3, 3, -4, 5, 0,
03714 3, 1, 13, 1, 3, -1, 5, 0,
03715 3, 1, 13, 2, 3, -2, 4, 0,
03716 3, 1, 13, 1, 3, -1, 6, 0,
03717 4, 2, 10, -1, 13, 3, 3, -3, 5, 0,
03718 4, 2, 10, -1, 13, 6, 2, -7, 3, 0,
03719 2, 1, 12, 1, 13, 2,
03720 4, 2, 10, -1, 13, 3, 3, -2, 5, 0,
03721 4, 2, 10, 1, 12, -1, 13, -2, 11, 0,
03722 2, 1, 10, 2, 12, 0,
03723 2, 1, 10, -2, 11, 0,
03724 3, 1, 13, 2, 2, -2, 3, 0,
03725 3, 1, 12, -1, 13, 2, 11, 0,
03726 4, 2, 10, -1, 13, 5, 2, -5, 3, 0,
03727 3, 1, 13, 2, 3, -3, 5, 0,
03728 2, 2, 10, -3, 13, 0,
03729 3, 1, 13, 2, 3, -2, 5, 0,
03730 3, 1, 13, 3, 2, -3, 3, 0,
03731 3, 1, 10, -1, 12, -2, 13, 0,
03732 4, 2, 10, -1, 13, 6, 2, -6, 3, 0,
03733 2, 2, 12, 1, 13, 1,
03734 3, 2, 10, -1, 13, -2, 11, 0,
03735 3, 1, 10, -1, 12, -2, 11, 0,
03736 3, 2, 10, 1, 13, -4, 11, 0,
03737 3, 1, 13, 4, 2, -4, 3, 0,
03738 4, 2, 10, -1, 13, 7, 2, -7, 3, 0,
03739 3, 2, 10, -1, 12, -3, 13, 1,
03740 2, 3, 12, 1, 13, 0,
03741 4, 2, 10, -1, 12, -1, 13, -2, 11, 0,
03742 3, 1, 13, 5, 2, -5, 3, 0,
03743 4, 2, 10, -1, 13, 8, 2, -8, 3, 0,
03744 3, 2, 10, -2, 12, -3, 13, 0,
03745 4, 2, 10, -1, 13, 9, 2, -9, 3, 0,
03746 3, 4, 10, -3, 12, -2, 13, 0,
03747 2, 2, 10, -4, 12, 0,
03748 3, 4, 10, -2, 12, -2, 13, 1,
03749 2, 6, 10, -4, 13, 0,
03750 3, 4, 10, -1, 12, -2, 11, 0,
03751 2, 2, 10, -3, 12, 1,
03752 3, 3, 10, -2, 12, -1, 13, 0,
03753 3, -2, 10, 3, 3, -2, 5, 0,
03754 3, 4, 10, -1, 12, -2, 13, 1,
03755 3, -2, 10, 3, 3, -3, 5, 0,
03756 2, 5, 10, -3, 13, 0,
03757 3, -2, 10, 4, 2, -4, 3, 0,
03758 3, -2, 10, 2, 2, -1, 3, 0,
03759 2, 4, 10, -2, 11, 0,
03760 2, 2, 10, -2, 12, 2,
03761 3, -2, 10, 3, 3, -2, 4, 0,
03762 3, -2, 10, 2, 3, -1, 5, 0,
03763 3, 3, 10, -1, 12, -1, 13, 1,
03764 3, -2, 10, 3, 2, -3, 3, 0,
03765 3, -2, 10, 2, 3, -2, 5, 0,
03766 2, 4, 10, -2, 13, 0,
03767 3, -2, 10, 2, 3, -3, 5, 0,
03768 2, -2, 10, 1, 2, 0,
03769 4, 2, 10, -1, 12, 2, 13, -2, 11, 0,
03770 3, -2, 10, 2, 2, -2, 3, 0,
03771 3, 3, 10, 1, 13, -2, 11, 0,
03772 3, 4, 10, 1, 12, -2, 11, 0,
03773 4, 2, 10, -1, 12, -1, 11, 1, 14, 0,
03774 4, -2, 10, -1, 13, 18, 2,-15, 3, 0,
03775 4, 2, 10, 3, 3, -8, 4, 3, 5, 0,
03776 2, 2, 10, -1, 12, 2,
03777 4, -2, 10, 5, 3, -8, 4, 3, 5, 0,
03778 4, 2, 10, -1, 13, 18, 2,-17, 3, 0,
03779 3, -2, 10, 1, 3, -1, 6, 0,
03780 3, -2, 10, 2, 3, -2, 4, 0,
03781 3, -2, 10, 1, 3, -1, 5, 0,
03782 2, 3, 10, -1, 13, 0,
03783 3, -2, 10, 3, 2, -4, 3, 0,
03784 3, -2, 10, 3, 3, -4, 4, 0,
03785 3, -2, 10, 1, 3, -2, 5, 0,
03786 3, 4, 10, 1, 12, -2, 13, 1,
03787 4, 2, 10, -1, 12, -2, 13, 2, 11, 0,
03788 3, -2, 10, 1, 2, -1, 3, 0,
03789 3, -2, 10, 2, 3, -3, 4, 0,
03790 3, 2, 10, 2, 13, -2, 11, 0,
03791 3, -2, 10, 2, 2, -3, 3, 0,
03792 2, 2, 12, -2, 13, 1,
03793 3, 2, 10, 2, 3, -4, 4, 0,
03794 3, 2, 10, 3, 2, -5, 3, 0,
03795 3, 1, 10, -1, 12, 1, 13, 1,
03796 3, -2, 13, 3, 2, -3, 3, 0,
03797 2, -2, 10, 1, 5, 0,
03798 3, 2, 10, 1, 3, -2, 4, 0,
03799 3, -2, 13, 2, 3, -2, 5, 0,
03800 3, 2, 10, -1, 11, 1, 14, 0,
03801 4, 4, 10, -2, 13, 2, 3, -3, 5, 0,
03802 3, -2, 10, 8, 2,-13, 3, 0,
03803 4, -2, 10, -1, 13, 18, 2,-16, 3, 1,
03804 4, -2, 10, 3, 2, -7, 3, 4, 4, 0,
03805 4, 2, 10, 4, 3, -8, 4, 3, 5, 1,
03806 1, 2, 10, 3,
03807 4, -2, 10, 4, 3, -8, 4, 3, 5, 1,
03808 4, 2, 10, 3, 2, -7, 3, 4, 4, 0,
03809 4, 2, 10, -1, 13, 18, 2,-16, 3, 1,
03810 3, 2, 10, 8, 2,-13, 3, 0,
03811 3, -2, 10, -1, 11, 1, 14, 0,
03812 4, 4, 10, -2, 13, 2, 3, -2, 5, 0,
03813 3, -2, 10, 1, 3, -2, 4, 0,
03814 2, 2, 10, 1, 5, 0,
03815 4, 4, 10, -2, 13, 3, 2, -3, 3, 0,
03816 3, 3, 10, 1, 12, -1, 13, 1,
03817 3, -2, 10, 3, 2, -5, 3, 0,
03818 3, -2, 10, 2, 3, -4, 4, 0,
03819 3, 4, 10, 2, 12, -2, 13, 0,
03820 3, 2, 10, 2, 2, -3, 3, 0,
03821 3, 2, 10, -2, 13, 2, 11, 0,
03822 3, 2, 10, 1, 2, -1, 3, 0,
03823 4, 2, 10, 1, 12, 2, 13, -2, 11, 0,
03824 2, 1, 12, -2, 13, 2,
03825 3, 2, 10, 1, 3, -2, 5, 0,
03826 3, -2, 13, 1, 3, -1, 5, 0,
03827 3, 2, 10, 3, 2, -4, 3, 0,
03828 2, 1, 10, 1, 13, 0,
03829 3, 2, 10, 1, 3, -1, 5, 0,
03830 3, 2, 10, 2, 3, -2, 4, 0,
03831 2, 2, 10, 1, 12, 2,
03832 2, 1, 12, -2, 11, 0,
03833 3, -2, 13, 1, 2, -1, 3, 0,
03834 3, 1, 10, -1, 13, 2, 11, 0,
03835 3, 2, 10, 2, 2, -2, 3, 0,
03836 3, 1, 10, 1, 12, -3, 13, 0,
03837 3, 2, 13, -1, 11, 1, 14, 0,
03838 3, 2, 10, 2, 3, -3, 5, 0,
03839 3, 2, 10, 6, 2, -8, 3, 0,
03840 3, -3, 13, 18, 2,-16, 3, 1,
03841 3, 2, 13, 2, 5, -5, 6, 0,
03842 4, 2, 13, 4, 3, -8, 4, 3, 5, 0,
03843 1, 2, 13, 0,
03844 4, -2, 13, 4, 3, -8, 4, 3, 5, 0,
03845 3, -2, 13, 2, 5, -5, 6, 0,
03846 3, 1, 13, 18, 2,-16, 3, 1,
03847 3, -2, 13, -1, 11, 1, 14, 0,
03848 3, 2, 10, 2, 3, -2, 5, 0,
03849 3, 2, 10, 3, 2, -3, 3, 0,
03850 3, 1, 10, 1, 12, 1, 13, 1,
03851 2, 2, 10, 2, 12, 1,
03852 2, 1, 11, 1, 14, 1,
03853 4, -1, 13, -2, 11, 18, 2,-16, 3, 0,
03854 1, 2, 11, 0,
03855 4, -1, 13, 2, 11, 18, 2,-16, 3, 0,
03856 2, -3, 11, 1, 14, 0,
03857 3, 2, 13, 1, 2, -1, 3, 0,
03858 3, 2, 10, 4, 2, -4, 3, 0,
03859 3, 2, 10, 1, 12, -4, 13, 0,
03860 2, 1, 10, -3, 13, 0,
03861 3, 2, 13, 1, 3, -1, 5, 0,
03862 2, 1, 12, 2, 13, 2,
03863 3, 1, 10, 2, 12, 1, 13, 0,
03864 3, 1, 10, -1, 13, -2, 11, 0,
03865 2, 1, 12, 2, 11, 1,
03866 3, 2, 10, 5, 2, -5, 3, 0,
03867 2, 2, 10, -4, 13, 0,
03868 3, 2, 10, 6, 2, -6, 3, 0,
03869 2, 2, 12, 2, 13, 0,
03870 3, 2, 10, -2, 13, -2, 11, 0,
03871 2, 2, 12, 2, 11, 0,
03872 2, 2, 10, -4, 11, 0,
03873 3, 2, 10, 7, 2, -7, 3, 0,
03874 3, 2, 10, -1, 12, -4, 13, 0,
03875 4, 2, 10, -1, 12, -2, 13, -2, 11, 0,
03876 3, 2, 10, 8, 2, -8, 3, 0,
03877 3, 2, 10, 9, 2, -9, 3, 0,
03878 3, 4, 10, -3, 12, -1, 13, 0,
03879 3, 6, 10, -1, 12, -3, 13, 0,
03880 3, 4, 10, -2, 12, -1, 13, 1,
03881 3, 5, 10, -1, 12, -2, 13, 0,
03882 2, 6, 10, -3, 13, 0,
03883 4, 4, 10, -1, 12, 1, 13, -2, 11, 0,
03884 3, 2, 10, -3, 12, 1, 13, 0,
03885 2, 3, 10, -2, 12, 0,
03886 3, 4, 10, -1, 12, -1, 13, 1,
03887 2, 5, 10, -2, 13, 0,
03888 3, 6, 10, 1, 12, -3, 13, 0,
03889 3, 4, 10, 1, 13, -2, 11, 0,
03890 3, 2, 10, -2, 12, 1, 13, 1,
03891 2, 3, 10, -1, 12, 0,
03892 4, -2, 10, -1, 13, 2, 3, -2, 5, 0,
03893 2, 4, 10, -1, 13, 0,
03894 4, 2, 10, -2, 12, -1, 13, 2, 11, 0,
03895 3, 4, 10, -3, 13, 2, 11, 0,
03896 4, -2, 10, -1, 13, 2, 2, -2, 3, 0,
03897 3, 2, 10, -1, 12, 1, 13, 2,
03898 4, -2, 10, -1, 13, 1, 3, -1, 5, 0,
03899 1, 3, 10, 0,
03900 3, 4, 10, 1, 12, -1, 13, 1,
03901 4, 2, 10, -1, 12, -1, 13, 2, 11, 1,
03902 4, -2, 10, -1, 13, 1, 2, -1, 3, 0,
03903 3, 2, 10, 3, 13, -2, 11, 0,
03904 2, 2, 12, -3, 13, 0,
03905 3, 1, 10, -1, 12, 2, 13, 0,
03906 4, 2, 10, 1, 13, -1, 11, 1, 14, 0,
03907 4, -2, 10, -2, 13, 18, 2,-16, 3, 0,
03908 5, 2, 10, 1, 13, 4, 3, -8, 4, 3, 5, 0,
03909 2, 2, 10, 1, 13, 1,
03910 5, -2, 10, -1, 13, 4, 3, -8, 4, 3, 5, 0,
03911 3, 2, 10, 18, 2,-16, 3, 0,
03912 4, -2, 10, -1, 13, -1, 11, 1, 14, 0,
03913 4, 4, 10, -1, 13, 2, 3, -2, 5, 0,
03914 4, 4, 10, -1, 13, 3, 2, -3, 3, 0,
03915 2, 3, 10, 1, 12, 1,
03916 3, 4, 10, 2, 12, -1, 13, 0,
03917 4, 2, 10, -1, 13, 1, 11, 1, 14, 0,
03918 3, 2, 10, -1, 13, 2, 11, 0,
03919 2, 1, 12, -3, 13, 1,
03920 2, 1, 10, 2, 13, 0,
03921 3, 2, 10, 1, 12, 1, 13, 1,
03922 3, 1, 12, -1, 13, -2, 11, 1,
03923 2, 1, 10, 2, 11, 0,
03924 4, 2, 10, 1, 12, -1, 13, 2, 11, 0,
03925 1, 3, 13, 0,
03926 4, 2, 10, 1, 13, 2, 3, -2, 5, 0,
03927 3, 1, 10, 1, 12, 2, 13, 0,
03928 3, 2, 10, 2, 12, 1, 13, 0,
03929 3, 1, 13, 1, 11, 1, 14, 0,
03930 2, 1, 13, 2, 11, 0,
03931 3, 1, 10, 1, 12, 2, 11, 0,
03932 4, 2, 10, 2, 12, -1, 13, 2, 11, 0,
03933 2, 1, 13, -4, 11, 0,
03934 2, 1, 10, -4, 13, 0,
03935 2, 1, 12, 3, 13, 1,
03936 3, 1, 12, 1, 13, 2, 11, 1,
03937 2, 2, 10, -5, 13, 0,
03938 3, 2, 10, -3, 13, -2, 11, 0,
03939 3, 2, 10, -1, 13, -4, 11, 0,
03940 3, 6, 10, -2, 12, -2, 13, 0,
03941 2, 4, 10, -3, 12, 0,
03942 3, 6, 10, -1, 12, -2, 13, 0,
03943 2, 4, 10, -2, 12, 1,
03944 2, 6, 10, -2, 13, 0,
03945 2, 4, 10, -1, 12, 1,
03946 2, 5, 10, -1, 13, 0,
03947 3, 6, 10, 1, 12, -2, 13, 0,
03948 4, 4, 10, -1, 12, -2, 13, 2, 11, 0,
03949 3, 4, 10, 2, 13, -2, 11, 0,
03950 3, 2, 10, -2, 12, 2, 13, 0,
03951 1, 4, 10, 0,
03952 3, 2, 10, -2, 12, 2, 11, 0,
03953 3, 4, 10, -2, 13, 2, 11, 0,
03954 3, 2, 10, -1, 12, 2, 13, 1,
03955 2, 3, 10, 1, 13, 0,
03956 2, 4, 10, 1, 12, 1,
03957 3, 2, 10, -1, 12, 2, 11, 1,
03958 3, 3, 10, -1, 13, 2, 11, 0,
03959 2, 2, 10, 2, 13, 0,
03960 3, 3, 10, 1, 12, 1, 13, 0,
03961 3, 2, 10, 1, 11, 1, 14, 0,
03962 2, 2, 10, 2, 11, 0,
03963 2, 1, 12, -4, 13, 0,
03964 2, 1, 10, 3, 13, 0,
03965 3, 2, 10, 1, 12, 2, 13, 1,
03966 3, 1, 12, -2, 13, -2, 11, 0,
03967 3, 1, 10, 1, 13, 2, 11, 0,
03968 3, 2, 10, 1, 12, 2, 11, 0,
03969 1, 4, 13, 0,
03970 3, 1, 10, 1, 12, 3, 13, 0,
03971 2, 2, 13, 2, 11, 0,
03972 4, 1, 10, 1, 12, 1, 13, 2, 11, 0,
03973 1, 4, 11, 0,
03974 2, 1, 12, 4, 13, 0,
03975 3, 1, 12, 2, 13, 2, 11, 0,
03976 3, 2, 10, -4, 13, -2, 11, 0,
03977 3, 6, 10, -2, 12, -1, 13, 0,
03978 2, 8, 10, -3, 13, 0,
03979 3, 6, 10, -1, 12, -1, 13, 0,
03980 3, 4, 10, -2, 12, 1, 13, 0,
03981 2, 6, 10, -1, 13, 0,
03982 3, 4, 10, -1, 12, 1, 13, 1,
03983 3, 6, 10, 1, 12, -1, 13, 0,
03984 4, 4, 10, -1, 12, -1, 13, 2, 11, 0,
03985 3, 2, 10, -2, 12, 3, 13, 0,
03986 2, 4, 10, 1, 13, 0,
03987 3, 4, 10, -1, 13, 2, 11, 0,
03988 3, 2, 10, -1, 12, 3, 13, 0,
03989 3, 4, 10, 1, 12, 1, 13, 0,
03990 4, 2, 10, -1, 12, 1, 13, 2, 11, 0,
03991 2, 2, 10, 3, 13, 0,
03992 3, 2, 10, 1, 13, 2, 11, 0,
03993 3, 2, 10, -1, 13, 4, 11, 0,
03994 3, 2, 10, 1, 12, 3, 13, 0,
03995 3, 1, 12, -3, 13, -2, 11, 0,
03996 3, 1, 10, 2, 13, 2, 11, 0,
03997 4, 2, 10, 1, 12, 1, 13, 2, 11, 0,
03998 1, 5, 13, 0,
03999 2, 3, 13, 2, 11, 0,
04000 2, 1, 13, 4, 11, 0,
04001 3, 1, 12, 3, 13, 2, 11, 0,
04002 2, 8, 10, -2, 13, 0,
04003 2, 6, 10, -1, 12, 0,
04004 1, 6, 10, 0,
04005 3, 6, 10, -2, 13, 2, 11, 0,
04006 3, 4, 10, -1, 12, 2, 13, 0,
04007 3, 4, 10, -1, 12, 2, 11, 0,
04008 2, 4, 10, 2, 13, 0,
04009 2, 4, 10, 2, 11, 0,
04010 3, 2, 10, -1, 12, 4, 13, 0,
04011 3, 4, 10, 1, 12, 2, 13, 0,
04012 4, 2, 10, -1, 12, 2, 13, 2, 11, 0,
04013 2, 2, 10, 4, 13, 0,
04014 3, 2, 10, 2, 13, 2, 11, 0,
04015 2, 2, 10, 4, 11, 0,
04016 1, 6, 13, 0,
04017 2, 4, 13, 2, 11, 0,
04018 2, 2, 13, 4, 11, 0,
04019 3, 6, 10, -1, 12, 1, 13, 0,
04020 2, 6, 10, 1, 13, 0,
04021 2, 4, 10, 3, 13, 0,
04022 3, 4, 10, 1, 13, 2, 11, 0,
04023 2, 2, 10, 5, 13, 0,
04024 3, 2, 10, 3, 13, 2, 11, 0,
04025 -1
04026 };
04027
04028 static long btabr[] = {-1};
04029 static long btabb[] = {-1};
04030 static long btabl[] = {
04031 -3, -4,
04032 4, -1856, 0, 8043,
04033 -9, -1082,
04034 -1, -310,
04035 -1, -522,
04036 -330, -1449, -853, 4656,
04037 -66, 7,
04038 -1, 9996928,
04039 -66, 6,
04040 23, 183,
04041 0, 173,
04042 0, -56,
04043 0, 50,
04044 0, -785,
04045 1, 51,
04046 0, -60,
04047 1, 11843, 0, -50754,
04048 0, 1834, 1, -7910,
04049 0, -48060,
04050 1, 56,
04051 0, 13141, -1, -56318,
04052 0, 2541,
04053 -1, -649,
04054 -133, 778,
04055 -46, 8,
04056 1, 1665737,
04057 -47, 7,
04058 0, 65,
04059 0, 45,
04060 0, -138,
04061 0, -1005,
04062 0, -2911,
04063 0, -47,
04064 0, 96,
04065 0, -394,
04066 2, 76,
04067 2, -17302, 0, 74337,
04068 0, -101,
04069 0, 58,
04070 0, -171,
04071 0, -77,
04072 0, -1283, 0, 2686,
04073 0, -55,
04074 0, 99,
04075 0, 55,
04076 0, 397,
04077 0, 540,
04078 0, 626,
04079 -1, -5188, 0, 10857,
04080 0, -216,
04081 -2, -123,
04082 0, 6337,
04083 2, 224,
04084 -152, -23472, -29, -74336, 0, 295775,
04085 -20, 149,
04086 -2, 84,
04087 9, 304,
04088 0, -3051,
04089 -70, -6,
04090 -57, 34,
04091 0, -638,
04092 0, -201,
04093 -73, 9,
04094 0, -100,
04095 -101, -8,
04096 0, -57,
04097 0, -207,
04098 -3, 80,
04099 -45, 45,
04100 -5, 102,
04101 -59, -23,
04102 52, 201,
04103 -48, 233, -220, 71,
04104 4, 2810, 0, 6236541,
04105 -61, 218, -216, 67,
04106 51, 201,
04107 -59, -23,
04108 -144, -837, -457, 3029,
04109 -45, 42,
04110 -15, 73,
04111 -6, -169,
04112 0, 135,
04113 -64, -7,
04114 0, -16245,
04115 0, -81,
04116 -74, -10,
04117 0, 702, 0, -3013,
04118 0, -5889,
04119 1, 141,
04120 58, 9598, 12, 30443, 1, -120946,
04121 -1, -84,
04122 -2, 11246, -1, -48391,
04123 0, 1393,
04124 0, 200,
04125 -136, -17,
04126 0, 558,
04127 -64, -8,
04128 0, -71,
04129 0, 317577,
04130 -28, 183,
04131 1, 219,
04132 0, 421,
04133 0, -133,
04134 501, -139,
04135 3, 354,
04136 -101, -13,
04137 74, 7,
04138 144, -84,
04139 59, -2,
04140 1, 64,
04141 -2931, 12559, -4641, 2638, -303, -2058,
04142 -13, -100, -123, -79,
04143 -19214, 6084, 1494, 26993, 15213, -82219,
04144 42, 52, 48, -101,
04145 -53, -4,
04146 4, 47,
04147 58, -131,
04148 46, 14,
04149 -21, -6,
04150 -1311, -8791, 10198, -4185, 2815, 5640,
04151 167, 422, -229, 83,
04152 3140, 39, 1221, 120, 96, -30,
04153 -1, 184612405,
04154 187, 416, -226, 81,
04155 -1985, -10083, 9983, -4464, 2807, 5643,
04156 -21, -9,
04157 113, -367,
04158 120, 580, -667, 27,
04159 8, 66,
04160 -56, -6,
04161 337, 95,
04162 -87, 3303,
04163 -1, 65,
04164 68, -374,
04165 0, -574,
04166 15, -94,
04167 0, -53,
04168 0, -1303,
04169 0, -236,
04170 283, 36,
04171 -1, -54,
04172 269, -35,
04173 0, -83,
04174 0, -52,
04175 0, 730, 0, -3129,
04176 0, 813,
04177 0, -4299,
04178 1, 59,
04179 -6, 5130, 1, 16239, -1, -64603,
04180 0, -80,
04181 91, 12,
04182 0, -561,
04183 133, -17,
04184 0, 250,
04185 -12, 71,
04186 0, 155664,
04187 82, -11,
04188 0, 106,
04189 0, -604,
04190 0, 21862,
04191 55, -7,
04192 0, -1514, 0, 6501,
04193 0, 906,
04194 0, -68,
04195 0, 241,
04196 0, 366,
04197 0, 70,
04198 0, -1382, 0, 5957,
04199 0, 113,
04200 0, -51,
04201 0, -55,
04202 0, 731,
04203 0, -264,
04204 0, 65788,
04205 1, -1504, 0, 3147,
04206 0, 217,
04207 0, -4105, 0, 17658,
04208 1, 69,
04209 0, -3518,
04210 0, -1767,
04211 -43, -7044, -10, -22304, 0, 88685,
04212 3, 91,
04213 0, -485,
04214 0, -57,
04215 -1, 333548,
04216 -24, 172,
04217 11, 544, 1, -1132,
04218 0, 353,
04219 0, -188,
04220 0, 53,
04221 0, 77,
04222 158, -887,
04223 35, 131,
04224 -54, 13,
04225 0, 1994821,
04226 -53, 14,
04227 36, 125,
04228 2, 56,
04229 0, -243,
04230 0, -364,
04231 -2, 1916, 0, -8227,
04232 0, 15700, -1, -67308,
04233 1, 66,
04234 0, -53686,
04235 1, 3058, 1, -13177,
04236 0, -72,
04237 0, -72,
04238 0, 61,
04239 0, 15812,
04240 0, 165,
04241 8, -96,
04242 318, 1341, 803, -4252,
04243 24, 193,
04244 1137, -226, 310, 622,
04245 -56, 30,
04246 -3, 10101666,
04247 -56, 30,
04248 1096, -225, 300, 600,
04249 -31, 409,
04250 -1, -507,
04251 0, -287,
04252 0, -1869, 0, 8026,
04253 1, 544, -1, -1133,
04254 0, 27984,
04255 0, -62,
04256 0, -249,
04257 0, 187,
04258 0, -1096,
04259 1, 53,
04260 2, 12388, 0, -53107,
04261 0, -322,
04262 0, -94,
04263 0, 15157,
04264 0, -582,
04265 0, 3291,
04266 0, 565,
04267 0, 106,
04268 0, 112,
04269 0, 306,
04270 0, 809,
04271 0, 130,
04272 0, -961, 0, 4149,
04273 0, 174,
04274 0, -105,
04275 0, 2196,
04276 0, 59,
04277 0, 36737,
04278 -1, -1832, 0, 3835,
04279 0, -139,
04280 0, 24138,
04281 0, 1325,
04282 1, 64,
04283 0, -361,
04284 0, -1162,
04285 -44, -6320, -10, -20003, 0, 79588,
04286 2, 80,
04287 0, -2059,
04288 0, -304,
04289 0, 21460,
04290 0, -166,
04291 0, -87,
04292 89, -493,
04293 32, 114,
04294 34, 510, 1, 1172616,
04295 31, 113,
04296 -1, 57,
04297 0, 214,
04298 0, -656,
04299 0, -646,
04300 0, 1850, 0, -7931,
04301 0, -6674,
04302 0, 2944, 0, -12641,
04303 0, 916,
04304 45, -255,
04305 16, 60,
04306 -1, 619116,
04307 16, 57,
04308 0, -58,
04309 0, 1045,
04310 0, -156,
04311 -15, 88,
04312 0, -62964,
04313 0, -126,
04314 0, 1490, 0, -6387,
04315 0, 119,
04316 0, 1338,
04317 0, -56,
04318 0, 204,
04319 0, 153,
04320 0, 940,
04321 0, 251,
04322 0, 312,
04323 0, 584,
04324 0, -786, 0, 3388,
04325 0, -52,
04326 0, 4733,
04327 0, 618,
04328 0, 29982,
04329 0, 101,
04330 0, -174,
04331 0, -2637, 0, 11345,
04332 0, -284,
04333 0, -524,
04334 0, -121,
04335 0, 1464,
04336 11, -60,
04337 -1, 151205,
04338 0, 139,
04339 0, -2448,
04340 0, -51,
04341 0, -768,
04342 0, -638,
04343 0, 552, 0, -2370,
04344 0, 70,
04345 0, 64,
04346 0, 57,
04347 0, 39840,
04348 0, 104,
04349 0, -10194,
04350 0, -635,
04351 0, 69,
04352 0, 113,
04353 0, 67,
04354 0, 96,
04355 0, 367,
04356 0, 134,
04357 0, 596,
04358 0, 63,
04359 0, 1622,
04360 0, 483,
04361 0, 72,
04362 0, 11917,
04363 0, -63,
04364 0, 1273,
04365 0, -66,
04366 0, -262,
04367 0, -97,
04368 0, 103,
04369 0, 15196,
04370 0, -1445,
04371 0, -66,
04372 0, -55,
04373 0, -323,
04374 0, 2632,
04375 0, -1179,
04376 0, 59,
04377 0, -56,
04378 0, 78,
04379 0, 65,
04380 0, 422,
04381 0, 309,
04382 0, 2125,
04383 0, -66,
04384 0, 124,
04385 0, -57,
04386 0, 1379,
04387 0, -304,
04388 0, 177,
04389 0, -118,
04390 0, 146,
04391 0, 283,
04392 0, 119,
04393 };
04394 static CHAR bargs[] = {
04395 0, 1,
04396 3, 1, 10, 1, 12, -1, 11, 1,
04397 4, 2, 10, 2, 12, -1, 13, -1, 11, 0,
04398 5, 2, 10, -1, 13, -1, 11, 3, 2, -3, 3, 0,
04399 5, 2, 10, -1, 13, -1, 11, 2, 3, -2, 5, 0,
04400 2, -1, 13, 1, 14, 1,
04401 5, -1, 13, 1, 11, 4, 3, -8, 4, 3, 5, 0,
04402 2, 1, 13, -1, 11, 0,
04403 5, 1, 13, -1, 11, 4, 3, -8, 4, 3, 5, 0,
04404 5, 2, 10, -1, 13, -1, 11, 2, 3, -3, 5, 0,
04405 4, 1, 10, 1, 12, -2, 13, 1, 11, 0,
04406 4, 1, 13, -1, 11, 1, 2, -1, 3, 0,
04407 5, 2, 10, -1, 13, -1, 11, 2, 2, -2, 3, 0,
04408 3, 1, 10, -2, 13, 1, 11, 0,
04409 4, 1, 13, -1, 11, 1, 3, -1, 5, 0,
04410 4, -1, 13, 1, 11, 1, 2, -1, 3, 0,
04411 3, 1, 12, 1, 13, -1, 11, 1,
04412 4, 2, 10, 1, 12, -1, 13, -1, 11, 1,
04413 2, 1, 10, -1, 11, 0,
04414 4, -1, 13, 1, 11, 1, 3, -1, 5, 0,
04415 3, 1, 12, -1, 13, 1, 11, 1,
04416 3, 2, 10, -3, 13, 1, 11, 0,
04417 3, 2, 12, 1, 13, -1, 11, 0,
04418 3, -2, 10, 1, 13, 1, 14, 0,
04419 6, -2, 10, 1, 13, 1, 11, 4, 3, -8, 4, 3, 5, 0,
04420 3, 2, 10, -1, 13, -1, 11, 0,
04421 6, 2, 10, -1, 13, -1, 11, 4, 3, -8, 4, 3, 5, 0,
04422 4, -1, 13, 1, 11, 2, 3, -2, 5, 0,
04423 4, -1, 13, 1, 11, 3, 2, -3, 3, 0,
04424 3, 1, 10, -1, 12, -1, 11, 0,
04425 3, 2, 12, -1, 13, 1, 11, 0,
04426 3, 2, 10, 1, 13, -3, 11, 0,
04427 5, -2, 10, 1, 13, 1, 11, 1, 2, -1, 3, 0,
04428 4, 2, 10, -1, 12, -3, 13, 1, 11, 0,
04429 3, 3, 10, -2, 13, -1, 11, 0,
04430 5, -2, 10, 1, 13, 1, 11, 1, 3, -1, 5, 0,
04431 4, 2, 10, -1, 12, -1, 13, -1, 11, 1,
04432 2, 3, 10, -3, 11, 0,
04433 5, -2, 10, 1, 13, 1, 11, 2, 2, -2, 3, 0,
04434 4, 2, 10, -1, 12, 1, 13, -3, 11, 0,
04435 3, 4, 10, -3, 13, -1, 11, 0,
04436 4, 2, 10, -2, 12, -1, 13, -1, 11, 1,
04437 3, 4, 10, -1, 13, -3, 11, 0,
04438 4, 2, 10, -3, 12, -1, 13, -1, 11, 0,
04439 3, 4, 10, -1, 12, -3, 11, 0,
04440 3, 2, 10, -3, 12, -1, 11, 0,
04441 4, 4, 10, -1, 12, -2, 13, -1, 11, 0,
04442 2, 4, 10, -3, 11, 0,
04443 3, 2, 10, -2, 12, -1, 11, 1,
04444 4, 3, 10, -1, 12, -1, 13, -1, 11, 0,
04445 4, -2, 10, 1, 11, 2, 3, -2, 5, 0,
04446 3, 4, 10, -2, 13, -1, 11, 0,
04447 4, -2, 10, 1, 11, 2, 2, -2, 3, 0,
04448 3, 2, 10, -1, 12, -1, 11, 2,
04449 3, -2, 10, 1, 12, 1, 14, 0,
04450 4, -2, 10, 1, 11, 2, 3, -2, 4, 0,
04451 4, -2, 10, 1, 11, 1, 3, -1, 5, 0,
04452 3, 3, 10, -1, 13, -1, 11, 0,
04453 4, -2, 10, 1, 11, 3, 2, -4, 3, 0,
04454 4, -2, 10, 1, 11, 1, 3, -2, 5, 0,
04455 4, 2, 10, -1, 12, -2, 13, 1, 11, 0,
04456 4, -2, 10, 1, 11, 1, 2, -1, 3, 0,
04457 2, -1, 10, 1, 2, 0,
04458 3, 2, 10, 2, 13, -3, 11, 0,
04459 4, -2, 10, 1, 11, 2, 2, -3, 3, 0,
04460 3, 2, 12, -2, 13, 1, 11, 0,
04461 4, 1, 10, -1, 12, 1, 13, -1, 11, 0,
04462 3, -2, 10, 1, 11, 1, 5, 0,
04463 4, 2, 10, -1, 11, 1, 3, -2, 4, 0,
04464 3, 2, 10, -2, 11, 1, 14, 0,
04465 4, -2, 10, 1, 11, 8, 2,-13, 3, 0,
04466 5, -2, 10, -1, 13, 1, 11, 18, 2,-16, 3, 0,
04467 5, 2, 10, -1, 11, 4, 3, -8, 4, 3, 5, 1,
04468 2, 2, 10, -1, 11, 1,
04469 5, -2, 10, 1, 11, 4, 3, -8, 4, 3, 5, 1,
04470 5, 2, 10, -1, 13, -1, 11, 18, 2,-16, 3, 0,
04471 4, 2, 10, -1, 11, 8, 2,-13, 3, 0,
04472 2, -2, 10, 1, 14, 1,
04473 4, -2, 10, 1, 11, 1, 3, -2, 4, 0,
04474 3, 2, 10, -1, 11, 1, 5, 0,
04475 2, 2, 12, -1, 11, 0,
04476 4, 3, 10, 1, 12, -1, 13, -1, 11, 0,
04477 4, 2, 10, -1, 11, 2, 2, -3, 3, 0,
04478 3, 2, 10, -2, 13, 1, 11, 0,
04479 4, 2, 10, -1, 11, 1, 2, -1, 3, 0,
04480 3, 1, 10, 1, 2, -2, 3, 0,
04481 3, 1, 12, -2, 13, 1, 11, 1,
04482 3, 1, 10, 1, 13, -1, 11, 0,
04483 4, 2, 10, -1, 11, 1, 3, -1, 5, 0,
04484 3, 2, 10, 1, 12, -1, 11, 2,
04485 3, -2, 10, -1, 12, 1, 14, 0,
04486 2, 1, 12, -1, 11, 1,
04487 3, 1, 10, -1, 13, 1, 11, 0,
04488 4, 2, 10, -1, 11, 2, 2, -2, 3, 0,
04489 3, 1, 10, 2, 2, -3, 3, 0,
04490 4, 2, 10, 1, 12, -2, 13, 1, 11, 0,
04491 3, -1, 10, 1, 2, -2, 3, 0,
04492 3, -1, 11, 1, 2, -1, 3, 0,
04493 2, 2, 13, -1, 11, 0,
04494 2, -2, 13, 1, 14, 0,
04495 4, 2, 10, -1, 11, 2, 3, -2, 5, 0,
04496 4, 2, 10, -1, 11, 3, 2, -3, 3, 0,
04497 4, 2, 10, 2, 12, -2, 13, -1, 11, 0,
04498 3, 1, 10, 1, 3, -2, 5, 0,
04499 4, 1, 10, 1, 12, 1, 13, -1, 11, 0,
04500 3, 1, 10, 3, 2, -4, 3, 0,
04501 3, 1, 10, 1, 3, -1, 5, 0,
04502 3, 1, 10, 1, 3, -2, 6, 0,
04503 3, 1, 10, 2, 3, -2, 4, 0,
04504 4, 1, 10, 1, 12, -1, 13, -1, 11, 0,
04505 3, 2, 10, 2, 12, -1, 11, 2,
04506 4, 1, 10, 1, 3, 2, 5, -5, 6, 1,
04507 1, 1, 14, 2,
04508 3, 1, 10, 8, 2,-12, 3, 1,
04509 5, -2, 10, 1, 13, -1, 11, 20, 2,-21, 3, 0,
04510 5, 2, 10, -2, 13, 1, 11, 2, 3, -3, 5, 0,
04511 3, 1, 10, 1, 3, 1, 6, 0,
04512 4, -1, 13, -1, 11, 26, 2,-29, 3, 0,
04513 3, -1, 11, 8, 2,-13, 3, 0,
04514 4, -1, 13, -1, 11, 18, 2,-16, 3, 2,
04515 4, -1, 13, 1, 11, 10, 2, -3, 3, 1,
04516 1, 1, 11, 3,
04517 4, -1, 13, -1, 11, 10, 2, -3, 3, 1,
04518 4, -1, 13, 1, 11, 18, 2,-16, 3, 2,
04519 3, 1, 11, 8, 2,-13, 3, 0,
04520 2, 1, 10, 2, 4, 0,
04521 4, 2, 10, -1, 11, 5, 2, -6, 3, 1,
04522 5, 2, 10, -2, 13, -1, 11, 2, 3, -3, 5, 0,
04523 5, -2, 10, 1, 13, 1, 11, 20, 2,-21, 3, 0,
04524 3, 1, 10, 1, 3, 1, 5, 0,
04525 2, -2, 11, 1, 14, 0,
04526 5, 2, 10, -2, 13, 1, 11, 2, 3, -2, 5, 0,
04527 3, 1, 10, 5, 2, -7, 3, 0,
04528 4, 1, 10, 1, 12, -1, 13, 1, 11, 0,
04529 3, 1, 10, 2, 2, -2, 3, 0,
04530 4, 2, 10, 2, 12, -2, 13, 1, 11, 0,
04531 2, 2, 13, -3, 11, 0,
04532 4, 2, 10, -1, 11, 4, 2, -4, 3, 0,
04533 3, 1, 10, 4, 2, -5, 3, 0,
04534 3, 1, 10, -3, 13, 1, 11, 0,
04535 2, 1, 10, 1, 2, 0,
04536 3, 1, 11, 1, 2, -1, 3, 0,
04537 4, 2, 10, -1, 11, 3, 3, -3, 5, 0,
04538 3, 1, 12, 2, 13, -1, 11, 1,
04539 4, 2, 10, 1, 12, -2, 13, -1, 11, 0,
04540 3, 1, 10, -1, 13, -1, 11, 0,
04541 3, 1, 11, 1, 3, -1, 5, 0,
04542 2, 1, 12, 1, 11, 2,
04543 4, 2, 10, -1, 11, 5, 2, -5, 3, 0,
04544 3, 1, 10, 5, 2, -6, 3, 0,
04545 3, 2, 10, 1, 12, -3, 11, 0,
04546 3, 1, 10, 2, 2, -1, 3, 0,
04547 3, 2, 10, -4, 13, 1, 11, 0,
04548 3, -2, 10, 2, 13, 1, 14, 0,
04549 3, 2, 10, -2, 13, -1, 11, 0,
04550 3, 1, 10, 3, 2, -2, 3, 0,
04551 4, 1, 10, -1, 12, -1, 13, -1, 11, 0,
04552 2, 2, 12, 1, 11, 0,
04553 2, 2, 10, -3, 11, 0,
04554 3, 1, 10, 4, 2, -3, 3, 0,
04555 4, 2, 10, -1, 12, -2, 13, -1, 11, 1,
04556 3, 2, 10, -1, 12, -3, 11, 0,
04557 3, 4, 10, -4, 13, -1, 11, 0,
04558 4, 2, 10, -2, 12, -2, 13, -1, 11, 0,
04559 4, 4, 10, -2, 12, -1, 13, -1, 11, 0,
04560 3, 6, 10, -3, 13, -1, 11, 0,
04561 4, 4, 10, -1, 12, -1, 13, -1, 11, 1,
04562 4, 2, 10, -3, 12, -1, 13, 1, 11, 0,
04563 3, 5, 10, -2, 13, -1, 11, 0,
04564 3, 4, 10, 1, 13, -3, 11, 0,
04565 4, 2, 10, -2, 12, 1, 13, -1, 11, 0,
04566 3, 3, 10, -1, 12, -1, 11, 0,
04567 3, 4, 10, -1, 13, -1, 11, 0,
04568 4, 2, 10, -2, 12, -1, 13, 1, 11, 1,
04569 3, 4, 10, -3, 13, 1, 11, 0,
04570 4, 2, 10, -1, 12, 1, 13, -1, 11, 1,
04571 5, -2, 10, 1, 13, -1, 11, 2, 2, -2, 3, 0,
04572 2, 3, 10, -1, 11, 0,
04573 4, 4, 10, 1, 12, -1, 13, -1, 11, 0,
04574 4, 2, 10, -1, 12, -1, 13, 1, 11, 2,
04575 5, -2, 10, 1, 13, -1, 11, 1, 3, -1, 5, 0,
04576 3, 3, 10, -2, 13, 1, 11, 0,
04577 5, -2, 10, 1, 13, -1, 11, 1, 2, -1, 3, 0,
04578 3, 2, 10, 1, 13, -1, 11, 0,
04579 3, -2, 10, -1, 13, 1, 14, 0,
04580 3, 2, 12, -1, 13, -1, 11, 1,
04581 3, 3, 10, 1, 12, -1, 11, 0,
04582 3, 1, 10, -1, 12, 1, 11, 0,
04583 4, -1, 13, -1, 11, 3, 2, -3, 3, 0,
04584 4, -1, 13, -1, 11, 2, 3, -2, 5, 0,
04585 3, 2, 10, -1, 13, 1, 14, 0,
04586 4, -2, 10, -1, 11, 18, 2,-16, 3, 0,
04587 6, 2, 10, -1, 13, 1, 11, 4, 3, -8, 4, 3, 5, 0,
04588 3, 2, 10, -1, 13, 1, 11, 0,
04589 6, -2, 10, 1, 13, -1, 11, 4, 3, -8, 4, 3, 5, 0,
04590 5, 2, 10, -2, 13, 1, 11, 18, 2,-16, 3, 0,
04591 4, -2, 10, 1, 13, -2, 11, 1, 14, 0,
04592 3, 1, 12, -3, 13, 1, 11, 0,
04593 3, 1, 10, 2, 13, -1, 11, 0,
04594 4, 2, 10, 1, 12, 1, 13, -1, 11, 1,
04595 3, 1, 12, -1, 13, -1, 11, 1,
04596 4, -1, 13, -1, 11, 1, 3, -1, 5, 0,
04597 2, 1, 10, 1, 11, 0,
04598 4, 2, 10, 1, 12, -1, 13, 1, 11, 1,
04599 3, 1, 12, 1, 13, -3, 11, 0,
04600 4, -1, 13, -1, 11, 1, 2, -1, 3, 0,
04601 5, 2, 10, -1, 13, 1, 11, 2, 2, -2, 3, 0,
04602 2, 3, 13, -1, 11, 0,
04603 4, 1, 10, 1, 12, -2, 13, -1, 11, 0,
04604 4, 2, 10, 2, 12, 1, 13, -1, 11, 0,
04605 2, 1, 13, 1, 14, 1,
04606 5, 2, 10, -1, 13, 1, 11, 2, 3, -3, 5, 0,
04607 4, -2, 13, -1, 11, 18, 2,-16, 3, 1,
04608 5, 1, 13, 1, 11, 4, 3, -8, 4, 3, 5, 0,
04609 2, 1, 13, 1, 11, 0,
04610 5, -1, 13, -1, 11, 4, 3, -8, 4, 3, 5, 0,
04611 3, 1, 11, 18, 2,-16, 3, 1,
04612 3, -1, 13, -2, 11, 1, 14, 0,
04613 5, 2, 10, -1, 13, 1, 11, 2, 3, -2, 5, 0,
04614 5, 2, 10, -1, 13, 1, 11, 3, 2, -3, 3, 0,
04615 3, 1, 10, 1, 12, 1, 11, 1,
04616 4, 2, 10, 2, 12, -1, 13, 1, 11, 1,
04617 2, 1, 13, -3, 11, 0,
04618 4, 1, 13, 1, 11, 1, 2, -1, 3, 0,
04619 3, 1, 12, 3, 13, -1, 11, 0,
04620 4, 2, 10, 1, 12, -3, 13, -1, 11, 0,
04621 3, 1, 10, -2, 13, -1, 11, 0,
04622 4, 1, 13, 1, 11, 1, 3, -1, 5, 0,
04623 3, 1, 12, 1, 13, 1, 11, 1,
04624 2, 1, 10, -3, 11, 0,
04625 3, 1, 12, -1, 13, 3, 11, 0,
04626 3, 2, 10, -3, 13, -1, 11, 0,
04627 3, 2, 12, 1, 13, 1, 11, 0,
04628 3, 2, 10, -1, 13, -3, 11, 0,
04629 4, 2, 10, -1, 12, -3, 13, -1, 11, 0,
04630 4, 2, 10, -1, 12, -1, 13, -3, 11, 0,
04631 4, 6, 10, -1, 12, -2, 13, -1, 11, 0,
04632 3, 4, 10, -2, 12, -1, 11, 0,
04633 3, 6, 10, -2, 13, -1, 11, 0,
04634 4, 4, 10, -2, 12, -2, 13, 1, 11, 0,
04635 3, 4, 10, -1, 12, -1, 11, 1,
04636 3, 2, 10, -3, 12, 1, 11, 0,
04637 3, 5, 10, -1, 13, -1, 11, 0,
04638 4, 4, 10, -1, 12, -2, 13, 1, 11, 0,
04639 4, 2, 10, -2, 12, 2, 13, -1, 11, 0,
04640 2, 4, 10, -1, 11, 0,
04641 3, 2, 10, -2, 12, 1, 11, 1,
04642 4, 3, 10, -1, 12, -1, 13, 1, 11, 0,
04643 3, 4, 10, -2, 13, 1, 11, 0,
04644 4, 2, 10, -1, 12, 2, 13, -1, 11, 0,
04645 4, -2, 10, -1, 11, 2, 2, -2, 3, 0,
04646 3, 3, 10, 1, 13, -1, 11, 0,
04647 3, 4, 10, 1, 12, -1, 11, 0,
04648 3, 2, 10, -1, 12, 1, 11, 2,
04649 4, -2, 10, -1, 11, 1, 3, -1, 5, 0,
04650 3, 3, 10, -1, 13, 1, 11, 0,
04651 4, 4, 10, 1, 12, -2, 13, 1, 11, 0,
04652 3, 2, 10, 2, 13, -1, 11, 0,
04653 3, 2, 12, -2, 13, -1, 11, 0,
04654 4, 1, 10, -1, 12, 1, 13, 1, 11, 0,
04655 2, 2, 10, 1, 14, 0,
04656 5, -2, 10, -1, 13, -1, 11, 18, 2,-16, 3, 0,
04657 2, 2, 10, 1, 11, 1,
04658 5, 2, 10, -1, 13, 1, 11, 18, 2,-16, 3, 0,
04659 3, -2, 10, -2, 11, 1, 14, 0,
04660 4, 3, 10, 1, 12, -1, 13, 1, 11, 0,
04661 3, 2, 10, -2, 13, 3, 11, 0,
04662 4, 2, 10, 1, 12, 2, 13, -1, 11, 0,
04663 3, 1, 12, -2, 13, -1, 11, 1,
04664 3, 1, 10, 1, 13, 1, 11, 0,
04665 3, 2, 10, 1, 12, 1, 11, 1,
04666 2, 4, 13, -1, 11, 0,
04667 2, 2, 13, 1, 14, 0,
04668 4, -3, 13, -1, 11, 18, 2,-16, 3, 0,
04669 2, 2, 13, 1, 11, 0,
04670 4, 1, 13, 1, 11, 18, 2,-16, 3, 0,
04671 4, 2, 10, 1, 11, 2, 3, -2, 5, 0,
04672 4, 1, 10, 1, 12, 1, 13, 1, 11, 0,
04673 3, 2, 10, 2, 12, 1, 11, 0,
04674 2, 2, 11, 1, 14, 0,
04675 1, 3, 11, 0,
04676 3, 1, 10, -3, 13, -1, 11, 0,
04677 3, 1, 12, 2, 13, 1, 11, 1,
04678 2, 1, 12, 3, 11, 0,
04679 3, 2, 10, -4, 13, -1, 11, 0,
04680 3, 2, 12, 2, 13, 1, 11, 0,
04681 3, 2, 10, -2, 13, -3, 11, 0,
04682 4, 6, 10, -1, 12, -1, 13, -1, 11, 0,
04683 3, 6, 10, -1, 13, -1, 11, 0,
04684 4, 4, 10, -2, 12, -1, 13, 1, 11, 0,
04685 3, 6, 10, -3, 13, 1, 11, 0,
04686 4, 4, 10, -1, 12, 1, 13, -1, 11, 0,
04687 4, 4, 10, -1, 12, -1, 13, 1, 11, 1,
04688 3, 5, 10, -2, 13, 1, 11, 0,
04689 3, 4, 10, 1, 13, -1, 11, 0,
04690 4, 2, 10, -2, 12, 1, 13, 1, 11, 0,
04691 3, 4, 10, -1, 13, 1, 11, 0,
04692 4, 2, 10, -1, 12, 3, 13, -1, 11, 0,
04693 4, 4, 10, 1, 12, 1, 13, -1, 11, 0,
04694 4, 2, 10, -1, 12, 1, 13, 1, 11, 1,
04695 2, 3, 10, 1, 11, 0,
04696 4, 4, 10, 1, 12, -1, 13, 1, 11, 0,
04697 4, 2, 10, -1, 12, -1, 13, 3, 11, 0,
04698 3, 2, 10, 3, 13, -1, 11, 0,
04699 3, 2, 10, 1, 13, 1, 14, 0,
04700 3, 2, 10, 1, 13, 1, 11, 0,
04701 3, 3, 10, 1, 12, 1, 11, 0,
04702 3, 2, 10, -1, 13, 3, 11, 0,
04703 4, 2, 10, 1, 12, 3, 13, -1, 11, 0,
04704 3, 1, 12, -3, 13, -1, 11, 0,
04705 3, 1, 10, 2, 13, 1, 11, 0,
04706 4, 2, 10, 1, 12, 1, 13, 1, 11, 1,
04707 3, 1, 12, -1, 13, -3, 11, 0,
04708 2, 1, 10, 3, 11, 0,
04709 2, 5, 13, -1, 11, 0,
04710 2, 3, 13, 1, 11, 0,
04711 4, 1, 10, 1, 12, 2, 13, 1, 11, 0,
04712 2, 1, 13, 3, 11, 0,
04713 3, 1, 12, 3, 13, 1, 11, 0,
04714 3, 1, 12, 1, 13, 3, 11, 0,
04715 3, 2, 10, -5, 13, -1, 11, 0,
04716 3, 6, 10, -1, 12, -1, 11, 0,
04717 4, 6, 10, -1, 12, -2, 13, 1, 11, 0,
04718 2, 6, 10, -1, 11, 0,
04719 3, 4, 10, -2, 12, 1, 11, 0,
04720 3, 6, 10, -2, 13, 1, 11, 0,
04721 4, 4, 10, -1, 12, 2, 13, -1, 11, 0,
04722 3, 4, 10, -1, 12, 1, 11, 0,
04723 3, 4, 10, 2, 13, -1, 11, 0,
04724 4, 2, 10, -2, 12, 2, 13, 1, 11, 0,
04725 2, 4, 10, 1, 11, 0,
04726 3, 4, 10, -2, 13, 3, 11, 0,
04727 4, 2, 10, -1, 12, 2, 13, 1, 11, 0,
04728 3, 3, 10, 1, 13, 1, 11, 0,
04729 3, 4, 10, 1, 12, 1, 11, 0,
04730 3, 2, 10, -1, 12, 3, 11, 0,
04731 3, 2, 10, 4, 13, -1, 11, 0,
04732 3, 2, 10, 2, 13, 1, 11, 0,
04733 2, 2, 10, 3, 11, 0,
04734 3, 1, 12, -4, 13, -1, 11, 0,
04735 3, 1, 10, 3, 13, 1, 11, 0,
04736 4, 2, 10, 1, 12, 2, 13, 1, 11, 0,
04737 2, 4, 13, 1, 11, 0,
04738 2, 2, 13, 3, 11, 0,
04739 1, 5, 11, 0,
04740 3, 1, 12, 4, 13, 1, 11, 0,
04741 4, 6, 10, -1, 12, -1, 13, 1, 11, 0,
04742 3, 6, 10, 1, 13, -1, 11, 0,
04743 3, 6, 10, -1, 13, 1, 11, 0,
04744 4, 4, 10, -1, 12, 1, 13, 1, 11, 0,
04745 3, 4, 10, 1, 13, 1, 11, 0,
04746 3, 4, 10, -1, 13, 3, 11, 0,
04747 4, 2, 10, -1, 12, 3, 13, 1, 11, 0,
04748 4, 4, 10, 1, 12, 1, 13, 1, 11, 0,
04749 3, 2, 10, 3, 13, 1, 11, 0,
04750 3, 2, 10, 1, 13, 3, 11, 0,
04751 2, 5, 13, 1, 11, 0,
04752 2, 3, 13, 3, 11, 0,
04753 2, 6, 10, 1, 11, 0,
04754 3, 4, 10, 2, 13, 1, 11, 0,
04755 3, 2, 10, 4, 13, 1, 11, 0,
04756 -1
04757 };
04758 struct plantbl moonlr = {
04759 { 3, 26, 29, 23, 5, 10, 0, 0, 0, 8, 4, 4, 6, 2, 0, 0, 0, 0,},
04760 3,
04761 lrargs,
04762 lrtabl,
04763 lrtabb,
04764 lrtabr,
04765 2.5735686895300000e-03,
04766 3.6525000000000000e+06,
04767 1.0000000000000000e-04,
04768 };
04769
04770 struct plantbl moonlat = {
04771 { 0, 26, 29, 8, 3, 5, 0, 0, 0, 6, 5, 3, 5, 1, 0, 0, 0, 0,},
04772 3,
04773 bargs,
04774 btabl,
04775 btabb,
04776 btabr,
04777 0.0000000000000000e+00,
04778 3.6525000000000000e+06,
04779 1.0000000000000000e-04,
04780 };
04781
04782
04783
04784
04785 static double mods3600(double x)
04786 {
04787 double y;
04788 y = x - 1296000. * floor( x/1296000.);
04789 return(y);
04790 }
04791
04792
04793
04794
04795 static void mean_elements (double JED)
04796 {
04797 double x, T, T2;
04798
04799
04800 T = (JED - MOSHIER_J2000) / 36525.0;
04801 T2 = T*T;
04802
04803
04804
04805
04806
04807 x = mods3600( 538101628.6889819 * T + 908103.213 );
04808 x += (6.39e-6 * T
04809 - 0.0192789) * T2;
04810 Args[0] = x;
04811
04812
04813 x = mods3600( 210664136.4335482 * T + 655127.236 );
04814 x += (-6.27e-6 * T
04815 + 0.0059381) * T2;
04816 Args[1] = x;
04817
04818
04819 x = mods3600( 129597742.283429 * T + 361679.198 );
04820 x += (-5.23e-6 * T
04821 - 2.04411e-2 ) * T2;
04822 Ea_arcsec = x;
04823 Args[2] = x;
04824
04825
04826 x = mods3600( 68905077.493988 * T + 1279558.751 );
04827 x += (-1.043e-5 * T
04828 + 0.0094264) * T2;
04829 Args[3] = x;
04830
04831
04832 x = mods3600( 10925660.377991 * T + 123665.420 );
04833 x += ((((-3.4e-10 * T
04834 + 5.91e-8) * T
04835 + 4.667e-6) * T
04836 + 5.706e-5) * T
04837 - 3.060378e-1)*T2;
04838 Args[4] = x;
04839
04840
04841 x = mods3600( 4399609.855372 * T + 180278.752 );
04842 x += (((( 8.3e-10 * T
04843 - 1.452e-7) * T
04844 - 1.1484e-5) * T
04845 - 1.6618e-4) * T
04846 + 7.561614E-1)*T2;
04847 Args[5] = x;
04848
04849
04850 x = mods3600( 1542481.193933 * T + 1130597.971 )
04851 + (0.00002156*T - 0.0175083)*T2;
04852 Args[6] = x;
04853
04854
04855 x = mods3600( 786550.320744 * T + 1095655.149 )
04856 + (-0.00000895*T + 0.0021103)*T2;
04857 Args[7] = x;
04858
04859
04860
04861 x = mods3600( 1.6029616009939659e+09 * T + 1.0722612202445078e+06 );
04862 x += (((((-3.207663637426e-013 * T
04863 + 2.555243317839e-011) * T
04864 + 2.560078201452e-009) * T
04865 - 3.702060118571e-005) * T
04866 + 6.9492746836058421e-03) * T
04867 - 6.7352202374457519e+00) * T2;
04868 Args[9] = x;
04869
04870
04871 x = mods3600( 1.7395272628437717e+09 * T + 3.3577951412884740e+05 );
04872 x += ((((( 4.474984866301e-013 * T
04873 + 4.189032191814e-011) * T
04874 - 2.790392351314e-009) * T
04875 - 2.165750777942e-006) * T
04876 - 7.5311878482337989e-04) * T
04877 - 1.3117809789650071e+01) * T2;
04878 NF_arcsec = x;
04879 Args[10] = x;
04880
04881
04882 x = mods3600(1.2959658102304320e+08 * T + 1.2871027407441526e+06);
04883 x += ((((((((
04884 1.62e-20 * T
04885 - 1.0390e-17 ) * T
04886 - 3.83508e-15 ) * T
04887 + 4.237343e-13 ) * T
04888 + 8.8555011e-11 ) * T
04889 - 4.77258489e-8 ) * T
04890 - 1.1297037031e-5 ) * T
04891 + 8.7473717367324703e-05) * T
04892 - 5.5281306421783094e-01) * T2;
04893 Args[11] = x;
04894
04895
04896 x = mods3600( 1.7179159228846793e+09 * T + 4.8586817465825332e+05 );
04897 x += (((((-1.755312760154e-012) * T
04898 + 3.452144225877e-011 * T
04899 - 2.506365935364e-008) * T
04900 - 2.536291235258e-004) * T
04901 + 5.2099641302735818e-02) * T
04902 + 3.1501359071894147e+01) * T2;
04903 Args[12] = x;
04904
04905
04906 x = mods3600( 1.7325643720442266e+09 * T + 7.8593980921052420e+05);
04907 x += ((((( 7.200592540556e-014 * T
04908 + 2.235210987108e-010) * T
04909 - 1.024222633731e-008) * T
04910 - 6.073960534117e-005) * T
04911 + 6.9017248528380490e-03) * T
04912 - 5.6550460027471399e+00) * T2;
04913 LP_equinox = x;
04914 Args[13] = x;
04915
04916
04917 x = ((((((((( -8.66e-20*T - 4.759e-17)*T
04918 + 2.424e-15)*T
04919 + 1.3095e-12)*T
04920 + 1.7451e-10)*T
04921 - 1.8055e-8)*T
04922 - 0.0000235316)*T
04923 + 0.000076)*T
04924 + 1.105414)*T
04925 + 5028.791959)*T;
04926
04927
04928
04929
04930 pA_precession = x;
04931
04932
04933
04934
04935
04936 Args[14] = mods3600( 4.48175409e7 * T + 8.060457e5 );
04937
04938
04939 Args[15] = mods3600( 5.36486787e6 * T - 391702.8 );
04940
04941 #if 0
04942
04943 Args[16] = mods3600( 1.7308227257e9 * T - 4.443583e5 );
04944 #endif
04945
04946 Args[17] = mods3600( 1.73573e6 * T );
04947 }
04948
04949
04950
04951
04952
04953 static int sscc (int k, double arg, int n)
04954 {
04955 double cu, su, cv, sv, s;
04956 int i;
04957
04958 s = STR * arg;
04959 su = sin (s);
04960 cu = cos (s);
04961 ss[k][0] = su;
04962 cc[k][0] = cu;
04963 sv = 2.0 * su * cu;
04964 cv = cu * cu - su * su;
04965 ss[k][1] = sv;
04966 cc[k][1] = cv;
04967 for (i = 2; i < n; i++)
04968 {
04969 s = su * cv + cu * sv;
04970 cv = cu * cv - su * sv;
04971 sv = s;
04972 ss[k][i] = sv;
04973 cc[k][i] = cv;
04974 }
04975 return (0);
04976 }
04977
04978
04979
04980
04981 static int g2plan (double J, struct plantbl *plan, double pobj[], int flag)
04982 {
04983 int i, j, k, m, k1, ip, np, nt;
04984
04985
04986 CHAR *p;
04987 long *pl, *pr;
04988 double su, cu, sv, cv;
04989 double t, sl, sr;
04990
04991 mean_elements (J);
04992
04993 if (flag)
04994 Args[13] -= pA_precession;
04995
04996 T = (J - MOSHIER_J2000) / plan->timescale;
04997
04998 for (i = 0; i < NARGS; i++)
04999 {
05000 if ((j = plan->max_harmonic[i]) > 0)
05001 {
05002 sscc (i, Args[i], j);
05003 }
05004 }
05005
05006
05007 p = plan->arg_tbl;
05008
05009 pl = plan->lon_tbl;
05010 pr = plan->rad_tbl;
05011 sl = 0.0;
05012 sr = 0.0;
05013
05014 for (;;)
05015 {
05016
05017
05018 np = *p++;
05019 if (np < 0)
05020 break;
05021 if (np == 0)
05022 {
05023 nt = *p++;
05024
05025 cu = *pl++;
05026 for (ip = 0; ip < nt; ip++)
05027 {
05028 cu = cu * T + *pl++;
05029 }
05030
05031 sl += cu;
05032
05033 cu = *pr++;
05034 for (ip = 0; ip < nt; ip++)
05035 {
05036 cu = cu * T + *pr++;
05037 }
05038 sr += cu;
05039 continue;
05040 }
05041 k1 = 0;
05042 cv = 0.0;
05043 sv = 0.0;
05044 for (ip = 0; ip < np; ip++)
05045 {
05046
05047 j = *p++;
05048
05049 m = *p++ - 1;
05050 if (j)
05051 {
05052 k = abs (j);
05053 k -= 1;
05054 su = ss[m][k];
05055 if (j < 0)
05056 su = -su;
05057 cu = cc[m][k];
05058 if (k1 == 0)
05059 {
05060 sv = su;
05061 cv = cu;
05062 k1 = 1;
05063 }
05064 else
05065 {
05066 t = su * cv + cu * sv;
05067 cv = cu * cv - su * sv;
05068 sv = t;
05069 }
05070 }
05071 }
05072
05073 nt = *p++;
05074
05075 cu = *pl++;
05076 su = *pl++;
05077 for (ip = 0; ip < nt; ip++)
05078 {
05079 cu = cu * T + *pl++;
05080 su = su * T + *pl++;
05081 }
05082 sl += cu * cv + su * sv;
05083
05084 cu = *pr++;
05085 su = *pr++;
05086 for (ip = 0; ip < nt; ip++)
05087 {
05088 cu = cu * T + *pr++;
05089 su = su * T + *pr++;
05090 }
05091 sr += cu * cv + su * sv;
05092 }
05093 t = plan->trunclvl;
05094 pobj[0] = t * sl;
05095 pobj[2] = t * sr;
05096 return (0);
05097 }
05098
05099
05100
05101
05102
05103
05104 static double g1plan (double J, struct plantbl *plan)
05105 {
05106 int i, j, k, m, k1, ip, np, nt;
05107
05108
05109 CHAR *p;
05110 long *pl;
05111 double su, cu, sv, cv;
05112 double t, sl;
05113
05114 T = (J - MOSHIER_J2000) / plan->timescale;
05115 mean_elements (J);
05116
05117 for (i = 0; i < NARGS; i++)
05118 {
05119 if ((j = plan->max_harmonic[i]) > 0)
05120 {
05121 sscc (i, Args[i], j);
05122 }
05123 }
05124
05125
05126 p = plan->arg_tbl;
05127
05128 pl = plan->lon_tbl;
05129 sl = 0.0;
05130
05131 for (;;)
05132 {
05133
05134
05135 np = *p++;
05136 if (np < 0)
05137 break;
05138 if (np == 0)
05139 {
05140 nt = *p++;
05141 cu = *pl++;
05142 for (ip = 0; ip < nt; ip++)
05143 {
05144 cu = cu * T + *pl++;
05145 }
05146
05147 sl += cu;
05148 continue;
05149 }
05150 k1 = 0;
05151 cv = 0.0;
05152 sv = 0.0;
05153 for (ip = 0; ip < np; ip++)
05154 {
05155
05156 j = *p++;
05157
05158 m = *p++ - 1;
05159 if (j)
05160 {
05161 k = abs (j);
05162 k -= 1;
05163 su = ss[m][k];
05164 if (j < 0)
05165 su = -su;
05166 cu = cc[m][k];
05167 if (k1 == 0)
05168 {
05169 sv = su;
05170 cv = cu;
05171 k1 = 1;
05172 }
05173 else
05174 {
05175 t = su * cv + cu * sv;
05176 cv = cu * cv - su * sv;
05177 sv = t;
05178 }
05179 }
05180 }
05181
05182 nt = *p++;
05183
05184 cu = *pl++;
05185 su = *pl++;
05186 for (ip = 0; ip < nt; ip++)
05187 {
05188 cu = cu * T + *pl++;
05189 su = su * T + *pl++;
05190 }
05191 sl += cu * cv + su * sv;
05192 }
05193 return (plan->trunclvl * sl);
05194 }
05195
05196
05197
05198
05199
05200
05201
05202
05203
05204 static int gecmoon (double J, struct plantbl *lrtab, struct plantbl *lattab, double pobj[])
05205 {
05206 double x;
05207
05208 g2plan (J, lrtab, pobj, 0);
05209 x = pobj[0];
05210 x += LP_equinox;
05211 if (x < -6.45e5)
05212 x += 1.296e6;
05213 if (x > 6.45e5)
05214 x -= 1.296e6;
05215 pobj[0] = STR * x;
05216 x = g1plan (J, lattab);
05217 pobj[1] = STR * x;
05218 pobj[2] = (STR * pobj[2] + 1.0) * lrtab->distance;
05219 return 0;
05220 }
05221
05222
05223
05224 static void moon_fast (double mjd, double *lam, double *bet,
05225 double *hp, double *msp, double *mdp);
05226
05227
05228
05229
05230
05231
05232
05233
05234
05235
05236
05237
05238 static void moon_fast (double mjd, double *lam, double *bet,
05239 double *hp, double *msp, double *mdp)
05240 {
05241 double t, t2;
05242 double ld;
05243 double ms;
05244 double md;
05245 double de;
05246 double f;
05247 double n;
05248 double a, sa, sn, b, sb, c, sc, e, e2, l, g, w1, w2;
05249 double m1, m2, m3, m4, m5, m6;
05250
05251 t = mjd/36525.;
05252 t2 = t*t;
05253
05254 m1 = mjd/27.32158213;
05255 m1 = 360.0*(m1-(long)m1);
05256 m2 = mjd/365.2596407;
05257 m2 = 360.0*(m2-(long)m2);
05258 m3 = mjd/27.55455094;
05259 m3 = 360.0*(m3-(long)m3);
05260 m4 = mjd/29.53058868;
05261 m4 = 360.0*(m4-(long)m4);
05262 m5 = mjd/27.21222039;
05263 m5 = 360.0*(m5-(long)m5);
05264 m6 = mjd/6798.363307;
05265 m6 = 360.0*(m6-(long)m6);
05266
05267 ld = 270.434164+m1-(.001133-.0000019*t)*t2;
05268 ms = 358.475833+m2-(.00015+.0000033*t)*t2;
05269 md = 296.104608+m3+(.009192+.0000144*t)*t2;
05270 de = 350.737486+m4-(.001436-.0000019*t)*t2;
05271 f = 11.250889+m5-(.003211+.0000003*t)*t2;
05272 n = 259.183275-m6+(.002078+.000022*t)*t2;
05273
05274 a = degrad(51.2+20.2*t);
05275 sa = sin(a);
05276 sn = sin(degrad(n));
05277 b = 346.56+(132.87-.0091731*t)*t;
05278 sb = .003964*sin(degrad(b));
05279 c = degrad(n+275.05-2.3*t);
05280 sc = sin(c);
05281 ld = ld+.000233*sa+sb+.001964*sn;
05282 ms = ms-.001778*sa;
05283 md = md+.000817*sa+sb+.002541*sn;
05284 f = f+sb-.024691*sn-.004328*sc;
05285 de = de+.002011*sa+sb+.001964*sn;
05286 e = 1-(.002495+7.52e-06*t)*t;
05287 e2 = e*e;
05288
05289 ld = degrad(ld);
05290 ms = degrad(ms);
05291 n = degrad(n);
05292 de = degrad(de);
05293 f = degrad(f);
05294 md = degrad(md);
05295
05296 l = 6.28875*sin(md)+1.27402*sin(2*de-md)+.658309*sin(2*de)+
05297 .213616*sin(2*md)-e*.185596*sin(ms)-.114336*sin(2*f)+
05298 .058793*sin(2*(de-md))+.057212*e*sin(2*de-ms-md)+
05299 .05332*sin(2*de+md)+.045874*e*sin(2*de-ms)+.041024*e*sin(md-ms);
05300 l = l-.034718*sin(de)-e*.030465*sin(ms+md)+.015326*sin(2*(de-f))-
05301 .012528*sin(2*f+md)-.01098*sin(2*f-md)+.010674*sin(4*de-md)+
05302 .010034*sin(3*md)+.008548*sin(4*de-2*md)-e*.00791*sin(ms-md+2*de)-
05303 e*.006783*sin(2*de+ms);
05304 l = l+.005162*sin(md-de)+e*.005*sin(ms+de)+.003862*sin(4*de)+
05305 e*.004049*sin(md-ms+2*de)+.003996*sin(2*(md+de))+
05306 .003665*sin(2*de-3*md)+e*.002695*sin(2*md-ms)+
05307 .002602*sin(md-2*(f+de))+e*.002396*sin(2*(de-md)-ms)-
05308 .002349*sin(md+de);
05309 l = l+e2*.002249*sin(2*(de-ms))-e*.002125*sin(2*md+ms)-
05310 e2*.002079*sin(2*ms)+e2*.002059*sin(2*(de-ms)-md)-
05311 .001773*sin(md+2*(de-f))-.001595*sin(2*(f+de))+
05312 e*.00122*sin(4*de-ms-md)-.00111*sin(2*(md+f))+.000892*sin(md-3*de);
05313 l = l-e*.000811*sin(ms+md+2*de)+e*.000761*sin(4*de-ms-2*md)+
05314 e2*.000704*sin(md-2*(ms+de))+e*.000693*sin(ms-2*(md-de))+
05315 e*.000598*sin(2*(de-f)-ms)+.00055*sin(md+4*de)+.000538*sin(4*md)+
05316 e*.000521*sin(4*de-ms)+.000486*sin(2*md-de);
05317 l = l+e2*.000717*sin(md-2*ms);
05318 *lam = ld+degrad(l);
05319 range (lam, 2*PI);
05320
05321 g = 5.12819*sin(f)+.280606*sin(md+f)+.277693*sin(md-f)+
05322 .173238*sin(2*de-f)+.055413*sin(2*de+f-md)+.046272*sin(2*de-f-md)+
05323 .032573*sin(2*de+f)+.017198*sin(2*md+f)+.009267*sin(2*de+md-f)+
05324 .008823*sin(2*md-f)+e*.008247*sin(2*de-ms-f);
05325 g = g+.004323*sin(2*(de-md)-f)+.0042*sin(2*de+f+md)+
05326 e*.003372*sin(f-ms-2*de)+e*.002472*sin(2*de+f-ms-md)+
05327 e*.002222*sin(2*de+f-ms)+e*.002072*sin(2*de-f-ms-md)+
05328 e*.001877*sin(f-ms+md)+.001828*sin(4*de-f-md)-e*.001803*sin(f+ms)-
05329 .00175*sin(3*f);
05330 g = g+e*.00157*sin(md-ms-f)-.001487*sin(f+de)-e*.001481*sin(f+ms+md)+
05331 e*.001417*sin(f-ms-md)+e*.00135*sin(f-ms)+.00133*sin(f-de)+
05332 .001106*sin(f+3*md)+.00102*sin(4*de-f)+.000833*sin(f+4*de-md)+
05333 .000781*sin(md-3*f)+.00067*sin(f+4*de-2*md);
05334 g = g+.000606*sin(2*de-3*f)+.000597*sin(2*(de+md)-f)+
05335 e*.000492*sin(2*de+md-ms-f)+.00045*sin(2*(md-de)-f)+
05336 .000439*sin(3*md-f)+.000423*sin(f+2*(de+md))+
05337 .000422*sin(2*de-f-3*md)-e*.000367*sin(ms+f+2*de-md)-
05338 e*.000353*sin(ms+f+2*de)+.000331*sin(f+4*de);
05339 g = g+e*.000317*sin(2*de+f-ms+md)+e2*.000306*sin(2*(de-ms)-f)-
05340 .000283*sin(md+3*f);
05341 w1 = .0004664*cos(n);
05342 w2 = .0000754*cos(c);
05343 *bet = degrad(g)*(1-w1-w2);
05344
05345 *hp = .950724+.051818*cos(md)+.009531*cos(2*de-md)+.007843*cos(2*de)+
05346 .002824*cos(2*md)+.000857*cos(2*de+md)+e*.000533*cos(2*de-ms)+
05347 e*.000401*cos(2*de-md-ms)+e*.00032*cos(md-ms)-.000271*cos(de)-
05348 e*.000264*cos(ms+md)-.000198*cos(2*f-md);
05349 *hp = *hp+.000173*cos(3*md)+.000167*cos(4*de-md)-e*.000111*cos(ms)+
05350 .000103*cos(4*de-2*md)-.000084*cos(2*md-2*de)-
05351 e*.000083*cos(2*de+ms)+.000079*cos(2*de+2*md)+.000072*cos(4*de)+
05352 e*.000064*cos(2*de-ms+md)-e*.000063*cos(2*de+ms-md)+
05353 e*.000041*cos(ms+de);
05354 *hp = *hp+e*.000035*cos(2*md-ms)-.000033*cos(3*md-2*de)-
05355 .00003*cos(md+de)-.000029*cos(2*(f-de))-e*.000029*cos(2*md+ms)+
05356 e2*.000026*cos(2*(de-ms))-.000023*cos(2*(f-de)+md)+
05357 e*.000019*cos(4*de-ms-md);
05358 *hp = degrad(*hp);
05359
05360 *msp = ms;
05361 *mdp = md;
05362 }
05363
05364
05365 #define EarthRadius 6378.16
05366
05367
05368
05369
05370
05371
05372
05373
05374
05375
05376
05377
05378
05379
05380
05381
05382
05383 void moon (double mjd, double *lam, double *bet, double *rho, double *msp, double *mdp)
05384 {
05385 double pobj[3], dt;
05386 double hp;
05387
05388 if (mjd >= MOSHIER_BEGIN && mjd <= MOSHIER_END) {
05389
05390 moon_fast (mjd, lam, bet, &hp, msp, mdp);
05391 *rho = EarthRadius/AUKM/sin(hp);
05392 dt = *rho * 5.7755183e-3;
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];
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 }
05407
05408 #define NUT_SCALE 1e4
05409 #define NUT_SERIES 106
05410 #define NUT_MAXMUL 4
05411 #define SECPERCIRC (3600.*360.)
05412
05413
05414 static double delaunay[5][4] = {
05415 {485866.733, 1717915922.633, 31.310, 0.064},
05416 {1287099.804, 129596581.224, -0.577, -0.012},
05417 {335778.877, 1739527263.137, -13.257, 0.011},
05418 {1072261.307, 1602961601.328, -6.891, 0.019},
05419 {450160.280, -6962890.539, 7.455, 0.008},
05420 };
05421
05422
05423 static short multarg[NUT_SERIES][5] = {
05424
05425 {0, 0, 0, 0, 1},
05426 {0, 0, 0, 0, 2},
05427 {-2, 0, 2, 0, 1},
05428 {2, 0, -2, 0, 0},
05429 {-2, 0, 2, 0, 2},
05430 {1, -1, 0, -1, 0},
05431 {0, -2, 2, -2, 1},
05432 {2, 0, -2, 0, 1},
05433 {0, 0, 2, -2, 2},
05434 {0, 1, 0, 0, 0},
05435 {0, 1, 2, -2, 2},
05436 {0, -1, 2, -2, 2},
05437 {0, 0, 2, -2, 1},
05438 {2, 0, 0, -2, 0},
05439 {0, 0, 2, -2, 0},
05440 {0, 2, 0, 0, 0},
05441 {0, 1, 0, 0, 1},
05442 {0, 2, 2, -2, 2},
05443 {0, -1, 0, 0, 1},
05444 {-2, 0, 0, 2, 1},
05445 {0, -1, 2, -2, 1},
05446 {2, 0, 0, -2, 1},
05447 {0, 1, 2, -2, 1},
05448 {1, 0, 0, -1, 0},
05449 {2, 1, 0, -2, 0},
05450 {0, 0, -2, 2, 1},
05451 {0, 1, -2, 2, 0},
05452 {0, 1, 0, 0, 2},
05453 {-1, 0, 0, 1, 1},
05454 {0, 1, 2, -2, 0},
05455 {0, 0, 2, 0, 2},
05456 {1, 0, 0, 0, 0},
05457 {0, 0, 2, 0, 1},
05458 {1, 0, 2, 0, 2},
05459 {1, 0, 0, -2, 0},
05460 {-1, 0, 2, 0, 2},
05461 {0, 0, 0, 2, 0},
05462 {1, 0, 0, 0, 1},
05463 {-1, 0, 0, 0, 1},
05464 {-1, 0, 2, 2, 2},
05465 {1, 0, 2, 0, 1},
05466 {0, 0, 2, 2, 2},
05467 {2, 0, 0, 0, 0},
05468 {1, 0, 2, -2, 2},
05469 {2, 0, 2, 0, 2},
05470 {0, 0, 2, 0, 0},
05471 {-1, 0, 2, 0, 1},
05472 {-1, 0, 0, 2, 1},
05473 {1, 0, 0, -2, 1},
05474 {-1, 0, 2, 2, 1},
05475 {1, 1, 0, -2, 0},
05476 {0, 1, 2, 0, 2},
05477 {0, -1, 2, 0, 2},
05478 {1, 0, 2, 2, 2},
05479 {1, 0, 0, 2, 0},
05480 {2, 0, 2, -2, 2},
05481 {0, 0, 0, 2, 1},
05482 {0, 0, 2, 2, 1},
05483 {1, 0, 2, -2, 1},
05484 {0, 0, 0, -2, 1},
05485 {1, -1, 0, 0, 0},
05486 {2, 0, 2, 0, 1},
05487 {0, 1, 0, -2, 0},
05488 {1, 0, -2, 0, 0},
05489 {0, 0, 0, 1, 0},
05490 {1, 1, 0, 0, 0},
05491 {1, 0, 2, 0, 0},
05492 {1, -1, 2, 0, 2},
05493 {-1, -1, 2, 2, 2},
05494 {-2, 0, 0, 0, 1},
05495 {3, 0, 2, 0, 2},
05496 {0, -1, 2, 2, 2},
05497 {1, 1, 2, 0, 2},
05498 {-1, 0, 2, -2, 1},
05499 {2, 0, 0, 0, 1},
05500 {1, 0, 0, 0, 2},
05501 {3, 0, 0, 0, 0},
05502 {0, 0, 2, 1, 2},
05503 {-1, 0, 0, 0, 2},
05504 {1, 0, 0, -4, 0},
05505 {-2, 0, 2, 2, 2},
05506 {-1, 0, 2, 4, 2},
05507 {2, 0, 0, -4, 0},
05508 {1, 1, 2, -2, 2},
05509 {1, 0, 2, 2, 1},
05510 {-2, 0, 2, 4, 2},
05511 {-1, 0, 4, 0, 2},
05512 {1, -1, 0, -2, 0},
05513 {2, 0, 2, -2, 1},
05514 {2, 0, 2, 2, 2},
05515 {1, 0, 0, 2, 1},
05516 {0, 0, 4, -2, 2},
05517 {3, 0, 2, -2, 2},
05518 {1, 0, 2, -2, 0},
05519 {0, 1, 2, 0, 1},
05520 {-1, -1, 0, 2, 1},
05521 {0, 0, -2, 0, 1},
05522 {0, 0, 2, -1, 2},
05523 {0, 1, 0, 2, 0},
05524 {1, 0, -2, -2, 0},
05525 {0, -1, 2, 0, 1},
05526 {1, 1, 0, -2, 1},
05527 {1, 0, -2, 2, 0},
05528 {2, 0, 0, 2, 0},
05529 {0, 0, 2, 4, 2},
05530 {0, 1, 0, 1, 0}
05531 };
05532
05533
05534
05535
05536 static long ampsecul[][5] = {
05537 {0 ,-171996 ,-1742 ,92025 ,89},
05538 {1 ,2062 ,2 ,-895 ,5},
05539 {8 ,-13187 ,-16 ,5736 ,-31},
05540 {9 ,1426 ,-34 ,54 ,-1},
05541 {10 ,-517 ,12 ,224 ,-6},
05542 {11 ,217 ,-5 ,-95 ,3},
05543 {12 ,129 ,1 ,-70 ,0},
05544 {15 ,17 ,-1 ,0 ,0},
05545 {17 ,-16 ,1 ,7 ,0},
05546 {30 ,-2274 ,-2 ,977 ,-5},
05547 {31 ,712 ,1 ,-7 ,0},
05548 {32 ,-386 ,-4 ,200 ,0},
05549 {33 ,-301 ,0 ,129 ,-1},
05550 {37 ,63 ,1 ,-33 ,0},
05551 {38 ,-58 ,-1 ,32 ,0},
05552 { -1, }
05553 };
05554
05555
05556
05557
05558
05559 static short ampconst[NUT_SERIES][2] = {
05560 {0,0},
05561 {0,0},
05562 {46,-24},
05563 {11,0},
05564 {-3,1},
05565 {-3,0},
05566 {-2,1},
05567 {1,0},
05568 {0,0},
05569 {0,0},
05570 {0,0},
05571 {0,0},
05572 {0,0},
05573 {48,1},
05574 {-22,0},
05575 {0,0},
05576 {-15,9},
05577 {0,0},
05578 {-12,6},
05579 {-6,3},
05580 {-5,3},
05581 {4,-2},
05582 {4,-2},
05583 {-4,0},
05584 {1,0},
05585 {1,0},
05586 {-1,0},
05587 {1,0},
05588 {1,0},
05589 {-1,0},
05590 {0,0},
05591 {0,0},
05592 {0,0},
05593 {0,0},
05594 {-158,-1},
05595 {123,-53},
05596 {63,-2},
05597 {0,0},
05598 {0,0},
05599 {-59,26},
05600 {-51,27},
05601 {-38,16},
05602 {29,-1},
05603 {29,-12},
05604 {-31,13},
05605 {26,-1},
05606 {21,-10},
05607 {16,-8},
05608 {-13,7},
05609 {-10,5},
05610 {-7,0},
05611 {7,-3},
05612 {-7,3},
05613 {-8,3},
05614 {6,0},
05615 {6,-3},
05616 {-6,3},
05617 {-7,3},
05618 {6,-3},
05619 {-5,3},
05620 {5,0},
05621 {-5,3},
05622 {-4,0},
05623 {4,0},
05624 {-4,0},
05625 {-3,0},
05626 {3,0},
05627 {-3,1},
05628 {-3,1},
05629 {-2,1},
05630 {-3,1},
05631 {-3,1},
05632 {2,-1},
05633 {-2,1},
05634 {2,-1},
05635 {-2,1},
05636 {2,0},
05637 {2,-1},
05638 {1,-1},
05639 {-1,0},
05640 {1,-1},
05641 {-2,1},
05642 {-1,0},
05643 {1,-1},
05644 {-1,1},
05645 {-1,1},
05646 {1,0},
05647 {1,0},
05648 {1,-1},
05649 {-1,0},
05650 {-1,0},
05651 {1,0},
05652 {1,0},
05653 {-1,0},
05654 {1,0},
05655 {1,0},
05656 {-1,0},
05657 {-1,0},
05658 {-1,0},
05659 {-1,0},
05660 {-1,0},
05661 {-1,0},
05662 {-1,0},
05663 {1,0},
05664 {-1,0},
05665 {1,0}
05666 };
05667
05668
05669
05670
05671 void nutation (double mjd, double *deps, double *dpsi)
05672
05673
05674
05675 {
05676 static double lastmjd = -10000, lastdeps, lastdpsi;
05677 double T, T2, T3, T10;
05678 double prec;
05679 int i, isecul;
05680 static double delcache[5][2*NUT_MAXMUL+1];
05681
05682
05683
05684
05685
05686 if (mjd == lastmjd) {
05687 *deps = lastdeps;
05688 *dpsi = lastdpsi;
05689 return;
05690 }
05691
05692 prec = 0.0;
05693
05694 #if 0
05695 prec =* deps;
05696 if (prec < 0.0 || prec > 1.0)
05697 prec = 1.0;
05698 #endif
05699
05700
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
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
05719 x /= SECPERCIRC;
05720 x -= floor(x);
05721 x *= 2.*PI;
05722
05723
05724 for (j = 0; j <= 2*NUT_MAXMUL; ++j)
05725 delcache[i][j] = (j - NUT_MAXMUL) * x;
05726 }
05727
05728
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
05736 ampsin = ampconst[i][0];
05737 ampcos = ampconst[i][1];
05738 } else {
05739
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
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 }
05765
05766
05767
05768
05769 void nut_eq (double mjd, double *ra, double *dec)
05770 {
05771 static double lastmjd = -10000;
05772 static double a[3][3];
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
05783
05784
05785
05786
05787
05788
05789
05790
05791
05792
05793
05794
05795
05796
05797
05798
05799
05800
05801
05802
05803
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);
05833 if (*ra < 0.) *ra += 2.*PI;
05834 }
05835
05836 void ta_par (double tha, double tdec, double phi, double ht, double *rho, double *aha, double *adec)
05837 {
05838 static double last_phi = 1000.0, last_ht = -1000.0, xobs, zobs;
05839 double x, y, z;
05840
05841
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
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 }
05860
05861 static void elongation (double lam, double bet, double lsn, double *el);
05862 static void unrefractLT15 (double pr, double tr, double aa, double *ta);
05863 static void unrefractGE15 (double pr, double tr, double aa, double *ta);
05864
05865 void unrefract (double pr, double tr, double aa, double *ta)
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
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 }
05887
05888 static void unrefractGE15 (double pr, double tr, double aa, double *ta)
05889 {
05890 double r;
05891
05892 r = 7.888888e-5*pr/((273+tr)*tan(aa));
05893 *ta = aa - r;
05894 }
05895
05896 static void unrefractLT15 (double pr, double tr, double aa, double *ta)
05897 {
05898 double aadeg = raddeg(aa);
05899 double r, a, b;
05900
05901 a = ((2e-5*aadeg+1.96e-2)*aadeg+1.594e-1)*pr;
05902 b = (273+tr)*((8.45e-2*aadeg+5.05e-1)*aadeg+1);
05903 r = degrad(a/b);
05904
05905 *ta = aa - r;
05906 }
05907
05908
05909
05910
05911
05912
05913
05914
05915
05916 void refract (double pr,double tr, double ta,double *aa)
05917 {
05918 #define MAXRERR degrad(0.1/3600.)
05919
05920 double d, t, t0, a;
05921
05922
05923
05924
05925 unrefract (pr, tr, ta, &t);
05926 d = 0.8*(ta - t);
05927 t0 = t;
05928 a = ta;
05929
05930
05931
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 }
05944
05945 #define ABERR_CONST (20.49552/3600./180.*PI)
05946 #define AB_ECL_EOD 0
05947 #define AB_EQ_EOD 1
05948
05949 static void ab_aux (double mjd, double *x, double *y, double lsn,int mode);
05950
05951
05952
05953
05954
05955
05956 void ab_ecl (double mjd, double lsn, double *lam, double *bet)
05957 {
05958 ab_aux(mjd, lam, bet, lsn, AB_ECL_EOD);
05959 }
05960
05961
05962
05963
05964
05965 void ab_eq (double mjd, double lsn, double *ra, double *dec)
05966 {
05967 ab_aux(mjd, ra, dec, lsn, AB_EQ_EOD);
05968 }
05969
05970
05971
05972
05973
05974
05975 static void ab_aux (double mjd, double *x, double *y, double lsn, int mode)
05976 {
05977 static double lastmjd = -10000;
05978 static double eexc;
05979 static double leperi;
05980 static char dirty = 1;
05981
05982 if (mjd != lastmjd) {
05983 double T;
05984
05985 T = (mjd - J2000)/36525.;
05986 eexc = 0.016708617 - (42.037e-6 + 0.1236e-6 * T) * T;
05987 leperi = degrad(102.93735 + (0.71953 + 0.00046 * T) * T);
05988 lastmjd = mjd;
05989 dirty = 1;
05990 }
05991
05992 switch (mode) {
05993 case AB_ECL_EOD:
05994 {
05995 double *lam = x, *bet = y;
05996 double dlsun, dlperi;
05997
05998 dlsun = lsn - *lam;
05999 dlperi = leperi - *lam;
06000
06001
06002 *lam -= ABERR_CONST/cos(*bet) * (cos(dlsun) -
06003 eexc*cos(dlperi));
06004 *bet -= ABERR_CONST*sin(*bet) * (sin(dlsun) -
06005 eexc*sin(dlperi));
06006 }
06007 break;
06008
06009 case AB_EQ_EOD:
06010 {
06011 double *ra = x, *dec = y;
06012 double sr, cr, sd, cd, sls, cls;
06013 static double cp, sp, ce, se;
06014 double dra, ddec;
06015
06016 if (dirty) {
06017 double eps;
06018
06019 cp = cos(leperi);
06020 sp = sin(leperi);
06021 obliquity(mjd, &eps);
06022 se = sin(eps);
06023 ce = cos(eps);
06024 dirty = 0;
06025 }
06026
06027 sr = sin(*ra);
06028 cr = cos(*ra);
06029 sd = sin(*dec);
06030 cd = cos(*dec);
06031 sls = sin(lsn);
06032 cls = cos(lsn);
06033
06034 dra = ABERR_CONST/cd * ( -(cr * cls * ce + sr * sls) +
06035 eexc * (cr * cp * ce + sr * sp));
06036
06037 ddec = se/ce * cd - sr * sd;
06038 ddec = ABERR_CONST * ( -(cls * ce * ddec + cr * sd * sls) +
06039 eexc * (cp * ce * ddec + cr * sd * sp) );
06040
06041 *ra += dra;
06042 range (ra, 2*PI);
06043 *dec += ddec;
06044 }
06045 break;
06046
06047 default:
06048
06049
06050 break;
06051
06052 }
06053 }
06054
06055
06056
06057
06058
06059
06060 static void ecleq_aux (int sw, double mjd, double x, double y,
06061 double *p, double *q);
06062
06063 #define EQtoECL 1
06064 #define ECLtoEQ (-1)
06065
06066
06067
06068
06069
06070
06071
06072
06073 void eq_ecl (double mjd, double ra, double dec, double *lat, double *lng)
06074 {
06075 ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
06076 }
06077
06078
06079
06080
06081
06082
06083
06084 void ecl_eq (double mjd, double lat, double lng, double *ra, double *dec)
06085 {
06086 ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
06087 }
06088
06089 static void ecleq_aux (int sw, double mjd, double x, double y, double *p,double *q)
06090
06091
06092
06093
06094
06095 {
06096 static double lastmjd = -10000;
06097 static double seps, ceps;
06098 double sx, cx, sy, cy, ty;
06099
06100 if (mjd != lastmjd) {
06101 double eps;
06102 obliquity (mjd, &eps);
06103 seps = sin(eps);
06104 ceps = cos(eps);
06105 lastmjd = mjd;
06106 }
06107
06108 sy = sin(y);
06109 cy = cos(y);
06110 if (fabs(cy)<1e-20) cy = 1e-20;
06111 ty = sy/cy;
06112 cx = cos(x);
06113 sx = sin(x);
06114 *q = asin((sy*ceps)-(cy*seps*sx*sw));
06115 *p = atan(((sx*ceps)+(ty*seps*sw))/cx);
06116 if (cx<0) *p += PI;
06117 range (p, 2*PI);
06118 }
06119
06120 double geoc_lat (double phi);
06121 static void aaha_aux (double lat, double x, double y, double *p, double *q);
06122
06123
06124
06125
06126
06127 void aa_hadec (double lat, double alt,double az, double *ha, double *dec)
06128 {
06129 aaha_aux (lat, az, alt, ha, dec);
06130 if (*ha > PI)
06131 *ha -= 2*PI;
06132 }
06133
06134
06135
06136
06137
06138 void hadec_aa (double lat, double ha, double dec, double *alt, double *az)
06139 {
06140 aaha_aux (lat, ha, dec, az, alt);
06141 }
06142
06143 #ifdef NEED_GEOC
06144
06145
06146
06147 double geoc_lat (double phi)
06148 {
06149 #define MAXLAT degrad(89.9999)
06150 return (fabs(phi)>MAXLAT ? phi : atan(tan(phi)/1.00674));
06151 }
06152 #endif
06153
06154
06155
06156
06157
06158 static void aaha_aux (double lat, double x, double y, double *p, double *q)
06159
06160
06161
06162 {
06163 static double last_lat = -3434, slat, clat;
06164 double cap, B;
06165
06166 if (lat != last_lat) {
06167 slat = sin(lat);
06168 clat = cos(lat);
06169 last_lat = lat;
06170 }
06171
06172 solve_sphere (-x, PI/2-y, slat, clat, &cap, &B);
06173 *p = B;
06174 *q = PI/2 - acos(cap);
06175 }
06176
06177
06178
06179
06180
06181
06182
06183
06184
06185
06186
06187
06188
06189
06190
06191
06192 void solve_sphere (double A, double b, double cc, double sc, double *cap, double *Bp)
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
06210 B = PI - A;
06211 } else if (cc < -.99999) {
06212
06213 B = A;
06214 } else {
06215
06216
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 }
06233
06234
06235 #define TABSTART 1620.0
06236 #define TABEND 2004.0
06237 #define TABSIZ 385
06238
06239
06240
06241
06242
06243 static short dt[TABSIZ] = {
06244
06245 12400, 11900, 11500, 11000, 10600, 10200, 9800, 9500, 9100, 8800,
06246 8500, 8200, 7900, 7700, 7400, 7200, 7000, 6700, 6500, 6300,
06247 6200, 6000, 5800, 5700, 5500, 5400, 5300, 5100, 5000, 4900,
06248 4800, 4700, 4600, 4500, 4400, 4300, 4200, 4100, 4000, 3800,
06249
06250 3700, 3600, 3500, 3400, 3300, 3200, 3100, 3000, 2800, 2700,
06251 2600, 2500, 2400, 2300, 2200, 2100, 2000, 1900, 1800, 1700,
06252 1600, 1500, 1400, 1400, 1300, 1200, 1200, 1100, 1100, 1000,
06253 1000, 1000, 900, 900, 900, 900, 900, 900, 900, 900,
06254
06255 900, 900, 900, 900, 900, 900, 900, 900, 1000, 1000,
06256 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1100, 1100, 1100,
06257 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100,
06258 1100, 1100, 1100, 1100, 1200, 1200, 1200, 1200, 1200, 1200,
06259
06260 1200, 1200, 1200, 1200, 1300, 1300, 1300, 1300, 1300, 1300,
06261 1300, 1400, 1400, 1400, 1400, 1400, 1400, 1400, 1500, 1500,
06262 1500, 1500, 1500, 1500, 1500, 1600, 1600, 1600, 1600, 1600,
06263 1600, 1600, 1600, 1600, 1600, 1700, 1700, 1700, 1700, 1700,
06264
06265 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700,
06266 1700, 1700, 1600, 1600, 1600, 1600, 1500, 1500, 1400, 1400,
06267
06268 1370, 1340, 1310, 1290, 1270, 1260, 1250, 1250, 1250, 1250,
06269 1250, 1250, 1250, 1250, 1250, 1250, 1250, 1240, 1230, 1220,
06270
06271 1200, 1170, 1140, 1110, 1060, 1020, 960, 910, 860, 800,
06272 750, 700, 660, 630, 600, 580, 570, 560, 560, 560,
06273 570, 580, 590, 610, 620, 630, 650, 660, 680, 690,
06274 710, 720, 730, 740, 750, 760, 770, 770, 780, 780,
06275
06276 788, 782, 754, 697, 640, 602, 541, 410, 292, 182,
06277 161, 10, -102, -128, -269, -324, -364, -454, -471, -511,
06278 -540, -542, -520, -546, -546, -579, -563, -564, -580, -566,
06279 -587, -601, -619, -664, -644, -647, -609, -576, -466, -374,
06280
06281 -272, -154, -2, 124, 264, 386, 537, 614, 775, 913,
06282 1046, 1153, 1336, 1465, 1601, 1720, 1824, 1906, 2025, 2095,
06283 2116, 2225, 2241, 2303, 2349, 2362, 2386, 2449, 2434, 2408,
06284 2402, 2400, 2387, 2395, 2386, 2393, 2373, 2392, 2396, 2402,
06285
06286 2433, 2483, 2530, 2570, 2624, 2677, 2728, 2778, 2825, 2871,
06287 2915, 2957, 2997, 3036, 3072, 3107, 3135, 3168, 3218, 3268,
06288 3315, 3359, 3400, 3447, 3503, 3573, 3654, 3743, 3829, 3920,
06289 4018, 4117, 4223, 4337, 4449, 4548, 4646, 4752, 4853, 4959,
06290
06291 5054, 5138, 5217, 5296, 5379, 5434, 5487, 5532, 5582, 5630,
06292 5686, 5757, 5831, 5912, 5998, 6078,
06293
06294 6163, 6230,
06295
06296 6296, 6420,
06297 6510, 6600, 6700, 6800, 6900
06298
06299
06300
06301
06302
06303 };
06304
06305
06306
06307
06308
06309 double deltat(double mjd)
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
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
06336
06337
06338 B = 0.01*(Y - 2000.0);
06339 ans = (23.58 * B + 100.3)*B + 101.6;
06340 } else {
06341
06342 B = 0.01*(Y - 2000.0) + 3.75;
06343 ans = 35.0 * B * B + 40.;
06344 }
06345 return(ans);
06346 }
06347
06348
06349
06350
06351
06352
06353
06354 double a=Y;
06355 p = floor(a);
06356 iy =(int)( p - TABSTART);
06357
06358
06359 ans = dt[iy];
06360 k = iy + 1;
06361 if( k >= TABSIZ )
06362 goto done;
06363
06364
06365
06366 p = Y - p;
06367
06368
06369
06370 ans += p*(dt[k] - dt[iy]);
06371 if( (iy-1 < 0) || (iy+2 >= TABSIZ) )
06372 goto done;
06373
06374
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
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
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
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
06411
06412
06413
06414
06415
06416
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 }
06425
06426
06427
06428
06429
06430
06431
06432
06433
06434
06435
06436
06437
06438
06439
06440
06441
06442