Logo Search packages:      
Sourcecode: nautic version File versions  Download package

starpos.cpp


#include <cmath>
#include "starpos.h"
#include <wx/arrimpl.cpp> // This is a magic incantation which must be done!
#define WILLIAMS 1
#define SIMON 0

#define TABSTART 1620
#define TABEND 2013
#define TABSIZ (TABEND - TABSTART + 1)
#define double double
#define COS cos
#define SIN sin
#define mods3600(x) ((x) - 1.296e6 * floor ((x)/1.296e6))
#define mod3600(x) ((x) - 1296000. * floor ((x)/1296000.))
#define NCON 89

/* Earth radii per au */
#define DISFAC1 2.3454780e4

/* cosine of 90 degrees 50 minutes: */
#define COSSUN -0.014543897651582657
/* cosine of 90 degrees 34 minutes: */
#define COSZEN -9.8900378587411476e-3

#define RTOH (12.0/PI)

double DTR = 1.7453292519943295769e-2;
double RTD = 5.7295779513082320877e1;
double RTS = 2.0626480624709635516e5; /* arc seconds per radian */
double STR = 4.8481368110953599359e-6; /* radians per arc second */
double PI = 3.14159265358979323846;


Starpos::Starpos()
{
}

double JD;
double  rearth[3];
double dtgiven = 0.0;

double  dradt;
double  ddecdt;

double jdeps = -1.0; /* Date for which obliquity was last computed */
double eps = 0.0; /* The computed obliquity in radians */
double coseps = 0.0; /* Cosine of the obliquity */
double sineps = 0.0; /* Sine of the obliquity */



#define  J2000  2451545.0     /* 2000 January 1.5 */
/*
         
    starpos.cpp star position implementation
    Copyright (C) 2011 Enas Giovanni <gio.enas@alice.it>
 
    This file is part of Nautic.

    Nautic is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Nautic is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Nautic.  If not, see <http://www.gnu.org/licenses/>.

 */

#define  B1950  2433282.423   /* 1950 January 0.923 Besselian epoch */
#define  J1900  2415020.0     /* 1900 January 0, 12h UT */

double  SE; /* earth-sun distance */
double  SO; /* object-sun distance */
double  EO; /* object-earth distance */

double TDT;
double UT;

int jdflag;
int prtflg = 1;

struct orbit *elobject;
static char starnam[80] = {'s','t','a','r','.','c','a','t','\0'};
static char orbnam[80] = {'o','r','b','i','t','.','c','a','t','\0'};
const char * pathname = "/usr/share/nautic/star";
//static int linenum = 1;




static double pAcof[] = {
#if 1
  /* Corrections to Williams (1994) introduced in DE403.  */
 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
 -0.235316, 0.076, 110.5414, 50287.91959
#else
 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
 -0.235316, 0.076, 110.5407, 50287.70000
#endif
 };

#if SIMON
/* Precession coefficients from Simon et al: */
static double pAcof[] = {
 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
 -0.235316, 0.07732, 111.2022, 50288.200 };
#endif
#if LASKAR
/* Precession coefficients taken from Laskar's paper: */
static double pAcof[] = {
 -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
 -0.235316, 0.07732, 111.1971, 50290.966 };
#endif
#if WILLIAMS
/* Pi from Williams' 1994 paper, in radians.  No change in DE403.  */
static double nodecof[] = {
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10,
-3.54e-9, -1.8103e-7,  1.26e-7,  7.436169e-5,
-0.04207794833,  3.052115282424};
/* pi from Williams' 1994 paper, in radians.  No change in DE403.  */
static double inclcof[] = {
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
-5.4000441e-11, 1.32115526e-9, -6.012e-7, -1.62442e-5,
 0.00227850649, 0.0 };
#endif
#if SIMON
static double nodecof[] = {
6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 1.9e-10,
-3.54e-9, -1.8103e-7, 2.579e-8, 7.4379679e-5,
-0.0420782900, 3.0521126906};

static double inclcof[] = {
1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11,
-5.4000441e-11, 1.32115526e-9, -5.99908e-7, -1.624383e-5,
 0.002278492868, 0.0 };
#endif

/* Node and inclination of the earth's orbit computed from
 * Laskar's data as done in Bretagnon and Francou's paper.
 * Units are radians.
 */


short dt[TABSIZ] = {
/* 1620.0 thru 1659.0 */
12400, 11900, 11500, 11000, 10600, 10200, 9800, 9500, 9100, 8800,
8500, 8200, 7900, 7700, 7400, 7200, 7000, 6700, 6500, 6300,
6200, 6000, 5800, 5700, 5500, 5400, 5300, 5100, 5000, 4900,
4800, 4700, 4600, 4500, 4400, 4300, 4200, 4100, 4000, 3800,
/* 1660.0 thru 1699.0 */
3700, 3600, 3500, 3400, 3300, 3200, 3100, 3000, 2800, 2700,
2600, 2500, 2400, 2300, 2200, 2100, 2000, 1900, 1800, 1700,
1600, 1500, 1400, 1400, 1300, 1200, 1200, 1100, 1100, 1000,
1000, 1000, 900, 900, 900, 900, 900, 900, 900, 900,
/* 1700.0 thru 1739.0 */
900, 900, 900, 900, 900, 900, 900, 900, 1000, 1000,
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1100, 1100, 1100,
1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100, 1100,
1100, 1100, 1100, 1100, 1200, 1200, 1200, 1200, 1200, 1200,
/* 1740.0 thru 1779.0 */
1200, 1200, 1200, 1200, 1300, 1300, 1300, 1300, 1300, 1300,
1300, 1400, 1400, 1400, 1400, 1400, 1400, 1400, 1500, 1500,
1500, 1500, 1500, 1500, 1500, 1600, 1600, 1600, 1600, 1600,
1600, 1600, 1600, 1600, 1600, 1700, 1700, 1700, 1700, 1700,
/* 1780.0 thru 1799.0 */
1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700,
1700, 1700, 1600, 1600, 1600, 1600, 1500, 1500, 1400, 1400,
/* 1800.0 thru 1819.0 */
1370, 1340, 1310, 1290, 1270, 1260, 1250, 1250, 1250, 1250,
1250, 1250, 1250, 1250, 1250, 1250, 1250, 1240, 1230, 1220,
/* 1820.0 thru 1859.0 */
1200, 1170, 1140, 1110, 1060, 1020, 960, 910, 860, 800,
750, 700, 660, 630, 600, 580, 570, 560, 560, 560,
570, 580, 590, 610, 620, 630, 650, 660, 680, 690,
710, 720, 730, 740, 750, 760, 770, 770, 780, 780,
/* 1860.0 thru 1899.0 */
788, 782, 754, 697, 640, 602, 541, 410, 292, 182,
161, 10, -102, -128, -269, -324, -364, -454, -471, -511,
-540, -542, -520, -546, -546, -579, -563, -564, -580, -566,
-587, -601, -619, -664, -644, -647, -609, -576, -466, -374,
/* 1900.0 thru 1939.0 */
-272, -154, -2, 124, 264, 386, 537, 614, 775, 913,
1046, 1153, 1336, 1465, 1601, 1720, 1824, 1906, 2025, 2095,
2116, 2225, 2241, 2303, 2349, 2362, 2386, 2449, 2434, 2408,
2402, 2400, 2387, 2395, 2386, 2393, 2373, 2392, 2396, 2402,
/* 1940.0 thru 1979.0 */
 2433, 2483, 2530, 2570, 2624, 2677, 2728, 2778, 2825, 2871,
 2915, 2957, 2997, 3036, 3072, 3107, 3135, 3168, 3218, 3268,
 3315, 3359, 3400, 3447, 3503, 3573, 3654, 3743, 3829, 3920,
 4018, 4117, 4223, 4337, 4449, 4548, 4646, 4752, 4853, 4959,
/* 1980.0 thru 2003.0 */
 5054, 5138, 5217, 5296, 5379, 5434, 5487, 5532, 5582, 5630,
 5686, 5757, 5831, 5912, 5998, 6078, 6163, 6230, 6297, 6347,
 6383, 6409, 6430, 6447,
/* Extrapolated values */
                         6456, 6600, 6700, 6800, 6900, 7000,
 7100, 7200, 7300, 7400,
};

struct plantbl moonlat = {
 14,
  {  0, 26, 29,  8,  3,  5,  0,  0,  0,  6,  5,  3,  5,  1,  0,  0,  0,  0,},
 3,
 argsl,
 tabll,
 tabbl,
 tabrl,
 0.0000000000000000e+00,
 3.6525000000000000e+06,
 1.0000000000000000e-04,
};

struct plantbl moonlr = {
 14,
  {  3, 26, 29, 23,  5, 10,  0,  0,  0,  8,  4,  4,  6,  2,  0,  0,  0,  0,},
 3,
 argsm,
 tablm,
 tabbm,
 tabrm,
 2.5735686895300000e-03,
 3.6525000000000000e+06,
 1.0000000000000000e-04,
};

struct plantbl ear404 = {
 12,
  {  2, 11, 14, 19,  6, 10,  2,  2,  0,  1,  1,  4,  0,  0,  0,  0,  0,  0,},
 3,
 args,
 tabl,
 tabb,
 tabr,
 1.0001398729597080e+00,
 3.6525000000000000e+06,
 1.0000000000000000e-04,
};

struct plantbl mar404 = {
 9,
  {  0,  5, 12, 24,  9,  7,  3,  2,  0,},
 5,
 argsmr,
 tablmr,
 tabbmr,
 tabrmr,
 1.5303348827100001e+00,
 3652500.0,
 1.0
};

struct plantbl ven404 = {
 9,
  {  5, 14, 13,  8,  4,  5,  1,  0,  0,},
 5,
 argsv,
 tablv,
 tabbv,
 tabrv,
 7.2332982000000001e-01,
 3652500.0,
 1.0
};

struct plantbl jup404 = {
 9,
  {  0,  0,  1,  0,  9, 16,  7,  5,  0,},
 6,
 argsj,
 tablj,
 tabbj,
 tabrj,
 5.2026032092000003e+00,
 3652500.0,
 1.0
};

struct plantbl sat404 = {
 9,
  {  0,  0,  1,  0,  8, 18,  9,  5,  0,},
 7,
 argsst,
 tablst,
 tabbst,
 tabrst,
 9.5575813548599999e+00,
 3652500.0,
 1.0
};

00299 struct orbit
      {
      char obname[16]; /* name of the object */
      double epoch;     /* epoch of orbital elements */
      double i;   /* inclination    */
      double W;   /* longitude of the ascending node */
      double w;   /* argument of the perihelion */
      double a;   /* mean distance (semimajor axis) */
      double dm;  /* daily motion */
      double ecc; /* eccentricity */
      double M;   /* mean anomaly */
      double equinox;   /* epoch of equinox and ecliptic */
      double mag; /* visual magnitude at 1AU from earth and sun */
      double sdiam;     /* equatorial semidiameter at 1au, arc seconds */
/* The following used by perterbation formulas: */
      struct plantbl *ptable;
      double L;   /* computed mean longitude */
      double r;   /* computed radius vector */
      double plat;      /* perturbation in ecliptic latitude */
      };

00320 struct star
      {
      char obname[32];  /* Object name (31 chars) */
      double epoch;           /* Epoch of coordinates */
      double ra;        /* Right Ascension, radians */
      double dec;       /* Declination, radians */
      double px;        /* Parallax, radians */
      double mura;            /* proper motion in R.A., rad/century */
      double mudec;           /* proper motion in Dec., rad/century */
      double v;         /* radial velocity, km/s */
      double equinox;         /* Epoch of equinox and ecliptic */
      double mag;       /* visual magnitude */
      };

      static double A[3] = {-1.62557e-6, -3.1919e-7, - 1.3843e-7};
static double Ad[3] = {1.244e-3, -1.579e-3, -6.60e-4};

/* Transformation matrix for unit direction vector,
 * and motion vector in arc seconds per century
 */
static double Mat[36] = {
0.9999256782, -0.0111820611, -4.8579477e-3,
      2.42395018e-6, -2.710663e-8, -1.177656e-8,
0.0111820610, 0.9999374784, -2.71765e-5,
      2.710663e-8, 2.42397878e-6, -6.587e-11,
4.8579479e-3, -2.71474e-5, 0.9999881997,
      1.177656e-8, -6.582e-11, 2.42410173e-6,
-5.51e-4, -0.238565, 0.435739,
      0.99994704, -0.01118251, -4.85767e-3,
0.238514, -2.667e-3, -8.541e-3,
      0.01118251, 0.99995883, -2.718e-5,
-0.435623, 0.012254, 2.117e-3,
      4.85767e-3, -2.714e-5, 1.00000956
};

char *months[12] = {
            (char*)"January",
            (char*)"February",
            (char*)"March",
            (char*)"April",
            (char*)"May",
            (char*)"June",
            (char*)"July",
            (char*)"August",
            (char*)"September",
            (char*)"October",
            (char*)"November",
            (char*)"December"
};

char *days[7] = {
      (char*)"Sunday",
      (char*)"Monday",
      (char*)"Tuesday",
      (char*)"Wednesday",
      (char*)"Thursday",
      (char*)"Friday",
      (char*)"Saturday"
};



struct orbit earth =
{
  "Earth          ",
  2446800.5,
  0.0,
  0.0,
  102.884,
  0.999999,
  0.985611,
  0.016713,
  1.1791,
  2446800.5,
  -3.86,
  0.0,
  &ear404,
  0.0,
  0.0,
  0.0
};

struct orbit venus = {
"Venus          ",
2446800.5,
3.3946,
76.561,
54.889,
0.723329,
1.60214,
0.006757,
9.0369,
2446800.5,
/* Note the calculated apparent visual magnitude for Venus
 * is not very accurate.
 */
-4.40,
8.34,
&ven404,
0.0,
0.0,
0.0
};

struct orbit mars = {
"Mars           ",
2446800.5,
1.8498,
49.457,
286.343,
1.523710,
0.524023,
0.093472,
53.1893,
2446800.5,
-1.52,
4.68,
&mar404,
0.0,
0.0,
0.0
};

struct orbit jupiter = {
"Jupiter        ",
2446800.5,
1.3051,
100.358,
275.129,
5.20265,
0.0830948,
0.048100,
344.5086,
2446800.5,
-9.40,
98.44,
&jup404,
0.0,
0.0,
0.0
};

struct orbit saturn = {
"Saturn         ",
2446800.5,
2.4858,
113.555,
337.969,
9.54050,
0.0334510,
0.052786,
159.6327,
2446800.5,
-8.88,
82.73,
&sat404,
0.0,
0.0,
0.0
};

static double freqs[] =
{
/* Arc sec per 10000 Julian years.  */
  53810162868.8982,
  21066413643.3548,
  12959774228.3429,
  6890507749.3988,
  1092566037.7991,
  439960985.5372,
  154248119.3933,
  78655032.0744,
  52272245.1795
};

static double phases[] =
{
/* Arc sec.  */
  252.25090552 * 3600.,
  181.97980085 * 3600.,
  100.46645683 * 3600.,
  355.43299958 * 3600.,
  34.35151874 * 3600.,
  50.07744430 * 3600.,
  314.05500511 * 3600.,
  304.34866548 * 3600.,
  860492.1546,
};

short nt[105*9] = {
 0, 0, 0, 0, 2, 2062, 2,-895, 5,
-2, 0, 2, 0, 1, 46, 0,-24, 0,
 2, 0,-2, 0, 0, 11, 0, 0, 0,
-2, 0, 2, 0, 2,-3, 0, 1, 0,
 1,-1, 0,-1, 0,-3, 0, 0, 0,
 0,-2, 2,-2, 1,-2, 0, 1, 0,
 2, 0,-2, 0, 1, 1, 0, 0, 0,
 0, 0, 2,-2, 2,-13187,-16, 5736,-31,
 0, 1, 0, 0, 0, 1426,-34, 54,-1,
 0, 1, 2,-2, 2,-517, 12, 224,-6,
 0,-1, 2,-2, 2, 217,-5,-95, 3,
 0, 0, 2,-2, 1, 129, 1,-70, 0,
 2, 0, 0,-2, 0, 48, 0, 1, 0,
 0, 0, 2,-2, 0,-22, 0, 0, 0,
 0, 2, 0, 0, 0, 17,-1, 0, 0,
 0, 1, 0, 0, 1,-15, 0, 9, 0,
 0, 2, 2,-2, 2,-16, 1, 7, 0,
 0,-1, 0, 0, 1,-12, 0, 6, 0,
-2, 0, 0, 2, 1,-6, 0, 3, 0,
 0,-1, 2,-2, 1,-5, 0, 3, 0,
 2, 0, 0,-2, 1, 4, 0,-2, 0,
 0, 1, 2,-2, 1, 4, 0,-2, 0,
 1, 0, 0,-1, 0,-4, 0, 0, 0,
 2, 1, 0,-2, 0, 1, 0, 0, 0,
 0, 0,-2, 2, 1, 1, 0, 0, 0,
 0, 1,-2, 2, 0,-1, 0, 0, 0,
 0, 1, 0, 0, 2, 1, 0, 0, 0,
-1, 0, 0, 1, 1, 1, 0, 0, 0,
 0, 1, 2,-2, 0,-1, 0, 0, 0,
 0, 0, 2, 0, 2,-2274,-2, 977,-5,
 1, 0, 0, 0, 0, 712, 1,-7, 0,
 0, 0, 2, 0, 1,-386,-4, 200, 0,
 1, 0, 2, 0, 2,-301, 0, 129,-1,
 1, 0, 0,-2, 0,-158, 0,-1, 0,
-1, 0, 2, 0, 2, 123, 0,-53, 0,
 0, 0, 0, 2, 0, 63, 0,-2, 0,
 1, 0, 0, 0, 1, 63, 1,-33, 0,
-1, 0, 0, 0, 1,-58,-1, 32, 0,
-1, 0, 2, 2, 2,-59, 0, 26, 0,
 1, 0, 2, 0, 1,-51, 0, 27, 0,
 0, 0, 2, 2, 2,-38, 0, 16, 0,
 2, 0, 0, 0, 0, 29, 0,-1, 0,
 1, 0, 2,-2, 2, 29, 0,-12, 0,
 2, 0, 2, 0, 2,-31, 0, 13, 0,
 0, 0, 2, 0, 0, 26, 0,-1, 0,
-1, 0, 2, 0, 1, 21, 0,-10, 0,
-1, 0, 0, 2, 1, 16, 0,-8, 0,
 1, 0, 0,-2, 1,-13, 0, 7, 0,
-1, 0, 2, 2, 1,-10, 0, 5, 0,
 1, 1, 0,-2, 0,-7, 0, 0, 0,
 0, 1, 2, 0, 2, 7, 0,-3, 0,
 0,-1, 2, 0, 2,-7, 0, 3, 0,
 1, 0, 2, 2, 2,-8, 0, 3, 0,
 1, 0, 0, 2, 0, 6, 0, 0, 0,
 2, 0, 2,-2, 2, 6, 0,-3, 0,
 0, 0, 0, 2, 1,-6, 0, 3, 0,
 0, 0, 2, 2, 1,-7, 0, 3, 0,
 1, 0, 2,-2, 1, 6, 0,-3, 0,
 0, 0, 0,-2, 1,-5, 0, 3, 0,
 1,-1, 0, 0, 0, 5, 0, 0, 0,
 2, 0, 2, 0, 1,-5, 0, 3, 0,
 0, 1, 0,-2, 0,-4, 0, 0, 0,
 1, 0,-2, 0, 0, 4, 0, 0, 0,
 0, 0, 0, 1, 0,-4, 0, 0, 0,
 1, 1, 0, 0, 0,-3, 0, 0, 0,
 1, 0, 2, 0, 0, 3, 0, 0, 0,
 1,-1, 2, 0, 2,-3, 0, 1, 0,
-1,-1, 2, 2, 2,-3, 0, 1, 0,
-2, 0, 0, 0, 1,-2, 0, 1, 0,
 3, 0, 2, 0, 2,-3, 0, 1, 0,
 0,-1, 2, 2, 2,-3, 0, 1, 0,
 1, 1, 2, 0, 2, 2, 0,-1, 0,
-1, 0, 2,-2, 1,-2, 0, 1, 0,
 2, 0, 0, 0, 1, 2, 0,-1, 0,
 1, 0, 0, 0, 2,-2, 0, 1, 0,
 3, 0, 0, 0, 0, 2, 0, 0, 0,
 0, 0, 2, 1, 2, 2, 0,-1, 0,
-1, 0, 0, 0, 2, 1, 0,-1, 0,
 1, 0, 0,-4, 0,-1, 0, 0, 0,
-2, 0, 2, 2, 2, 1, 0,-1, 0,
-1, 0, 2, 4, 2,-2, 0, 1, 0,
 2, 0, 0,-4, 0,-1, 0, 0, 0,
 1, 1, 2,-2, 2, 1, 0,-1, 0,
 1, 0, 2, 2, 1,-1, 0, 1, 0,
-2, 0, 2, 4, 2,-1, 0, 1, 0,
-1, 0, 4, 0, 2, 1, 0, 0, 0,
 1,-1, 0,-2, 0, 1, 0, 0, 0,
 2, 0, 2,-2, 1, 1, 0,-1, 0,
 2, 0, 2, 2, 2,-1, 0, 0, 0,
 1, 0, 0, 2, 1,-1, 0, 0, 0,
 0, 0, 4,-2, 2, 1, 0, 0, 0,
 3, 0, 2,-2, 2, 1, 0, 0, 0,
 1, 0, 2,-2, 0,-1, 0, 0, 0,
 0, 1, 2, 0, 1, 1, 0, 0, 0,
-1,-1, 0, 2, 1, 1, 0, 0, 0,
 0, 0,-2, 0, 1,-1, 0, 0, 0,
 0, 0, 2,-1, 2,-1, 0, 0, 0,
 0, 1, 0, 2, 0,-1, 0, 0, 0,
 1, 0,-2,-2, 0,-1, 0, 0, 0,
 0,-1, 2, 0, 1,-1, 0, 0, 0,
 1, 1, 0,-2, 1,-1, 0, 0, 0,
 1, 0,-2, 2, 0,-1, 0, 0, 0,
 2, 0, 0, 2, 0, 1, 0, 0, 0,
 0, 0, 2, 4, 2,-1, 0, 0, 0,
 0, 1, 0, 1, 0, 1, 0, 0, 0,
};

