Initial commit. This version doesn't do much but evaluate VSOP87/Smart97 series.

1 files changed, 195 insertions, 0 deletions

diff --git a/org/madore/ephem/Comput.java b/org/madore/ephem/Comput.java
new file mode 100644
index 0000000..dee77f0
--- /dev/null
+++ b/org/madore/ephem/Comput.java
@@ -0,0 +1,195 @@
+package org.madore.ephem;
+
+import java.util.List;
+import java.util.ArrayList;
+import java.util.Arrays;
+
+public final class Comput {
+
+    public static final double ipow(double x, int n) {
+	if ( n == 0 )
+	    return 1;
+	else if ( n == 1 )
+	    return x;
+	else if ( n < 0 ) {
+	    x = 1./x;
+	    n = -n;
+	}
+	double y = x;
+	int k = 1;
+	double z = 1;
+	int m = 0;
+	while ( true ) {
+	    if ( ( n & k ) != 0 ) {
+		z *= y;
+		m += k;
+	    }
+	    if ( m < n ) {
+		y *= y;
+		k *= 2;
+	    } else
+		return z;
+	}
+    }
+
+    private static Double[] boxDoubleArray(double[] c) {
+	Double[] boxed = new Double[c.length];
+	for ( int i=0 ; i<boxed.length ; i++ )
+	    boxed[i] = c[i];
+	return boxed;
+    }
+
+    public static interface Function {
+	public double v(double t);
+    }
+
+    public static final class Constant implements Function {
+	double val;
+	public Constant(double val) {
+	    this.val = val;
+	}
+	public double v(double t) {
+	    return val;
+	}
+    }
+
+    public static final class Identity implements Function {
+	public double v(double t) {
+	    return t;
+	}
+    }
+
+    public static class SumFunction implements Function {
+	List<? extends Function> tms;
+	public SumFunction(List<? extends Function> tms) {
+	    this.tms = tms;
+	}
+	public SumFunction(Function[] tms) {
+	    this(Arrays.asList(tms));
+	}
+	public double v(double t) {
+	    double z = 0;
+	    for ( int i=0 ; i<tms.size() ; i++ )
+		z += tms.get(i).v(t);
+	    return z;
+	}	
+    }
+
+    public static class LinCombFunction implements Function {
+	List<Double> c;
+	List<Function> tms;
+	public LinCombFunction(List<Double> c, List<Function> tms) {
+	    if ( c.size() != tms.size() )
+		throw new IllegalArgumentException();
+	    this.c = c;
+	    this.tms = tms;
+	}
+	public LinCombFunction(List<Double> c, Function[] tms) {
+	    this(c, Arrays.asList(tms));
+	}
+	public LinCombFunction(Double[] c, List<Function> tms) {
+	    this(Arrays.asList(c), tms);
+	}
+	public LinCombFunction(Double[] c, Function[] tms) {
+	    this(Arrays.asList(c), Arrays.asList(tms));
+	}
+	public LinCombFunction(double[] c, List<Function> tms) {
+	    this(boxDoubleArray(c), tms);
+	}
+	public LinCombFunction(double[] c, Function[] tms) {
+	    this(boxDoubleArray(c), tms);
+	}
+	public double v(double t) {
+	    double z = 0;
+	    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 {
+	int deg;
+	Function fn;
+	public TimesPower(int deg, Function fn) {
+	    this.deg = deg;
+	    this.fn = fn;
+	}
+	public double v(double t) {
+	    return ipow(t,deg) * fn.v(t);
+	}
+    }
+
+    public static final class Polynomial implements Function {
+	List<Double> c;
+	public Polynomial(List<Double> c) {
+	    this.c = c;
+	}
+	public Polynomial(Double[] c) {
+	    this(Arrays.asList(c));
+	}
+	public Polynomial(double[] c) {
+	    this(boxDoubleArray(c));
+	}
+	public double v(double t) {
+	    double z = 0;
+	    double y = 1;
+	    for ( int i=0 ; i<c.size() ; i++ ) {
+		z += c.get(i) * y;
+		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))
+	int deg;
+	double c;
+	double phi;	
+	double om;
+	PoissonTerm(int deg, double c, double phi, double om) {
+	    this.c = c;
+	    this.deg = deg;
+	    this.phi = phi;
+	    this.om = om;
+	}
+	public static PoissonTerm cphi(int deg, double c, double phi, double om) {
+	    return new PoissonTerm(deg, c, phi, om);
+	}
+	public static PoissonTerm ab(int deg, double a, double b, double om) {
+	    double c = Math.sqrt(a*a + b*b);
+	    double phi = Math.atan2(-b, a);
+	    return new PoissonTerm(deg, c, phi, om);
+	}
+	public double v(double t) {
+	    return c * ipow(t, deg) * Math.cos(om*t + phi);
+	}
+    }
+
+    public static class SumPoissonTerms extends SumFunction
+	implements Function {
+	List<? extends PoissonTerm> tms;
+	public SumPoissonTerms(List<? extends PoissonTerm> tms) {
+	    super(tms);
+	    this.tms = tms;
+	}
+	public SumPoissonTerms(PoissonTerm[] tms) {
+	    this(Arrays.asList(tms));
+	}
+	public double v(double t) {
+	    double z = 0;
+	    final int npow = 6;
+	    List<Double> tpow = new ArrayList<Double>(npow);
+	    for ( int d=0 ; d<npow ; d++ )
+		tpow.add(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);
+		z += tm.c * tp * Math.cos(tm.om*t + tm.phi);
+	    }
+	    return z;
+	}	
+    }
+
+}