Changeset 71

Timestamp:

May 28, 2010 6:15:43 PM (13 years ago)

Author:

toby

Message:

Update orthogonalization routine w/unit-tests

Location:

trunk

Files:

• GSASIIlattice.py (modified) (6 diffs)
• testinp/gencelltests.py (modified) (2 diffs)
• unit_tests.py (modified) (2 diffs)

trunk/GSASIIlattice.py

@@ +1 @@
+'''Perform lattice-related computations'''
+
 import numpy as np
 import numpy.linalg as nl
@@ -107 +109 @@
     return ainv,binv,cinv,alpinv,betinv,gaminv
 
-# not tested yet (broken?)
 def cell2AB(cell):
     '''Computes orthogonalization matrix from unit cell constants
@@ -115 +116 @@
        B (inverse) for crystal to Cartesian transformation B*x = X
     '''
+    G,g = cell2Gmat(cell) 
+    cellstar = Gmat2cell(G)
+    A = np.zeros(shape=(3,3))
+    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
+    A[0][0] = cell[0]                # a
+    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
+    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
+    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
+    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
+    A[2][2] = 1/cellstar[2]         # 1/c*
+    B = nl.inv(A)
+    return A,B
+    # bob's code:
     G,g = cell2Gmat(cell)       #reciprocal & real metric tensors
     cosAlpStar = G[2][1]/np.sqrt(G[1][1]*G[2][2])
@@ -157 +171 @@
        [  3.36690552e-18,   3.36690552e-18,   2.77777778e-02]]), (0.32991443953692895, 0.32991443953692895, 0.16666666666666669, 90.0, 90.0, 60.000000000000021), 63.652867178156257, 0.015710211406520427],
 ]
+CoordTestData = [
+# cell, ((frac, ortho),...)
+  ((4,4,4,90,90,90,), [
+ ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
+ ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
+ ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
+ ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
+ ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
+ ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
+ ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
+]),
+# cell, ((frac, ortho),...)
+  ((4.1,5.2,6.3,100,80,130,), [
+ ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
+ ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
+ ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
+ ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
+ ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
+ ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
+ ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
+]),
+# cell, ((frac, ortho),...)
+  ((3.5,3.5,6,90,90,120,), [
+ ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
+ ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
+ ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
+ ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
+ ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
+ ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
+ ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
+]),
+]
 
 def test0():
@@ -202 +248 @@
         assert np.allclose(rcell,trcell),msg
 
+def test6():
+    msg = 'test cell2AB'
+    for (cell,coordlist) in CoordTestData:
+        A,B = cell2AB(cell)
+        for (frac,ortho) in coordlist:
+            to = np.inner(A,frac)
+            tf = np.inner(B,to)
+            assert np.allclose(ortho,to), msg
+            assert np.allclose(frac,tf), msg
+
 if __name__ == '__main__':
     test0()
@@ -209 +265 @@
     test4()
     test5()
+    test6()
     print "OK"

trunk/testinp/gencelltests.py

@@ -6 +6 @@
 
 fp = sys.stdout
-fp.write("# output from uctbx computed on platform %s on %s\n" %
+fp.write("# output from uctbx (cctbx) computed on platform %s on %s\n" %
          (sys.platform, datetime.date.today()) )
-fp.write("import numpy as np\n")
+fp.write("#import numpy as np\n")
 fp.write("array = np.array\n\n")
 fp.write("CellTestData = [\n")
@@ -40 +40 @@
 fp.write("]\n")
 
+fp.write("CoordTestData = [\n")
+for cell in [
+    (4,4,4,90,90,90),
+    (4.1,5.2,6.3,100,80,130),
+    (3.5,3.5,6,90,90,120),
+    ]:
+    fp.write("# cell, ((frac, ortho),...)\n")
+    fp.write("  ((%s,%s,%s,%s,%s,%s,), [\n" % cell)
+    for frac in [
+        (0.1,0.,0.),
+        (0.,0.1,0.),
+        (0.,0.,0.1),
+        (0.1,0.2,0.3),
+        (0.2,0.3,0.1),
+        (0.3,0.2,0.1),
+        (0.5,0.5,0.5),
+        ]:
+        fp.write(" (%s,%s),\n" % (frac, uc.unit_cell(cell).orthogonalize(frac)))
+    fp.write("]),\n")
+fp.write("]\n")
+    
+sys.exit()

trunk/unit_tests.py

@@ -21 +21 @@
 def test_GSASIIlattice5():
     G2l.test5()
-
-
-
+def test_GSASIIlattice6():
+    G2l.test6()
 
 if __name__ == '__main__':
@@ -36 +35 @@
     G2l.test4()
     G2l.test5()
+    G2l.test6()
     print "OK"