/* Greek letters
 */

#define NGREEK 24
const char *greek[NGREEK] = {
"alpha",
"beta",
"gamma",
"delta",
"epsilon",
"zeta",
"eta",
"theta",
"iota",
"kappa",
"lambda",
"mu",
"nu",
"xi",
"omicron",
"pi",
"rho",
"sigma",
"tau",
"upsilon",
"phi",
"chi",
"psi",
"omega",
};
const char *constel[NCON] = {
"And Andromedae",
"Ant Antliae",
"Aps Apodis",
"Aql Aquilae",
"Aqr Aquarii",
"Ari Arietis",
"Ara Arae",
"Aur Aurigae",
"Boo Bootis",
"Cae Caeli",
"Cam Camelopardalis",
"Can Cancri",           /* also abbreviated Cnc */
"Cap Capricorni",
"Car Carinae",
"Cas Cassiopeiae",
"Cen Centauri",
"Cep Cephei",
"Cet Ceti",
"Cha Chamaeleontis",
"Cir Circini",
"CMa Canis Majoris",
"CMi Canis Minoris",
"Cnc Cancri",
"Col Columbae",
"Com Comae Berenices",
"CrA Coronae Austrinae",
"CrB Coronae Borealis",
"Crt Crateris",
"Cru Crucis",
"Crv Corvi",
"CVn Canum Venaticorum",
"Cyg Cygni",
"Del Delphini",
"Dor Doradus",
"Dra Draconis",
"Equ Equulei",
"Eri Eridani",
"For Fornacis",
"Gem Geminorum",
"Gru Gruis",
"Her Herculis",
"Hor Horologii",
"Hya Hydrae",
"Hyi Hydri",
"Ind Indi",
"Lac Lacertae",
"Leo Leonis",
"Lep Leporis",
"Lib Librae",
"LMi Leonis Minoris",
"Lup Lupi",
"Lyn Lyncis",
"Lyr Lyrae",
"Men Mensae",
"Mic Microscopii",
"Mon Monocerotis",
"Mus Muscae",
"Nor Normae",
"Oct Octantis",
"Oph Ophiuchi",
"Ori Orionis",
"Pav Pavonis",
"Peg Pegasi",
"Per Persei",
"Phe Phoenicis",
"Pic Pictoris",
"PsA Piscis Austrini",
"Psc Piscium",
"Pup Puppis",
"Pyx Pyxidis",
"Ret Reticuli",
"Scl Sculptoris",
"Sco Scorpii",
"Sct Scuti",
"Ser Serpentis",
"Sex Sextantis",
"Sge Sagittae",
"Sgr Sagittarii",
"Tau Tauri",
"Tel Telescopii",
"TrA Trianguli Australis",
"Tri Trianguli",
"Tuc Tucanae",
"UMa Ursae Majoris",
"UMi Ursae Minoris",
"Vel Velorum",
"Vir Virginis",
"Vol Volantis",
"Vul Vulpeculae",
};

#define NBNDRIES 357
long bndries[4*NBNDRIES] = {
     0,  86400, 316800,  84,
 28800,  52200, 311400,  84,
 75600,  82800, 310200,  84,
 64800,  75600, 309600,  84,
     0,  28800, 306000,  16,
 33000,  38400, 295200,  10,
     0,  18000, 288000,  16,
 38400,  52200, 288000,  10,
 63000,  64800, 288000,  84,
 72600,  75600, 288000,  34,
     0,  12630, 277200,  16,
 41400,  48900, 277200,  10,
 59520,  63000, 270000,  84,
 72600,  74400, 270000,  16,
 28680,  33000, 264600,  10,
 33000,  40800, 264600,  34,
 46800,  59520, 252000,  84,
 11160,  12300, 244800,  14,
 73500,  74400, 241200,  34,
 40800,  43200, 239400,  34,
     0,   1200, 237600,  16,
 50400,  56400, 237600,  84,
 84900,  86400, 237600,  16,
 43200,  48600, 230400,  34,
 48600,  51900, 226800,  34,
 83400,  84900, 226800,  16,
 21960,  25200, 223200,  10,
 72000,  73500, 221400,  34,
 73932,  74160, 219300,  16,
 25200,  28680, 216000,  10,
 28680,  30300, 216000,  83,
 71160,  72000, 214200,  34,
 72000,  73932, 214200,  16,
 82320,  83400, 212700,  16,
     0,   8760, 210600,  14,
 69900,  71160, 208800,  34,
  6120,   6870, 207000,  14,
  8760,  11160, 205200,  14,
 11160,  11400, 205200,  10,
 80340,  82320, 202500,  16,
 18000,  21960, 201600,  10,
 50520,  51900, 199800,  83,
 51900,  69900, 199800,  34,
 11400,  12000, 198000,  10,
 79680,  80340, 198000,  16,
 74160,  79080, 197400,  16,
     0,   6120, 194400,  14,
 21960,  23400, 194400,  51,
 43500,  48600, 190800,  83,
 54900,  56700, 190800,  34,
 79080,  79680, 189900,  16,
 12000,  18000, 189000,  10,
 82320,  84000, 189000,  14,
 56700,  61200, 185400,  34,
  7350,   9060, 181800,  63,
 61200,  65640, 181800,  34,
     0,   4920, 180000,  14,
  4920,   6000, 180000,  63,
 23400,  24480, 180000,  51,
 84000,  86400, 180000,  14,
 48600,  50520, 174600,  83,
     0,   4020, 172800,  14,
 84900,  86400, 172800,  14,
 65430,  65640, 171000,  40,
 65640,  68700, 171000,  34,
 68700,  69000, 171000,  31,
  6000,   7350, 169200,  63,
 30300,  33000, 169200,  83,
   600,   3120, 165600,  14,
 43200,  43500, 162000,  83,
 24480,  26520, 160200,  51,
 78870,  79080, 158400,  31,
 78750,  78870, 157500,  31,
 69000,  69840, 156600,  31,
 33000,  36600, 151200,  83,
 36600,  38820, 144000,  83,
 55560,  56700, 144000,   8,
 56700,  58800, 144000,  40,
 33300,  34500, 143100,  51,
     0,   9060, 132300,   0,
  9060,   9240, 132300,  63,
 69690,  69840, 131400,  52,
 16200,  16890, 129600,  63,
 78240,  78750, 129600,  31,
 78750,  79200, 129600,  45,
 23520,  26520, 127800,   7,
 26520,  27900, 127800,  51,
     0,   7200, 126000,   0,
 79200,  82140, 126000,  45,
 82140,  82320, 124200,  45,
 82320,  84600, 124200,   0,
  9240,   9780, 122400,  63,
 38820,  39600, 122400,  83,
 43200,  44400, 122400,  30,
 27900,  33300, 120600,  51,
 33300,  35580, 120600,  49,
  2580,   5070, 118800,   0,
 54660,  55560, 118800,   8,
 84600,  85500, 115500,   0,
 44400,  47700, 115200,  30,
 85500,  86400, 112800,   0,
 50250,  50520, 110700,  30,
  8700,   9780, 110400,  81,
  9780,  16200, 110400,  63,
 16200,  17100, 108000,   7,
 65430,  69690, 108000,  52,
 39600,  43200, 104400,  83,
 70800,  75300, 104400,  31,
 17100,  21180, 102600,   7,
 35580,  37800, 102600,  49,
 47700,  50250, 102600,  30,
     0,    240, 100800,   0,
  5070,   6000, 100800,  81,
 21180,  23520, 100800,   7,
 28380,  28800, 100800,  38,
 75300,  78240, 100800,  31,
 69330,  70800,  99000,  31,
  6900,   8700,  98100,  81,
 58200,  58800,  97200,  26,
 54300,  54660,  93600,   8,
 54660,  58200,  93600,  26,
 66120,  67920,  93600,  52,
 38700,  39600,  91800,  49,
 67920,  69330,  91800,  52,
  6000,   6900,  90000,  81,
  2580,   3060,  85500,  67,
 37800,  38700,  84600,  49,
 76500,  77100,  84600,  88,
 20520,  21180,  82200,  78,
   240,    510,  79200,   0,
 57300,  57720,  79200,  74,
 21180,  22380,  77400,  38,
 71400,  72900,  76500,  88,
 67920,  69300,  75900,  88,
   510,   3060,  75600,   0,
 72900,  74040,  73800,  88,
 28110,  28380,  72000,  38,
 74040,  76500,  70200,  88,
 69300,  71400,  69000,  88,
 11820,  12120,  68400,   5,
 67920,  68400,  66600,  76,
 20520,  20760,  64800,  60,
 22380,  22710,  63000,  38,
 68400,  71400,  58200,  76,
 17880,  19200,  57600,  78,
 57300,  57900,  57600,  40,
 71400,  72900,  56700,  76,
 16620,  17880,  55800,  78,
 19200,  20160,  55800,  78,
 46200,  48600,  54000,  24,
 62100,  65700,  51600,  40,
 42720,  46200,  50400,  24,
 27000,  28110,  48600,  38,
 60300,  62100,  46200,  40,
     0,    510,  45000,  62,
 20160,  20760,  45000,  78,
 25200,  27000,  45000,  38,
 76020,  76800,  45000,  62,
 22710,  24960,  43200,  38,
 65700,  67920,  43200,  40,
 75150,  75780,  42600,  32,
 75780,  76020,  42600,  62,
 41460,  42720,  39600,  46,
 22470,  22710,  36000,  60,
 24960,  25200,  36000,  38,
 28110,  28530,  36000,  22,
 85800,  86400,  36000,  62,
  6000,  11820,  35700,   5,
 72510,  73080,  30600,  32,
 48600,  54300,  28800,   8,
 81900,  85800,  27000,  62,
 28530,  33300,  25200,  22,
 33300,  38700,  25200,  46,
 65700,  67184,  22500,  59,
 67184,  67920,  22500,   3,
 75000,  75150,  21600,  32,
 25200,  25260,  19800,  21,
 65700,  66330,  16200,  74,
 57900,  60300,  14400,  40,
 65700,  66330,  10800,  59,
 77280,  78000,   9900,  62,
     0,   7200,   7200,  67,
 66900,  67920,   7200,  74,
 73080,  75000,   7200,  32,
 75000,  76800,   7200,  35,
 76800,  77280,   7200,  62,
 79200,  81900,   7200,  62,
 78000,  79200,   6300,  62,
 25260,  25920,   5400,  21,
 12900,  16620,      0,  78,
 16620,  16800,      0,  60,
 25920,  29100,      0,  21,
 52800,  54300,      0,  86,
 64200,  65700,      0,  59,
  9540,  11820,  -6300,  17,
 11820,  12900,  -6300,  78,
 54300,  58560, -11700,  74,
 16800,  18300, -14400,  60,
 21000,  22470, -14400,  60,
 64200,  64680, -14400,  74,
 65700,  66900, -14400,  74,
 66900,  67920, -14400,   3,
 81900,  85800, -14400,  67,
 38700,  41460, -21600,  46,
 41460,  42600, -21600,  86,
     0,   1200, -25200,  67,
 85800,  86400, -25200,  67,
 51300,  52800, -28800,  86,
 57300,  58560, -28800,  59,
 72000,  73920, -32400,   3,
 76800,  78720, -32400,   4,
 61800,  64680, -36000,  59,
 21000,  29100, -39600,  55,
 17700,  18300, -39600,  36,
 18300,  21000, -39600,  60,
 29100,  30120, -39600,  42,
 34500,  38700, -39600,  75,
 42600,  46200, -39600,  86,
 63300,  63600, -42000,  59,
 67920,  72000, -43320,   3,
 17400,  17700, -52200,  36,
 73920,  76800, -54000,   4,
 61800,  65700, -57600,  74,
 65700,  67920, -57600,  73,
 30120,  30900, -61200,  42,
 58560,  58950, -65700,  59,
 30900,  32700, -68400,  42,
 38700,  39000, -68400,  27,
 58560,  58950, -69300,  59,
 56400,  57300, -72000,  48,
 45300,  46200, -79200,  29,
 46200,  51300, -79200,  86,
 32700,  35100, -86400,  42,
  6000,   9540, -87780,  17,
  9540,  13500, -87780,  36,
 39000,  42600, -88200,  27,
 42600,  45300, -88200,  29,
 51300,  53700, -88200,  48,
 58560,  60300, -88500,  59,
     0,   6000, -91800,  17,
 76800,  78720, -91800,  12,
 78720,  85800, -91800,   4,
 85800,  86400, -91800,  17,
 35100,  36900, -95400,  42,
 16920,  17400, -98100,  36,
 17400,  22020, -98100,  47,
 72000,  76800, -100800,  12,
 36900,  38100, -105000,  42,
 45300,  53700, -106200,  42,
 53700,  56400, -106200,  48,
 56400,  57600, -106200,  72,
 16500,  16920, -108000,  36,
 60300,  63360, -108000,  59,
 63360,  64200, -108000,  77,
 38100,  39000, -112200,  42,
 22020,  26520, -118800,  20,
 44100,  45300, -118800,  42,
 39000,  44100, -126000,  42,
 12600,  13500, -129600,  37,
 30120,  33720, -132300,  69,
 15360,  16500, -133200,  36,
 64200,  69000, -133200,  77,
 76800,  82800, -133200,  66,
 82800,  84000, -133200,  71,
 10800,  12600, -142500,  37,
 33720,  39600, -143100,   1,
     0,   6000, -144000,  71,
  6000,  10800, -144000,  37,
 13920,  15360, -144000,  36,
 84000,  86400, -144000,  71,
 51000,  53700, -151200,  15,
 56400,  57600, -151200,  50,
 57600,  59115, -151200,  72,
 17400,  18000, -154800,   9,
 18000,  23700, -154800,  23,
 28800,  30120, -154800,  68,
 12300,  13920, -158400,  36,
 59115,  64200, -163800,  72,
 64200,  69000, -163800,  25,
 69000,  73200, -163800,  77,
 73200,  76800, -163800,  54,
 10800,  12300, -165600,  36,
 16200,  17400, -167400,   9,
 55200,  56400, -172800,  50,
     0,   8400, -173400,  64,
  9600,  10800, -176400,  36,
 14700,  15360, -176400,  41,
 15360,  16200, -176400,   9,
 76800,  79200, -180000,  39,
 21600,  28800, -182700,  68,
 28800,  29400, -182700,  85,
  8700,   9600, -183600,  36,
 13800,  14700, -183600,  41,
     0,   6600, -185400,  64,
 21600,  22200, -189000,  13,
 29400,  30420, -190800,  85,
 12600,  13800, -191400,  41,
 13800,  14400, -191400,  33,
     0,   5700, -192600,  64,
  7800,   8700, -194400,  36,
 16200,  18000, -194400,  65,
 54180,  55200, -194400,  50,
 30420,  31800, -196200,  85,
 22200,  23400, -198000,  13,
 42600,  46200, -198000,  15,
 51000,  54180, -198000,  50,
 54180,  55200, -198000,  57,
 14400,  15600, -203400,  33,
 31800,  39600, -203400,  85,
 39600,  40500, -203400,  15,
 63000,  64800, -205200,   6,
 64800,  73200, -205200,  79,
 79200,  84000, -205200,  39,
 11520,  12600, -207000,  41,
 18000,  19800, -207000,  65,
 23400,  24600, -208800,  13,
     0,   4800, -210600,  64,
  4800,   7800, -210600,  36,
 84000,  86400, -210600,  64,
 15600,  16500, -212400,  33,
 55200,  59115, -216000,  57,
 73200,  76800, -216000,  44,
 19800,  21600, -219600,  65,
 54600,  55200, -219600,  19,
 59115,  59700, -219600,   6,
 53700,  54600, -228900,  19,
 59700,  60300, -228900,   6,
 21600,  24600, -230400,  65,
 24600,  32520, -230400,  13,
 40500,  42600, -230400,  15,
 42600,  46200, -230400,  28,
 46200,  52320, -230400,  15,
 48600,  49200, -234000,  19,
 60300,  60600, -234000,   6,
  7800,  11520, -243000,  41,
 11520,  16500, -243000,  70,
 53100,  53700, -243000,  19,
 60600,  63000, -243000,   6,
 63000,  64800, -243000,  61,
 79200,  84000, -243000,  82,
 16500,  23700, -252000,  33,
 49200,  53100, -252000,  19,
 53100,  61200, -252000,  80,
     0,   4800, -270000,  82,
 12600,  16500, -270000,  43,
 23700,  32520, -270000,  87,
 32520,  40500, -270000,  13,
 40500,  49200, -270000,  56,
 64800,  76800, -270000,  61,
 76800,  84000, -270000,  44,
 84000,  86400, -270000,  82,
  2700,   4800, -273600,  82,
     0,  12600, -297000,  43,
 27600,  49200, -297000,  18,
 49200,  64800, -297000,   2,
 12600,  27600, -306000,  53,
     0,  86400, -324000,  58,
};


long cyear = 1986L;
static int month = 1;
static double day = 1.0;
short yerend = 0;

double  pq; /* cosine of sun-object-earth angle */
double  ep; /* -cosine of sun-earth-object angle */
double  qe; /* cosine of earth-sun-object angle */
/* correction vector, saved for display  */
double dp[3];
double jvearth = -1.0;
double vearth[3] = {0.0};
double Clightaud; /* C in au/day  */
double emrat = 81.300585;

