1 files changed, 73 insertions, 0 deletions
diff --git a/org/madore/ephem/Nutation.java b/org/madore/ephem/Nutation.java
new file mode 100644
index 0000000..ecc08a4
--- /dev/null
+++ b/org/madore/ephem/Nutation.java
@@ -0,0 +1,73 @@
+package org.madore.ephem;
+
+import java.util.List;
+import java.util.ArrayList;
+import java.util.EnumMap;
+import java.io.InputStream;
+import java.io.InputStreamReader;
+import java.io.BufferedReader;
+
+public final class Nutation {
+
+    public static enum Variable {
+	PSI("psi"), EPSILON("epsilon");
+	final String name;
+	Variable(String name) { this.name = name; }
+    }
+
+    // The argument is in Julian centuries (value is in radians).
+    public static Comput.Polynomial[] lambdas = new Comput.Polynomial[] {
+	new Comput.Polynomial(134.96340251*Comput.degree, 1717915923.2178*Comput.arcsecond, 31.8792*Comput.arcsecond +0.051635*Comput.arcsecond, -0.00024470*Comput.arcsecond), // Delaunay l
+	new Comput.Polynomial(357.52910918*Comput.degree, 129596581.0481*Comput.arcsecond, -0.5532*Comput.arcsecond +0.000136*Comput.arcsecond, -0.00001149*Comput.arcsecond), // Delaunay l′
+	new Comput.Polynomial(93.27209062*Comput.degree, 1739527262.8478*Comput.arcsecond, -12.7512*Comput.arcsecond -0.001037*Comput.arcsecond, 0.00000417*Comput.arcsecond), // Delaunay F
+	new Comput.Polynomial(297.85019547*Comput.degree, 1602961601.2090*Comput.arcsecond, -6.3706*Comput.arcsecond +0.006593*Comput.arcsecond, -0.00003169*Comput.arcsecond), // Delaunay D
+	new Comput.Polynomial(125.04455501*Comput.degree, -6962890.5431*Comput.arcsecond, 7.4722*Comput.arcsecond +0.007702*Comput.arcsecond, -0.00005939*Comput.arcsecond), // Delaunay Ω
+	new Comput.Polynomial(4.402608842, 2608.7903141574), // L(Mercury)
+	new Comput.Polynomial(3.176146697, 1021.3285546211), // L(Venus)
+	new Comput.Polynomial(1.753470314, 628.3075849991), // L(Earth)
+	new Comput.Polynomial(6.203480913, 334.0612426700), // L(Mars)
+	new Comput.Polynomial(0.599546497, 52.9690962641), // L(Jupiter)
+	new Comput.Polynomial(0.874016757, 21.3299104960), // L(Saturn)
+	new Comput.Polynomial(5.481293872, 7.4781598567), // L(Uranus)
+	new Comput.Polynomial(5.311886287, 3.8133035638), // L(Neptune)
+	new Comput.Polynomial(0, 0.02438175, 0.00000538691) // General precession
+    };
+
+    private static EnumMap<Variable, Comput.SumPoisson2Terms> data = new EnumMap<Variable, Comput.SumPoisson2Terms>(Variable.class);
+
+    public static Comput.SumPoisson2Terms getFunc(Variable v) {
+	if ( ! data.containsKey(v) ) {
+	    List<Comput.Poisson2Term> series = new ArrayList<Comput.Poisson2Term>();
+	    try {
+		InputStream str = Nutation.class.getResourceAsStream("nutation.dat");
+		BufferedReader in = new BufferedReader(new InputStreamReader(str, "utf-8"));
+		String s;
+		while ( ( s = in.readLine() ) != null ) {
+		    String[] fields = s.split("\t");
+		    if ( ! fields[1].equals(v.name) )
+			continue;
+		    int deg = Integer.parseInt(fields[2]);
+		    String[] subfields = fields[3].split(",");
+		    double[] om = new double[14];
+		    for ( int j=0 ; j<14 ; j++ )
+			om[j] = Double.parseDouble(subfields[j]);
+		    double b = Double.parseDouble(fields[4]);
+		    double a = Double.parseDouble(fields[5]);
+		    series.add(Comput.Poisson2Term.ab(deg, a, b, om, lambdas));
+		}
+	    } catch ( Exception e ) {
+		throw new RuntimeException(e);
+	    }
+	    data.put(v, new Comput.SumPoisson2Terms(series));
+	}
+	return data.get(v);
+    }
+
+    public static Frames.RotationMatrix matrix(double t) {
+	Frames.RotationMatrix mat1 = Frames.RotationMatrix.rotx(Precession.epsilon.v(t)*Comput.arcsecond);
+	Frames.RotationMatrix mat2 = Frames.RotationMatrix.rotz(-getFunc(Variable.PSI).v(t)*Comput.arcsecond).apply(mat1);
+	Frames.RotationMatrix mat3 = Frames.RotationMatrix.rotx(-(Precession.epsilon.v(t)+getFunc(Variable.EPSILON).v(t))*Comput.arcsecond).apply(mat2);
+	return mat3;
+    }
+
+}