continued with a3 calc while on the way to petrozavodsk

tools/tascalc/TasCalc.java

@@ -10,6 +10,10 @@ public class TasCalc
     // calculated with scipy, see tascalc.py
     protected static final double E_to_k2 = 0.482596406464;
 
+    //protected static final double a3_offs = 0.;             // for sics: a3 is angle between ki and orient1
+    //protected static final double a3_offs = np.pi/2.;       // for takin: Q along orient1 => a3:=a4/2
+    protected static final double a3_offs = Math.PI;          // for nomad: a3 is angle between ki and orient1 + 180 deg
+
 
     public TasCalc()
     {
@@ -52,7 +56,7 @@ public class TasCalc
             }
         };
 
-        return A;
+        return Calc.transpose(A);
     }
 
 
@@ -186,5 +190,63 @@ public class TasCalc
     {
         return (ki*ki - kf*kf) / E_to_k2;
     }
+
+
+    /**
+     * calculate a3 & a4 angles
+     */
+    public static double[] get_a3a4(double ki, double kf, double[] Q_rlu, 
+        double[] orient_rlu, double[] orient_up_rlu, double[][] B, double sense_sample)
+        throws Exception
+    {
+        // metric
+        double[][] metric = Calc.get_metric(B);
+
+        // length of Q
+        double Qlen = Math.sqrt(Calc.dot(Q_rlu, Q_rlu, metric));
+
+        // angle xi between Q and orientation reflex
+        double xi = Calc.angle(Q_rlu, orient_rlu, metric);
+
+        // sign of xi
+        if(Calc.dot(Calc.cross(orient_rlu, Q_rlu, B), orient_up_rlu, metric) < 0.)
+            xi = -xi;
+
+        // distance to plane
+        double up_len = Math.sqrt(Calc.dot(orient_up_rlu, orient_up_rlu, metric));
+        double dist_Q_plane = Calc.dot(Q_rlu, orient_up_rlu, metric) / up_len;
+
+        // angle psi enclosed by ki and Q
+        double psi = get_psi(ki, kf, Qlen, sense_sample);
+
+        double a3 = - psi - xi + a3_offs;
+        double a4 = get_a4(ki, kf, Qlen);
+
+        //System.out.println("xi = " + xi/Math.PI*180. + ", psi = " + psi/Math.PI*180. + ", offs = " + a3_offs);
+        return new double[]{a3, a4, dist_Q_plane};
+    }
+
+
+    /**
+     * get Q position
+     */
+    public static double[] get_hkl(double ki, double kf, double a3, double Qlen, 
+        double[] orient_rlu, double[] orient_up_rlu, double[][] B, double sense_sample)
+        throws Exception
+    {
+        double[][] B_inv = Calc.inv(B);
+
+        // angle enclosed by ki and Q
+        double psi = TasCalc.get_psi(ki, kf, Qlen, sense_sample);
+
+        // angle between Q and orientation reflex
+        double xi = - a3 + a3_offs - psi;
+
+        double[] Q_lab = Calc.rotate(Calc.dot(B, orient_up_rlu), Calc.dot(B, Calc.mul(orient_rlu, Qlen)), xi);
+        Q_lab = Calc.mul(Q_lab, Qlen / Calc.norm_2(Q_lab));
+        double[] Q_rlu = Calc.dot(B_inv, Q_lab);
+
+        return Q_rlu;
+    }
     // ------------------------------------------------------------------------
 }

tools/tascalc/TestTasCalc.java

@@ -14,12 +14,8 @@ import org.junit.Assert;
 
 public class TestTasCalc extends TestCase
 {
-    private TasCalc m_tas;
-
-
     public void setUp()
     {   
-        m_tas = new TasCalc();
     }
 
 
@@ -154,4 +150,32 @@ public class TestTasCalc extends TestCase
         double E_2 = TasCalc.get_E(ki, kf);
         assertEquals("Wrong E!", E, E_2, 1e-3);
     }
+
+
+    public void testAngles()
+        throws Exception
+    {
+        double kf = 2.662;
+        double E = 2.;
+        double ki = TasCalc.get_ki(kf, E);
+        double[] Q_rlu = new double[]{1., 2., 2.};
+        double[] orient_rlu = new double[]{1., 0., 0.};
+        double[] orient_up_rlu = new double[]{-1./3., -2./3., 2./3.};
+
+        double[] lattice = new double[]{5., 5., 5.};
+        double[] angles = new double[]{90./180.*Math.PI, 90./180.*Math.PI, 60./180.*Math.PI};
+
+        double[][] B = TasCalc.get_B(lattice, angles);
+        double[] a3a4 = TasCalc.get_a3a4(ki, kf, Q_rlu, orient_rlu, orient_up_rlu, B, 1.);
+
+        assertEquals("Wrong distance to scattering plane!", 0., a3a4[2], 1e-3);
+        assertEquals("Wrong a4!", 80.457, a3a4[1]/Math.PI*180., 1e-3);
+        assertEquals("Wrong a3!", 42.389, a3a4[0]/Math.PI*180., 1e-3);
+
+
+        double[][] metric = Calc.get_metric(B);
+        double Qlen = Math.sqrt(Calc.dot(Q_rlu, Q_rlu, metric));
+        double[] Qhkl = TasCalc.get_hkl(ki, kf, a3a4[0], Qlen, orient_rlu, orient_up_rlu, B, 1.);
+        Assert.assertArrayEquals("Wrong Q position!", new double[]{1., 2., 2.}, Qhkl, 1e-4);
+    }
 }

tools/tascalc/tascalc.py

@@ -187,7 +187,7 @@ def get_a3a4(ki, kf, Q_rlu, orient_rlu, orient_up_rlu, B, sense_sample=1.):
 	a3 = - psi - xi + a3_offs
 	a4 = get_a4(ki, kf, Qlen)
 
-	#print("xi = " + str(xi/np.pi*180.) + ", psi = " + str(psi/np.pi*180.))
+	#print("xi = " + str(xi/np.pi*180.) + ", psi = " + str(psi/np.pi*180.) + ", offs = " + str(a3_offs/np.pi*180.))
 	return [a3, a4, dist_Q_plane]