static double Args[NARGS];
static double LP_equinox;
static double NF_arcsec;
static double Ea_arcsec;
static double pA_precession;
static double ss[NARGS][31];
static double cc[NARGS][31];

double jdnut = -1.0;    /* time to which the nutation applies */
double nutl;      /* nutation in longitude (radians) */
double nuto;      /* nutation in obliquity (radians) */
double j;   /* Julian date and fraction */
      /* East longitude of observer, degrees */
//double tlong = 0;//-71.13;  /* Cambridge, Massachusetts */
//double tlat = 0;//42.38; /* geocentric */
//double glat = 0;//42.17; /* geodetic */

double r_trnsit;
double r_rise;
double r_set;
int f_trnsit;

double t_rise;
double t_trnsit;
double t_set;
//int objnum = 88;                  /* I.D. number of object */
double obpolar[3];
double eapolar[3];
double trho = 0.9985;
double aearth = 6378137.;  /* Radius of the earth, in meters.  */
double au = 1.49597870691e8; /* Astronomical unit, in kilometers.  */
static double DISFAC = 2.3454780e4;
double atpress; /* millibars */
double attemp; /* degrees C */
static double flat = 298.257222;
static double height = 0.0;
double Clight = 2.99792458e5;  /* Speed of light, km/sec  */

const char *intfmt = "%d";
const char *lngfmt = "%ld";
const char *dblfmt = "%lf";
const char *strfmt = "%s";
double Rearth;
struct star * starobject;
static double djd = 1.0;
static int ntab = 1;

/* Space for star description read from a disc file.
 */
struct star fstar;

/* Space for orbit read from a disc file. Entering 99 for the
 * planet number yields a prompt for a file name containg ASCII
 * strings specifying the elements.
 */
struct orbit forbit;

static int hours = 0;
static int minutes = 0;
static double seconds = 0.0;
bool flag = false;
bool fPlanet = false;
static double  ra;      /* Right Ascension */
static double  dec;     /* Declination */

double robject[3] = {0.0};

/*////////////////////////////////funzioni////////////////////////////////////////////////////*/

int Starpos::trnsit(double J, double lha, double dec)
{
      double x, y, z, N, D, TPI;
      double lhay, cosdec, sindec, coslat, sinlat;

      f_trnsit = 0;
      TPI = 2.0*PI;
      /* Initialize to no-event flag value. */
      r_rise = -10.0;
      r_set = -10.0;
      /* observer's geodetic latitude, in radians */
      x = glat * DTR;
      coslat = cos(x);
      sinlat = sin(x);

      cosdec = cos(dec);
      sindec = sin(dec);

      /*  x = (J - floor(J-0.5) - 0.5) * TPI; */
       x = floor(J) - 0.5;
       x = (J - x) * TPI;
      /* adjust local hour angle */
      y = lha;
      /* printf ("%.7f,", lha); */
      while( y < -PI )
            y += TPI;
      while( y > PI )
            y -= TPI;
      lhay = y;
      y =  y/( -dradt/TPI + 1.00273790934);
      r_trnsit = x - y;
      /* printf ("rt %.7f ", r_trnsit); */
      /* Ordinarily never print here.  */
      if( prtflg > 1 )
            {
            //printf( "approx local meridian transit" );
            hms( r_trnsit );

            //printf( "UT  (date + %.5f)\n", r_trnsit/TPI );
            }
      if( (coslat == 0.0) || (cosdec == 0.0) )
            goto norise;

      /* The time at which the upper limb of the body meets the
       * horizon depends on the body's angular diameter.
       */
      switch( objnum )
            {
      /* Sun */
            case 0: N = COSSUN; break;

      /* Moon, elevation = -34' - semidiameter + parallax
       * semidiameter = 0.272453 * parallax + 0.0799"
       */
            case 3:
                  N = 1.0/(DISFAC*obpolar[2]);
                  D = asin( N ); /* the parallax */
                  N =  - 9.890199094634534e-3
                        + (1.0 - 0.2725076)*D
                        - 3.874e-7;
                  N = sin(N);
                  break;

      /* Other object */
            default: N = COSZEN; break;
            }
      y = (N - sinlat*sindec)/(coslat*cosdec);

      if( (y < 1.0) && (y > -1.0) )
            {
            f_trnsit = 1;
      /* Derivative of y with respect to declination
       * times rate of change of declination:
       */
            z = -ddecdt*(sinlat + COSZEN*sindec);
            z /= TPI*coslat*cosdec*cosdec;
      /* Derivative of acos(y): */
            z /= sqrt( 1.0 - y*y);
            y = acos(y);
            D = -dradt/TPI + 1.00273790934;
            r_rise = x - (lhay + y)*(1.0 + z)/D;
            r_set = x - (lhay - y)*(1.0 - z)/D;
            /* Ordinarily never print here.  */
                  if( prtflg > 1 )
                  {
                  //printf( "rises approx " );
                  hms(r_rise);
                  //printf( "UT,  sets approx " );
                  hms(r_set);
                  //printf( "UT\n" );
                  }
            }
      norise:           ;
      return(0);
}


void Starpos::iter_func(double t, int (* func)())

{
  int prtsave;

  prtsave = prtflg;
  prtflg = 0;
  JD = t;
  update(); /* Set UT and TDT */
  kepler( TDT, &earth, rearth, eapolar );
  (*func)();
  prtflg = prtsave;
}

double Starpos::sidrlt(double j,double tlong )
//double j; /* Julian date and fraction */
//double tlong;   /* East longitude of observer, degrees */
{
double jd0;     /* Julian day at midnight Universal Time */
double secs;   /* Time of day, UT seconds since UT midnight */
double eqeq, gmst, jd, T0, msday;
/*long il;*/

  /* Julian day at given UT */
jd = j;
jd0 = floor(jd);
secs = j - jd0;
if( secs < 0.5 )
      {
      jd0 -= 0.5;
      secs += 0.5;
      }
else
      {
      jd0 += 0.5;
      secs -= 0.5;
      }
secs *= 86400.0;

  /* Julian centuries from standard epoch J2000.0 */
/* T = (jd - J2000)/36525.0; */
  /* Same but at 0h Universal Time of date */
T0 = (jd0 - J2000)/36525.0;

/* The equation of the equinoxes is the nutation in longitude
 * times the cosine of the obliquity of the ecliptic.
 * We already have routines for these.
 */
nutlo(TDT);
epsiln(TDT);
/* nutl is in radians; convert to seconds of time
 * at 240 seconds per degree
 */
eqeq = 240.0 * RTD * nutl * coseps;
  /* Greenwich Mean Sidereal Time at 0h UT of date */
#if 1
#if 0
  /* J. G. Williams, "Contributions to the Earth's obliquity rate, precession,
     and nutation,"  Astronomical Journal 108, p. 711 (1994). */
gmst = (((-2.0e-6*T0 - 3.e-7)*T0 + 9.27695e-2)*T0 + 8640184.7928613)*T0
       + 24110.54841;
/* UT days per sidereal day at date T0  */
msday = (((-(4. * 2.0e-6)*T0 - (3. * 3.e-7))*T0 + (2. * 9.27695e-2))*T0
          + 8640184.7928613)/(86400.*36525.) + 1.0;
#else
/* Corrections to Williams (1994) introduced in DE403.  */
gmst = (((-2.0e-6*T0 - 3.e-7)*T0 + 9.27701e-2)*T0 + 8640184.7942063)*T0
       + 24110.54841;
msday = (((-(4. * 2.0e-6)*T0 - (3. * 3.e-7))*T0 + (2. * 9.27701e-2))*T0
          + 8640184.7942063)/(86400.*36525.) + 1.0;
#endif
#else
/* This is the 1976 IAU formula. */
gmst = (( -6.2e-6*T0 + 9.3104e-2)*T0 + 8640184.812866)*T0 + 24110.54841;
/* mean solar days per sidereal day at date T0  */
msday = 1.0 + ((-1.86e-5*T0 + 0.186208)*T0 + 8640184.812866)/(86400.*36525.);
#endif
  /* Local apparent sidereal time at given UT */
gmst = gmst + msday*secs + eqeq + 240.0*tlong;
  /* Sidereal seconds modulo 1 sidereal day */
gmst = gmst - 86400.0 * floor( gmst/86400.0 );
/*
 * il = gmst/86400.0;
 * gmst = gmst - 86400.0 * il;
 * if( gmst < 0.0 )
 *    gmst += 86400.0;
 */
return( gmst );
}

int Starpos::diurab(double last,double *ra,double *dec )
//double last;    /* local apparent sidereal time, radians */
//double *ra;     /* right ascension, radians */
//double *dec;    /* declination, radians */
{
double lha, coslha, sinlha, cosdec, sindec;
double coslat, N, D;
//trho =0;

lha = last - *ra;
coslha = cos(lha);
sinlha = sin(lha);
cosdec = cos(*dec);
sindec = sin(*dec);
coslat = cos( DTR*tlat );

if( cosdec != 0.0 )
      N = 1.5472e-6*trho*coslat*coslha/cosdec;
else
      N = 0.0;
*ra += N;

D = 1.5472e-6*trho*coslat*sinlha*sindec;
*dec += D;
/*
if( prtflg )
      {
printf( "diurnal aberration dRA %.3fs dDec %.2f\"\n", RTS*N/15.0, RTS*D );
      }*/
return(0);
}

int Starpos::diurpx(double last,double *ra,double *dec,double dist )
//double last;  /* local apparent sidereal time, radians */
//double *ra;  /* right ascension, radians */
//double *dec; /* declination, radians */
//double dist; /* Earth - object distance, au */
{
      double cosdec, sindec, coslat, sinlat;
      double p[3], dp[3], x, y, z, D;
      
      /* Don't bother with this unless the equatorial horizontal parallax
       * is at least 0.005"
       */
      if( dist > 1758.8 )
            return(-1);

      DISFAC = au / (0.001 * aearth);
      cosdec = cos(*dec);
      sindec = sin(*dec);

      /* Observer's astronomical latitude
       */
      x = tlat * DTR;
      coslat = cos(x);
      sinlat = sin(x);

      /* Convert to equatorial rectangular coordinates
       * in which unit distance = earth radius
       */
      D = dist * DISFAC;
      p[0] = D*cosdec*cos(*ra);
      p[1] = D*cosdec*sin(*ra);
      p[2] = D*sindec;

      dp[0] = -trho*coslat*cos(last);
      dp[1] = -trho*coslat*sin(last);
      dp[2] = -trho*sinlat;

      x = p[0] + dp[0];
      y = p[1] + dp[1];
      z = p[2] + dp[2];
      D = x*x + y*y + z*z;
      D = sqrt(D);      /* topocentric distance */

      /* recompute ra and dec */
      *ra = zatan2(x,y);
      *dec = asin(z/D);
      showcor((char*) "diurnal parallax", p, dp );
      return(0);
}

double Starpos::refrac(double alt)
//double alt;     /* altitude in degrees */
{
double y, y0, D0, N, D, P, Q;
int i;

if( (alt < -2.0) || (alt >= 90.0) )
      return(0.0);

/* For high altitude angle, AA page B61
 * Accuracy "usually about 0.1' ".
 */
if( alt > 15.0 )
      {
      D = 0.00452*atpress/((273.0+attemp)*tan( DTR*alt ));
      return(D);
      }

/* Formula for low altitude is from the Almanac for Computers.
 * It gives the correction for observed altitude, so has
 * to be inverted numerically to get the observed from the true.
 * Accuracy about 0.2' for -20C < T < +40C and 970mb < P < 1050mb.
 */

/* Start iteration assuming correction = 0
 */
y = alt;
D = 0.0;
/* Invert Almanac for Computers formula numerically
 */
P = (atpress - 80.0)/930.0;
Q = 4.8e-3 * (attemp - 10.0);
y0 = y;
D0 = D;

for( i=0; i<4; i++ )
      {
      N = y + (7.31/(y+4.4));
      N = 1.0/tan(DTR*N);
      D = N*P/(60.0 + Q * (N + 39.0));
      N = y - y0;
      y0 = D - D0 - N; /* denominator of derivative */

      if( (N != 0.0) && (y0 != 0.0) )
/* Newton iteration with numerically estimated derivative */
            N = y - N*(alt + D - y)/y0;
      else
/* Can't do it on first pass */
            N = alt + D;

      y0 = y;
      D0 = D;
      y = N;
      }
return( D );
}

int Starpos::altaz(double pol[],double J )

{
      double dec, cosdec, sindec, lha, coslha, sinlha;
      double ra, dist, last, alt, az, coslat, sinlat;
      double N, D, x, y, z, TPI;

      
      ra = pol[0];
      dec = pol[1];
      dist = pol[2];
      TPI = 2.0*PI;

      /* local apparent sidereal time, seconds converted to radians
       */
      last = sidrlt( J, tlong ) * DTR/240.0;
      lha = last - ra; /* local hour angle, radians */
      if( prtflg )
            {
            //printf( "Local apparent sidereal time " );
            hms( last );
            //printf( "\n" );
            }
      /* Display rate at which ra and dec are changing
       */
      /*
       *if( prtflg )
       *    {
       *    x = RTS/24.0;
       *    N = x*dradt;
       *    D = x*ddecdt;
       *    if( N != 0.0 )
       *          printf( "dRA/dt %.2f\"/h, dDec/dt %.2f\"/h\n", N, D );
       *    }
       */

      diurab( last, &ra, &dec );
      /* Do rise, set, and transit times
         trnsit.c takes diurnal parallax into account,
         but not diurnal aberration.  */
      lha = last - ra;
      trnsit( J, lha, dec );

      /* Diurnal parallax
       */
      diurpx( last, &ra, &dec, dist );

      /* Diurnal aberration
       */
      /*diurab( last, &ra, &dec );*/

      /* Convert ra and dec to altitude and azimuth
       */
      cosdec = cos(dec);
      sindec = sin(dec);
      lha = last - ra;
      coslha = cos(lha);
      sinlha = sin(lha);

      /* Use the geodetic latitude for altitude and azimuth */
      x = DTR * glat;
      coslat = cos(x);
      sinlat = sin(x);

      N = -cosdec * sinlha;
      D =  sindec * coslat  -  cosdec * coslha * sinlat;
      az = RTD * zatan2( D, N );
      alt = sindec * sinlat  +  cosdec * coslha * coslat;   
      alt = RTD * asin(alt);


      /* Correction for atmospheric refraction
       * unit = degrees
       */
      D = refrac( alt );
      alt += D;

      /* Convert back to R.A. and Dec.
       */
      y = sin(DTR*alt);
      x = cos(DTR*alt);
      z = cos(DTR*az);
      sinlha = -x * sin(DTR*az);
      coslha = y*coslat - x*z*sinlat;
      sindec = y*sinlat + x*z*coslat;
      lha = zatan2( coslha, sinlha );


      y = ra; /* save previous values, before refrac() */
      z = dec;
      dec = asin( sindec );
      ra = last - lha;
      //dec = dec + 0.000029088;//aggiunti da me
      //ra = ra + 0.000203616;//aggiunti da me
      y = ra - y; /* change in ra */
      while( y < -PI )
            y += TPI;
      while( y > PI )
            y -= TPI;
      y = RTS*y/15.0;
      z = RTS*(dec - z);
      if( prtflg )
            {
            printf( "atmospheric refraction %.3f deg  dRA %.3fs dDec %.2f\"\n", D, y, z );
            printf( "Topocentric:  Altitude %.3f deg, ", alt );
            printf( "Azimuth %.3f deg\n", az );

            printf( "Topocentric: R.A." );
            hms( ra );
            printf( " Dec." );
            dms( dec );
            printf( "\n" );
            }
      return(0);
}

/* Iterative computation of rise, transit, and set times.  */

int Starpos::iter_trnsit( int (* func)() )
{
  double JDsave, TDTsave, UTsave;
  double date, date_save, date_trnsit, t0, t1;
  double rise1, set1, trnsit1, loopctr, retry;
  double TPI;
  int prtsave;

  loopctr = 0;
  prtsave = prtflg;
  TPI = PI + PI;
  JDsave = JD;
  TDTsave = TDT;
  UTsave = UT;
  date = floor(UT - 0.5) + 0.5;
  retry = 0;
  date_save = date;
  t1 = date_save;

  /* Find transit time. */
  do
    {
      date = date_save;
      t0 = t1;
      iter_func(t0, func);
      t1 = date + r_trnsit / TPI;
      if (++loopctr > 10)
      {
        printf ("? Transit time did not converge.\n");
        goto no_trnsit;
      }

  /* Reject transit if it is more than half a day from entered date.  */
  if (retry == 0)
    {
      if ((UTsave - t1) > 0.5)
      {
        date_save += 1.0;
        t1 = t_trnsit + 1.0;
        retry = 1;
        /*    goto do_retry; */
      }
      if ((UTsave - t1) < -0.5)
      {
        date_save -= 1.0;
        t1 = t_trnsit - 1.0;
        retry = 1;
        /*    goto do_retry; */
      }
    }
    }
  while (fabs(t1 - t0) > .0001);

  t_trnsit = t1;
  trnsit1 = r_trnsit;
  set1 = r_set;


  if (f_trnsit == 0)
      goto prtrnsit;

  /* Set current date to be that of the transit just found.  */
  date_trnsit = date;
  t1 = date + r_rise / TPI;
  /* Choose rising no later than transit.  */
  if (t1 >= t_trnsit)
    {
      date -= 1.0;
      t1 = date + r_rise / TPI;
    }
  loopctr = 0;
  do
    {
      t0 = t1;
      iter_func(t0, func);
      /* Skip out if no event found.  Don't report rise or set.  */
      if ((f_trnsit == 0) || (++loopctr > 10))
      {
        if (loopctr > 10)
          {
            printf ("? Rise time did not converge.\n");
          }
        f_trnsit = 0;
        goto prtrnsit;
      }
      t1 = date + r_rise / TPI;
      if (t1 > t_trnsit)
      {
        date -= 1;
        t1 = date + r_rise / TPI;
      }
    }
  while (fabs(t1 - t0) > .0001);
  rise1 = r_rise;
  t_rise = t1;

  /* Set current date to be that of the transit.  */
  date = date_trnsit;
  r_set = set1;
  /* Choose setting no earlier than transit.  */
  t1 = date + r_set / TPI;
  if (t1 <= t_trnsit)
    {
      date += 1.0;
      t1 = date + r_set / TPI;
    }
  loopctr = 0;
  do
    {
      t0 = t1;
      iter_func(t0, func);
      if ((f_trnsit == 0) || (++loopctr > 10))
      {
        if (loopctr > 10)
          {
            printf ("? Set time did not converge.\n");
          }
        f_trnsit = 0;
        goto prtrnsit;
      }
      t1 = date + r_set / TPI;
      if (t1 < t_trnsit)
      {
        date += 1.0;
        t1 = date + r_set / TPI;
      }
    }
  while (fabs(t1 - t0) > .0001);

  t_set = t1;
  r_trnsit = trnsit1;
  r_rise = rise1;

prtrnsit:
  prtflg = 1;
  printf( "local meridian transit " );
  UT = t_trnsit;
  jtocal (t_trnsit);
  if (f_trnsit != 0)
    {
      printf( "rises " );
      UT = t_rise;
      jtocal (t_rise);
      printf( "sets " );
      UT = t_set;
      jtocal (t_set);
      t0 = t_set - t_rise;
      if ((t0 > 0.0) && (t0 < 1.0))
      printf("Visible hours %.4f\n", 24.0 * t0);
      if ((fabs(JDsave - t_rise) > 0.5) && (fabs(JDsave - t_trnsit) > 0.5)
        && (fabs(JDsave - t_set) > 0.5))
      printf("? wrong event date.\n");
    }
no_trnsit:
  JD = JDsave;
  TDT = TDTsave;
  UT = UTsave;
  /* Reset to original input date entry.  */
  prtflg = 0;
  update();
  prtflg = prtsave;
  return 0;
}

