Changeset 74


Ignore:
Timestamp:
Jun 2, 2010 12:27:50 PM (12 years ago)
Author:
vondreel
Message:

rearrange routines to be in blocks for indexing & image processing - these can now be put into separate modules

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIcomp.py

    r72 r74  
    3535    return '%d:%2d:%.2f'%(H,M,S)
    3636
    37 def ValEsd(value,esd=0,nTZ=False):
     37def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
    3838    # returns value(esd) string; nTZ=True for no trailing zeros
    3939    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
     
    5959
    6060       
    61        
     61#GSASII peak fitting routine: Thompson, Cox & Hastings; Finger, Cox & Jephcoat model       
    6262
    6363def DoPeakFit(peaks,background,limits,inst,data):
     
    302302    return rdsq
    303303   
    304 def scaleAbyV(A,V):
    305     v = calc_V(A)
    306     scale = math.exp(math.log(v/V)/3.)**2
    307     for i in range(6):
    308         A[i] *= scale
    309    
    310 def ranaxis(dmin,dmax):
    311     import random as rand
    312     return rand.random()*(dmax-dmin)+dmin
    313    
    314 def ran2axis(k,N):
    315     import random as rand
    316     T = 1.5+0.49*k/N
    317 #    B = 0.99-0.49*k/N
    318 #    B = 0.99-0.049*k/N
    319     B = 0.99-0.149*k/N
    320     R = (T-B)*rand.random()+B
    321     return R
    322    
    323 def ranNaxis(k,N):
    324     import random as rand
    325     T = 1.0+1.0*k/N
    326     B = 1.0-1.0*k/N
    327     R = (T-B)*rand.random()+B
    328     return R
    329    
    330 def ranAbyV(Bravais,dmin,dmax,V):
    331     cell = [0,0,0,0,0,0]
    332     bad = True
    333     while bad:
    334         bad = False
    335         cell = rancell(Bravais,dmin,dmax)
    336         G,g = cell2Gmat(cell)
    337         A = Gmat2A(G)
    338         if calc_rVsq(A) < 1:
    339             scaleAbyV(A,V)
    340             cell = A2cell(A)
    341             for i in range(3):
    342                 bad |= cell[i] < dmin
    343     return A
    344    
    345 def ranAbyR(Bravais,A,k,N,ranFunc):
    346     R = ranFunc(k,N)
    347     if Bravais in [0,1,2]:          #cubic - not used
    348         A[0] = A[1] = A[2] = A[0]*R
    349         A[3] = A[4] = A[5] = 0.
    350     elif Bravais in [3,4]:          #hexagonal/trigonal
    351         A[0] = A[1] = A[3] = A[0]*R
    352         A[2] *= R
    353         A[4] = A[5] = 0.       
    354     elif Bravais in [5,6]:          #tetragonal
    355         A[0] = A[1] = A[0]*R
    356         A[2] *= R
    357         A[3] = A[4] = A[5] = 0.       
    358     elif Bravais in [7,8,9,10]:     #orthorhombic
    359         A[0] *= R
    360         A[1] *= R
    361         A[2] *= R
    362         A[3] = A[4] = A[5] = 0.       
    363     elif Bravais in [11,12]:        #monoclinic
    364         A[0] *= R
    365         A[1] *= R
    366         A[2] *= R
    367         A[4] *= R
    368         A[3] = A[5] = 0.       
    369     else:                           #triclinic
    370         A[0] *= R
    371         A[1] *= R
    372         A[2] *= R
    373         A[3] *= R
    374         A[4] *= R
    375         A[5] *= R
    376     return A
    377    
    378 def rancell(Bravais,dmin,dmax):
    379     if Bravais in [0,1,2]:          #cubic
    380         a = b = c = ranaxis(dmin,dmax)
    381         alp = bet = gam = 90
    382     elif Bravais in [3,4]:          #hexagonal/trigonal
    383         a = b = ranaxis(dmin,dmax)
    384         c = ranaxis(dmin,dmax)
    385         alp = bet =  90
    386         gam = 120
    387     elif Bravais in [5,6]:          #tetragonal
    388         a = b = ranaxis(dmin,dmax)
    389         c = ranaxis(dmin,dmax)
    390         alp = bet = gam = 90
    391     elif Bravais in [7,8,9,10]:       #orthorhombic - F,I,C,P - a<b<c convention
    392         abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
    393         abc.sort()
    394         a = abc[0]
    395         b = abc[1]
    396         c = abc[2]
    397         alp = bet = gam = 90
    398     elif Bravais in [11,12]:        #monoclinic - C,P - a<c convention
    399         ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
    400         ac.sort()
    401         a = ac[0]
    402         b = ranaxis(dmin,dmax)
    403         c = ac[1]
    404         alp = gam = 90
    405         bet = ranaxis(90.,130.)
    406     else:                           #triclinic - a<b<c convention
    407         abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
    408         abc.sort()
    409         a = abc[0]
    410         b = abc[1]
    411         c = abc[2]
    412         r = 0.5*b/c
    413         alp = ranaxis(acosd(r),acosd(-r))
    414         r = 0.5*a/c
    415         bet = ranaxis(acosd(r),acosd(-r))
    416         r = 0.5*a/b
    417         gam = ranaxis(acosd(r),acosd(-r)) 
    418     return [a,b,c,alp,bet,gam]
    419    
    420 def makeMat(Angle,Axis):
    421     #Make rotation matrix from Angle in degrees,Axis =0 for rotation about x, =1 for about y, etc.
    422     cs = cosd(Angle)
    423     ss = sind(Angle)
    424     M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
    425     return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
    426                    
    427304def MaxIndex(dmin,A):
    428305    #finds maximum allowed hkl for given A within dmin
     
    456333    return X
    457334   
    458 def sortM20(cells):
    459     #cells is M20,X20,Bravais,a,b,c,alp,bet,gam
    460     #sort highest M20 1st
    461     T = []
    462     for i,M in enumerate(cells):
    463         T.append((M[0],i))
    464     D = dict(zip(T,cells))
    465     T.sort()
    466     T.reverse()
    467     X = []
    468     for key in T:
    469         X.append(D[key])
    470     return X
    471                
    472335 
    473336def GenHBravais(dmin,Bravais,A):
    474337    '''Generate the positionally unique powder diffraction reflections
    475     for a lattice and Bravais type
    476     dmin - minimum d-spacing
    477     Bravais in range(14) to indicate Bravais lattice; 0-2 cubic, 3,4 - hexagonal/trigonal,
    478     5,6 - tetragonal, 7-10 - orthorhombic, 11,12 - monoclinic, 13 - triclinic
    479     A - as defined in calc_rDsq
    480     returns HKL = [h,k,l,d,0] sorted so d largest first'''
     338    for a lattice and Bravais type'''
     339# dmin - minimum d-spacing
     340# Bravais in range(14) to indicate Bravais lattice; 0-2 cubic, 3,4 - hexagonal/trigonal,
     341# 5,6 - tetragonal, 7-10 - orthorhombic, 11,12 - monoclinic, 13 - triclinic
     342# A - as defined in calc_rDsq
     343# returns HKL = [h,k,l,d,0] sorted so d largest first
    481344    import math
    482345    if Bravais in [9,11]:
     
    566429    return sortHKLd(HKL,True,False)
    567430   
    568 def SwapXY(x,y):
    569     return [y,x]
    570    
    571431def SwapIndx(Axis,H):
    572432    if Axis in [1,-1]:
     
    597457    '''Generate the crystallographically unique powder diffraction reflections
    598458    for a lattice and Bravais type
    599     dmin - minimum d-spacing
    600     Laue - Laue group symbol = '-1','2/m','mmmm','4/m','6/m','4/mmm','6/mmm',
    601                              '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
    602     Cent - lattice centering = 'P','A','B','C','I','F'
    603     Axis - code for unique monoclinic axis = 'a','b','c'
    604     A - 6 terms as defined in calc_rDsq
    605     returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique
    606     part of reciprocal space ignoring anomalous dispersion
    607459    '''
     460# dmin - minimum d-spacing
     461# Laue - Laue group symbol = '-1','2/m','mmmm','4/m','6/m','4/mmm','6/mmm',
     462#                            '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
     463# Cent - lattice centering = 'P','A','B','C','I','F'
     464# Axis - code for unique monoclinic axis = 'a','b','c'
     465# A - 6 terms as defined in calc_rDsq
     466# returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique
     467# part of reciprocal space ignoring anomalous dispersion
    608468    import math
    609469    Hmax = MaxIndex(dmin,A)
     
    695555    return sortHKLd(HKL,True,True)
    696556   
     557#GSASII cell indexing program: variation on that of A. Coehlo
     558#   includes cell refinement from peal positions (not zero as yet)
     559   
     560def scaleAbyV(A,V):
     561    v = calc_V(A)
     562    scale = math.exp(math.log(v/V)/3.)**2
     563    for i in range(6):
     564        A[i] *= scale
     565   
     566def ranaxis(dmin,dmax):
     567    import random as rand
     568    return rand.random()*(dmax-dmin)+dmin
     569   
     570def ran2axis(k,N):
     571    import random as rand
     572    T = 1.5+0.49*k/N
     573#    B = 0.99-0.49*k/N
     574#    B = 0.99-0.049*k/N
     575    B = 0.99-0.149*k/N
     576    R = (T-B)*rand.random()+B
     577    return R
     578   
     579def ranNaxis(k,N):
     580    import random as rand
     581    T = 1.0+1.0*k/N
     582    B = 1.0-1.0*k/N
     583    R = (T-B)*rand.random()+B
     584    return R
     585   
     586def ranAbyV(Bravais,dmin,dmax,V):
     587    cell = [0,0,0,0,0,0]
     588    bad = True
     589    while bad:
     590        bad = False
     591        cell = rancell(Bravais,dmin,dmax)
     592        G,g = cell2Gmat(cell)
     593        A = Gmat2A(G)
     594        if calc_rVsq(A) < 1:
     595            scaleAbyV(A,V)
     596            cell = A2cell(A)
     597            for i in range(3):
     598                bad |= cell[i] < dmin
     599    return A
     600   
     601def ranAbyR(Bravais,A,k,N,ranFunc):
     602    R = ranFunc(k,N)
     603    if Bravais in [0,1,2]:          #cubic - not used
     604        A[0] = A[1] = A[2] = A[0]*R
     605        A[3] = A[4] = A[5] = 0.
     606    elif Bravais in [3,4]:          #hexagonal/trigonal
     607        A[0] = A[1] = A[3] = A[0]*R
     608        A[2] *= R
     609        A[4] = A[5] = 0.       
     610    elif Bravais in [5,6]:          #tetragonal
     611        A[0] = A[1] = A[0]*R
     612        A[2] *= R
     613        A[3] = A[4] = A[5] = 0.       
     614    elif Bravais in [7,8,9,10]:     #orthorhombic
     615        A[0] *= R
     616        A[1] *= R
     617        A[2] *= R
     618        A[3] = A[4] = A[5] = 0.       
     619    elif Bravais in [11,12]:        #monoclinic
     620        A[0] *= R
     621        A[1] *= R
     622        A[2] *= R
     623        A[4] *= R
     624        A[3] = A[5] = 0.       
     625    else:                           #triclinic
     626        A[0] *= R
     627        A[1] *= R
     628        A[2] *= R
     629        A[3] *= R
     630        A[4] *= R
     631        A[5] *= R
     632    return A
     633   
     634def rancell(Bravais,dmin,dmax):
     635    if Bravais in [0,1,2]:          #cubic
     636        a = b = c = ranaxis(dmin,dmax)
     637        alp = bet = gam = 90
     638    elif Bravais in [3,4]:          #hexagonal/trigonal
     639        a = b = ranaxis(dmin,dmax)
     640        c = ranaxis(dmin,dmax)
     641        alp = bet =  90
     642        gam = 120
     643    elif Bravais in [5,6]:          #tetragonal
     644        a = b = ranaxis(dmin,dmax)
     645        c = ranaxis(dmin,dmax)
     646        alp = bet = gam = 90
     647    elif Bravais in [7,8,9,10]:       #orthorhombic - F,I,C,P - a<b<c convention
     648        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
     649        abc.sort()
     650        a = abc[0]
     651        b = abc[1]
     652        c = abc[2]
     653        alp = bet = gam = 90
     654    elif Bravais in [11,12]:        #monoclinic - C,P - a<c convention
     655        ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
     656        ac.sort()
     657        a = ac[0]
     658        b = ranaxis(dmin,dmax)
     659        c = ac[1]
     660        alp = gam = 90
     661        bet = ranaxis(90.,130.)
     662    else:                           #triclinic - a<b<c convention
     663        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
     664        abc.sort()
     665        a = abc[0]
     666        b = abc[1]
     667        c = abc[2]
     668        r = 0.5*b/c
     669        alp = ranaxis(acosd(r),acosd(-r))
     670        r = 0.5*a/c
     671        bet = ranaxis(acosd(r),acosd(-r))
     672        r = 0.5*a/b
     673        gam = ranaxis(acosd(r),acosd(-r)) 
     674    return [a,b,c,alp,bet,gam]
     675   
    697676def calc_M20(peaks,HKL):
    698677    diff = 0
     
    723702    return M20,X20
    724703   
     704def sortM20(cells):
     705    #cells is M20,X20,Bravais,a,b,c,alp,bet,gam
     706    #sort highest M20 1st
     707    T = []
     708    for i,M in enumerate(cells):
     709        T.append((M[0],i))
     710    D = dict(zip(T,cells))
     711    T.sort()
     712    T.reverse()
     713    X = []
     714    for key in T:
     715        X.append(D[key])
     716    return X
     717               
    725718def IndexPeaks(peaks,HKL):
    726719    import bisect
     
    11291122        return False,0,0
    11301123       
     1124#GSASII image calculations: ellipse fitting & image integration       
     1125       
     1126def makeMat(Angle,Axis):
     1127    #Make rotation matrix from Angle in degrees,Axis =0 for rotation about x, =1 for about y, etc.
     1128    cs = cosd(Angle)
     1129    ss = sind(Angle)
     1130    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
     1131    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
     1132                   
     1133def SwapXY(x,y):
     1134    return [y,x]
     1135   
    11311136def FitRing(ring):
    11321137    Err,parms = FitCircle(ring)
     
    16011606    finally:
    16021607        dlg.Destroy()
    1603        
    1604 def test():
    1605     cell = [5,5,5,90,90,90]
    1606     A = cell2A(cell)
    1607     assert ( calc_V(A) == 125. )
    1608    
    1609     print 'test passed'
    1610 
    1611 print __name__
    1612 if __name__ == '__main__':
    1613     test()
    1614    
     1608           
    16151609       
Note: See TracChangeset for help on using the changeset viewer.