Implement more general Poisson series.

1 files changed, 77 insertions, 14 deletions

diff --git a/org/madore/ephem/Comput.java b/org/madore/ephem/Comput.java
index 214d5e8..77e305b 100644
--- a/org/madore/ephem/Comput.java
+++ b/org/madore/ephem/Comput.java
@@ -1,7 +1,7 @@
 package org.madore.ephem;
 
 import java.util.List;
-import java.util.ArrayList;
+import java.util.HashMap;
 import java.util.Arrays;
 
 public final class Comput {
@@ -75,7 +75,7 @@ public final class Comput {
 	    for ( int i=0 ; i<tms.size() ; i++ )
 		z += tms.get(i).v(t);
 	    return z;
-	}	
+	}
     }
 
     public static class LinCombFunction implements Function {
@@ -107,7 +107,7 @@ public final class Comput {
 	    for ( int i=0 ; i<tms.size() ; i++ )
 		z += c.get(i) * tms.get(i).v(t);
 	    return z;
-	}	
+	}
     }
 
     public static final class TimesPower implements Function {
@@ -141,20 +141,83 @@ public final class Comput {
 		y *= t;
 	    }
 	    return z;
-	}	
+	}
     }
 
-    public static final class PoissonTerm implements Function {
-	// Returns c*t^d*cos(phi+om*t) = t^d*(a*cos(om*t) + b*sin(om*t))
+    public static class Poisson2Term implements Function {
+	// Returns c*t^d*cos(phi+sum(om[i]*fn[i](t)))
 	int deg;
 	double c;
-	double phi;	
-	double om;
-	PoissonTerm(int deg, double c, double phi, double om) {
+	double phi;
+	double[] om;
+	Function[] fn;
+	Poisson2Term(int deg, double c, double phi, double[] om, Function[] fn) {
+	    if ( om.length != fn.length )
+		throw new IllegalArgumentException();
 	    this.c = c;
 	    this.deg = deg;
 	    this.phi = phi;
 	    this.om = om;
+	    this.fn = fn;
+	}
+	public static Poisson2Term cphi(int deg, double c, double phi, double[] om, Function[] fn) {
+	    return new Poisson2Term(deg, c, phi, om, fn);
+	}
+	public static Poisson2Term ab(int deg, double a, double b, double om[], Function[] fn) {
+	    double c = Math.sqrt(a*a + b*b);
+	    double phi = Math.atan2(-b, a);
+	    return new Poisson2Term(deg, c, phi, om, fn);
+	}
+	public double v(double t) {
+	    double arg = 0;
+	    for ( int i=0 ; i<om.length ; i++ )
+		arg += om[i]*fn[i].v(t);
+	    return c * ipow(t, deg) * Math.cos(arg + phi);
+	}
+    }
+
+    public static class SumPoisson2Terms extends SumFunction
+	implements Function {
+	List<? extends Poisson2Term> tms;
+	public SumPoisson2Terms(List<? extends Poisson2Term> tms) {
+	    super(tms);
+	    this.tms = tms;
+	}
+	public SumPoisson2Terms(Poisson2Term[] tms) {
+	    this(Arrays.asList(tms));
+	}
+	public double v(double t) {
+	    double z = 0;
+	    final int npow = 6;
+	    double[] tpow = new double[npow];
+	    for ( int d=0 ; d<npow ; d++ )
+		tpow[d] = ipow(t,d);
+	    HashMap<Function,Double> fncache = new HashMap<Function,Double>();
+	    for ( int i=0 ; i<tms.size() ; i++ ) {
+		Poisson2Term tm = tms.get(i);
+		int deg = tm.deg;
+		double tp = deg<npow ? tpow[deg] : ipow(t,deg);
+		double arg = 0;
+		for ( int j=0 ; j<tm.om.length ; j++ ) {
+		    if ( ! fncache.containsKey(tm.fn[j]) ) {
+			double v = tm.fn[j].v(t);
+			fncache.put(tm.fn[j], v);
+		    }
+		    arg += tm.om[j]*fncache.get(tm.fn[j]);
+		}
+		z += tm.c * tp * Math.cos(arg + tm.phi);
+	    }
+	    return z;
+	}
+    }
+
+    public static final class PoissonTerm extends Poisson2Term implements Function {
+	// Returns c*t^d*cos(phi+om*t) = t^d*(a*cos(om*t) + b*sin(om*t))
+	double om;
+	private static Function identity = new Identity();
+	PoissonTerm(int deg, double c, double phi, double om) {
+	    super(deg, c, phi, new double[] { om }, new Function[] { identity });
+	    this.om = om;
 	}
 	public static PoissonTerm cphi(int deg, double c, double phi, double om) {
 	    return new PoissonTerm(deg, c, phi, om);
@@ -169,7 +232,7 @@ public final class Comput {
 	}
     }
 
-    public static class SumPoissonTerms extends SumFunction
+    public static class SumPoissonTerms extends SumPoisson2Terms
 	implements Function {
 	List<? extends PoissonTerm> tms;
 	public SumPoissonTerms(List<? extends PoissonTerm> tms) {
@@ -182,17 +245,17 @@ public final class Comput {
 	public double v(double t) {
 	    double z = 0;
 	    final int npow = 6;
-	    List<Double> tpow = new ArrayList<Double>(npow);
+	    double[] tpow = new double[npow];
 	    for ( int d=0 ; d<npow ; d++ )
-		tpow.add(ipow(t,d));
+		tpow[d] = ipow(t,d);
 	    for ( int i=0 ; i<tms.size() ; i++ ) {
 		PoissonTerm tm = tms.get(i);
 		int deg = tm.deg;
-		double tp = deg<npow ? tpow.get(deg) : ipow(t,deg);
+		double tp = deg<npow ? tpow[deg] : ipow(t,deg);
 		z += tm.c * tp * Math.cos(tm.om*t + tm.phi);
 	    }
 	    return z;
-	}	
+	}
     }
 
 }