void Starpos:: mean_elements (double J)
{
  double x, T, T2;

  /* Time variables.  T is in Julian centuries.  */
  T = (J - 2451545.0) / 36525.0;
  T2 = T*T;

  /* Mean longitudes of planets (Simon et al, 1994)
     .047" subtracted from constant term for offset to DE403 origin. */

  /* Mercury */
  x = mods3600( 538101628.6889819 * T + 908103.213 );
  x += (6.39e-6 * T
       - 0.0192789) * T2;
  Args[0] = STR * x;

  /* Venus */
  x = mods3600( 210664136.4335482 * T + 655127.236 );
  x += (-6.27e-6  * T
       + 0.0059381) * T2;
  Args[1] = STR * x;

  /* Earth  */
  x = mods3600( 129597742.283429 * T + 361679.198 );
  x += (-5.23e-6 * T
       - 2.04411e-2 ) * T2;
  Ea_arcsec = x;
  Args[2] = STR * x;

  /* Mars */
  x = mods3600(  68905077.493988 * T +  1279558.751 );
  x += (-1.043e-5 * T
       + 0.0094264) * T2;
  Args[3] = STR * x;

  /* Jupiter */
  x = mods3600( 10925660.377991 * T + 123665.420 );
  x += ((((-3.4e-10 * T
          + 5.91e-8) * T
         + 4.667e-6) * T
        + 5.706e-5) * T
         - 3.060378e-1)*T2;
  Args[4] = STR * x;

   /* Saturn */
  x = mods3600( 4399609.855372 * T + 180278.752 );
  x += (((( 8.3e-10 * T
        - 1.452e-7) * T
        - 1.1484e-5) * T
         - 1.6618e-4) * T
       + 7.561614E-1)*T2;
  Args[5] = STR * x;

  /* Uranus */
  x = mods3600( 1542481.193933 * T + 1130597.971 )
       + (0.00002156*T - 0.0175083)*T2;
  Args[6] = STR * x;

  /* Neptune */
  x = mods3600( 786550.320744 * T + 1095655.149 )
       + (-0.00000895*T + 0.0021103)*T2;
  Args[7] = STR * x;

  /* Copied from cmoon.c, DE404 version.  */
  /* Mean elongation of moon = D */
  x = mods3600( 1.6029616009939659e+09 * T + 1.0722612202445078e+06 );
  x += (((((-3.207663637426e-013 * T
             + 2.555243317839e-011) * T
            + 2.560078201452e-009) * T
           - 3.702060118571e-005) * T
            + 6.9492746836058421e-03) * T /* D, t^3 */
           - 6.7352202374457519e+00) * T2; /* D, t^2 */
  Args[9] = STR * x;

  /* Mean distance of moon from its ascending node = F */
  x = mods3600( 1.7395272628437717e+09 * T + 3.3577951412884740e+05 );
  x += ((((( 4.474984866301e-013 * T
             + 4.189032191814e-011) * T
             - 2.790392351314e-009) * T
            - 2.165750777942e-006) * T
            - 7.5311878482337989e-04) * T /* F, t^3 */
           - 1.3117809789650071e+01) * T2; /* F, t^2 */
  NF_arcsec = x;
  Args[10] = STR * x;

/* Mean anomaly of sun = l' (J. Laskar) */
  x = mods3600(1.2959658102304320e+08 * T + 1.2871027407441526e+06);
  x += ((((((((
             1.62e-20 * T
             - 1.0390e-17 ) * T
            - 3.83508e-15 ) * T
           + 4.237343e-13 ) * T
          + 8.8555011e-11 ) * T
         - 4.77258489e-8 ) * T
        - 1.1297037031e-5 ) * T
       + 8.7473717367324703e-05) * T
      - 5.5281306421783094e-01) * T2;
  Args[11] = STR * x;

  /* Mean anomaly of moon = l */
  x = mods3600( 1.7179159228846793e+09 * T + 4.8586817465825332e+05 );
  x += (((((-1.755312760154e-012 * T
            + 3.452144225877e-011) * T
            - 2.506365935364e-008) * T
             - 2.536291235258e-004) * T
              + 5.2099641302735818e-02) * T /* l, t^3 */
             + 3.1501359071894147e+01) * T2; /* l, t^2 */
  Args[12] = STR * x;

  /* Mean longitude of moon, re mean ecliptic and equinox of date = L  */
  x = mods3600( 1.7325643720442266e+09 * T + 7.8593980921052420e+05);
  x += ((((( 7.200592540556e-014 * T
           + 2.235210987108e-010) * T
          - 1.024222633731e-008) * T
         - 6.073960534117e-005) * T
       + 6.9017248528380490e-03) * T /* L, t^3 */
      - 5.6550460027471399e+00) * T2; /* L, t^2 */
  LP_equinox = x;
  Args[13] = STR * x;

  /* Precession of the equinox  */
 x = ((((((((( -8.66e-20*T - 4.759e-17)*T
           + 2.424e-15)*T
           + 1.3095e-12)*T
           + 1.7451e-10)*T
           - 1.8055e-8)*T
           - 0.0000235316)*T
           + 0.000076)*T
           + 1.105414)*T
           + 5028.791959)*T;
  /* Moon's longitude re fixed J2000 equinox.  */
 /*
   Args[13] -= x;
 */
   pA_precession = STR * x;

  /* Free librations.  */
  /* longitudinal libration 2.891725 years */
  x = mods3600( 4.48175409e7 * T + 8.060457e5 );
  Args[14] = STR * x;
  /* libration P, 24.2 years */
  x = mods3600(  5.36486787e6 * T - 391702.8 );
  Args[15] = STR * x;

#if 0
  Args[16] = 0.0;
#endif
  /* libration W, 74.7 years. */
 x = mods3600( 1.73573e6 * T );
 Args[17] = STR * x;
}


int Starpos::sscc (int k,double arg,int n)
{
  double cu, su, cv, sv, s;
  int i;

  su = sin(arg);
  cu = cos(arg);
  ss[k][0] = su;        /* sin(L) */
  cc[k][0] = cu;        /* cos(L) */
  sv = 2.0 * su * cu;
  cv = cu * cu - su * su;
  ss[k][1] = sv;        /* sin(2L) */
  cc[k][1] = cv;
  for (i = 2; i < n; i++)
    {
      s = su * cv + cu * sv;
      cv = cu * cv - su * sv;
      sv = s;
      ss[k][i] = sv;          /* sin( i+1 L ) */
      cc[k][i] = cv;
    }
  return (0);
}

int Starpos::nutlo(double J)

{
      double f, g, T, T2, T10;
      double MM, MS, FF, DD, OM;
      double cu, su, cv, sv, sw;
      double C, D;
      int i, j, k, k1, m;
      short *p;

      if( jdnut == J )
            return(0);
      jdnut = J;

      /* Julian centuries from 2000 January 1.5,
       * barycentric dynamical time
       */
      T = (J-2451545.0)/36525.0;
      T2 = T * T;
      T10 = T / 10.0;

      /* Fundamental arguments in the FK5 reference system.  */

      /* longitude of the mean ascending node of the lunar orbit
       * on the ecliptic, measured from the mean equinox of date
       */
      OM = (mod3600 (-6962890.539 * T + 450160.280) + (0.008 * T + 7.455) * T2)
            * STR;

      /* mean longitude of the Sun minus the
       * mean longitude of the Sun's perigee
       */
      MS = (mod3600 (129596581.224 * T + 1287099.804) - (0.012 * T + 0.577) * T2)
            * STR;

      /* mean longitude of the Moon minus the
       * mean longitude of the Moon's perigee
       */
      MM = (mod3600 (1717915922.633 * T + 485866.733) + (0.064 * T + 31.310) * T2)
            * STR;

      /* mean longitude of the Moon minus the
       * mean longitude of the Moon's node
       */
      FF = (mod3600 (1739527263.137 * T + 335778.877) + (0.011 * T - 13.257) * T2)
            * STR;

      /* mean elongation of the Moon from the Sun.
       */
      DD = (mod3600 (1602961601.328 * T + 1072261.307) + (0.019 * T - 6.891) * T2)
            * STR;

      /* Calculate sin( i*MM ), etc. for needed multiple angles
       */
      sscc( 0, MM, 3 );
      sscc( 1, MS, 2 );
      sscc( 2, FF, 4 );
      sscc( 3, DD, 4 );
      sscc( 4, OM, 2 );

      C = 0.0;
      D = 0.0;
      p = &nt[0]; /* point to start of table */

      for( i=0; i<105; i++ )
            {
      /* argument of sine and cosine */
            k1 = 0;
            cv = 0.0;
            sv = 0.0;
            for( m=0; m<5; m++ )
                  {
                  j = *p++;
                  if( j )
                        {
                        k = j;
                        if( j < 0 )
                              k = -k;
                        su = ss[m][k-1]; /* sin(k*angle) */
                        if( j < 0 )
                              su = -su;
                        cu = cc[m][k-1];
                        if( k1 == 0 )
                              { /* set first angle */
                              sv = su;
                              cv = cu;
                              k1 = 1;
                              }
                        else
                              { /* combine angles */
                              sw = su*cv + cu*sv;
                              cv = cu*cv - su*sv;
                              sv = sw;
                              }
                        }
                  }
      /* longitude coefficient */
            f  = *p++;
            if( (k = *p++) != 0 )
                  f += T10 * k;

      /* obliquity coefficient */
            g = *p++;
            if( (k = *p++) != 0 )
                  g += T10 * k;

      /* accumulate the terms */
            C += f * sv;
            D += g * cv;
            }
      /* first terms, not in table: */
      C += (-1742.*T10 - 171996.)*ss[4][0];     /* sin(OM) */
      D += (   89.*T10 +  92025.)*cc[4][0];     /* cos(OM) */
      /*
      printf( "nutation: in longitude %.3f\", in obliquity %.3f\"\n", C, D );
      */
      /* Save answers, expressed in radians */
      nutl = 0.0001 * STR * C;
      nuto = 0.0001 * STR * D;
      return(0);
}


int Starpos::nutate(double J,double p[] )
{
      double ce, se, cl, sl, sino, f;
      double dp[3], p1[3];
      int i;

      nutlo(J); /* be sure we calculated nutl and nuto */
      epsiln(J); /* and also the obliquity of date */

      f = eps + nuto;
      ce = cos( f );
      se = sin( f );
      sino = sin(nuto);
      cl = cos( nutl );
      sl = sin( nutl );

      /* Apply adjustment
       * to equatorial rectangular coordinates of object.
       *
       * This is a composite of three rotations: rotate about x axis
       * to ecliptic of date; rotate about new z axis by the nutation
       * in longitude; rotate about new x axis back to equator of date
       * plus nutation in obliquity.
       */
      p1[0] =   cl*p[0]
                  - sl*coseps*p[1]
                  - sl*sineps*p[2];

      p1[1] =   sl*ce*p[0]
                  + ( cl*coseps*ce + sineps*se )*p[1]
                  - ( sino + (1.0-cl)*sineps*ce )*p[2];

      p1[2] =   sl*se*p[0]
                  + ( sino + (cl-1.0)*se*coseps )*p[1]
                  + ( cl*sineps*se + coseps*ce )*p[2];

      for( i=0; i<3; i++ )
            dp[i] = p1[i] - p[i];
      showcor((char*) "nutation", p, dp );

      for( i=0; i<3; i++ )
            p[i] = p1[i];
      return(0);
}


int Starpos::g3plan (double J, plantbl *plan, double pobj[],int objnum)
    
{

        int i, j, k, m, n, k1, ip, np, nt;
      #ifdef _MSC_VER
        char FAR *p;
        long FAR *pl;
        long FAR *pb;
        long FAR *pr;
      #else
        /* On some systems such as Silicon Graphics, char is unsigned
             by default.  */
      #ifdef vax
        /* VAX CC rejects "signed char."  */
        char *p;
      #else
      #ifdef __STDC__
         char *p;
      #else
        char *p;
      #endif
      #endif
        long *pl, *pb, *pr;
      #endif
        double su, cu, sv, cv;
        double T, t, sl, sb, sr;

        mean_elements (J);
      #if 0
        /* For librations, moon's longitude is sidereal.  */
        if (flag)
            Args[13] -= pA_precession;
      #endif
        T = (J - J2000) / plan->timescale;
        n = plan->maxargs;
        /* Calculate sin( i*MM ), etc. for needed multiple angles.  */
        for (i = 0; i < n; i++)
            {
              if ((j = plan->max_harmonic[i]) > 0)
            {
              sscc (i, Args[i], j);
            }
            }

        /* Point to start of table of arguments. */

        p = plan->arg_tbl;
        /* Point to tabulated cosine and sine amplitudes.  */
      #ifdef _MSC_VER
        pl = (long *) plan->lon_tbl;
        pb = (long *) plan->lat_tbl;
        pr = (long *) plan->rad_tbl;
      #else
        pl = (long *) plan->lon_tbl;
        pb = (long *) plan->lat_tbl;
        pr = (long *) plan->rad_tbl;
      #endif
        sl = 0.0;
        sb = 0.0;
        sr = 0.0;

        for (;;)
            {
              /* argument of sine and cosine */
              /* Number of periodic arguments. */
              np = *p++;
              if (np < 0)
            break;
              if (np == 0)
            {                 /* It is a polynomial term.  */
              nt = *p++;
              /* "Longitude" polynomial (phi). */
              cu = *pl++;
              for (ip = 0; ip < nt; ip++)
                  {
                    cu = cu * T + *pl++;
                  }
              /*    sl +=  mods3600 (cu); */
              sl += cu;
              /* "Latitude" polynomial (theta). */
              cu = *pb++;
              for (ip = 0; ip < nt; ip++)
                  {
                    cu = cu * T + *pb++;
                  }
              sb += cu;
              /* Radius polynomial (psi). */
              cu = *pr++;
              for (ip = 0; ip < nt; ip++)
                  {
                    cu = cu * T + *pr++;
                  }
              sr += cu;
              continue;
            }
              k1 = 0;
              cv = 0.0;
              sv = 0.0;
              for (ip = 0; ip < np; ip++)
            {
              /* What harmonic.  */
              j = *p++;
              /* Which planet.  */
              m = *p++ - 1;
              if (j)
                  {
      /*          k = abs (j); */
                  if (j < 0)
                        k = -j;
                  else
                        k = j;
                    k -= 1;
                    su = ss[m][k];  /* sin(k*angle) */
                    if (j < 0)
                  su = -su;
                    cu = cc[m][k];
                    if (k1 == 0)
                  {           /* set first angle */
                    sv = su;
                    cv = cu;
                    k1 = 1;
                  }
                    else
                  {           /* combine angles */
                    t = su * cv + cu * sv;
                    cv = cu * cv - su * sv;
                    sv = t;
                  }
                  }
            }
              /* Highest power of T.  */
              nt = *p++;
              /* Longitude. */
              cu = *pl++;
              su = *pl++;
              for (ip = 0; ip < nt; ip++)
            {
              cu = cu * T + *pl++;
              su = su * T + *pl++;
            }
              sl += cu * cv + su * sv;
              /* Latitiude. */
              cu = *pb++;
              su = *pb++;
              for (ip = 0; ip < nt; ip++)
            {
              cu = cu * T + *pb++;
              su = su * T + *pb++;
            }
              sb += cu * cv + su * sv;
              /* Radius. */
              cu = *pr++;
              su = *pr++;
              for (ip = 0; ip < nt; ip++)
            {
              cu = cu * T + *pr++;
              su = su * T + *pr++;
            }
              sr += cu * cv + su * sv;
            }
        t = plan->trunclvl;
        pobj[0] = Args[objnum - 1] + STR * t * sl;
        pobj[1] = STR * t * sb;
        pobj[2] = plan->distance * (1.0 + STR * t * sr);
        return (0);
}


int Starpos:: gplan (double J, plantbl *plan,double pobj[])
{
  register double su, cu, sv, cv, T;
  double t, sl, sb, sr;
  int i, j, k, m, n, k1, ip, np, nt;
#ifdef _MSC_VER
  char FAR *p;
  double FAR *pl;
  double FAR *pb;
  double FAR *pr;
#else
  /* On some systems such as Silicon Graphics, char is unsigned
     by default.  */
#ifdef vax
  /* VAX CC rejects "signed char."  */
  char *p;
#else
#ifdef __STDC__
   char *p;
#else
  char *p;
#endif
#endif
  double *pl, *pb, *pr;
#endif

  T = (J - J2000) / plan->timescale;
  n = plan->maxargs;
  /* Calculate sin( i*MM ), etc. for needed multiple angles.  */
  for (i = 0; i < n; i++)
    {
      if ((j = plan->max_harmonic[i]) > 0)
      {
        sr = (mods3600 (freqs[i] * T) + phases[i]) * STR;
        sscc (i, sr, j);
      }
    }

  /* Point to start of table of arguments. */
  p = plan->arg_tbl;
  /* Point to tabulated cosine and sine amplitudes.  */
#ifdef _MSC_VER
  pl = (double FAR *) plan->lon_tbl;
  pb = (double FAR *) plan->lat_tbl;
  pr = (double FAR *) plan->rad_tbl;
#else
  pl = (double *) plan->lon_tbl;
  pb = (double *) plan->lat_tbl;
  pr = (double *) plan->rad_tbl;
#endif
  sl = 0.0;
  sb = 0.0;
  sr = 0.0;

  for (;;)
    {
      /* argument of sine and cosine */
      /* Number of periodic arguments. */
      np = *p++;
      if (np < 0)
      break;
      if (np == 0)
      {                 /* It is a polynomial term.  */
        nt = *p++;
        /* Longitude polynomial. */
        cu = *pl++;
        for (ip = 0; ip < nt; ip++)
          {
            cu = cu * T + *pl++;
          }
        sl +=  mods3600 (cu);
        /* Latitude polynomial. */
        cu = *pb++;
        for (ip = 0; ip < nt; ip++)
          {
            cu = cu * T + *pb++;
          }
        sb += cu;
        /* Radius polynomial. */
        cu = *pr++;
        for (ip = 0; ip < nt; ip++)
          {
            cu = cu * T + *pr++;
          }
        sr += cu;
        continue;
      }
      k1 = 0;
      cv = 0.0;
      sv = 0.0;
      for (ip = 0; ip < np; ip++)
      {
        /* What harmonic.  */
        j = *p++;
        /* Which planet.  */
        m = *p++ - 1;
        if (j)
          {
            k = j;
            if (j < 0)
            k = -k;
            k -= 1;
            su = ss[m][k];    /* sin(k*angle) */
            if (j < 0)
            su = -su;
            cu = cc[m][k];
            if (k1 == 0)
            {           /* set first angle */
              sv = su;
              cv = cu;
              k1 = 1;
            }
            else
            {           /* combine angles */
              t = su * cv + cu * sv;
              cv = cu * cv - su * sv;
              sv = t;
            }
          }
      }
      /* Highest power of T.  */
      nt = *p++;
      /* Longitude. */
      cu = *pl++;
      su = *pl++;
      for (ip = 0; ip < nt; ip++)
      {
        cu = cu * T + *pl++;
        su = su * T + *pl++;
      }
      sl += cu * cv + su * sv;
      /* Latitiude. */
      cu = *pb++;
      su = *pb++;
      for (ip = 0; ip < nt; ip++)
      {
        cu = cu * T + *pb++;
        su = su * T + *pb++;
      }
      sb += cu * cv + su * sv;
      /* Radius. */
      cu = *pr++;
      su = *pr++;
      for (ip = 0; ip < nt; ip++)
      {
        cu = cu * T + *pr++;
        su = su * T + *pr++;
      }
      sr += cu * cv + su * sv;
    }
  pobj[0] = STR * sl;
  pobj[1] = STR * sb;
  pobj[2] = STR * plan->distance * sr + plan->distance;
  return (0);
}

