Changeset 71


Ignore:
Timestamp:
May 28, 2010 6:15:43 PM (12 years ago)
Author:
toby
Message:

Update orthogonalization routine w/unit-tests

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIlattice.py

    r70 r71  
     1'''Perform lattice-related computations'''
     2
    13import numpy as np
    24import numpy.linalg as nl
     
    107109    return ainv,binv,cinv,alpinv,betinv,gaminv
    108110
    109 # not tested yet (broken?)
    110111def cell2AB(cell):
    111112    '''Computes orthogonalization matrix from unit cell constants
     
    115116       B (inverse) for crystal to Cartesian transformation B*x = X
    116117    '''
     118    G,g = cell2Gmat(cell)
     119    cellstar = Gmat2cell(G)
     120    A = np.zeros(shape=(3,3))
     121    # from Giacovazzo (Fundamentals 2nd Ed.) p.75
     122    A[0][0] = cell[0]                # a
     123    A[0][1] = cell[1]*cosd(cell[5])  # b cos(gamma)
     124    A[0][2] = cell[2]*cosd(cell[4])  # c cos(beta)
     125    A[1][1] = cell[1]*sind(cell[5])  # b sin(gamma)
     126    A[1][2] = -cell[2]*cosd(cellstar[3])*sind(cell[4]) # - c cos(alpha*) sin(beta)
     127    A[2][2] = 1/cellstar[2]         # 1/c*
     128    B = nl.inv(A)
     129    return A,B
     130    # bob's code:
    117131    G,g = cell2Gmat(cell)       #reciprocal & real metric tensors
    118132    cosAlpStar = G[2][1]/np.sqrt(G[1][1]*G[2][2])
     
    157171       [  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],
    158172]
     173CoordTestData = [
     174# cell, ((frac, ortho),...)
     175  ((4,4,4,90,90,90,), [
     176 ((0.10000000000000001, 0.0, 0.0),(0.40000000000000002, 0.0, 0.0)),
     177 ((0.0, 0.10000000000000001, 0.0),(2.4492935982947065e-17, 0.40000000000000002, 0.0)),
     178 ((0.0, 0.0, 0.10000000000000001),(2.4492935982947065e-17, -2.4492935982947065e-17, 0.40000000000000002)),
     179 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.40000000000000013, 0.79999999999999993, 1.2)),
     180 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.80000000000000016, 1.2, 0.40000000000000002)),
     181 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(1.2, 0.80000000000000004, 0.40000000000000002)),
     182 ((0.5, 0.5, 0.5),(2.0, 1.9999999999999998, 2.0)),
     183]),
     184# cell, ((frac, ortho),...)
     185  ((4.1,5.2,6.3,100,80,130,), [
     186 ((0.10000000000000001, 0.0, 0.0),(0.40999999999999998, 0.0, 0.0)),
     187 ((0.0, 0.10000000000000001, 0.0),(-0.33424955703700043, 0.39834311042186865, 0.0)),
     188 ((0.0, 0.0, 0.10000000000000001),(0.10939835193016617, -0.051013289294572106, 0.6183281045774256)),
     189 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(0.069695941716497567, 0.64364635296002093, 1.8549843137322766)),
     190 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(-0.073350319180835066, 1.1440160419710339, 0.6183281045774256)),
     191 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.67089923785616512, 0.74567293154916525, 0.6183281045774256)),
     192 ((0.5, 0.5, 0.5),(0.92574397446582857, 1.7366491056364828, 3.0916405228871278)),
     193]),
     194# cell, ((frac, ortho),...)
     195  ((3.5,3.5,6,90,90,120,), [
     196 ((0.10000000000000001, 0.0, 0.0),(0.35000000000000003, 0.0, 0.0)),
     197 ((0.0, 0.10000000000000001, 0.0),(-0.17499999999999993, 0.3031088913245536, 0.0)),
     198 ((0.0, 0.0, 0.10000000000000001),(3.6739403974420595e-17, -3.6739403974420595e-17, 0.60000000000000009)),
     199 ((0.10000000000000001, 0.20000000000000001, 0.29999999999999999),(2.7675166561703527e-16, 0.60621778264910708, 1.7999999999999998)),
     200 ((0.20000000000000001, 0.29999999999999999, 0.10000000000000001),(0.17500000000000041, 0.90932667397366063, 0.60000000000000009)),
     201 ((0.29999999999999999, 0.20000000000000001, 0.10000000000000001),(0.70000000000000018, 0.6062177826491072, 0.60000000000000009)),
     202 ((0.5, 0.5, 0.5),(0.87500000000000067, 1.5155444566227676, 3.0)),
     203]),
     204]
    159205
    160206def test0():
     
    202248        assert np.allclose(rcell,trcell),msg
    203249
     250def test6():
     251    msg = 'test cell2AB'
     252    for (cell,coordlist) in CoordTestData:
     253        A,B = cell2AB(cell)
     254        for (frac,ortho) in coordlist:
     255            to = np.inner(A,frac)
     256            tf = np.inner(B,to)
     257            assert np.allclose(ortho,to), msg
     258            assert np.allclose(frac,tf), msg
     259
    204260if __name__ == '__main__':
    205261    test0()
     
    209265    test4()
    210266    test5()
     267    test6()
    211268    print "OK"
  • trunk/testinp/gencelltests.py

    r70 r71  
    66
    77fp = sys.stdout
    8 fp.write("# output from uctbx computed on platform %s on %s\n" %
     8fp.write("# output from uctbx (cctbx) computed on platform %s on %s\n" %
    99         (sys.platform, datetime.date.today()) )
    10 fp.write("import numpy as np\n")
     10fp.write("#import numpy as np\n")
    1111fp.write("array = np.array\n\n")
    1212fp.write("CellTestData = [\n")
     
    4040fp.write("]\n")
    4141
     42fp.write("CoordTestData = [\n")
     43for cell in [
     44    (4,4,4,90,90,90),
     45    (4.1,5.2,6.3,100,80,130),
     46    (3.5,3.5,6,90,90,120),
     47    ]:
     48    fp.write("# cell, ((frac, ortho),...)\n")
     49    fp.write("  ((%s,%s,%s,%s,%s,%s,), [\n" % cell)
     50    for frac in [
     51        (0.1,0.,0.),
     52        (0.,0.1,0.),
     53        (0.,0.,0.1),
     54        (0.1,0.2,0.3),
     55        (0.2,0.3,0.1),
     56        (0.3,0.2,0.1),
     57        (0.5,0.5,0.5),
     58        ]:
     59        fp.write(" (%s,%s),\n" % (frac, uc.unit_cell(cell).orthogonalize(frac)))
     60    fp.write("]),\n")
     61fp.write("]\n")
     62   
     63sys.exit()
  • trunk/unit_tests.py

    r70 r71  
    2121def test_GSASIIlattice5():
    2222    G2l.test5()
    23 
    24 
    25 
     23def test_GSASIIlattice6():
     24    G2l.test6()
    2625
    2726if __name__ == '__main__':
     
    3635    G2l.test4()
    3736    G2l.test5()
     37    G2l.test6()
    3838    print "OK"
Note: See TracChangeset for help on using the changeset viewer.