Changeset 32740

Timestamp:

Nov 22, 2011, 8:24:35 AM (14 years ago)

Author:

eugene

Message:

update for more accurate par factor (still low-ish precision)

Location:

trunk/Ohana/src/relastro

Files:

• include/relastro.h (modified) (1 diff)
• src/ParFactor.c (modified) (2 diffs)

trunk/Ohana/src/relastro/include/relastro.h

@@ -315 +315 @@
 int MeasFilterTest(Measure *measure, int applySigmaLim);
 
-int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
-int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean);
+int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius);
+int ParFactor (double *pR, double *pD, double RA, double Dec, double Time, double Tmean);
 int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB);
 int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts);

trunk/Ohana/src/relastro/src/ParFactor.c

@@ -24 +24 @@
 # endif
 
+# if (0)
 /* code borrowed from Skycalc : fix this stuff XXX */
 /* Low precision formulae for the sun, from Almanac p. C24 (1990) */
@@ -40 +41 @@
   return TRUE;
 }
+# endif
+
+/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
+int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius) {
+
+# define J2000 2451545.       /* Julian date at standard epoch */
+
+  double n = jd - J2000;              // day number
+  double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration)
+  double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly
+
+  *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees
+  *beta = 0.0;                                     // approx latitude
+  *epsilon = (23.439 - 0.0000004 * n);             // obliquity of ecliptic in degrees
+  *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU
+  return TRUE;
+}
 
 /* given RA, DEC, Time, calculate the parallax factor */
-int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean) {
+int ParFactor (double *pR, double *pD, double RA, double DEC, double Time, double Tmean) {
 
-  double jd;
-  double L, B, E, e, s, r, d;
+  double jd, lambda, beta, epsilon, Radius;
 
   /* 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);
+  jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean));
+  sun_ecliptic (jd, &lambda, &beta, &epsilon, &Radius);
 
-  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);
+  double lambda_rad = lambda*RAD_DEG;
+  double epsilon_rad = epsilon*RAD_DEG;
+  double RA_rad = RA*RAD_DEG;
+  double DEC_rad = DEC*RAD_DEG;
+
+  double x = Radius*cos(lambda_rad);
+  double y = Radius*cos(epsilon_rad)*sin(lambda_rad);
+  double z = Radius*sin(epsilon_rad)*sin(lambda_rad);
+
+  // original terms:
+  // *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);
+
+  // convert e->eps, s->lam, etc
+  // *pR =  +(cos(eps)*sin(lam)*cos(ra) - cos(lam)*sin(ra));
+  // *pD =  -(cos(eps)*sin(lam)*sin(ra) + cos(lam)*cos(ra))*sin(dec) + sin(eps)*sin(lam)*cos(dec);
+
+  // convert to x,y,z
+  // *pR = +(y*cos(ra) - x*sin(ra));
+  // *pD = -(y*sin(ra) + x*cos(ra))*sin(dec) + z*cos(dec);
+
+  // NOTE: seems to be identical to the old values, except I now include the varying solar distance
+  *pR = +(y*cos(RA_rad) - x*sin(RA_rad));
+  *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad);
+
   return TRUE;
 }