double Starpos::g1plan (double J, plantbl *plan)
{
  int i, j, k, m, k1, ip, np, nt;
#ifdef _MSC_VER
  char FAR *p;
  long FAR *pl;
#else
  /* On some systems such as Silicon Graphics, char is unsigned
     by default.  */
#ifdef vax
  /* VAX CC rejects "signed char."  */
  char *p;
#else
#ifdef __STDC__
   char *p;
#else
  char *p;
#endif
#endif
  long *pl;
#endif
  double su, cu, sv, cv;
  double T, t, sl;

  T = (J - J2000) / plan->timescale;
  mean_elements (J);
  /* Calculate sin( i*MM ), etc. for needed multiple angles.  */
  for (i = 0; i < NARGS; i++)
    {
      if ((j = plan->max_harmonic[i]) > 0)
      {
        sscc (i, Args[i], j);
      }
    }

  /* Point to start of table of arguments. */
  p = plan->arg_tbl;
  /* Point to tabulated cosine and sine amplitudes.  */
#ifdef _MSC_VER
  pl = (long FAR *) plan->lon_tbl;
#else
  pl = (long *) plan->lon_tbl;
#endif
  sl = 0.0;

  for (;;)
    {
      /* argument of sine and cosine */
      /* Number of periodic arguments. */
      np = *p++;
      if (np < 0)
      break;
      if (np == 0)
      {                 /* It is a polynomial term.  */
        nt = *p++;
        cu = *pl++;
        for (ip = 0; ip < nt; ip++)
          {
            cu = cu * T + *pl++;
          }
        /*    sl +=  mods3600 (cu); */
        sl += cu;
        continue;
      }
      k1 = 0;
      cv = 0.0;
      sv = 0.0;
      for (ip = 0; ip < np; ip++)
      {
        /* What harmonic.  */
        j = *p++;
        /* Which planet.  */
        m = *p++ - 1;
        if (j)
          {
/*          k = abs (j); */
            if (j < 0)
                  k = -j;
            else
                  k = j;
            k -= 1;
            su = ss[m][k];    /* sin(k*angle) */
            if (j < 0)
            su = -su;
            cu = cc[m][k];
            if (k1 == 0)
            {           /* set first angle */
              sv = su;
              cv = cu;
              k1 = 1;
            }
            else
            {           /* combine angles */
              t = su * cv + cu * sv;
              cv = cu * cv - su * sv;
              sv = t;
            }
          }
      }
      /* Highest power of T.  */
      nt = *p++;
      /* Cosine and sine coefficients.  */
      cu = *pl++;
      su = *pl++;
      for (ip = 0; ip < nt; ip++)
      {
        cu = cu * T + *pl++;
        su = su * T + *pl++;
      }
      sl += cu * cv + su * sv;
    }
  return (plan->trunclvl * sl);
}



int Starpos::g2plan (double J, plantbl *plan,double pobj[])
{
  int i, j, k, m, n, k1, ip, np, nt;
#if _MSC_VER
  char FAR *p;
  long FAR *pl;
  long FAR *pr;
#else
  /* On some systems such as Silicon Graphics, char is unsigned
     by default.  */
#ifdef vax
  /* VAX CC rejects "signed char."  */
  char *p;
#else
#ifdef __STDC__
   char *p;
#else
  char *p;
#endif
#endif
  long *pl, *pr;
#endif
  double su, cu, sv, cv;
  double T, t, sl, sr;

  mean_elements (J);
#if 0
  /* For librations, moon's longitude is sidereal.  */
  if (flag)
    Args[13] -= pA_precession;
#endif
  T = (J - J2000) / plan->timescale;
  n = plan->maxargs;
  /* Calculate sin( i*MM ), etc. for needed multiple angles.  */
  for (i = 0; i < n; i++)
    {
      if ((j = plan->max_harmonic[i]) > 0)
      {
        sscc (i, Args[i], j);
      }
    }

  /* Point to start of table of arguments. */
  p = plan->arg_tbl;
  /* Point to tabulated cosine and sine amplitudes.  */
#ifdef _MSC_VER
  pl = (long FAR *) plan->lon_tbl;
  pr = (long FAR *) plan->rad_tbl;
#else
  pl = (long *) plan->lon_tbl;
  pr = (long *) plan->rad_tbl;
#endif
  sl = 0.0;
  sr = 0.0;

  for (;;)
    {
      /* argument of sine and cosine */
      /* Number of periodic arguments. */
      np = *p++;
      if (np < 0)
      break;
      if (np == 0)
      {                 /* It is a polynomial term.  */
        nt = *p++;
        /* Longitude polynomial. */
        cu = *pl++;
        for (ip = 0; ip < nt; ip++)
          {
            cu = cu * T + *pl++;
          }
        /*    sl +=  mods3600 (cu); */
        sl += cu;
        /* Radius polynomial. */
        cu = *pr++;
        for (ip = 0; ip < nt; ip++)
          {
            cu = cu * T + *pr++;
          }
        sr += cu;
        continue;
      }
      k1 = 0;
      cv = 0.0;
      sv = 0.0;
      for (ip = 0; ip < np; ip++)
      {
        /* What harmonic.  */
        j = *p++;
        /* Which planet.  */
        m = *p++ - 1;
        if (j)
          {
/*          k = abs (j); */
            if (j < 0)
                  k = -j;
            else
                  k = j;
            k -= 1;
            su = ss[m][k];    /* sin(k*angle) */
            if (j < 0)
            su = -su;
            cu = cc[m][k];
            if (k1 == 0)
            {           /* set first angle */
              sv = su;
              cv = cu;
              k1 = 1;
            }
            else
            {           /* combine angles */
              t = su * cv + cu * sv;
              cv = cu * cv - su * sv;
              sv = t;
            }
          }
      }
      /* Highest power of T.  */
      nt = *p++;
      /* Longitude. */
      cu = *pl++;
      su = *pl++;
      for (ip = 0; ip < nt; ip++)
      {
        cu = cu * T + *pl++;
        su = su * T + *pl++;
      }
      sl += cu * cv + su * sv;
      /* Radius. */
      cu = *pr++;
      su = *pr++;
      for (ip = 0; ip < nt; ip++)
      {
        cu = cu * T + *pr++;
        su = su * T + *pr++;
      }
      sr += cu * cv + su * sv;
    }
  t = plan->trunclvl;
  pobj[0] = t * sl;
  pobj[2] = t * sr;
  return (0);
}


int Starpos::gmoon (double J,double rect[],double pol[])
{
  double x, cosB, sinB, cosL, sinL;

  g2plan (J, &moonlr, pol);
  x = pol[0];
  x += LP_equinox;
  if (x < -6.48e5)
    x += 1.296e6;
  if (x > 6.48e5)
    x -= 1.296e6;
  pol[0] = STR * x;
  x = g1plan (J, &moonlat);
  pol[1] = STR * x;
  x = (1.0 + STR * pol[2]) * moonlr.distance;
  pol[2] = x;
  /* Convert ecliptic polar to equatorial rectangular coordinates.  */
  epsiln(J);
  cosB = cos(pol[1]);
  sinB = sin(pol[1]);
  cosL = cos(pol[0]);
  sinL = sin(pol[0]);
  rect[0] = cosB * cosL * x;
  rect[1] = (coseps * cosB * sinL - sineps * sinB) * x;
  rect[2] = (sineps * cosB * sinL + coseps * sinB) * x;
  return 0;
}

int Starpos::embofs(double J,double ea[],double * pr )
{
      double pm[3], polm[3];
      double a, b;
      int i;

      /* Compute the vector Moon - Earth.  */
      gmoon( J, pm, polm );

      /* Precess the lunar position
       * to ecliptic and equinox of J2000.0
       */
      precess( pm, J, 1 );

      /* Adjust the coordinates of the Earth
       */
      a = 1.0 / (emrat +  1.0);
      b = 0.0;
      for( i=0; i<3; i++ )
            {
            ea[i] = ea[i] - a * pm[i];
            b = b + ea[i] * ea[i];
            }
      /* Sun-Earth distance.  */
      *pr = sqrt(b);
      return(0);
}

/* Reduce x modulo 2 pi
 */
#define TPI (2.0*PI)
double Starpos:: modtp(double x)
{
      double y;

      y = floor( x/TPI );
      y = x - y * TPI;
      while( y < 0.0 )
            y += TPI;
      while( y >= TPI )
            y -= TPI;
      return(y);
}

/* Reduce x modulo 360 degrees
 */
double Starpos:: mod360(double x)

{
      long k;
      double y;

      k = (long) (x/360.0);
      y = x  -  k * 360.0;
      while( y < 0.0 )
            y += 360.0;
      while( y > 360.0 )
            y -= 360.0;
      return(y);
}



int Starpos::kepler(double J, struct orbit *e,double rect[],double polar[])

{
      double alat, E, M, W, temp;
      double epoch, inclination, ascnode, argperih;
      double meandistance, dailymotion, eccent, meananomaly;
      double r, coso, sino, cosa;

      /* Call program to compute position, if one is supplied.  */
      if( e->ptable )
            {
            if( e == &earth )
                  g3plan (J, e->ptable, polar, 3);
            else
                  gplan (J, e->ptable, polar);
            E = polar[0]; /* longitude */
            e->L = E;
            W = polar[1]; /* latitude */
            r = polar[2]; /* radius */
            e->r = r;
            e->epoch = J;
            e->equinox = J2000;
            goto kepdon;
            }

      /* Decant the parameters from the data structure
       */
      epoch = e->epoch;
      inclination = e->i;
      ascnode = e->W * DTR;
      argperih = e->w;
      meandistance = e->a; /* semimajor axis */
      dailymotion = e->dm;
      eccent = e->ecc;
      meananomaly = e->M;
      /* Check for parabolic orbit. */
      if( eccent == 1.0 )
            {
      /* meandistance = perihelion distance, q
       * epoch = perihelion passage date
       */
            temp = meandistance * sqrt(meandistance);
            W = (J - epoch ) * 0.0364911624 / temp;
      /* The constant above is 3 k / sqrt(2),
       * k = Gaussian gravitational constant = 0.01720209895 . */
            E = 0.0;
            M = 1.0;
            while( fabs(M) > 1.0e-11 )
                  {
                  temp = E * E;
                  temp = (2.0 * E * temp + W)/( 3.0 * (1.0 + temp));
                  M = temp - E;
                  if( temp != 0.0 )
                        M /= temp;
                  E = temp;
                  }
            r = meandistance * (1.0 + E * E );
            M = atan( E );
            M = 2.0 * M;
            alat = M + DTR*argperih;
            goto parabcon;
            }
      if( eccent > 1.0 )
            {
      /* The equation of the hyperbola in polar coordinates r, theta
       * is r = a(e^2 - 1)/(1 + e cos(theta))
       * so the perihelion distance q = a(e-1),
       * the "mean distance"  a = q/(e-1).
       */
            meandistance = meandistance/(eccent - 1.0);
            temp = meandistance * sqrt(meandistance);
            W = (J - epoch ) * 0.01720209895 / temp;
      /* solve M = -E + e sinh E */
            E = W/(eccent - 1.0);
            M = 1.0;
            while( fabs(M) > 1.0e-11 )
                  {
                  M = -E + eccent * sinh(E) - W;
                  E += M/(1.0 - eccent * cosh(E));
                  }
            r = meandistance * (-1.0 + eccent * cosh(E));
            temp = (eccent + 1.0)/(eccent - 1.0);
            M = sqrt(temp) * tanh( 0.5*E );
            M = 2.0 * atan(M);
            alat = M + DTR*argperih;
            goto parabcon;
            }
      /* Calculate the daily motion, if it is not given.
       */
      if( dailymotion == 0.0 )
            {
      /* The constant is 180 k / pi, k = Gaussian gravitational constant.
       * Assumes object in heliocentric orbit is massless.
       */
            dailymotion = 0.9856076686/(e->a*sqrt(e->a));
            }
      dailymotion *= J - epoch;
      /* M is proportional to the area swept out by the radius
       * vector of a circular orbit during the time between
       * perihelion passage and Julian date J.
       * It is the mean anomaly at time J.
       */
      M = DTR*( meananomaly + dailymotion );
      M = modtp(M);
      /* If mean longitude was calculated, adjust it also
       * for motion since epoch of elements.
       */
      if( e->L )
            {
            e->L += dailymotion;
            e->L = mod360( e->L );
            }

      /* By Kepler's second law, M must be equal to
       * the area swept out in the same time by an
       * elliptical orbit of same total area.
       * Integrate the ellipse expressed in polar coordinates
       *     r = a(1-e^2)/(1 + e cosW)
       * with respect to the angle W to get an expression for the
       * area swept out by the radius vector.  The area is given
       * by the mean anomaly; the angle is solved numerically.
       *
       * The answer is obtained in two steps.  We first solve
       * Kepler's equation
       *    M = E - eccent*sin(E)
       * for the eccentric anomaly E.  Then there is a
       * closed form solution for W in terms of E.
       */

      E = M; /* Initial guess is same as circular orbit. */
      temp = 1.0;
      do
            {
      /* The approximate area swept out in the ellipse */
            temp = E - eccent * sin(E)
      /* ...minus the area swept out in the circle */
                  - M;
      /* ...should be zero.  Use the derivative of the error
       * to converge to solution by Newton's method.
       */
            E -= temp/(1.0 - eccent*cos(E));
            }
      while( fabs(temp) > 1.0e-11 );

      /* The exact formula for the area in the ellipse is
       *    2.0*atan(c2*tan(0.5*W)) - c1*eccent*sin(W)/(1+e*cos(W))
       * where
       *    c1 = sqrt( 1.0 - eccent*eccent )
       *    c2 = sqrt( (1.0-eccent)/(1.0+eccent) ).
       * Substituting the following value of W
       * yields the exact solution.
       */
      temp = sqrt( (1.0+eccent)/(1.0-eccent) );
      W = 2.0 * atan( temp * tan(0.5*E) );

      /* The true anomaly.
       */
      W = modtp(W);

      meananomaly *= DTR;
      /* Orbital longitude measured from node
       * (argument of latitude)
       */
      if( e->L )
            alat = (e->L)*DTR + W - meananomaly - ascnode;
      else
            alat = W + DTR*argperih; /* mean longitude not given */

      /* From the equation of the ellipse, get the
       * radius from central focus to the object.
       */
      r = meandistance*(1.0-eccent*eccent)/(1.0+eccent*cos(W));

      parabcon:

      /* The heliocentric ecliptic longitude of the object
       * is given by
       *   tan( longitude - ascnode )  =  cos( inclination ) * tan( alat ).
       */
      coso = cos( alat );
      sino = sin( alat );
      inclination *= DTR;
      W = sino * cos( inclination );
      E = zatan2( coso, W ) + ascnode;

      /* The ecliptic latitude of the object
       */
      W = sino * sin( inclination );
      W = asin(W);

      kepdon:

      /* Convert to rectangular coordinates,
       * using the perturbed latitude.
       */
      rect[2] = r * sin(W);
      cosa = cos(W);
      rect[1] = r * cosa * sin(E);
      rect[0] = r * cosa * cos(E);

      /* Convert from heliocentric ecliptic rectangular
       * to heliocentric equatorial rectangular coordinates
       * by rotating eps radians about the x axis.
       */
      epsiln( e->equinox );
      W = coseps*rect[1] - sineps*rect[2];
      M = sineps*rect[1] + coseps*rect[2];
      rect[1] = W;
      rect[2] = M;

      /* Precess the position
       * to ecliptic and equinox of J2000.0
       * if not already there.
       */
      precess( rect, e->equinox, 1 );

      /* If earth, adjust from earth-moon barycenter to earth
       * by AA page E2.
       */
      if( e == &earth )
            {
            embofs( J, rect, &r ); /* see below */
            }

      /* Rotate back into the ecliptic.  */
      epsiln( J2000 );
      W = coseps*rect[1] + sineps*rect[2];
      M = -sineps*rect[1] + coseps*rect[2];

      /* Convert to polar coordinates */
      E = zatan2( rect[0], W );
      W = asin( M/r );

      /* Output the polar cooordinates
       */
      polar[0] = E; /* longitude */
      polar[1] = W; /* latitude */
      polar[2] = r; /* radius */

      return(0);
}




/*//////////caluculate velociti vector/////////*/
int Starpos::velearth(double J )

{
      double e[3], p[3], t;
      int i;
      #if DEBUG
      double x[3], A, q;
      #endif

      if( J == jvearth )
            return(0);

      jvearth = J;

      /* calculate heliocentric position of the earth
       * as of a short time ago.
       */
      t = 0.005;
      kepler( TDT-t, &earth, e, p );

      for( i=0; i<3; i++ )
            vearth[i] = (rearth[i] - e[i])/t;

      #if DEBUG
      /* Generate display for comparison with Almanac values. */
      for( i=0; i<3; i++ )
            {
            q = vearth[i];
            A += q*q;
            x[i] = q;
            }
      A = sqrt(A);
      precess( x, TDT, 1 );
      printf( "Vearth %.6e, X %.6f, Y %.6f, Z %.6f\n", A, x[0], x[1], x[2] );
      #endif
      return(0);
}

/*//////////////////annual abberration///////////////////////*/
int Starpos:: annuab(double p[] )
//double p[]; /* unit vector pointing from earth to object */
{
      double A, B, C;
      double betai, pV;
      double x[3], V[3];
      int i;

      /* Calculate the velocity of the earth (see vearth.c).
       */
      velearth( TDT );
      betai = 0.0;
      pV = 0.0;
      for( i=0; i<3; i++ )
            {
            A = vearth[i]/Clightaud;
            V[i] = A;
            betai += A*A;
            pV += p[i] * A;
            }
      /* Make the adjustment for aberration.
       */
      betai = sqrt( 1.0 - betai );
      C = 1.0 + pV;
      A = betai/C;
      B = (1.0  +  pV/(1.0 + betai))/C;

      for( i=0; i<3; i++ )
            {
            C = A * p[i]  +  B * V[i];
            x[i] = C;
            dp[i] = C - p[i];
            }

      showcor((char*) "annual aberration", p, dp );
      for( i=0; i<3; i++ )
            p[i] = x[i];
      return(0);
}

/* jd from gregorian calendar */
double Starpos::caltoj(long year,int month,double day )

{
      long y, a, b, c, e, m;
      double J;


      /* The origin should be chosen to be a century year
       * that is also a leap year.  We pick 4801 B.C.
       */
      y = year + 4800;
      if( year < 0 )
            {
            y += 1;
            }

      /* The following magic arithmetic calculates a sequence
       * whose successive terms differ by the correct number of
       * days per calendar month.  It starts at 122 = March; January
       * and February come after December.
       */
      m = month;
      if( m <= 2 )
            {
            m += 12;
            y -= 1;
            }
      e = (306 * (m+1))/10;

      a = y/100;  /* number of centuries */
      if( year <= 1582L )
            {
            if( year == 1582L )
                  {
                  if( month < 10 )
                        goto julius;
                  if( month > 10)
                        goto gregor;
                  if( day >= 15 )
                        goto gregor;
                  }
      julius:
            //printf( " Julian Calendar assumed.\n" );
            b = -38;
            }
      else
            { /* -number of century years that are not leap years */
      gregor:
            b = (a/4) - a;
            }

      c = (36525L * y)/100; /* Julian calendar years and leap years */

      /* Add up these terms, plus offset from J 0 to 1 Jan 4801 B.C.
       * Also fudge for the 122 days from the month algorithm.
       */
      J = b + c + e + day - 32167.5;
      return( J );
}


