# include "relastro.h"

# if (0)
/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
/* ra and dec are returned as decimal hours and decimal degrees. */
void lpsun (double jd, double *ra, double *dec) {

  double n, L, g, lambda,epsilon,alpha,delta,x,y,z;

  n = jd - J2000;
  L = 280.460 + 0.9856474 * n;
  g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN;
  lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN;
  epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN;

  // this is the conversion from ecliptic to celestial coords
  x = cos(lambda);
  y = cos(epsilon)*sin(lambda);
  z = sin(epsilon)*sin(lambda);

  *ra = (atan_circ(x,y))*HRS_IN_RADIAN;
  *dec = (asin(z))*DEG_IN_RADIAN;
}
# endif

/* code borrowed from Skycalc : fix this stuff XXX */
/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) {

  double n, L, g;

# define J2000 2451545.       /* Julian date at standard epoch */

  n = jd - J2000;
  L = 280.460 + 0.9856474 * n;
  g = (357.528 + 0.9856003 * n)*RAD_DEG;
  *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees
  *beta = 0.0;					  // approx latitude
  *epsilon = (23.439 - 0.0000004 * n);		  // obliquity of ecliptic in degrees
  return TRUE;
}

/* given RA, DEC, Time, calculate the parallax factor */
int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean) {

  double jd;
  double L, B, E, e, s, r, d;

  /* given a time T in UNIX seconds, determine the solar longitude S */

  jd = ohana_sec_to_jd (365.25*86400.0*(T + Tmean));
  sun_ecliptic (jd, &L, &B, &E);

  e = E * RAD_DEG;
  s = L * RAD_DEG;
  r = R * RAD_DEG;
  d = D * RAD_DEG;
  
  *pR =  +(cos(e)*sin(s)*cos(r) - cos(s)*sin(r));
  *pD =  -(cos(e)*sin(s)*sin(r) + cos(s)*cos(r))*sin(d) + sin(e)*sin(s)*cos(d);
  return TRUE;
}
