switch to full version of lea406

c/Makefile

@@ -1,9 +1,13 @@
 CC = gcc
 CFLAGS = -Wall -O2
 LIBS = -lm
-SRCS = lunarcal.c astro.c lea406.c vsop.c nutation.c julian.c lunarcalbase.c
+SRCS = lunarcal.c lunarcalbase.c lunarcalbase.h \
+	astro.c astro.h vsop.c nutation.c julian.c \
+	lea406-full.c lea406-full.h 
 OBJS = $(SRCS:.c=.o)
-TESTSRCS = testastro.c astro.c lea406.c vsop.c nutation.c julian.c lunarcalbase.c
+TESTSRCS = testastro.c \
+	astro.c astro.h vsop.c nutation.c julian.c \
+	lea406-full.c lea406-full.h 
 TESTOBJS = $(TESTSRCS:.c=.o)
 LUNARCAL = lunarcal
 TESTASTRO = testastro

c/astro.c

@@ -42,6 +42,7 @@ double rootbysecand(double (*f)(double , double),
         x1 = x2;
         //printf("debug in rootbysecand: iter = %d angle = %.9f \n", i, fx1);
     }
+    printf("debug in rootbysecand: not found after %d iterations \n", i);
     return -1;
 }
 

c/astro.h

@@ -2,6 +2,7 @@
  * header for astro functions
  */
 
+#include <stdio.h>
 #define PI        3.14159265358979323846
 #define TWOPI     6.28318530717958647693
 #define RAD2DEG  57.295779513082320876798

c/lea406-full.c

@@ -0,0 +1,60 @@
+/*
+ copyright 2014, Chen Wei <weichen302@gmail.com>
+ version 0.0.3
+Implement astronomical algorithms for finding solar terms and moon phases.
+
+Truncated LEA-406 for calculate Moon's apparent longitude;
+
+Reference:
+    LEA-406: S. M. Kudryavtsev (2007) "Long-term harmonic development of
+             lunar ephemeris", Astronomy and Astrophysics 471, 1069-1075
+*/
+
+#include <stdio.h>
+#include <math.h>
+#include "astro.h"
+#include "lea406-full.h"
+
+/*
+ * LEA-406 Moon Solution
+ *
+ * Reference:
+ *     Long-term harmonic development of lunar ephemeris.
+ *         Kudryavtsev S.M.  <Astron. Astrophys. 471, 1069 (2007)>
+ *
+ */
+
+/* compute moon ecliptic longitude using lea406 */
+double lea406(double jd, int ignorenutation) {
+    double t, tm, tm2;
+    t = (jd - J2000) / 36525.0;
+    tm = t / 10.0;
+    tm2 = tm * tm;
+
+    double V, arg;
+    V = FRM[0] + (((FRM[4] * t + FRM[3]) * t + FRM[2]) * t + FRM[1]) * t;
+
+    int i;
+    for (i = 0; i < LEA406TERMS; i++) {
+        arg = (M_ARG[i][0] + t * (M_ARG[i][1] + M_ARG[i][2] * t)) * ASEC2RAD;
+        V +=    M_AP[i][0] * sin(arg + M_AP[i][3] * DEG2RAD)
+              + M_AP[i][1] * sin(arg + M_AP[i][4] * DEG2RAD) * tm
+              + M_AP[i][2] * sin(arg + M_AP[i][5] * DEG2RAD) * tm2;
+    }
+
+    V *= ASEC2RAD;
+
+    if (!ignorenutation) {
+        V += nutation(jd);
+        /* printf("debug lea406,  nutation been adjusted"); */
+    }
+    return V;
+}
+
+/* calculate the apparent position of the Moon, it is an alias to the
+ * lea406 function
+ */
+double apparentmoon(double jd, int ignorenutation)
+{
+    return lea406(jd, ignorenutation);
+}

c/lunarcalbase.c

@@ -129,14 +129,15 @@ int get_cached_lc(struct lunarcal *p[], int year)
     }
 
     len = mark_month_day(p);
-    for (i = 0; i < len; i++) {
-        if (cachep > CACHESIZE)
-            cachep = 0;
-        cached_lcs[cachep][i] = p[i];
-    }
 
     /* add to cache */
+    if (cachep >= CACHESIZE)
+        cachep = 0;
+    for (i = 0; i < len; i++) {
+        cached_lcs[cachep][i] = p[i];
+    }
     cached_year[cachep++] = year;
+
     return len;
 }
 

c/lunarcalbase.h

@@ -1,7 +1,7 @@
 #define MAX_SOLARTERMS 27
 #define MAX_NEWMOONS 15
 #define MAX_DAYS 450
-#define CACHESIZE 5
+#define CACHESIZE 3
 
 struct lunarcal {
     double jd;
@@ -26,5 +26,3 @@ struct lunarcal *lcalloc(double jd);
 void print_lunarcal(struct lunarcal *p[], int len);
 
 int get_cache_index(int year);
-
-

c/testastro.c

@@ -64,14 +64,16 @@ void testnewmoon_solarterm(void)
     int n;
     char isodt[30];
     int i;
-    for (n = 0; n < 5; n++) {
+    for (n = 0; n < 50; n++) {
         jd = g2jd(year, 1, 1.0);
         findnewmoons(newmoons, NMCOUNT, jd);
         year += 1;
-    }
-    for (i = 0; i < NMCOUNT; i++) {
-        jdftime(isodt, newmoons[i], "%y-%m-%d %H:%M:%S", 8.0, 1);
-        printf("found newmoon: %s %.8f\n", isodt, newmoons[i]);
+
+        for (i = 0; i < NMCOUNT; i++) {
+            jdftime(isodt, newmoons[i], "%y-%m-%d %H:%M:%S", 8.0, 1);
+            printf("found newmoon: %s %.8f\n", isodt, newmoons[i]);
+        }
+
     }
 
     double angle;
@@ -168,7 +170,7 @@ void verify_apparent_sun_moon(void)
 
     lensun  = parsejplhorizon("jpl_sun.txt",  jplsun);
     lenmoon = parsejplhorizon("jpl_moon.txt", jplmoon);
-    step = 1;
+    step = 2;
     i = 0;
     count = 0;
     delta_sun_n = 0;