/*
 * Convert change in rectangular coordinatates to change
 * in right ascension and declination.
 * For changes greater than about 0.1 degree, the
 * coordinates are converted directly to R.A. and Dec.
 * and the results subtracted.  For small changes,
 * the change is calculated to first order by differentiating
 *   tan(R.A.) = y/x
 * to obtain
 *    dR.A./cos**2(R.A.) = dy/x  -  y dx/x**2
 * where
 *    cos**2(R.A.)  =  1/(1 + (y/x)**2).
 *
 * The change in declination arcsin(z/R) is
 *   d asin(u) = du/sqrt(1-u**2)
 *   where u = z/R.
 *
 * p0 is the initial object - earth vector and
 * p1 is the vector after motion or aberration.
 *
 */

int Starpos::deltap(double p0[],double p1[],double* dr,double * dd )

{
      double dp[3], A, B, P, Q, x, y, z;
      int i;

      P = 0.0;
      Q = 0.0;
      z = 0.0;
      for( i=0; i<3; i++ )
            {
            x = p0[i];
            y = p1[i];
            P += x * x;
            Q += y * y;
            y = y - x;
            dp[i] = y;
            z += y*y;
            }

      A = sqrt(P);
      B = sqrt(Q);

      if( (A < 1.e-7) || (B < 1.e-7) || (z/(P+Q)) > 5.e-7 )
            {
            P = zatan2( p0[0], p0[1] );
            Q = zatan2( p1[0], p1[1] );
            Q = Q - P;
            while( Q < -PI )
                  Q += 2.0*PI;
            while( Q > PI )
                  Q -= 2.0*PI;
            *dr = Q;
            P = asin( p0[2]/A );
            Q = asin( p1[2]/B );
            *dd = Q - P;
            return(0);
            }


      x = p0[0];
      y = p0[1];
      if( x == 0.0 )
            {
            *dr = 1.0e38;
            }
      else
            {
            Q = y/x;
            Q = (dp[1]  -  dp[0]*y/x)/(x * (1.0 + Q*Q));
            *dr = Q;
            }

      x = p0[2]/A;
      P = sqrt( 1.0 - x*x );
      *dd = (p1[2]/B - x)/P;
      return(0);
}

/* Display magnitude of correction vector
 * in arc seconds
 */
int Starpos::showcor(char* strng,double p[],double dp[] )

{
      double p1[3], dr, dd;
      int i;

      if( prtflg == 0 )
            return(0);

      for( i=0; i<3; i++ )
            p1[i] = p[i] + dp[i];

      deltap( p, p1, &dr, &dd );
      printf( "%s dRA %.3fs dDec %.2f\"\n", strng, RTS*dr/15.0, RTS*dd );
      return(0);
}


int Starpos::relativity(double p[],double q[],double e[] )
/* unit vector from earth to object:

double p[];
/* heliocentric ecliptic rectangular coordinates
 * of earth and object:

double q[], e[];*/
{
      double C;
      int i;

      C = 1.974e-8/(SE*(1.0+qe));
      for( i=0; i<3; i++ )
            {
            dp[i] = C*(pq*e[i]/SE - ep*q[i]/SO);
            p[i] += dp[i];
            }
      if( prtflg )
            //printf( "elongation from sun %.2f degrees, ", acos( -ep )/DTR );
      showcor((char*) "light defl.", p, dp );
      return(0);
}


/* Radians to degrees, minutes, seconds
 */
int Starpos::dms(double x )

{
      double s;
      int d, m;

      s = x * RTD;
      if(fPlanet)
      {
            Dec = s;
      }

      if( s < 0.0 )
            {
            printf( " -" );
            s = -s;
            }
      else
            printf( "  " );
      d = (int) s;
      s -= d;
      s *= 60;
      m = (int) s;
      s -= m;
      s *= 60;
      printf( "%3dd %02d\' %05.2f\"  ", d, m, s );
      return(0);
}


/* Display Right Ascension and Declination
 * from input equatorial rectangular unit vector.
 * Output vector pol[] contains R.A., Dec., and radius.
 */
int Starpos:: showrd(char* msg,double p[],double pol[] )
{
      double x, y, r;
      int i;

      r = 0.0;
      for( i=0; i<3; i++ )
            {
            x = p[i];
            r += x * x;
            }
      r = sqrt(r);

      x =  zatan2( p[0], p[1] );
      pol[0] = x;

      y = asin( p[2]/r );
      pol[1] = y;
      pol[2] = r;

      if(flag)
      {
            RA = x * RTOH * 15;
            RA = RA - 360;
            RA = fabs(RA);
            Dec = y * RTD;
      }
      if (prtflg != 0)
        {
            printf( "%s  R.A. ", msg );
            hms( x );
            printf( "Dec. " );
            dms( y );
            printf( "\n" );
            constellation = whatconstel (x, y);
        }
      return(0);
}


int Starpos::angles(double p[],double q[],double e[] )
{
      double a, b, s;
      int i;

      EO = 0.0;
      SE = 0.0;
      SO = 0.0;
      pq = 0.0;
      ep = 0.0;
      qe = 0.0;
      for( i=0; i<3; i++ )
            {
            a = e[i];
            b = q[i];
            s = p[i];
            EO += s * s;
            SE += a * a;
            SO += b * b;
            pq += s * b;
            ep += a * s;
            qe += b * a;
            }
      EO = sqrt(EO); /* Distance between Earth and object */
      SO = sqrt(SO); /* Sun - object */
      SE = sqrt(SE); /* Sun - earth */
      /* Avoid fatality: if object equals sun, SO is zero.  */
      if( SO > 1.0e-12 )
            {
            pq /= EO*SO;      /* cosine of sun-object-earth */
            qe /= SO*SE;      /* cosine of earth-sun-object */
            }
      ep /= SE*EO;      /* -cosine of sun-earth-object */
      return(0);
}


int Starpos::hms( double x )

{
      int h, m;
      long sint, sfrac;
      double s;

      s = x * RTOH;
      if(fPlanet == true)
      {
            RA = s * 15;
            RA = RA - 360;
            RA = fabs(RA);
            RA = fmod(RA,360);
      }
      
      if( s < 0.0 )
            s += 24.0;
      h = (int) s;
      s -= h;
      s *= 60;
      m = (int) s;
      s -= m;
      s *= 60;
      /* Handle shillings and pence roundoff. */
      sfrac = (long) (1000.0 * s + 0.5);
      if( sfrac >= 60000L )
        {
            sfrac -= 60000L;
            m += 1;
            if( m >= 60 )
              {
            m -= 60;
            h += 1;
              }
        }
      sint = sfrac / 1000;
      sfrac -= sint * 1000;
      printf( "%3dh %02dm %02ld.%03lds  ", h, m, sint, sfrac );
      return(0);
}


int Starpos::epsiln(double J)
//double J; /* Julian date input */
{
      double T;

      if( J == jdeps )
            return(0);
      T = (J - 2451545.0)/36525.0;

      #if WILLIAMS
      /* DE403 values. */
            T /= 10.0;
            eps = ((((((((( 2.45e-10*T + 5.79e-9)*T + 2.787e-7)*T
                  + 7.12e-7)*T - 3.905e-5)*T - 2.4967e-3)*T
            - 5.138e-3)*T + 1.9989)*T - 0.0175)*T - 468.33960)*T
            + 84381.406173;
      #else
      /* This expansion is from the AA.
       * Note the official 1976 IAU number is 23d 26' 21.448", but
       * the JPL numerical integration found 21.4119".
       */
      #if SIMON
            T /= 10.0;
            eps = ((((((((( 2.45e-10*T + 5.79e-9)*T + 2.787e-7)*T
                  + 7.12e-7)*T - 3.905e-5)*T - 2.4967e-3)*T
            - 5.138e-3)*T + 1.9989)*T - 0.0152)*T - 468.0927)*T
            + 84381.412;
      #else
      if( fabs(T) < 2.0 )
            eps = ((1.813e-3*T - 5.9e-4)*T - 46.8150)*T + 84381.448;

      /* This expansion is from Laskar, cited above.
       * Bretagnon and Simon say, in Planetary Programs and Tables, that it
       * is accurate to 0.1" over a span of 6000 years. Laskar estimates the
       * precision to be 0.01" after 1000 years and a few seconds of arc
       * after 10000 years.
       */
      else
            {
            eps = ((((((((( 2.45e-10*T + 5.79e-9)*T + 2.787e-7)*T
                  + 7.12e-7)*T - 3.905e-5)*T - 2.4967e-3)*T
            - 5.138e-3)*T + 1.99925)*T - 0.0155)*T - 468.093)*T
            + 84381.448;
            }
      #endif /* not SIMON */
      #endif /* not WILLIAMS */
      eps *= STR;
      coseps = cos( eps );
      sineps = sin( eps );
      jdeps = J;
      return(0);
}


int Starpos:: jtocal(double J )

{
      int month, day;
      long year, a, c, d, x, y, jd;
      int BC;
      double dd;

      if( J < 1721425.5 ) /* January 1.0, 1 A.D. */
            BC = 1;
      else
            BC = 0;

      jd = (long) (J + 0.5); /* round Julian date up to integer */

      /* Find the number of Gregorian centuries
       * since March 1, 4801 B.C.
       */
      a = (100*jd + 3204500L)/3652425L;

      /* Transform to Julian calendar by adding in Gregorian century years
       * that are not leap years.
       * Subtract 97 days to shift origin of JD to March 1.
       * Add 122 days for magic arithmetic algorithm.
       * Add four years to ensure the first leap year is detected.
       */
      c = jd + 1486;
      if( jd >= 2299160.5 )
            c += a - a/4;
      else
            c += 38;
      /* Offset 122 days, which is where the magic arithmetic
       * month formula sequence starts (March 1 = 4 * 30.6 = 122.4).
       */
      d = (100*c - 12210L)/36525L;
      /* Days in that many whole Julian years */
      x = (36525L * d)/100L;

      /* Find month and day. */
      y = ((c-x)*100L)/3061L;
      day = (int) (c - x - ((306L*y)/10L));
      month = (int) (y - 1);
      if( y > 13 )
            month -= 12;

      /* Get the year right. */
      year = d - 4715;
      if( month > 2 )
            year -= 1;

      /* Day of the week. */
      a = (jd + 1) % 7;

      /* Fractional part of day. */
      dd = day + J - jd + 0.5;

      /* post the year. */
      cyear = year;

      if( BC )
            {
            year = -year + 1;
            cyear = -year;
            if( prtflg )
                  printf( "%ld B.C. ", year );
            }
      else
            {
            if( prtflg )
                  printf( "%ld ", year );
            }
      day = (int) dd;
      if( prtflg )
            //printf( "%s %d %s", months[month-1], day, days[(int) a] );

      /* Flag last or first day of year */
      if( ((month == 1) && (day == 1))
            || ((month == 12) && (day == 31)) )
            yerend = 1;
      else
            yerend = 0;

      /* Display fraction of calendar day
       * as clock time.
       */
      a = (long) dd;
      dd = dd - a;
      if( prtflg )
            {
            hms( 2.0*PI*dd );
/*          if( J == TDT )
                  //printf( "TDT\n" ); /* Indicate Terrestrial Dynamical Time
            else if( J == UT )
                  //printf( "UT\n" ); /* Universal Time
            else
                  //printf( "\n" );*/
            }
      return(0);
}









int Starpos::fk4fk5(double p[],double m[], struct star *el )
{
      double a, b, c;
      double *u, *v;
      double R[6];
      int i, j;

      printf( "Converting to FK5 system\n" );

      /* Note the direction vector and motion vector
       * are already supplied by rstar.c.
       */
      a = 0.0;
      b = 0.0;
      for( i=0; i<3; i++ )
            {
            m[i] *= RTS; /* motion must be in arc seconds per century */
            a += A[i] * p[i];
            b += Ad[i] * p[i];
            }
      /* Remove E terms of aberration from FK4
       */
      for( i=0; i<3; i++ )
            {
            R[i] = p[i] - A[i] + a * p[i];
            R[i+3] = m[i] - Ad[i] + b * p[i];
            }

      /* Perform matrix multiplication
       */
      v = &Mat[0];
      for( i=0; i<6; i++ )
            {
            a = 0.0;
            u = &R[0];
            for( j=0; j<6; j++ )
                  a += *u++ * *v++;
            if( i < 3 )
                  p[i] = a;
            else
                  m[i-3] = a;
            }

      /* Transform the answers into J2000 catalogue entries
       * in radian measure.
       */

      b = p[0]*p[0] + p[1]*p[1];
      a = b + p[2]*p[2];
      c = a;
      a = sqrt(a);

      el->ra = zatan2( p[0], p[1] );
      el->dec = asin( p[2]/a );

      /* Note motion converted back to radians per (Julian) century */
      el->mura = (p[0]*m[1] - p[1]*m[0])/(RTS*b);
      el->mudec = (m[2]*b - p[2]*(p[0]*m[0] + p[1]*m[1]) )/(RTS*c*sqrt(b));

      if( el->px > 0.0 )
            {
            c = 0.0;
            for( i=0; i<3; i++ )
                  c += p[i] * m[i];

      /* divide by RTS to deconvert m (and therefore c)
       * from arc seconds back to radians
       */
            el->v = c/(21.094952663 * el->px * RTS * a);
            }
      el->px = el->px/a;      /* a is dimensionless */
      el->epoch = J2000;

      #if DEBUG
      /* Display the computed J2000 catalogue entries
       */
      hms( el->ra );
      dms( el->dec );

      u = (double *)&el->px;
      for( i=0; i<3; i++ )
            printf( " %.4f ", *u++ * RTS );
      printf( "\n" );
      #endif
      return(0);
}



/* Precession coefficients from Simon et al: */





/* Subroutine arguments:
 *
 * R = rectangular equatorial coordinate vector to be precessed.
 *     The result is written back into the input vector.
 * J = Julian date
 * direction =
 *      Precess from J to J2000: direction = 1
 *      Precess from J2000 to J: direction = -1
 * Note that if you want to precess from J1 to J2, you would
 * first go from J1 to J2000, then call the program again
 * to go from J2000 to J2.
 */

int Starpos::precess(double R[],double J,int direction)

{
      double A, B, T, pA, W, z;
      double x[3];
      double *p;
      int i;
      #if IAU
      double sinth, costh, sinZ, cosZ, sinz, cosz, Z, TH;
      #endif

      if( J == J2000 )
            return(0);
      /* Each precession angle is specified by a polynomial in
       * T = Julian centuries from J2000.0.  See AA page B18.
       */
      T = (J - J2000)/36525.0;

      #if IAU
      /* Use IAU formula only for a few centuries, if at all.  */
      if( FABS(T) > Two )
            goto laskar;

      Z =  (( 0.017998*T + 0.30188)*T + 2306.2181)*T*STR;
      z =  (( 0.018203*T + 1.09468)*T + 2306.2181)*T*STR;
      TH = ((-0.041833*T - 0.42665)*T + 2004.3109)*T*STR;

      sinth = SIN(TH);
      costh = COS(TH);
      sinZ = SIN(Z);
      cosZ = COS(Z);
      sinz = SIN(z);
      cosz = COS(z);
      A = cosZ*costh;
      B = sinZ*costh;

      if( direction < 0 )
            { /* From J2000.0 to J */
            x[0] =    (A*cosz - sinZ*sinz)*R[0]
                        - (B*cosz + cosZ*sinz)*R[1]
                                      - sinth*cosz*R[2];

            x[1] =    (A*sinz + sinZ*cosz)*R[0]
                        - (B*sinz - cosZ*cosz)*R[1]
                                      - sinth*sinz*R[2];

            x[2] =              cosZ*sinth*R[0]
                                      - sinZ*sinth*R[1]
                                             + costh*R[2];
            }
      else
            { /* From J to J2000.0 */
            x[0] =    (A*cosz - sinZ*sinz)*R[0]
                        + (A*sinz + sinZ*cosz)*R[1]
                                      + cosZ*sinth*R[2];

            x[1] =   -(B*cosz + cosZ*sinz)*R[0]
                        - (B*sinz - cosZ*cosz)*R[1]
                                      - sinZ*sinth*R[2];

            x[2] =             -sinth*cosz*R[0]
                                      - sinth*sinz*R[1]
                                             + costh*R[2];
            }
      goto done;

      laskar:
      #endif /* IAU */

      /* Implementation by elementary rotations using Laskar's expansions.
       * First rotate about the x axis from the initial equator
       * to the ecliptic. (The input is equatorial.)
       */
      if( direction == 1 )
            epsiln( J ); /* To J2000 */
      else
            epsiln( J2000 ); /* From J2000 */
      x[0] = R[0];
      z = coseps*R[1] + sineps*R[2];
      x[2] = -sineps*R[1] + coseps*R[2];
      x[1] = z;

      /* Precession in longitude
       */
      T /= 10.0; /* thousands of years */
      p = pAcof;
      pA = *p++;
      for( i=0; i<9; i++ )
            pA = pA * T + *p++;
      pA *= STR * T;

      /* Node of the moving ecliptic on the J2000 ecliptic.
       */
      p = nodecof;
      W = *p++;
      for( i=0; i<10; i++ )
            W = W * T + *p++;

      /* Rotate about z axis to the node.
       */
      if( direction == 1 )
            z = W + pA;
      else
            z = W;
      B = COS(z);
      A = SIN(z);
      z = B * x[0] + A * x[1];
      x[1] = -A * x[0] + B * x[1];
      x[0] = z;

      /* Rotate about new x axis by the inclination of the moving
       * ecliptic on the J2000 ecliptic.
       */
      p = inclcof;
      z = *p++;
      for( i=0; i<10; i++ )
            z = z * T + *p++;
      if( direction == 1 )
            z = -z;
      B = COS(z);
      A = SIN(z);
      z = B * x[1] + A * x[2];
      x[2] = -A * x[1] + B * x[2];
      x[1] = z;

      /* Rotate about new z axis back from the node.
       */
      if( direction == 1 )
            z = -W;
      else
            z = -W - pA;
      B = COS(z);
      A = SIN(z);
      z = B * x[0] + A * x[1];
      x[1] = -A * x[0] + B * x[1];
      x[0] = z;

      /* Rotate about x axis to final equator.
       */
      if( direction == 1 )
            epsiln( J2000 );
      else
            epsiln( J );
      z = coseps * x[1] - sineps * x[2];
      x[2] = sineps * x[1] + coseps * x[2];
      x[1] = z;

      #if IAU
      done:
      #endif

      for( i=0; i<3; i++ )
            R[i] = x[i];
      return(0);
}

double Starpos::deltat(double Y)

