diff options
Diffstat (limited to 'org/madore')
-rw-r--r-- | org/madore/ephem/Frames.java | 82 | ||||
-rw-r--r-- | org/madore/ephem/Nutation.java | 8 | ||||
-rw-r--r-- | org/madore/ephem/Precession.java | 10 |
3 files changed, 46 insertions, 54 deletions
diff --git a/org/madore/ephem/Frames.java b/org/madore/ephem/Frames.java index e85412c..a8a4667 100644 --- a/org/madore/ephem/Frames.java +++ b/org/madore/ephem/Frames.java @@ -10,75 +10,67 @@ public final class Frames { } public static final class Vector { - double[] crd; - public Vector(double... crd) { - if ( crd.length != 3 ) + final double[] v; // v.length == 3 + public Vector(double... v) { + if ( v.length != 3 ) throw new IllegalArgumentException("Vector constructor expects 3 coordinates"); - this.crd = crd; + this.v = v; } - public Vector(Double... crd) { - this(unboxDoubleArray(crd)); + public Vector(Double... v) { + this(unboxDoubleArray(v)); } - public double dotprod(Vector v) { + public double dotprod(Vector w) { double res = 0; for ( int i=0 ; i<3 ; i++ ) - res += crd[i]*v.crd[i]; + res += v[i]*w.v[i]; return res; } public double sqnorm() { double res = 0; for ( int i=0 ; i<3 ; i++ ) - res += crd[i]*crd[i]; + res += v[i]*v[i]; return res; } } - public static final class RotationMatrix { - double[][] mat; - public RotationMatrix(double[][] mat) { - if ( mat.length != 3 ) - throw new IllegalArgumentException("RotationMatrix constructor expects 3 lines"); - for ( int i=0 ; i<3 ; i++ ) - if ( mat[i].length != 3 ) - throw new IllegalArgumentException("RotationMatrix constructor expects 3 columns"); - this.mat = mat; - } - private static double[][] matOfVectors(Vector[] vecs) { - double[][] mat = new double[vecs.length][]; - for ( int i=0 ; i<mat.length ; i++ ) - mat[i] = vecs[i].crd; - return mat; + public static final class Rotation { + final double[] q; // q.length == 4 + public Rotation(double... q) { + if ( q.length != 4 ) + throw new IllegalArgumentException("Rotation constructor expects 4 coordinates"); + this.q = q; } - public RotationMatrix(Vector... vecs) { - this(matOfVectors(vecs)); + public Rotation(Double... q) { + this(unboxDoubleArray(q)); } public Vector apply(double[] v) { - double out[] = new double[3]; - for ( int i=0 ; i<3 ; i++ ) - for ( int j=0 ; j<3 ; j++ ) - out[i] += v[j]*mat[j][i]; + if ( v.length != 3 ) + throw new IllegalArgumentException("Rotation.apply() expects 3 coordinates"); + double[] out = new double[3]; + out[0] = (q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3])*v[0] + (2*q[1]*q[2] + 2*q[0]*q[3])*v[1] + (-2*q[0]*q[2] + 2*q[1]*q[3])*v[2]; + out[1] = (2*q[1]*q[2] - 2*q[0]*q[3])*v[0] + (q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3])*v[1] + (2*q[0]*q[1] + 2*q[2]*q[3])*v[2]; + out[2] = (2*q[0]*q[2] + 2*q[1]*q[3])*v[0] + (-2*q[0]*q[1] + 2*q[2]*q[3])*v[1] + (q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3])*v[2]; return new Vector(out); } public Vector apply(Vector v) { - return apply(v.crd); + return apply(v.v); } - public RotationMatrix apply(RotationMatrix mv) { - return new RotationMatrix(this.apply(mv.mat[0]), this.apply(mv.mat[1]), this.apply(mv.mat[2])); + public Rotation apply(Rotation h) { + double[] outq = new double[4]; + outq[0] = q[0]*h.q[0] - q[1]*h.q[1] - q[2]*h.q[2] - q[3]*h.q[3]; + outq[1] = q[1]*h.q[0] + q[0]*h.q[1] + q[3]*h.q[2] - q[2]*h.q[3]; + outq[2] = q[2]*h.q[0] - q[3]*h.q[1] + q[0]*h.q[2] + q[1]*h.q[3]; + outq[3] = q[3]*h.q[0] + q[2]*h.q[1] - q[1]*h.q[2] + q[0]*h.q[3]; + return new Rotation(outq); } - public static RotationMatrix rotx(double theta) { - return new RotationMatrix(new Vector(1, 0, 0), - new Vector(0, Math.cos(theta), -Math.sin(theta)), - new Vector(0, Math.sin(theta), Math.cos(theta))); + public static Rotation rotx(double theta) { + return new Rotation(Math.cos(theta/2), Math.sin(theta/2), 0, 0); } - public static RotationMatrix roty(double theta) { - return new RotationMatrix(new Vector(Math.cos(theta), 0, Math.sin(theta)), - new Vector(0, 1, 0), - new Vector(-Math.sin(theta), 0, Math.cos(theta))); + public static Rotation roty(double theta) { + return new Rotation(Math.cos(theta/2), 0, Math.sin(theta/2), 0); } - public static RotationMatrix rotz(double theta) { - return new RotationMatrix(new Vector(Math.cos(theta), -Math.sin(theta), 0), - new Vector(Math.sin(theta), Math.cos(theta), 0), - new Vector(0, 0, 1)); + public static Rotation rotz(double theta) { + return new Rotation(Math.cos(theta/2), 0, 0, Math.sin(theta/2)); } } diff --git a/org/madore/ephem/Nutation.java b/org/madore/ephem/Nutation.java index dcfdc30..bc6e721 100644 --- a/org/madore/ephem/Nutation.java +++ b/org/madore/ephem/Nutation.java @@ -63,10 +63,10 @@ public final class Nutation { 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); + public static Frames.Rotation matrix(double t) { + Frames.Rotation mat1 = Frames.Rotation.rotx(Precession.epsilon.v(t)*Comput.arcsecond); + Frames.Rotation mat2 = Frames.Rotation.rotz(-getFunc(Variable.PSI).v(t)*Comput.arcsecond).apply(mat1); + Frames.Rotation mat3 = Frames.Rotation.rotx(-(Precession.epsilon.v(t)+getFunc(Variable.EPSILON).v(t))*Comput.arcsecond).apply(mat2); return mat3; } diff --git a/org/madore/ephem/Precession.java b/org/madore/ephem/Precession.java index 8b556d5..af8a1dc 100644 --- a/org/madore/ephem/Precession.java +++ b/org/madore/ephem/Precession.java @@ -15,11 +15,11 @@ public final class Precession { public static Comput.Polynomial eqorg = new Comput.Polynomial(0.014506, 4612.156534, 1.3915817, -0.00000044, -0.000029956, -0.0000000368); - public static Frames.RotationMatrix matrix(double t) { - Frames.RotationMatrix mat1 = Frames.RotationMatrix.rotx(epsilon0*Comput.arcsecond); - Frames.RotationMatrix mat2 = Frames.RotationMatrix.rotz(-psi.v(t)*Comput.arcsecond).apply(mat1); - Frames.RotationMatrix mat3 = Frames.RotationMatrix.rotx(-omega.v(t)*Comput.arcsecond).apply(mat2); - Frames.RotationMatrix mat4 = Frames.RotationMatrix.rotz(chi.v(t)*Comput.arcsecond).apply(mat3); + public static Frames.Rotation matrix(double t) { + Frames.Rotation mat1 = Frames.Rotation.rotx(epsilon0*Comput.arcsecond); + Frames.Rotation mat2 = Frames.Rotation.rotz(-psi.v(t)*Comput.arcsecond).apply(mat1); + Frames.Rotation mat3 = Frames.Rotation.rotx(-omega.v(t)*Comput.arcsecond).apply(mat2); + Frames.Rotation mat4 = Frames.Rotation.rotz(chi.v(t)*Comput.arcsecond).apply(mat3); return mat4; } |