Changeset 78


Ignore:
Timestamp:
Jun 3, 2010 1:22:05 PM (12 years ago)
Author:
vondreel
Message:

split IndexPeaks?, etc. from GSASIIcomp.py & put in new GSASIIindex.py

Location:
trunk
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASII.py

    r68 r78  
    11041104        info = wx.AboutDialogInfo()
    11051105        info.Name = 'GSAS-II'
    1106         info.Version = '0.0.1'
     1106        info.Version = __version__
    11071107        info.Copyright = '''
    11081108Robert B. Von Dreele
  • trunk/GSASIIcomp.py

    r77 r78  
    2828npatand = lambda x: 180.*np.arctan(x)/np.pi
    2929npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
    30 
    31 def sec2HMS(sec):
    32     H = int(sec/3600)
    33     M = int(sec/60-H*60)
    34     S = sec-3600*H-60*M
    35     return '%d:%2d:%.2f'%(H,M,S)
    3630
    3731def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
     
    289283    return True,smin,Rwp,runtime,GoOn
    290284   
    291 #reflection generation routines
    292 #for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
    293 #           cell - a,b,c,alp,bet,gam in A & deg
    294    
    295 def calc_rDsq(H,A):
    296     rdsq = H[0]*H[0]*A[0]+H[1]*H[1]*A[1]+H[2]*H[2]*A[2]+H[0]*H[1]*A[3]+H[0]*H[2]*A[4]+H[1]*H[2]*A[5]
    297     return rdsq
    298    
    299 def calc_rDsqZ(H,A,Z,tth,lam):
    300     rpd = math.pi/180.
    301     rdsq = calc_rDsq(H,A)+Z*math.sin(tth*rpd)*2.0*rpd/(lam*lam)
    302     return rdsq
    303        
    304 def MaxIndex(dmin,A):
    305     #finds maximum allowed hkl for given A within dmin
    306     Hmax = [0,0,0]
    307     try:
    308         cell = G2lat.A2cell(A)
    309     except:
    310         cell = [1,1,1,90,90,90]
    311     for i in range(3):
    312         Hmax[i] = int(round(cell[i]/dmin))
    313     return Hmax
    314    
    315 def sortHKLd(HKLd,ifreverse,ifdup):
    316     #HKLd is a list of [h,k,l,d,...]; ifreverse=True for largest d first
    317     #ifdup = True if duplicate d-spacings allowed
    318     T = []
    319     for i,H in enumerate(HKLd):
    320         if ifdup:
    321             T.append((H[3],i))
    322         else:
    323             T.append(H[3])           
    324     D = dict(zip(T,HKLd))
    325     T.sort()
    326     if ifreverse:
    327         T.reverse()
    328     X = []
    329     okey = ''
    330     for key in T:
    331         if key != okey: X.append(D[key])    #remove duplicate d-spacings
    332         okey = key
    333     return X
    334    
    335 def SwapIndx(Axis,H):
    336     if Axis in [1,-1]:
    337         return H
    338     elif Axis in [2,-3]:
    339         return [H[1],H[2],H[0]]
    340     else:
    341         return [H[2],H[0],H[1]]
    342        
    343 def CentCheck(Cent,H):
    344     h,k,l = H
    345     if Cent == 'A' and (k+l)%2:
    346         return False
    347     elif Cent == 'B' and (h+l)%2:
    348         return False
    349     elif Cent == 'C' and (h+k)%2:
    350         return False
    351     elif Cent == 'I' and (h+k+l)%2:
    352         return False
    353     elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
    354         return False
    355     elif Cent == 'R' and (-h+k+l)%3:
    356         return False
    357     else:
    358         return True
    359                                    
    360 def GenHBravais(dmin,Bravais,A):
    361     '''Generate the positionally unique powder diffraction reflections
    362     for a lattice and Bravais type'''
    363 # dmin - minimum d-spacing
    364 # Bravais in range(14) to indicate Bravais lattice; 0-2 cubic, 3,4 - hexagonal/trigonal,
    365 # 5,6 - tetragonal, 7-10 - orthorhombic, 11,12 - monoclinic, 13 - triclinic
    366 # A - as defined in calc_rDsq
    367 # returns HKL = [h,k,l,d,0] sorted so d largest first
    368     import math
    369     if Bravais in [9,11]:
    370         Cent = 'C'
    371     elif Bravais in [1,5,8]:
    372         Cent = 'I'
    373     elif Bravais in [0,7]:
    374         Cent = 'F'
    375     elif Bravais in [3]:
    376         Cent = 'R'
    377     else:
    378         Cent = 'P'
    379     Hmax = MaxIndex(dmin,A)
    380     dminsq = 1./(dmin**2)
    381     HKL = []
    382     if Bravais == 13:                       #triclinic
    383         for l in range(-Hmax[2],Hmax[2]+1):
    384             for k in range(-Hmax[1],Hmax[1]+1):
    385                 hmin = 0
    386                 if (k < 0): hmin = 1
    387                 if (k ==0 and l < 0): hmin = 1
    388                 for h in range(hmin,Hmax[0]+1):
    389                     H=[h,k,l]
    390                     rdsq = calc_rDsq(H,A)
    391                     if 0 < rdsq <= dminsq:
    392                         HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
    393     elif Bravais in [11,12]:                #monoclinic - b unique
    394         Hmax = SwapIndx(2,Hmax)
    395         for h in range(Hmax[0]+1):
    396             for k in range(-Hmax[1],Hmax[1]+1):
    397                 lmin = 0
    398                 if k < 0:lmin = 1
    399                 for l in range(lmin,Hmax[2]+1):
    400                     [h,k,l] = SwapIndx(-2,[h,k,l])
    401                     H = []
    402                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    403                     if H:
    404                         rdsq = calc_rDsq(H,A)
    405                         if 0 < rdsq <= dminsq:
    406                             HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
    407                     [h,k,l] = SwapIndx(2,[h,k,l])
    408     elif Bravais in [7,8,9,10]:            #orthorhombic
    409         for h in range(Hmax[0]+1):
    410             for k in range(Hmax[1]+1):
    411                 for l in range(Hmax[2]+1):
    412                     H = []
    413                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    414                     if H:
    415                         rdsq = calc_rDsq(H,A)
    416                         if 0 < rdsq <= dminsq:
    417                             HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
    418     elif Bravais in [5,6]:                  #tetragonal
    419         for l in range(Hmax[2]+1):
    420             for k in range(Hmax[1]+1):
    421                 for h in range(k,Hmax[0]+1):
    422                     H = []
    423                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    424                     if H:
    425                         rdsq = calc_rDsq(H,A)
    426                         if 0 < rdsq <= dminsq:
    427                             HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
    428     elif Bravais in [3,4]:
    429         lmin = 0
    430         if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
    431         for l in range(lmin,Hmax[2]+1):
    432             for k in range(Hmax[1]+1):
    433                 hmin = k
    434                 if l < 0: hmin += 1
    435                 for h in range(hmin,Hmax[0]+1):
    436                     H = []
    437                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    438                     if H:
    439                         rdsq = calc_rDsq(H,A)
    440                         if 0 < rdsq <= dminsq:
    441                             HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
    442 
    443     else:                                   #cubic
    444         for l in range(Hmax[2]+1):
    445             for k in range(l,Hmax[1]+1):
    446                 for h in range(k,Hmax[0]+1):
    447                     H = []
    448                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    449                     if H:
    450                         rdsq = calc_rDsq(H,A)
    451                         if 0 < rdsq <= dminsq:
    452                             HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
    453     return sortHKLd(HKL,True,False)
    454    
    455 def GenHLaue(dmin,Laue,Cent,Axis,A):
    456     '''Generate the crystallographically unique powder diffraction reflections
    457     for a lattice and Bravais type
    458     '''
    459 # dmin - minimum d-spacing
    460 # Laue - Laue group symbol = '-1','2/m','mmmm','4/m','6/m','4/mmm','6/mmm',
    461 #                            '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
    462 # Cent - lattice centering = 'P','A','B','C','I','F'
    463 # Axis - code for unique monoclinic axis = 'a','b','c'
    464 # A - 6 terms as defined in calc_rDsq
    465 # returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique
    466 # part of reciprocal space ignoring anomalous dispersion
    467     import math
    468     Hmax = MaxIndex(dmin,A)
    469     dminsq = 1./(dmin**2)
    470     HKL = []
    471     if Laue == '-1':                       #triclinic
    472         for l in range(-Hmax[2],Hmax[2]+1):
    473             for k in range(-Hmax[1],Hmax[1]+1):
    474                 hmin = 0
    475                 if (k < 0) or (k ==0 and l < 0): hmin = 1
    476                 for h in range(hmin,Hmax[0]+1):
    477                     H = []
    478                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    479                     rdsq = calc_rDsq(H,A)
    480                     if 0 < rdsq <= dminsq:
    481                         HKL.append([h,k,l,1/math.sqrt(rdsq)])
    482     elif Laue == '2/m':                #monoclinic
    483         Hmax = SwapIndx(Axis,Hmax)
    484         for h in range(Hmax[0]+1):
    485             for k in range(-Hmax[1],Hmax[1]+1):
    486                 lmin = 0
    487                 if k < 0:lmin = 1
    488                 for l in range(lmin,Hmax[2]+1):
    489                     [h,k,l] = SwapIndx(-Axis,[h,k,l])
    490                     H = []
    491                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    492                     if H:
    493                         rdsq = calc_rDsq(H,A)
    494                         if 0 < rdsq <= dminsq:
    495                             HKL.append([h,k,l,1/math.sqrt(rdsq)])
    496                     [h,k,l] = SwapIndx(Axis,[h,k,l])
    497     elif Laue in ['mmm','4/m','6/m']:            #orthorhombic
    498         for l in range(Hmax[2]+1):
    499             for h in range(Hmax[0]+1):
    500                 kmin = 1
    501                 if Laue == 'mmm' or h ==0: kmin = 0
    502                 for k in range(kmin,Hmax[2]+1):
    503                     H = []
    504                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    505                     if H:
    506                         rdsq = calc_rDsq(H,A)
    507                         if 0 < rdsq <= dminsq:
    508                             HKL.append([h,k,l,1/math.sqrt(rdsq)])
    509     elif Laue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
    510         for l in range(Hmax[2]+1):
    511             for h in range(Hmax[0]+1):
    512                 for k in range(h+1):
    513                     H = []
    514                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    515                     if H:
    516                         rdsq = calc_rDsq(H,A)
    517                         if 0 < rdsq <= dminsq:
    518                             HKL.append([h,k,l,1/math.sqrt(rdsq)])
    519     elif Laue in ['3m1','31m','3','3R','3mR']:                  #trigonals
    520         for l in range(-Hmax[2],Hmax[2]+1):
    521             hmin = 0
    522             if l < 0: hmin = 1
    523             for h in range(hmin,Hmax[0]+1):
    524                 if Laue in ['3R','3']:
    525                     kmax = h
    526                     kmin = -int((h-1)/2.)
    527                 else:
    528                     kmin = 0
    529                     kmax = h
    530                     if Laue in ['3m1','3mR'] and l < 0: kmax = h-1
    531                     if Laue == '31m' and l < 0: kmin = 1
    532                 for k in range(kmin,kmax+1):
    533                     H = []
    534                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    535                     if H:
    536                         rdsq = calc_rDsq(H,A)
    537                         if 0 < rdsq <= dminsq:
    538                             HKL.append([h,k,l,1/math.sqrt(rdsq)])
    539     else:                                   #cubic
    540         for h in range(Hmax[0]+1):
    541             for k in range(h+1):
    542                 lmin = 0
    543                 lmax = k
    544                 if Laue =='m3':
    545                     lmax = h-1
    546                     if h == k: lmax += 1
    547                 for l in range(lmin,lmax+1):
    548                     H = []
    549                     if CentCheck(Cent,[h,k,l]): H=[h,k,l]
    550                     if H:
    551                         rdsq = calc_rDsq(H,A)
    552                         if 0 < rdsq <= dminsq:
    553                             HKL.append([h,k,l,1/math.sqrt(rdsq)])
    554     return sortHKLd(HKL,True,True)
    555    
    556 #GSASII cell indexing program: variation on that of A. Coehlo
    557 #   includes cell refinement from peak positions (not zero as yet)
    558    
    559 def scaleAbyV(A,V):
    560     v = G2lat.calc_V(A)
    561     scale = math.exp(math.log(v/V)/3.)**2
    562     for i in range(6):
    563         A[i] *= scale
    564    
    565 def ranaxis(dmin,dmax):
    566     import random as rand
    567     return rand.random()*(dmax-dmin)+dmin
    568    
    569 def ran2axis(k,N):
    570     import random as rand
    571     T = 1.5+0.49*k/N
    572 #    B = 0.99-0.49*k/N
    573 #    B = 0.99-0.049*k/N
    574     B = 0.99-0.149*k/N
    575     R = (T-B)*rand.random()+B
    576     return R
    577    
    578 #def ranNaxis(k,N):
    579 #    import random as rand
    580 #    T = 1.0+1.0*k/N
    581 #    B = 1.0-1.0*k/N
    582 #    R = (T-B)*rand.random()+B
    583 #    return R
    584    
    585 def ranAbyV(Bravais,dmin,dmax,V):
    586     cell = [0,0,0,0,0,0]
    587     bad = True
    588     while bad:
    589         bad = False
    590         cell = rancell(Bravais,dmin,dmax)
    591         G,g = G2lat.cell2Gmat(cell)
    592         A = G2lat.Gmat2A(G)
    593         if G2lat.calc_rVsq(A) < 1:
    594             scaleAbyV(A,V)
    595             cell = G2lat.A2cell(A)
    596             for i in range(3):
    597                 bad |= cell[i] < dmin
    598     return A
    599    
    600 def ranAbyR(Bravais,A,k,N,ranFunc):
    601     R = ranFunc(k,N)
    602     if Bravais in [0,1,2]:          #cubic - not used
    603         A[0] = A[1] = A[2] = A[0]*R
    604         A[3] = A[4] = A[5] = 0.
    605     elif Bravais in [3,4]:          #hexagonal/trigonal
    606         A[0] = A[1] = A[3] = A[0]*R
    607         A[2] *= R
    608         A[4] = A[5] = 0.       
    609     elif Bravais in [5,6]:          #tetragonal
    610         A[0] = A[1] = A[0]*R
    611         A[2] *= R
    612         A[3] = A[4] = A[5] = 0.       
    613     elif Bravais in [7,8,9,10]:     #orthorhombic
    614         A[0] *= R
    615         A[1] *= R
    616         A[2] *= R
    617         A[3] = A[4] = A[5] = 0.       
    618     elif Bravais in [11,12]:        #monoclinic
    619         A[0] *= R
    620         A[1] *= R
    621         A[2] *= R
    622         A[4] *= R
    623         A[3] = A[5] = 0.       
    624     else:                           #triclinic
    625         A[0] *= R
    626         A[1] *= R
    627         A[2] *= R
    628         A[3] *= R
    629         A[4] *= R
    630         A[5] *= R
    631     return A
    632    
    633 def rancell(Bravais,dmin,dmax):
    634     if Bravais in [0,1,2]:          #cubic
    635         a = b = c = ranaxis(dmin,dmax)
    636         alp = bet = gam = 90
    637     elif Bravais in [3,4]:          #hexagonal/trigonal
    638         a = b = ranaxis(dmin,dmax)
    639         c = ranaxis(dmin,dmax)
    640         alp = bet =  90
    641         gam = 120
    642     elif Bravais in [5,6]:          #tetragonal
    643         a = b = ranaxis(dmin,dmax)
    644         c = ranaxis(dmin,dmax)
    645         alp = bet = gam = 90
    646     elif Bravais in [7,8,9,10]:       #orthorhombic - F,I,C,P - a<b<c convention
    647         abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
    648         abc.sort()
    649         a = abc[0]
    650         b = abc[1]
    651         c = abc[2]
    652         alp = bet = gam = 90
    653     elif Bravais in [11,12]:        #monoclinic - C,P - a<c convention
    654         ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
    655         ac.sort()
    656         a = ac[0]
    657         b = ranaxis(dmin,dmax)
    658         c = ac[1]
    659         alp = gam = 90
    660         bet = ranaxis(90.,130.)
    661     else:                           #triclinic - a<b<c convention
    662         abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
    663         abc.sort()
    664         a = abc[0]
    665         b = abc[1]
    666         c = abc[2]
    667         r = 0.5*b/c
    668         alp = ranaxis(acosd(r),acosd(-r))
    669         r = 0.5*a/c
    670         bet = ranaxis(acosd(r),acosd(-r))
    671         r = 0.5*a/b
    672         gam = ranaxis(acosd(r),acosd(-r)) 
    673     return [a,b,c,alp,bet,gam]
    674    
    675 def calc_M20(peaks,HKL):
    676     diff = 0
    677     X20 = 0
    678     for Nobs20,peak in enumerate(peaks):
    679         if peak[3]:
    680             Qobs = 1.0/peak[7]**2
    681             Qcalc = 1.0/peak[8]**2
    682             diff += abs(Qobs-Qcalc)
    683         elif peak[2]:
    684             X20 += 1
    685         if Nobs20 == 19:
    686             d20 = peak[7]
    687             break
    688     else:
    689         d20 = peak[7]
    690         Nobs20 = len(peaks)
    691     for N20,hkl in enumerate(HKL):
    692         if hkl[3] < d20:
    693             break               
    694     eta = diff/Nobs20
    695     Q20 = 1.0/d20**2
    696     if diff:
    697         M20 = Q20/(2.0*diff)
    698     else:
    699         M20 = 0
    700     M20 /= (1.+X20)
    701     return M20,X20
    702    
    703 def sortM20(cells):
    704     #cells is M20,X20,Bravais,a,b,c,alp,bet,gam
    705     #sort highest M20 1st
    706     T = []
    707     for i,M in enumerate(cells):
    708         T.append((M[0],i))
    709     D = dict(zip(T,cells))
    710     T.sort()
    711     T.reverse()
    712     X = []
    713     for key in T:
    714         X.append(D[key])
    715     return X
    716                
    717 def IndexPeaks(peaks,HKL):
    718     import bisect
    719     hklds = [1000.0]                                    #make bounded list of available d-spacings
    720     N = len(HKL)
    721     if N == 0: return False
    722     for hkl in HKL:
    723         hklds.append(hkl[3])
    724     hklds.append(0.0)
    725     hklds.sort()                                        # ascending sort - upper bound at end
    726     hklmax = [0,0,0]
    727     for ipk,peak in enumerate(peaks):
    728         if peak[2]:
    729             i = bisect.bisect_right(hklds,peak[7])          # find peak position in hkl list
    730             dm = peak[7]-hklds[i-1]                         # peak to neighbor hkls in list
    731             dp = hklds[i]-peak[7]
    732             pos = N-i                                       # reverse the order
    733             if dp > dm: pos += 1                            # closer to upper than lower
    734             hkl = HKL[pos]                                 # put in hkl
    735             if hkl[4] >= 0:                                 # peak already assigned - test if this one better
    736                 opeak = peaks[hkl[4]]
    737                 dold = abs(opeak[7]-hkl[3])
    738                 dnew = min(dm,dp)
    739                 if dold > dnew:                             # new better - zero out old
    740                     for j in range(3):
    741                         opeak[j+4] = 0
    742                     opeak[8] = 0.
    743                 else:                                       # old better - do nothing
    744                     continue               
    745             hkl[4] = ipk
    746             for j in range(3):
    747                 peak[j+4] = hkl[j]
    748             peak[8] = hkl[3]                                # fill in d-calc
    749     for peak in peaks:
    750         peak[3] = False
    751         if peak[2]:
    752             if peak[8] > 0.:
    753                 for j in range(3):
    754                     if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
    755                 peak[3] = True
    756     if hklmax[0]*hklmax[1]*hklmax[2] > 0:
    757         return True
    758     else:
    759         return False
    760    
    761 def FitHKL(ibrav,peaks,A,wtP):
    762     def ShiftTest(a,b):
    763         if b < -0.1*a:
    764             b = -0.0001*a
    765         return b
    766     smin = 0.
    767     first = True
    768     for peak in peaks:
    769         if peak[2] and peak[3]:
    770             h,k,l = H = peak[4:7]
    771             Qo = 1./peak[7]**2
    772             Qc = calc_rDsq(H,A)
    773             try:
    774                 peak[8] = 1./math.sqrt(Qc)
    775             except:
    776                 print G2lat.A2invcell(A)
    777             delt = Qo-Qc
    778             smin += delt**2
    779             dp = []
    780             if ibrav in [0,1,2]:            #m3m
    781                 dp.append(h*h+k*k+l*l)
    782             elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
    783                 dp.append(h*h+k*k+h*k)
    784                 dp.append(l*l)
    785             elif ibrav in [5,6]:            #4/mmm
    786                 dp.append(h*h+k*k)
    787                 dp.append(l*l)
    788             elif ibrav in [7,8,9,10]:       #mmm
    789                 dp.append(h*h)
    790                 dp.append(k*k)
    791                 dp.append(l*l)
    792             elif ibrav in [11,12]:          #2/m
    793                 dp.append(h*h)
    794                 dp.append(k*k)
    795                 dp.append(l*l)
    796                 dp.append(h*l)
    797             else:                           #1
    798 #    derivatives for a0*h^2+a1*k^2+a2*l^2+a3*h*k+a4*h*l+a5*k*l
    799                 dp.append(h*h)
    800                 dp.append(k*k)
    801                 dp.append(l*l)
    802                 dp.append(h*k)
    803                 dp.append(h*l)
    804                 dp.append(k*l)
    805             if first:
    806                 first = False
    807                 M = len(dp)
    808                 B = np.zeros(shape=(M,M))
    809                 V = np.zeros(shape=(M))
    810             dwt = peak[7]**wtP
    811             B,V = pyp.buildmv(delt*dwt,dwt,M,dp,B,V)
    812     if nl.det(B) > 0.0:
    813         try:
    814             b = nl.solve(B,V)
    815             B = nl.inv(B)
    816             sig = np.diag(B)
    817         except SingularMatrix:
    818             return False,0
    819         if ibrav in [0,1,2]:            #m3m
    820             A[0] += ShiftTest(A[0],b[0])
    821             A[1] = A[2] = A[0]
    822         elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
    823             A[0] += ShiftTest(A[0],b[0])
    824             A[1] = A[3] = A[0]
    825             A[2] += ShiftTest(A[2],b[1])
    826         elif ibrav in [5,6]:            #4/mmm
    827             A[0] += ShiftTest(A[0],b[0])
    828             A[1] = A[0]
    829             A[2] += ShiftTest(A[2],b[1])
    830         elif ibrav in [7,8,9,10]:       #mmm
    831             for i in range(3):
    832                 A[i] += ShiftTest(A[i],b[i])
    833         elif ibrav in [11,12]:          #2/m
    834             for i in range(3):
    835                 A[i] += ShiftTest(A[i],b[i])
    836             A[4] += ShiftTest(A[4],b[3])
    837             A[4] = min(1.4*math.sqrt(A[0]*A[2]),A[4])   #min beta star = 45
    838         else:                           #1
    839             oldV = math.sqrt(1./G2lat.calc_rVsq(A))
    840             oldA = A[:]
    841             for i in range(6):
    842                 A[i] += b[i]*0.2
    843             A[3] = min(1.1*math.sqrt(max(0,A[1]*A[2])),A[3])
    844             A[4] = min(1.1*math.sqrt(max(0,A[0]*A[2])),A[4])
    845             A[5] = min(1.1*math.sqrt(max(0,A[0]*A[1])),A[5])
    846             ratio = math.sqrt(1./G2lat.calc_rVsq(A))/oldV
    847             if 0.9 > ratio or ratio > 1.1:
    848                 A = oldA
    849 #    else:
    850 #        print B
    851 #        print V
    852 #        for peak in peaks:
    853 #            print peak
    854     return True,smin
    855        
    856 def FitHKLZ(ibrav,peaks,Z,A):
    857     return A,Z
    858    
    859 def rotOrthoA(A):
    860     return [A[1],A[2],A[0],0,0,0]
    861    
    862 def swapMonoA(A):
    863     return [A[2],A[1],A[0],0,A[4],0]
    864    
    865 def oddPeak(indx,peaks):
    866     noOdd = True
    867     for peak in peaks:
    868         H = peak[4:7]
    869         if H[indx] % 2:
    870             noOdd = False
    871     return noOdd
    872    
    873 def halfCell(ibrav,A,peaks):
    874     if ibrav in [0,1,2]:
    875         if oddPeak(0,peaks):
    876             A[0] *= 2
    877             A[1] = A[2] = A[0]
    878     elif ibrav in [3,4,5,6]:
    879         if oddPeak(0,peaks):
    880             A[0] *= 2
    881             A[1] = A[0]
    882         if oddPeak(2,peaks):
    883             A[2] *=2
    884     else:
    885         if oddPeak(0,peaks):
    886             A[0] *=2
    887         if oddPeak(1,peaks):
    888             A[1] *=2
    889         if oddPeak(2,peaks):
    890             A[2] *=2
    891     return A
    892    
    893 def getDmin(peaks):
    894     return peaks[-1][7]
    895    
    896 def getDmax(peaks):
    897     return peaks[0][7]
    898    
    899 def refinePeaks(peaks,ibrav,A):
    900     dmin = getDmin(peaks)
    901     smin = 1.0e10
    902     pwr = 6
    903     maxTries = 10
    904     if ibrav == 13:
    905         pwr = 4
    906         maxTries = 10
    907     OK = False
    908     tries = 0
    909     HKL = GenHBravais(dmin,ibrav,A)
    910     while IndexPeaks(peaks,HKL):
    911         Pwr = pwr - 2*(tries % 2)
    912         HKL = []
    913         tries += 1
    914         osmin = smin
    915         oldA = A
    916         OK,smin = FitHKL(ibrav,peaks,A,Pwr)
    917         for a in A[:3]:
    918             if a < 0:
    919                 A = oldA
    920                 OK = False
    921                 break
    922         if OK:
    923             HKL = GenHBravais(dmin,ibrav,A)
    924         if len(HKL) == 0: break                         #absurd cell obtained!
    925         rat = (osmin-smin)/smin
    926         if abs(rat) < 1.0e-5 or not OK: break
    927         if tries > maxTries: break
    928     if OK:
    929         OK,smin = FitHKL(ibrav,peaks,A,4)
    930     M20,X20 = calc_M20(peaks,HKL)
    931     return len(HKL),M20,X20
    932        
    933 def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
    934 # dlg & ncMax are used for wx progress bar
    935 # A != 0 find the best A near input A,
    936 # A = 0 for random cell, volume normalized to V1;
    937 # returns number of generated hkls, M20, X20 & A for best found
    938     mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10]
    939     dmin = getDmin(peaks)-0.05
    940     amin = 2.5
    941     amax = 5.*getDmax(peaks)
    942     Asave = []
    943     GoOn = True
    944     if A:
    945         HKL = GenHBravais(dmin,ibrav,A[:])
    946         if len(HKL) > mHKL[ibrav]:
    947             IndexPeaks(peaks,HKL)
    948             Asave.append([calc_M20(peaks,HKL),A[:]])
    949     tries = 0
    950     while tries < Ntries:
    951         if A:
    952             Anew = ranAbyR(ibrav,A[:],tries+1,Ntries,ran2axis)
    953             if ibrav in [11,12,13]:
    954                 Anew = ranAbyR(ibrav,A[:],tries/10+1,Ntries,ran2axis)
    955         else:
    956             Anew = ranAbyV(ibrav,amin,amax,V1)
    957         HKL = GenHBravais(dmin,ibrav,Anew)
    958         if len(HKL) > mHKL[ibrav] and IndexPeaks(peaks,HKL):
    959             Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
    960             Asave.append([calc_M20(peaks,HKL),Anew[:]])
    961             if ibrav == 9:                          #C-centered orthorhombic
    962                 for i in range(2):
    963                     Anew = rotOrthoA(Anew[:])
    964                     Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
    965                     HKL = GenHBravais(dmin,ibrav,Anew)
    966                     IndexPeaks(peaks,HKL)
    967                     Asave.append([calc_M20(peaks,HKL),Anew[:]])
    968             elif ibrav == 11:                      #C-centered monoclinic
    969                 Anew = swapMonoA(Anew[:])
    970                 Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
    971                 HKL = GenHBravais(dmin,ibrav,Anew)
    972                 IndexPeaks(peaks,HKL)
    973                 Asave.append([calc_M20(peaks,HKL),Anew[:]])
    974         else:
    975             break
    976         Nc = len(HKL)
    977         if Nc >= ncMax:
    978             GoOn = False
    979         elif dlg:
    980             GoOn = dlg.Update(Nc)[0]
    981             if not GoOn:
    982                 break
    983         tries += 1   
    984     X = sortM20(Asave)
    985     if X:
    986         Lhkl,M20,X20 = refinePeaks(peaks,ibrav,X[0][1])
    987         return GoOn,Lhkl,M20,X20,X[0][1]
    988     else:
    989         return GoOn,0,0,0,0
    990        
    991 def monoCellReduce(ibrav,A):
    992     a,b,c,alp,bet,gam = G2lat.A2cell(A)
    993     G,g = G2lat.A2Gmat(A)
    994     if ibrav in [11]:
    995         u = [0,0,-1]
    996         v = [1,0,2]
    997         anew = math.sqrt(np.dot(np.dot(v,g),v))
    998         if anew < a:
    999             cang = np.dot(np.dot(u,g),v)/(anew*c)
    1000             beta = acosd(-abs(cang))
    1001             A = G2lat.cell2A([anew,b,c,90,beta,90])
    1002     else:
    1003         u = [-1,0,0]
    1004         v = [1,0,1]
    1005         cnew = math.sqrt(np.dot(np.dot(v,g),v))
    1006         if cnew < c:
    1007             cang = np.dot(np.dot(u,g),v)/(a*cnew)
    1008             beta = acosd(-abs(cang))
    1009             A = G2lat.cell2A([a,b,cnew,90,beta,90])
    1010     return A
    1011 
    1012 def DoIndexPeaks(peaks,inst,controls,bravais):
    1013    
    1014     def peakDspace(peaks,A):
    1015         for peak in peaks:
    1016             if peak[3]:
    1017                 dsq = calc_rDsq(peak[4:7],A)
    1018                 if dsq > 0:
    1019                     peak[8] = 1./math.sqrt(dsq)
    1020         return
    1021     delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
    1022     amin = 2.5
    1023     amax = 5.0*getDmax(peaks)
    1024     dmin = getDmin(peaks)-delt
    1025     bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
    1026         'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C',
    1027         'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic']
    1028     tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
    1029     N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,  50,50, 200]
    1030     N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,   2,2,   4]
    1031     Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,   2,2,   4]
    1032     Nobs = len(peaks)
    1033     wave = inst[1]
    1034     if len(inst) > 10:
    1035         zero = inst[3]
    1036     else:
    1037         zero = inst[2]
    1038     print "%s %8.5f %6.3f" % ('wavelength, zero =',wave,zero)
    1039     print "%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax)
    1040     ifzero,maxzero,ncno = controls[:3]
    1041     ncMax = Nobs*ncno
    1042     print "%s %d %s %d %s %d" % ('change zero =',ifzero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs)
    1043     cells = []
    1044     for ibrav in range(14):
    1045         begin = time.time()
    1046         if bravais[ibrav]:
    1047             print 'cell search for ',bravaisNames[ibrav]
    1048             print '      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test'
    1049             V1 = controls[3]
    1050             bestM20 = 0
    1051             topM20 = 0
    1052             cycle = 0
    1053             while cycle < 5:
    1054                 dlg = wx.ProgressDialog("Generated reflections",tries[cycle]+" cell search for "+bravaisNames[ibrav],ncMax,
    1055                     style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
    1056                 screenSize = wx.DisplaySize()
    1057                 Size = dlg.GetSize()
    1058                 dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
    1059                 try:
    1060                     GoOn = True
    1061                     while GoOn:                                                 #Loop over increment of volume
    1062                         N2 = 0
    1063                         while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
    1064                             if ibrav > 2:
    1065                                 if not N2:
    1066                                     GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
    1067                                 if A:
    1068                                     GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
    1069                             else:
    1070                                 GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
    1071                             if Nc >= ncMax:
    1072                                 GoOn = False
    1073                                 break
    1074                             elif 3*Nc < Nobs:
    1075                                 N2 = 10
    1076                                 break
    1077                             else:
    1078                                 if not GoOn:
    1079                                     break
    1080                                 if M20 > 1.0:
    1081                                     bestM20 = max(bestM20,M20)
    1082                                     A = halfCell(ibrav,A[:],peaks)
    1083                                     if ibrav in [12]:
    1084                                         A = monoCellReduce(ibrav,A[:])
    1085                                     HKL = GenHBravais(dmin,ibrav,A)
    1086                                     IndexPeaks(peaks,HKL)
    1087                                     a,b,c,alp,bet,gam = G2lat.A2cell(A)
    1088                                     V = G2lat.calc_V(A)
    1089                                     print "%10.3f %3d %3d %10.5f %10.5f %10.5f %10.3f %10.3f %10.3f %10.2f %10.2f" % (M20,X20,Nc,a,b,c,alp,bet,gam,V,V1)
    1090                                     if M20 >= 2.0:
    1091                                         cells.append([M20,X20,ibrav,a,b,c,alp,bet,gam,V,False])
    1092                             if not GoOn:
    1093                                 break
    1094                             N2 += 1
    1095                         if ibrav < 11:
    1096                             V1 *= 1.1
    1097                         elif ibrav in range(11,14):
    1098                             V1 *= 1.05
    1099                         if not GoOn:
    1100                             if bestM20 > topM20:
    1101                                 topM20 = bestM20
    1102                                 if cells:
    1103                                     V1 = cells[0][9]
    1104                                 else:
    1105                                     V1 = 25
    1106                                 ncMax += Nobs
    1107                                 cycle += 1
    1108                                 print 'Restart search, new Max Nc = ',ncMax
    1109                             else:
    1110                                 cycle = 10
    1111                 finally:
    1112                     dlg.Destroy()
    1113             print '%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
    1114                 ', elapsed time = ',sec2HMS(time.time()-begin))
    1115            
    1116     if cells:
    1117         cells = sortM20(cells)
    1118         cells[0][-1] = True
    1119         return True,dmin,cells
    1120     else:
    1121         return False,0,0
    1122        
    1123285#GSASII image calculations: ellipse fitting & image integration       
    1124286       
     
    1130292    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
    1131293                   
    1132 def SwapXY(x,y):
    1133     return [y,x]
    1134    
    1135294def FitRing(ring):
    1136295    Err,parms = FitCircle(ring)
     
    1184343        sr2 = math.sqrt(f/b)
    1185344        if sr1 > sr2:
    1186             sr1,sr2 = SwapXY(sr1,sr2)
     345            sr1,sr2 = sr2,sr1
    1187346            phi -= 90.
    1188347            if phi < -90.:
  • trunk/GSASIIgrid.py

    r75 r78  
    99import GSASIIcomp as G2cmp
    1010import GSASIIlattice as G2lat
     11import GSASIIindex as G2indx
    1112import GSASIIspc as G2spc
    1213import GSASIIElem as G2elem
     
    372373        inst = self.PatternTree.GetItemPyData(GetPatternTreeItemId(self,PatternId, 'Instrument Parameters'))
    373374        data = self.PatternTree.GetItemPyData(PatternId)[1]
    374         OK,smin,Rwp,runtime,GoOn = G2cmp.DoPeakFit(peaks,background,limits,inst,data)
     375        OK,smin,Rwp,runtime,GoOn = G2indx.DoPeakFit(peaks,background,limits,inst,data)
    375376        UpdatePeakGrid(self,peaks)
    376377        G2plt.PlotPatterns(self)
     
    408409        while GoOn:
    409410            osmin = smin
    410             OK,smin,Rwp,runtime,GoOn = G2cmp.DoPeakFit(peaks,background,limits,inst,data)
     411            OK,smin,Rwp,runtime,GoOn = G2indx.DoPeakFit(peaks,background,limits,inst,data)
    411412            UpdatePeakGrid(self,peaks)
    412413            if not OK:
     
    759760                    ibrav = cell[2]
    760761                    A = G2lat.cell2A(cell[3:9])
    761                     self.HKL = G2cmp.GenHBravais(dmin,ibrav,A)
    762                     G2cmp.IndexPeaks(data,self.HKL)
     762                    self.HKL = G2lat.GenHBravais(dmin,ibrav,A)
     763                    G2indx.IndexPeaks(data,self.HKL)
    763764                    for hkl in self.HKL:
    764765                        hkl.append(2.0*asind(inst[1]/(2.*hkl[3])))             
     
    815816        print controls[5]
    816817        ibrav = bravaisSymb.index(controls[5])
    817         dmin = G2cmp.getDmin(peaks)-0.005
    818         Lhkl,M20,X20 = G2cmp.refinePeaks(peaks,ibrav,A)
     818        dmin = G2indx.getDmin(peaks)-0.005
     819        Lhkl,M20,X20 = G2indx.refinePeaks(peaks,ibrav,A)
    819820        controls[6:12] = G2lat.A2cell(A)
    820821        controls[12] = G2lat.calc_V(A)
    821822        data = [controls,bravais,cells,dmin]
    822823        self.PatternTree.SetItemPyData(GetPatternTreeItemId(self,PatternId, 'Unit Cells List'),data)
    823         self.HKL = G2cmp.GenHBravais(dmin,ibrav,A)
     824        self.HKL = G2lat.GenHBravais(dmin,ibrav,A)
    824825        UpdateUnitCellsGrid(self,data)
    825826        print "%s%10.3f" % ('refinement M20 = ',M20)
     
    852853        self.dataFrame.IndexPeaks.Enable(False)
    853854        self.dataFrame.CopyCell.Enable(False)
    854         OK,dmin,cells = G2cmp.DoIndexPeaks(peaks,inst,controls,bravais)
     855        OK,dmin,cells = G2indx.DoIndexPeaks(peaks,inst,controls,bravais)
    855856        if OK:
    856857            data = [controls,bravais,cells,dmin]
     
    859860            bestCell = cells[0]
    860861            if bestCell[0] > 10.:
    861                 self.HKL = G2cmp.GenHBravais(dmin,bestCell[2],G2lat.cell2A(bestCell[3:9]))
     862                self.HKL = G2lat.GenHBravais(dmin,bestCell[2],G2lat.cell2A(bestCell[3:9]))
    862863                for hkl in self.HKL:
    863864                    hkl.append(2.0*asind(inst[1]/(2.*hkl[3])))             
     
    897898                ibrav = cells[r][2]
    898899                A = G2lat.cell2A(cells[r][3:9])
    899                 self.HKL = G2cmp.GenHBravais(dmin,ibrav,A)
     900                self.HKL = G2lat.GenHBravais(dmin,ibrav,A)
    900901                for hkl in self.HKL:
    901902                    hkl.append(2.0*asind(inst[1]/(2.*hkl[3])))
     
    10261027                if cell[-1]:
    10271028                    A = G2lat.cell2A(cell[3:9])
    1028                     self.HKL = G2cmp.GenHBravais(dmin,cell[2],A)
     1029                    self.HKL = G2lat.GenHBravais(dmin,cell[2],A)
    10291030                    for hkl in self.HKL:
    10301031                        hkl.append(2.0*asind(inst[1]/(2.*hkl[3])))
     
    12641265        Utth = float(self.OuterTth.GetValue())
    12651266        if Ltth > Utth:
    1266             G2cmp.SwapXY(Ltth,Utth)
     1267            Ltth,Utth = Utth,Ltth
    12671268        data['IOtth'] = [Ltth,Utth]
    12681269        self.InnerTth.SetValue("%8.2f" % (Ltth))
  • trunk/GSASIIlattice.py

    r76 r78  
    1313acosd = lambda x: 180.*np.arccos(x)/np.pi
    1414rdsq2d = lambda x,p: round(1.0/np.sqrt(x),p)
     15
     16def sec2HMS(sec):
     17    H = int(sec/3600)
     18    M = int(sec/60-H*60)
     19    S = sec-3600*H-60*M
     20    return '%d:%2d:%.2f'%(H,M,S)
    1521
    1622def fillgmat(cell):
     
    129135    return A,B
    130136
    131 
     137#reflection generation routines
     138#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
     139#           cell - a,b,c,alp,bet,gam in A & deg
     140   
     141def calc_rDsq(H,A):
     142    rdsq = H[0]*H[0]*A[0]+H[1]*H[1]*A[1]+H[2]*H[2]*A[2]+H[0]*H[1]*A[3]+H[0]*H[2]*A[4]+H[1]*H[2]*A[5]
     143    return rdsq
     144   
     145def calc_rDsqZ(H,A,Z,tth,lam):
     146    rpd = math.pi/180.
     147    rdsq = calc_rDsq(H,A)+Z*math.sin(tth*rpd)*2.0*rpd/(lam*lam)
     148    return rdsq
     149       
     150def MaxIndex(dmin,A):
     151    #finds maximum allowed hkl for given A within dmin
     152    Hmax = [0,0,0]
     153    try:
     154        cell = A2cell(A)
     155    except:
     156        cell = [1,1,1,90,90,90]
     157    for i in range(3):
     158        Hmax[i] = int(round(cell[i]/dmin))
     159    return Hmax
     160   
     161def sortHKLd(HKLd,ifreverse,ifdup):
     162    #HKLd is a list of [h,k,l,d,...]; ifreverse=True for largest d first
     163    #ifdup = True if duplicate d-spacings allowed
     164    T = []
     165    for i,H in enumerate(HKLd):
     166        if ifdup:
     167            T.append((H[3],i))
     168        else:
     169            T.append(H[3])           
     170    D = dict(zip(T,HKLd))
     171    T.sort()
     172    if ifreverse:
     173        T.reverse()
     174    X = []
     175    okey = ''
     176    for key in T:
     177        if key != okey: X.append(D[key])    #remove duplicate d-spacings
     178        okey = key
     179    return X
     180   
     181def SwapIndx(Axis,H):
     182    if Axis in [1,-1]:
     183        return H
     184    elif Axis in [2,-3]:
     185        return [H[1],H[2],H[0]]
     186    else:
     187        return [H[2],H[0],H[1]]
     188       
     189def CentCheck(Cent,H):
     190    h,k,l = H
     191    if Cent == 'A' and (k+l)%2:
     192        return False
     193    elif Cent == 'B' and (h+l)%2:
     194        return False
     195    elif Cent == 'C' and (h+k)%2:
     196        return False
     197    elif Cent == 'I' and (h+k+l)%2:
     198        return False
     199    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
     200        return False
     201    elif Cent == 'R' and (-h+k+l)%3:
     202        return False
     203    else:
     204        return True
     205                                   
     206def GenHBravais(dmin,Bravais,A):
     207    '''Generate the positionally unique powder diffraction reflections
     208    for a lattice and Bravais type'''
     209# dmin - minimum d-spacing
     210# Bravais in range(14) to indicate Bravais lattice; 0-2 cubic, 3,4 - hexagonal/trigonal,
     211# 5,6 - tetragonal, 7-10 - orthorhombic, 11,12 - monoclinic, 13 - triclinic
     212# A - as defined in calc_rDsq
     213# returns HKL = [h,k,l,d,0] sorted so d largest first
     214    import math
     215    if Bravais in [9,11]:
     216        Cent = 'C'
     217    elif Bravais in [1,5,8]:
     218        Cent = 'I'
     219    elif Bravais in [0,7]:
     220        Cent = 'F'
     221    elif Bravais in [3]:
     222        Cent = 'R'
     223    else:
     224        Cent = 'P'
     225    Hmax = MaxIndex(dmin,A)
     226    dminsq = 1./(dmin**2)
     227    HKL = []
     228    if Bravais == 13:                       #triclinic
     229        for l in range(-Hmax[2],Hmax[2]+1):
     230            for k in range(-Hmax[1],Hmax[1]+1):
     231                hmin = 0
     232                if (k < 0): hmin = 1
     233                if (k ==0 and l < 0): hmin = 1
     234                for h in range(hmin,Hmax[0]+1):
     235                    H=[h,k,l]
     236                    rdsq = calc_rDsq(H,A)
     237                    if 0 < rdsq <= dminsq:
     238                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
     239    elif Bravais in [11,12]:                #monoclinic - b unique
     240        Hmax = SwapIndx(2,Hmax)
     241        for h in range(Hmax[0]+1):
     242            for k in range(-Hmax[1],Hmax[1]+1):
     243                lmin = 0
     244                if k < 0:lmin = 1
     245                for l in range(lmin,Hmax[2]+1):
     246                    [h,k,l] = SwapIndx(-2,[h,k,l])
     247                    H = []
     248                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     249                    if H:
     250                        rdsq = calc_rDsq(H,A)
     251                        if 0 < rdsq <= dminsq:
     252                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
     253                    [h,k,l] = SwapIndx(2,[h,k,l])
     254    elif Bravais in [7,8,9,10]:            #orthorhombic
     255        for h in range(Hmax[0]+1):
     256            for k in range(Hmax[1]+1):
     257                for l in range(Hmax[2]+1):
     258                    H = []
     259                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     260                    if H:
     261                        rdsq = calc_rDsq(H,A)
     262                        if 0 < rdsq <= dminsq:
     263                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
     264    elif Bravais in [5,6]:                  #tetragonal
     265        for l in range(Hmax[2]+1):
     266            for k in range(Hmax[1]+1):
     267                for h in range(k,Hmax[0]+1):
     268                    H = []
     269                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     270                    if H:
     271                        rdsq = calc_rDsq(H,A)
     272                        if 0 < rdsq <= dminsq:
     273                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
     274    elif Bravais in [3,4]:
     275        lmin = 0
     276        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
     277        for l in range(lmin,Hmax[2]+1):
     278            for k in range(Hmax[1]+1):
     279                hmin = k
     280                if l < 0: hmin += 1
     281                for h in range(hmin,Hmax[0]+1):
     282                    H = []
     283                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     284                    if H:
     285                        rdsq = calc_rDsq(H,A)
     286                        if 0 < rdsq <= dminsq:
     287                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
     288
     289    else:                                   #cubic
     290        for l in range(Hmax[2]+1):
     291            for k in range(l,Hmax[1]+1):
     292                for h in range(k,Hmax[0]+1):
     293                    H = []
     294                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     295                    if H:
     296                        rdsq = calc_rDsq(H,A)
     297                        if 0 < rdsq <= dminsq:
     298                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
     299    return sortHKLd(HKL,True,False)
     300   
     301def GenHLaue(dmin,Laue,Cent,Axis,A):
     302    '''Generate the crystallographically unique powder diffraction reflections
     303    for a lattice and Bravais type
     304    '''
     305# dmin - minimum d-spacing
     306# Laue - Laue group symbol = '-1','2/m','mmmm','4/m','6/m','4/mmm','6/mmm',
     307#                            '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
     308# Cent - lattice centering = 'P','A','B','C','I','F'
     309# Axis - code for unique monoclinic axis = 'a','b','c'
     310# A - 6 terms as defined in calc_rDsq
     311# returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique
     312# part of reciprocal space ignoring anomalous dispersion
     313    import math
     314    Hmax = MaxIndex(dmin,A)
     315    dminsq = 1./(dmin**2)
     316    HKL = []
     317    if Laue == '-1':                       #triclinic
     318        for l in range(-Hmax[2],Hmax[2]+1):
     319            for k in range(-Hmax[1],Hmax[1]+1):
     320                hmin = 0
     321                if (k < 0) or (k ==0 and l < 0): hmin = 1
     322                for h in range(hmin,Hmax[0]+1):
     323                    H = []
     324                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     325                    rdsq = calc_rDsq(H,A)
     326                    if 0 < rdsq <= dminsq:
     327                        HKL.append([h,k,l,1/math.sqrt(rdsq)])
     328    elif Laue == '2/m':                #monoclinic
     329        Hmax = SwapIndx(Axis,Hmax)
     330        for h in range(Hmax[0]+1):
     331            for k in range(-Hmax[1],Hmax[1]+1):
     332                lmin = 0
     333                if k < 0:lmin = 1
     334                for l in range(lmin,Hmax[2]+1):
     335                    [h,k,l] = SwapIndx(-Axis,[h,k,l])
     336                    H = []
     337                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     338                    if H:
     339                        rdsq = calc_rDsq(H,A)
     340                        if 0 < rdsq <= dminsq:
     341                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
     342                    [h,k,l] = SwapIndx(Axis,[h,k,l])
     343    elif Laue in ['mmm','4/m','6/m']:            #orthorhombic
     344        for l in range(Hmax[2]+1):
     345            for h in range(Hmax[0]+1):
     346                kmin = 1
     347                if Laue == 'mmm' or h ==0: kmin = 0
     348                for k in range(kmin,Hmax[2]+1):
     349                    H = []
     350                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     351                    if H:
     352                        rdsq = calc_rDsq(H,A)
     353                        if 0 < rdsq <= dminsq:
     354                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
     355    elif Laue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
     356        for l in range(Hmax[2]+1):
     357            for h in range(Hmax[0]+1):
     358                for k in range(h+1):
     359                    H = []
     360                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     361                    if H:
     362                        rdsq = calc_rDsq(H,A)
     363                        if 0 < rdsq <= dminsq:
     364                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
     365    elif Laue in ['3m1','31m','3','3R','3mR']:                  #trigonals
     366        for l in range(-Hmax[2],Hmax[2]+1):
     367            hmin = 0
     368            if l < 0: hmin = 1
     369            for h in range(hmin,Hmax[0]+1):
     370                if Laue in ['3R','3']:
     371                    kmax = h
     372                    kmin = -int((h-1)/2.)
     373                else:
     374                    kmin = 0
     375                    kmax = h
     376                    if Laue in ['3m1','3mR'] and l < 0: kmax = h-1
     377                    if Laue == '31m' and l < 0: kmin = 1
     378                for k in range(kmin,kmax+1):
     379                    H = []
     380                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     381                    if H:
     382                        rdsq = calc_rDsq(H,A)
     383                        if 0 < rdsq <= dminsq:
     384                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
     385    else:                                   #cubic
     386        for h in range(Hmax[0]+1):
     387            for k in range(h+1):
     388                lmin = 0
     389                lmax = k
     390                if Laue =='m3':
     391                    lmax = h-1
     392                    if h == k: lmax += 1
     393                for l in range(lmin,lmax+1):
     394                    H = []
     395                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
     396                    if H:
     397                        rdsq = calc_rDsq(H,A)
     398                        if 0 < rdsq <= dminsq:
     399                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
     400    return sortHKLd(HKL,True,True)
     401   
    132402# output from uctbx computed on platform darwin on 2010-05-28
    133403array = np.array
Note: See TracChangeset for help on using the changeset viewer.