{
      double ans, p, B;
      int d[6];
      int i, iy, k;

      if( dtgiven != 0.0 )
            return( dtgiven );

      if( Y > TABEND )
            {
      #if 0
      /* Morrison, L. V. and F. R. Stephenson, "Sun and Planetary System"
       * vol 96,73 eds. W. Fricke, G. Teleki, Reidel, Dordrecht (1982)
       */
            B = 0.01*(Y-1800.0) - 0.1;
            ans = -15.0 + 32.5*B*B;
            return(ans);
      #else
      /* Extrapolate forward by a second-degree curve that agrees with
         the most recent data in value and slope, and vaguely fits
         over the past century.  This idea communicated by Paul Muller,
         who says NASA used to do something like it.  */
            B = Y - TABEND;
            /* slope */
            p = dt[TABSIZ-1] - dt[TABSIZ-2];
            /* square term */
            ans = (dt[TABSIZ - 101] - (dt[TABSIZ - 1] - 100.0 * p)) * 1e-4;
                  ans = 0.01 * (dt[TABSIZ-1] + p * B + ans * B * B);
            if( prtflg )
              printf("[extrapolated deltaT] ");
            return(ans);
      #endif
            }

      if( Y < TABSTART )
            {
            if( Y >= 948.0 )
                  {
      /* Stephenson and Morrison, stated domain is 948 to 1600:
       * 25.5(centuries from 1800)^2 - 1.9159(centuries from 1955)^2
       */
                  B = 0.01*(Y - 2000.0);
                  ans = (23.58 * B + 100.3)*B + 101.6;
                  }
            else
                  {
      /* Borkowski */
                  B = 0.01*(Y - 2000.0)  +  3.75;
                  ans = 35.0 * B * B  +  40.;
                  }
            return(ans);
            }

      /* Besselian interpolation from tabulated values.
       * See AA page K11.
       */

      /* Index into the table.
       */
      p = floor(Y);
      iy = (int) (p - TABSTART);
      /* Zeroth order estimate is value at start of year
       */
      ans = dt[iy];
      k = iy + 1;
      if( k >= TABSIZ )
            goto done; /* No data, can't go on. */

      /* The fraction of tabulation interval
       */
      p = Y - p;

      /* First order interpolated value
       */
      ans += p*(dt[k] - dt[iy]);
      if( (iy-1 < 0) || (iy+2 >= TABSIZ) )
            goto done; /* can't do second differences */

      /* Make table of first differences
       */
      k = iy - 2;
      for( i=0; i<5; i++ )
            {
            if( (k < 0) || (k+1 >= TABSIZ) )
                  {
                  d[i] = 0;
                  }
            else
                  d[i] = dt[k+1] - dt[k];
            k += 1;
            }

      /* Compute second differences
       */
      for( i=0; i<4; i++ )
            d[i] = d[i+1] - d[i];
      B = 0.25*p*(p-1.0);
      ans += B*(d[1] + d[2]);
      #if DEMO
      printf( "B %.4lf, ans %.4lf\n", B, ans );
      #endif
      if( iy+2 >= TABSIZ )
            goto done;

      /* Compute third differences
       */
      for( i=0; i<3; i++ )
            d[i] = d[i+1] - d[i];
      B = 2.0*B/3.0;
      ans += (p-0.5)*B*d[1];
      #if DEMO
      printf( "B %.4lf, ans %.4lf\n", B*(p-0.5), ans );
      #endif
      if( (iy-2 < 0) || (iy+3 > TABSIZ) )
            goto done;

      /* Compute fourth differences
       */
      for( i=0; i<2; i++ )
            d[i] = d[i+1] - d[i];
      B = 0.125*B*(p+1.0)*(p-2.0);
      ans += B*(d[0] + d[1]);
      #if DEMO
      printf( "B %.4lf, ans %.4lf\n", B, ans );
      #endif

      done:
      /* Astronomical Almanac table is corrected by adding the expression
       *     -0.000091 (ndot + 26)(year-1955)^2  seconds
       * to entries prior to 1955 (AA page K8), where ndot is the secular
       * tidal term in the mean motion of the Moon.
       *
       * Entries after 1955 are referred to atomic time standards and
       * are not affected by errors in Lunar or planetary theory.
       */
      ans *= 0.01;
      if( Y < 1955.0 )
            {
            B = (Y - 1955.0);
      #if 1
            ans += -0.000091 * (-25.8 + 26.0) * B * B;
      #else
            ans += -0.000091 * (-23.8946 + 26.0) * B * B;
      #endif
            }
      return( ans );
}

int Starpos::update()
{
      double T;

      /* Convert Julian date to Julian years re J2000.0.
       */
      T = 2000.0 + (JD - J2000)/365.25;

      switch( jdflag )
            {
            case 0:
                  TDT = JD;
                  UT = JD;
                  break;
            case 1:
                  TDT = JD;
                  UT = TDT - deltat(T)/86400.0;
                  jtocal(UT); /* display the other date */
                  break;
            case 2:
                  UT = JD;
                  TDT = UT + deltat(T)/86400.0;
                  jtocal(TDT);
                  break;
            }
      jtocal(JD);
      return(0);
}

int Starpos::dostar()
{
      
        rstar( (struct star *) elobject );
        return 0;
}

 int Starpos::mainStar(int Year, int month, int hour, int day, int min, int sec)
{
      int i;
      kinit();
      char * buff;
      //loop:
      objnum = 88;
      //prtflg = 1;
      //printf( "Enter starting date of tabulation\n" );
      JD = caltoj(Year,month,day); /* date */
      double t = (3600.0*hours + 60.0*minutes + seconds)/86400.0;;      /* time of day */
      JD += t;
      update(); /* find UT and ET */

      if(objnum ==88)
      {
              elobject = (struct orbit *)&fstar;
              i = getstar( (struct star *) elobject );

      }

      if( elobject == (struct orbit *)&fstar )
            showcname( &elobject->obname[0] );
      else if( elobject )
            buff = &elobject->obname[0];



      /* Always calculate heliocentric position of the earth */
            kepler( TDT, &earth, rearth, eapolar );

            if( objnum ==88)
                  {
                   dostar();
                  }

      }


int Starpos::rstar(struct star *el )

{
      double p[3], q[3], e[3], m[3], temp[3], polar[3];
      double T, vpi, epoch;
      double cosdec, sindec, cosra, sinra;
      int i;
      double log();


      /* Convert from RA and Dec to equatorial rectangular direction
       */
      loop:
      cosdec = cos( el->dec );
      sindec = sin( el->dec );
      cosra = cos( el->ra );
      sinra = sin( el->ra );
      q[0] = cosra * cosdec;
      q[1] = sinra * cosdec;
      q[2] = sindec;



      /* space motion */
      vpi = 21.094952663 * el->v * el->px;
      m[0] =    -el->mura * cosdec * sinra
            - el->mudec * sindec * cosra
            +       vpi * q[0];

      m[1] =     el->mura * cosdec * cosra
            - el->mudec * sindec * sinra
            +       vpi * q[1];

      m[2] =    el->mudec * cosdec +   vpi * q[2];

      epoch = el->epoch;

      /* Convert FK4 to FK5 catalogue */

      if( epoch == B1950 )
            {
            fk4fk5( q, m, el );
            goto loop;
            }


      for( i=0; i<3; i++ )
            e[i] = rearth[i];

      /* precess the earth to the star epoch */
      precess( e, epoch, -1 );

      /* Correct for proper motion and parallax
       */
      T = (TDT - epoch)/36525.0;
      for( i=0; i<3; i++ )
            {
            p[i] = q[i]  +  T * m[i]  -  el->px * e[i];
            }

      /* precess the star to J2000 */
      precess( p, epoch, 1 );
      /* reset the earth to J2000 */
      for( i=0; i<3; i++ )
            e[i] = rearth[i];

      /* Find Euclidean vectors between earth, object, and the sun
       * angles( p, q, e );
       */
      angles( p, p, e );//changed q for p

      /* Find unit vector from earth in direction of object
       */
      for( i=0; i<3; i++ )
            {
            p[i] /= EO;
            temp[i] = p[i];
            }

      if( prtflg )
            {
      //    printf( "approx. visual magnitude %.1f\n", el->mag );
      /* Report astrometric position
       */
            showrd((char*) "Astrometric J2000.0:", p, polar );

      /* Also in 1950 coordinates
       */
            precess( temp, B1950, -1 );
            showrd((char*) "Astrometric B1950.0:", temp, polar );

      /* For equinox of date: */
            for( i=0; i<3; i++ )
                  temp[i] = p[i];
            precess( temp, TDT, -1 );
            
            showrd((char*) "Astrometric of date:", temp, polar );
            }

      /* Correct position for light deflection
       * relativity( p, q, e );
       */
      relativity( p, p, e );


      /* Correct for annual aberration
       */
      annuab( p );

      /* Precession of the equinox and ecliptic
       * from J2000.0 to ephemeris date
       */
      precess( p, TDT, -1 );

      /* Ajust for nutation
       * at current ecliptic.
       */
      epsiln( TDT );
      nutate( TDT, p );
      
      /* Display the final apparent R.A. and Dec.
       * for equinox of date.
       */
      flag = true;
      showrd((char*) "    Apparent:", p, polar );
      flag = false;

      /* Go do topocentric reductions.
       */
      dradt = 0.0;
      ddecdt = 0.0;
      polar[2] = 1.0e38; /* make it ignore diurnal parallax */
      altaz( polar, UT );
      return(0);
}

int Starpos::kinit()
{
double a, b, fl, co, si, u;
FILE *f;// *fopen();
char s[84];

//printf( "\n\tSteve Moshier's Ephemeris Program v5.6\n\n" );
//printf( "Planetary and lunar positions approximate DE404.\n" );

f = fopen( "aa.ini", "r" );
if( f )
      {
      char* res =fgets( s, 80, f );
      sscanf( s, "%lf", &tlong );
      res =fgets( s, 80, f );
      sscanf( s, "%lf", &glat );
      res =fgets( s, 80, f );
      sscanf( s, "%lf", &height );
      u = glat * DTR;

/* Reduction from geodetic latitude to geocentric latitude
 * AA page K5
 */
      co = cos(u);
      si = sin(u);
      fl = 1.0 - 1.0/flat;
      fl = fl*fl;
      si = si*si;
      u = 1.0/sqrt( co*co + fl*si );
      a = aearth*u + height;
      b = aearth*fl*u  +  height;
      trho = sqrt( a*a*co*co + b*b*si );
      tlat = RTD * acos( a*co/trho );
      if( glat < 0.0 )
            tlat = -tlat;
      trho /= aearth;

/* Reduction from geodetic latitude to geocentric latitude///////////////////////////////////
 * AA page K5


      tlat = glat
            - 0.19242861 * sin(2.0*u)
            + 0.00032314 * sin(4.0*u)
            - 0.00000072 * sin(6.0*u);

      trho =    0.998327073
            + 0.001676438 * cos(2.0*u)
            - 0.000003519 * cos(4.0*u)
            + 0.000000008 * cos(6.0*u);
      trho += height/6378160.;/////////////////////////////*/

      //printf( "Terrestrial east longitude %.4f deg\n", tlong );
      //printf( "geocentric latitude %.4f deg\n", tlat );
      //printf( "Earth radius %.5f\n", trho );

      res =fgets( s, 80, f );
      sscanf( s, "%lf", &attemp );
      //printf( "temperature %.1f C\n", attemp );
      res =fgets( s, 80, f );
      sscanf( s, "%lf", &atpress );
      //printf( "pressure %.0f mb\n", atpress );
      res =fgets( s, 80, f );
      sscanf( s, "%d", &jdflag );
/*    switch( jdflag )
            {
            case 0: printf("TDT and UT assumed equal.\n");
                  break;
            case 1: printf("Input time is TDT.\n" );
                  break;
            case 2: printf("Input time is UT.\n" );
                  break;
            default: printf("Illegal jdflag\n" );
            exit(0);
            }*/
      res =fgets( s, 80, f );
      sscanf( s, "%lf", &dtgiven );
      if( dtgiven != 0.0 )
            //printf( "Using deltaT = %.2fs.\n", dtgiven );
      fclose(f);
      }
Clightaud = 86400.0 * Clight / au;
/* Radius of the earth in au
   Thanks to Min He <Min.He@businessobjects.com> for pointing out
res   this needs to be initialized early.  */
Rearth = 0.001 * aearth / au;
return(0);
}



/* Program to read in a file containing orbital parameters
 */
extern struct orbit earth;

int Starpos::getorbit(struct orbit *el)

{
FILE *f;
char s1[128], s2[128], *u, *v;
int i;


getnum((char *) "Name of orbit catalogue file: ", orbnam, strfmt );
f = fincat( orbnam, 2, s1, s2 );
if( f == 0 )
      return(-1); /* failure flag */


//printf( "%s\n", s1 );
//printf( "%s\n", s2 );

/* Read in ASCII floating point numbers
 */
sscanf( s1, "%lf %lf %lf %lf %lf %lf",
      &el->epoch, &el->i, &el->W, &el->w, &el->a, &el->dm );

sscanf( s2, "%lf %lf %lf %lf %lf %15s", &el->ecc, &el->M,
      &el->equinox, &el->mag, &el->sdiam, &el->obname[0] );

el->obname[15] = '\0';

/* Clear out the rest of the data structure
 */
el->ptable = 0;
el->L = 0.0;
el->r = 0.0;
el->plat = 0.0;
if( strcmp( &el->obname[0], "Earth" )  )
      {
      return(0);
      }
else
      {
      u = (char *)&earth;
      v = (char *)el;
      for( i=0; i < (int) sizeof(struct orbit); i++ )
            *u++ = *v++;
      printf( "Read in earth orbit\n" );
      return(1);
      }
}

/* Open catalogue and find line number
 */
FILE *Starpos::fincat(char * name,int n,char * str1,char * str2 )

//int n;   /* number of lines per catalogue entry */

{
      strcpy(name,pathname);
      int i;
      FILE *f;// *fopen();
      char *res;

      f = fopen(name, "r" );
      if( f == 0 )
            {
            printf( "Can't find file %s\n", name );
            return(0); /* failure flag */
            }

      //getnum((char*) "Line number", &linenum, intfmt );
      //linenum = 1; //prova di input

      if( linenum <= 0 )
            goto failure;

      for( i=0; i<linenum; i++ )
            {

            res =fgets( str1, 126, f );
            if( *str1 == '-' )
                  goto endf;
            if( n > 1 )
                  {
                  res=fgets( str2, 126, f );
                  if( *str2 == '-' )
                        goto endf;
                  }
            }
      fclose(f);
      return( f );


      endf:
      printf( "End of file reached.\n" );
      failure:
      fclose(f);
      return(0);
}

int Starpos:: getnum(char* msg,void* num,const char* format )

{
      char s[40];

      printf( "%s (", msg );
      if( format == strfmt )
            printf( format, (char *) num );
      else if( format == dblfmt )
            printf( format, *(double *)num );
      else if( format == intfmt )
            printf( format, *(int *)num );
      else if( format == lngfmt )
            printf( format, *(long *)num );
      else
            printf( "Illegal input format\n"  );
      printf( ") ? ");
      //char * res = gets(s);
      if( s[0] != '\0' )
            sscanf( s, format, num );
      return(0);
}

int Starpos:: getstar(struct star *el)

{
      int sign;
      char s[128];
      double rh, rm,rs, dd, dm, ds, x, z;

      FILE *f;
      char *p;
      char *str;
      char buffer[10];
      int i;

      f = fincat( starnam, 1, s, (char *)0 );
      if( f == 0 )
            return(-1); /* failure flag */

      /* Read in the ASCII string data and name of the object*/
      str = strtok(s," ");
      el->epoch = atof(str);
      str = strtok(NULL," ");
      rh = atof(str);
      str = strtok(NULL," ");
      rm = atof(str);
      str = strtok(NULL," ");
      rs = atof(str);
      str = strtok(NULL," ");
      dd = atof(str);
      str = strtok(NULL," ");
      dm = atof(str);
      str = strtok(NULL," ");
      ds = atof(str);
      str = strtok(NULL," ");
      el->mura = atof(str);
      str = strtok(NULL," ");
      el->mudec = atof(str);
      str = strtok(NULL," ");
      el->v = atof(str);
      str = strtok(NULL," ");
      el->px = atof(str);
      str = strtok(NULL," ");
      el->mag = atof(str);
      str = strtok(NULL," ");
      strcpy(el->obname,str);
      constname =el->obname;

/*    sscanf(s,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %s",&el->epoch, &rh, &rm, &rs, &dd, &dm, &ds,
                  &el->mura, &el->mudec, &el->v, &el->px, &el->mag, &el->obname[0] );*/
      
      x = el->epoch;
      if( x == 2000.0 )
            x = J2000;
      else if( x == 1950.0 )
            x = B1950;
      else if( x == 1900.0 )
            x = J1900;
      else
            x = J2000  +  365.25 * (x - 2000.0);
      el->epoch = x;

      /* read the right ascension */
      el->ra = 2.0 * PI * (3600.0*rh + 60.0*rm + rs)/86400.0;
      
      /* read the declination */
      sign = 1;
      if( (dd < 0.0) || (dm < 0.0) || (ds < 0.0) )
            sign = -1;
      z = (3600.0*fabs(dd) + 60.0*fabs(dm) + fabs(ds))/RTS;
      if( dd == 0.0 )
            {
      /* Scan the text for possible minus sign in front of declination 0 */
            p = s;
      /* skip over 4 fields */
            for( i=0; i<4; i++ )
                  {
                  while( *p++ == ' ' )
                        ;
                  while( *p++ != ' ' )
                        ;
                  }
            while( *p++ == ' ' )
                  ;
            --p;
            if( *p == '-' )
                  sign = -1;
            }
      if( sign < 0 )
            z = -z;
      el->dec = z;

      #if DEBUG
/*    printf( "%.2f\n", el->epoch );
      printf( "%.0f %.0f %.3f\n", rh, rm, rs );
      printf( "%.8f\n", el->ra );
      printf( "%.0f %.0f %.3f\n", dd, dm, ds );
      printf( "%.8f\n", el->dec );
      printf( "d %.3f mua %.3f mud %.3f v %.3f\n",
            el->px, el->mura, el->mudec, el->v );*/
      #endif

      el->mura *= 15.0/RTS;   /* s/century -> "/century -> rad/century */
      el->mudec /= RTS;
      z = el->px;
      if( z < 1.0 )
            {
            if( z <= 0.0 )
                  el->px = 0.0;
            else
                  el->px = STR * z;  /* assume px in arc seconds */
            }
      else
            {
            el->px = 1.0/(RTS * z); /* parsecs -> radians */
            }
      return(0);
}

int Starpos::showcname( char *in )

{
      char *g, *p, *q;
      char ans[80];
      int i;


      p = in;
      q = ans;


      skipwh(p);
      if( isnumber(p) )
            {
            while( isnumber(p) )
                  *q++ = *p++;
            }
      skipwh(p);
      *q++ = ' ';

      if( islow(p) )
            {
            for( i=0; i<NGREEK; i++ )
                  {
                  g =(char *) greek[i];
                  if( (*p == *g) && (*(p+1) == *(g+1)) )
                        break;
                  }
            if( i < NGREEK )
                  {
                  while( *g != '\0' )
                        *q++ = *g++;
                  }
            while( islow(p) )
                  ++p;
            }
      skipwh(p);
      /* copy things like "-a" until uppercase letter found */
      while( (*p != '\0') && !isup(p) )
            *q++ = *p++;

      *q++ = ' ';

      if( isup(p) )
            {
      /* Check the list of constellation names */
            for( i=0; i<NCON; i++ )
                  {
                  g =(char *)constel[i];
                  if( (*p == *g) && ( *(p+1) == *(g+1) )
                        && ( *(p+2) == *(g+2) ) )
                        break;
                  }
      /* Get the name found */
            if( i < NCON )
                  {
                  g += 4;
                  while( *g != '\0' )
                        *q++ = *g++;
                  p += 3;
                  }
            }
      skipwh(p);
      *q++ = ' ';
      while( *p )
            *q++ = *p++;
      *q++ = '\0';
      /* convert all '_' characters to spaces */
      q = ans;
      while( *q != '\0' )
            {
            if( *q == '_' )
                  *q = ' ';
            ++q;
            }
      //printf( "\n              %s\n", ans );
      return(0);
}




int Starpos:: islow(char* p)

{
      if( (*p >= 'a') && (*p <= 'z') )
            return(1);
      else
            return(0);
}

int Starpos:: isnumber(char *p)

{
      if( (*p >= '0') && (*p <= '9') )
            return(1);
      else
            return(0);
}




int Starpos::skipwh(char *p)
{
      while( ((*p == ' ') || (*p == '\t') || (*p == '_'))
            && (*p != '\0') && (*p != '\n') && (*p != '\r') )
                  ++p;
      return(0);
}

int Starpos:: isup(char *p)

{
      if( (*p >= 'A') && (*p <= 'Z') )
            return(1);
      else
            return(0);
}

double Starpos:: zgetdate()
{
      double J;

      /* Get operator to type in a date.
       */
      getnum((char*) "Calendar date: Year", &cyear, lngfmt );

      if( (cyear > 53994L) || (cyear < -4713L) )
            {
            printf( "Year out of range.\n" );
            goto err;
            }
      if( cyear == 0 )
            {
            printf( "There is no year 0.\n" );
      err:  J = 0.0;
            goto pdate;
            }

      //getnum( (char*)"Month (1-12)", &month, intfmt);
      //getnum( (char*)"Day.fraction", &day, dblfmt );

      /* Find the Julian day. */
      J = caltoj(cyear,month,day);
      /*printf( "Julian day %.1f\n", J );*/

      pdate:
      /* Convert back to calendar date. */
      /* jtocal( J ); */
      return(J);
}

double Starpos:: gethms()
{
      double t;

      getnum((char*) "Time: Hours", &hours, intfmt );
      getnum( (char*)"Minutes", &minutes, intfmt );
      getnum( (char*)"Seconds", &seconds, dblfmt );
      t = (3600.0*hours + 60.0*minutes + seconds)/86400.0;
      return(t);
}

double Starpos::zatan2(double x,double y )
{
      double z, w;
      short code;


      code = 0;

      if( x < 0.0 )
            code = 2;
      if( y < 0.0 )
            code |= 1;

      if( x == 0.0 )
            {
            if( code & 1 )
                  return( 1.5*PI );
            if( y == 0.0 )
                  return( 0.0 );
            return( 0.5*PI );
            }

      if( y == 0.0 )
            {
            if( code & 2 )
                  return( PI );
            return( 0.0 );
            }


      switch( code )
            {
            default:
            case 0: w = 0.0; break;
            case 1: w = 2.0 * PI; break;
            case 2:
            case 3: w = PI; break;
            }

      z = atan( y/x );

      return( w + z );
}

int Starpos::mainPlanet(int Year,int month, int day, int hour, int minute, int second)
{
      int i;
      kinit();
      prtflg = 1;
      //printf( "Enter starting date of tabulation\n" );
      JD = caltoj(Year,month,day); /* date */
      double t = (3600.0*hour + 60.0*minute + second)/86400.0;;   /* time of day */
      JD += t;
      update(); /* find UT and ET */
      //printf( "Julian day %.7f\n", JD );

      /*getnum( "Enter interval between tabulations in days", &djd, dblfmt );
      getnum( "Number of tabulations to display", &ntab, intfmt );
      if( ntab <= 0 )
            ntab = 1;

      loop1:
      getnum( "Planet number 0-9 or 88 to read star, 99 to read orbit",
            &objnum, intfmt );*/
       ////////////////////da cambiare per selezionaere pianeta//////////////////
      switch(objnum)
            {
            case -1: exit(0);
            
            case 1: elobject = &venus; break;
            
            case 2: elobject = &mars; break;
            case 3: elobject = &jupiter; break;
            case 4: elobject = &saturn; break;

            
            }

      if( elobject == (struct orbit *)&fstar )
            showcname( &elobject->obname[0] );
      else if( elobject )
            bodyname = &elobject->obname[0];


      
      /* print Julian date */

            update();

      /* Always calculate heliocentric position of the earth */
            kepler( TDT, &earth, rearth, eapolar );

            switch( objnum )
                  {
                  //case 0: dosun();  break;
                  case 0: domoon();  break;
                  default: doplanet();  break;
                  }
            //printf( "\n" );
            JD += djd;


      #ifdef _MSC_VER
      return 0;
      #endif
}

int Starpos::domoon()
{
      int i, prtsav;
      double ra0, dec0;
      double x, y, z, lon0;
      double pp[3], qq[3], pe[3], re[3], moonpp[3], moonpol[3];

      /* Geometric equatorial coordinates of the earth.  */
      for (i = 0; i < 3; i++)
            re[i] = rearth[i];

      /* Run the orbit calculation twice, at two different times,
       * in order to find the rate of change of R.A. and Dec.
       */

      /* Calculate for 0.001 day ago
       */
      prtsav = prtflg;
      prtflg = 0; /* disable display */
      moonll(TDT-0.001, moonpp, moonpol);
      ra0 = ra;
      dec0 = dec;
      lon0 = moonpol[0];
      prtflg = prtsav;

      /* Calculate for present instant.
       */
      moonll(TDT, moonpp, moonpol);

      if(prtflg)
        printf("Geometric lon %.3f deg, lat %.3f deg, rad %.4e au\n",
             RTD * obpolar[0], RTD * obpolar[1], obpolar[2]);
      prtflg = 0;

      /* The rates of change.  These are used by altaz() to
       * correct the time of rising, transit, and setting.
       */
      dradt = ra - ra0;
      if (dradt >= PI)
        dradt = dradt - 2.0 * PI;
      if (dradt <= -PI)
        dradt = dradt + 2.0 * PI;
      dradt = 1000.0 * dradt;
      ddecdt = 1000.0*(dec-dec0);

      /* Rate of change in longitude, degrees per day
       * used for phase of the moon
       */
      lon0 = 1000.0*RTD*(moonpol[0] - lon0);

      /* Get apparent coordinates for the earth.  */
      z = re[0] * re[0] + re[1] * re[1] + re[2] * re[2];
      z = sqrt(z);
      for (i = 0; i < 3; i++)
        re[i] /= z;
      annuab( re );  /* aberration of light.  */
      /* pe[0] -= STR * (20.496/(RTS*pe[2])); */
      precess( re, TDT, -1 );
      nutate( TDT, re );
      for (i = 0; i < 3; i++)
        re[i] *= z;
      lonlat( re, TDT, pe, 0 );
      prtflg = prtsav; /* reenable display */

      /* Find sun-moon-earth angles */
      for( i=0; i<3; i++ )
            qq[i] = re[i] + moonpp[i];
      angles( moonpp, qq, re );

      /* Display answers
       */
      if( prtflg )
            {

            mlong =     RTD * moonpol[0];//latitude
            mlat = RTD * moonpol[1];//longitude
            printf( "Distance %.3f Earth-radii\n", moonpol[2]/Rearth );
            printf( "Horizontal parallax" );
            x = Rearth/moonpol[2];
            prlx = asin(x) * RTD;//parallax
            dms( asin(x) ) ;
            printf( "Semidiameter" );
            x = 0.272453 * x  +  0.0799/RTS; /* AA page L6 */
            sdiam = x * RTD;//semidiameter
            dms( x );

            x = RTD * acos(-ep);
      /*    x = 180.0 - RTD * arcdot (re, pp); */
            printf( "\nElongation from sun %.2f deg,", x );
            x = 0.5 * (1.0 + pq);
            printf( "  Illuminated fraction %.2f\n", x );

      /* Find phase of the Moon by comparing Moon's longitude
       * with Earth's longitude.
       *
       * The number of days before or past indicated phase is
       * estimated by assuming the true longitudes change linearly
       * with time.  These rates are estimated for the date, but
       * do not stay constant.  The error can exceed 0.15 day in 4 days.
       */
            x = moonpol[0] - pe[0];
            x = modtp( x ) * RTD;   /* difference in longitude */
            i = (int) (x/90); /* number of quarters */
            x = (x - i*90.0); /* phase angle mod 90 degrees */

      /* days per degree of phase angle */
            z = moonpol[2]/(12.3685 * 0.00257357);

            if( x > 45.0 )
                  {
                  y = -(x - 90.0)*z;
                  if( y > 1.0 )
                        printf( "Phase %.1f days before ", y );
                  else
                        printf( "Phase %.2f days before ", y );
                  i = (i+1) & 3;
                  }
            else
                  {
                  y = x*z;
                  if( y > 1.0 )
                        printf( "Phase %.1f days past ", y );
                  else
                        printf( "Phase %.2f days past ", y );
                  }


            switch(i)
                  {
                  case 0: printf( "Full Moon\n" ); break;
                  case 1: printf( "Third Quarter\n" ); break;
                  case 2: printf( "New Moon\n" ); break;
                  case 3: printf( "First Quarter\n" ); break;
                  }
            fPlanet =true;
            printf( "    Apparent:  R.A." );
            hms(ra);
            printf( "Declination" );
            dms(dec);
            printf( "\n" );
            fPlanet = false;
            } /* if prtflg */

      /* Compute and display topocentric position (altaz.c)
       */
      pp[0] = ra;
      pp[1] = dec;
      pp[2] = moonpol[2];     

      altaz( pp, UT );

      constellation = whatconstel(ra,dec);/////////////////////////////aggiunto da me////////////////////

      return(0);
}

/* Calculate apparent latitude, longitude, and horizontal parallax
 * of the Moon at Julian date J.
 */
int Starpos::moonll(double J,double rect[],double pol[])

{
      double cosB, sinB, cosL, sinL, y, z;
      double qq[3], pp[3];
      int i;

      /* Compute obliquity of the ecliptic, coseps, and sineps.  */
      epsiln( J );
      /* Get geometric coordinates of the Moon.  */
      gmoon (J, rect, pol);
      /* Post the geometric ecliptic longitude and latitude, in radians,
       * and the radius in au.
       */
      obpolar[0] = pol[0];
      obpolar[1] = pol[1];
      obpolar[2] = pol[2];

      /* Light time correction to longitude,
       * about 0.7".
       */
      pol[0] -= 0.0118 * DTR * Rearth / pol[2];

      /* convert to equatorial system of date */
      cosB = cos(pol[1]);
      sinB = sin(pol[1]);
      cosL = cos(pol[0]);
      sinL = sin(pol[0]);
      rect[0] = cosB*cosL;
      rect[1] = coseps*cosB*sinL - sineps*sinB;
      rect[2] = sineps*cosB*sinL + coseps*sinB;

      /* Rotate to J2000. */
      precess( rect, TDT, 1 );


      /* Find Euclidean vectors and angles between earth, object, and the sun
       */
      for( i=0; i<3; i++ )
            {
            pp[i] = rect[i] * pol[2];
            qq[i] = rearth[i] + pp[i];
            }
      angles( pp, qq, rearth );

      /* Make rect a unit vector.  */
      /* for (i = 0; i < 3; i++) */
      /*  rect[i] /= EO; */

      /* Correct position for light deflection.
         (Ignore.)  */
      /* relativity( rect, qq, rearth ); */

      /* Aberration of light.
         The Astronomical Almanac (Section D, Daily Polynomial Coefficients)
         seems to omit this, even though the reference ephemeris is inertial.  */
      /* annuab (rect); */

      /* Precess to date.  */
      precess( rect, TDT, -1 );

      /* Correct for nutation at date TDT.
       */
      nutate( TDT, rect );

      /* Apparent geocentric right ascension and declination.  */
      ra = zatan2(rect[0],rect[1]);
      dec = asin(rect[2]);

      /* For apparent ecliptic coordinates, rotate from the true
         equator into the ecliptic of date.  */
      cosL = cos(eps+nuto);
      sinL  = sin(eps+nuto);
      y = cosL * rect[1] + sinL * rect[2];
      z = -sinL * rect[1] + cosL * rect[2];
      pol[0] = zatan2( rect[0], y );
      pol[1] = asin(z);

      /* Restore earth-moon distance.  */
      for( i=0; i<3; i++ )
        rect[i] *= EO;

      return(0);
}

int Starpos::doplanet()
{
  /* calculate heliocentric position of the object */
  kepler( TDT, elobject, robject, obpolar );
  /* apply correction factors and print apparent place */
  reduce( elobject, robject, rearth );
  return 0;
}

int Starpos::reduce( struct orbit *elemnt,double q[],double e[] )
//struct orbit *elemnt; /* orbital elements of q */
//double q[], e[];      /* heliocentric coordinates */
{
      double p[3], temp[3], polar[3];
      double a, b, s;
      int i;
      //double sqrt(), asin(), log();

      /* Save the geometric coordinates at TDT
       */
      for( i=0; i<3; i++ )
            temp[i] = q[i];

      /* Display ecliptic longitude and latitude, precessed to equinox
            of date.  */
      if( prtflg )
            lonlat( q, TDT, polar, 1 );

      /* Adjust for light time (planetary aberration)
      */
      lightt( elemnt, q, e );

      /* Find Euclidean vectors between earth, object, and the sun
       */
      for( i=0; i<3; i++ )
            p[i] = q[i] - e[i];

      angles( p, q, e );

      if( prtflg )
            {
            a = 0.0;
            for( i=0; i<3; i++ )
                  {
                  b = temp[i] - e[i];
                  a += b * b;
                  }
            a = sqrt(a);
            //printf( "true geocentric distance %.7f au    ", a );
            planet_distance = a;/* was EO */
            printf( "equatorial diameter %.2f\"\n", 2.0*elemnt->sdiam/EO );


      /* Calculate visual magnitude.
       * "Visual" refers to the spectrum of visible light.
       * Phase = 0.5(1+pq) = geometric fraction of disc illuminated.
       * where pq = cos( sun-object-earth angle )
       * The magnitude is
       *    V(1,0) + 2.5 log10( SE^2 SO^2 / Phase)
       * where V(1,0) = elemnt->mag is the magnitude at 1au from
       * both earth and sun and 100% illumination.
       */
            a = 0.5 * (1.0 + pq);
      /* Fudge the phase for light leakage in magnitude estimation.
       * Note this phase term estimate does not reflect reality well.
       * Calculated magnitudes of Mercury and Venus are inaccurate.
       */
            b = 0.5 * (1.01 + 0.99*pq);
            s = elemnt->mag + 2.1715 * log( EO*SO ) - 1.085*log(b);
            printf( "approx. visual magnitude %.1f, phase %.3f\n", s, a );
            }

      /* Find unit vector from earth in direction of object
       */
      for( i=0; i<3; i++ )
            {
            p[i] /= EO;
            temp[i] = p[i];
            }

      if( prtflg )
            {
      /* Report astrometric position
       */
            showrd((char *) "Astrometric J2000.0:", p, polar );

      /* Also in 1950 coordinates
       */
            precess( temp, B1950, -1 );
            showrd((char *) "Astrometric B1950.0:", temp, polar );
            }

      /* Correct position for light deflection
       */
      relativity( p, q, e );

      /* Correct for annual aberration
       */
      annuab( p );

      /* Precession of the equinox and ecliptic
       * from J2000.0 to ephemeris date
       */
      precess( p, TDT, -1 );

      /* Ajust for nutation
       * at current ecliptic.
       */
      epsiln( TDT );
      nutate( TDT, p );


      /* Display the final apparent R.A. and Dec.
       * for equinox of date.
       */

/*    if( prtflg)
      {
            constellation = whatconstel (ra, dec);
      }*/

      fPlanet = true;
      showrd((char *) "  Apparent:", p, polar );
      fPlanet = false;

      /* Geocentric ecliptic longitude and latitude.  */
      if( prtflg )
        {
            printf ("Apparent geocentric ");
            for( i=0; i<3; i++ )
            p[i] *= EO;
            lonlat( p, TDT, temp, 0 );
        }
      /* Go do topocentric reductions.
       */
      polar[2] = EO;
      altaz( polar, UT );
      return(0);
}

int Starpos::lonlat(double pp[],double J,double polar[],int ofdate )
//int ofdate; /* 1 means precess from J2000 to date J.  */
{
      double s[3], x, y, z, yy, zz, r;
      int i;

      /* Make local copy of position vector
       * and calculate radius.
       */
      r = 0.0;
      for( i=0; i<3; i++ )
            {
            x = pp[i];
            s[i] = x;
            r += x * x;
            }
      r = sqrt(r);

      /* Precess to equinox of date J
       */
      if( ofdate )
            precess( s, J, -1 );

      /* Convert from equatorial to ecliptic coordinates
       */
      epsiln(J);
      yy = s[1];
      zz = s[2];
      x  = s[0];
      y  =  coseps * yy  +  sineps * zz;
      z  = -sineps * yy  +  coseps * zz;

      yy = zatan2( x, y );
      zz = asin( z/r );

      polar[0] = yy;
      polar[1] = zz;
      polar[2] = r;

      if( prtflg == 0 )
            return(0);

      printf( "ecliptic long" );
      dms( yy );
      plong = yy * RTD;//////////////////longitude
      printf( " lat" );
      dms( zz );
      plat = zz *RTD; //////////////latitude
      printf( " rad %.6E\n", r );
      return(0);
}

/* Correction for light time from object to earth
 * including gravitational retardation due to the Sun.
 * AA page B36.
 */
int Starpos::lightt( struct orbit *elemnt,double q[],double e[] )
//double e[], q[];      /* rectangular position vectors */
//struct orbit *elemnt; /* orbital elements of object q */
{
      double p[3], p0[3], ptemp[3];
      double P, Q, E, t, x, y;
      int i, k;


      /* save initial q-e vector for display */
      for( i=0; i<3; i++ )
            {
            p0[i] = q[i] - e[i];
            }

      E = 0.0;
      for( i=0; i<3; i++ )
            E += e[i]*e[i];
      E = sqrt(E);

      for( k=0; k<2; k++ )
            {
            P = 0.0;
            Q = 0.0;
            for( i=0; i<3; i++ )
                  {
                  y = q[i];
                  x = y - e[i];
                  p[i] = x;
                  Q += y * y;
                  P += x * x;
                  }
            P = sqrt(P);
            Q = sqrt(Q);
      /* Note the following blows up if object equals sun. */
            t = (P + 1.97e-8 * log( (E+P+Q)/(E-P+Q) ) )/173.1446327;
            kepler( TDT-t, elemnt, q, ptemp );
            }

      if( prtflg )
            printf( "light time %.4fm,  ", 1440.0*t );

      /* Final object-earth vector and the amount by which it changed.
       */
      for( i=0; i<3; i++ )
            {
            x = q[i] - e[i];
            p[i] = x;
            dp[i] = x - p0[i];
            }
      showcor((char *) "aberration", p0, dp );

      /* Calculate dRA/dt and dDec/dt.
       * The desired correction of apparent coordinates is relative
       * to the equinox of date, but the coordinates here are
       * for J2000.  This introduces a slight error.
       *
       * Estimate object-earth vector t days ago.  We have
       * p(?) = q(J-t) - e(J), and must adjust to
       * p(J-t)  =  q(J-t) - e(J-t)  =  q(J-t) - (e(J) - Vearth * t)
       *         =  p(?) + Vearth * t.
       */
      velearth(TDT);
      for( i=0; i<3; i++ )
            p[i] += vearth[i]*t;

      deltap( p, p0, &dradt, &ddecdt );  /* see dms.c */
      dradt /= t;
      ddecdt /= t;
      return(0);
}

const char * Starpos::whatconstel(double ra, double dec)// (double pp[],double epoch)

{
        int i, k;
       /* double ra, dec, d;
      /*  double p[3];

        for (i = 0; i < 3; i++)
            p[i] = pp[i];

        /* Precess from given epoch to J2000.  */
       /* precess (p, epoch, 1);
        /* Precess from J2000 to Besselian epoch 1875.0.  */
      /*  precess (p, 2405889.25855, -1);
      /*  d = p[0] * p[0] + p[1] * p[1] + p[2] * p[2];
        d = sqrt (d);
        ra = zatan2 (p[1], p[0]) * (RTD * 3600. / 15.);
        if (ra < 0.0)
            ra += 86400.0;
        dec = asin (p[2] / d) * (RTD * 3600.);*/

         ra = ra * RTD * 3600. / 15.;
        if (ra < 0.0)
            ra += 86400.0;
        dec = dec * RTD * 3600.;

            /* FIND CONSTELLATION SUCH THAT THE DECLINATION ENTERED IS HIGHER THAN
               THE LOWER BOUNDARY OF THE CONSTELLATION WHEN THE UPPER AND LOWER
               RIGHT ASCENSIONS FOR THE CONSTELLATION BOUND THE ENTERED RIGHT
               ASCENSION
            */
        for (i = 0; i < NBNDRIES; i++)
            {
              k = i << 2;
              if (ra >= bndries[k] && ra < bndries[k+1] && dec > bndries[k+2])
            {
              k = bndries[k+3];
              return (constel[k]);
            }
            }
        return ("?? constellation not found");
}




Starpos::~Starpos()
{
}


Generated by  Doxygen 1.6.0   Back to index