source: trunk/GSASIIcomp.py @ 53

Last change on this file since 53 was 53, checked in by vondreel, 13 years ago

finish image integration
problems still with motion, picking pts., etc. on patterns
File size: 57.6 KB
Line 
1#GSASII computational module
2import sys
3import math
4import wx
5import time
6import numpy as np
7import numpy.linalg as nl
8import os.path as ospath
9# determine a binary path pased on the host OS and the python version, path is relative to
10# location of this file
11if sys.platform == "win32":
12    bindir = 'binwin%d.%d' % sys.version_info[0:2]
13elif sys.platform == "darwin":
14    bindir = 'binmac%d.%d' % sys.version_info[0:2]
15else:
16    bindir = 'bin'
17
18if ospath.exists(ospath.join(sys.path[0],bindir)) and ospath.join(sys.path[0],bindir) not in sys.path: 
19    sys.path.insert(0,ospath.join(sys.path[0],bindir))
20
21try: 
22    import pypowder as pyp
23except:
24    # create an app to display the error, since we are still loading routines at this point
25    app = wx.App()
26    app.MainLoop()
27    msg = wx.MessageDialog(None, message="Unable to load the GSAS powder computation module, pypowder",
28                     caption="Import Error",
29                     style=wx.ICON_ERROR | wx.OK | wx.STAY_ON_TOP)
30    msg.ShowModal()
31    # this error is non-recoverable, so just quit
32    exit()
33import GSASIIplot as G2plt
34
35# trig functions in degrees
36sind = lambda x: math.sin(x*math.pi/180.)
37asind = lambda x: 180.*math.asin(x)/math.pi
38tand = lambda x: math.tan(x*math.pi/180.)
39atand = lambda x: 180.*math.atan(x)/math.pi
40atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
41cosd = lambda x: math.cos(x*math.pi/180.)
42acosd = lambda x: 180.*math.acos(x)/math.pi
43rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
44#numpy versions
45npsind = lambda x: np.sin(x*np.pi/180.)
46npasind = lambda x: 180.*np.arcsin(x)/math.pi
47npcosd = lambda x: np.cos(x*math.pi/180.)
48nptand = lambda x: np.tan(x*math.pi/180.)
49npatand = lambda x: 180.*np.arctan(x)/np.pi
50npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
51
52def sec2HMS(sec):
53    H = int(sec/3600)
54    M = int(sec/60-H*60)
55    S = sec-3600*H-60*M
56    return '%d:%2d:%.2f'%(H,M,S)
57
58def ValEsd(value,esd=0,nTZ=False):
59    # returns value(esd) string; nTZ=True for no trailing zeros
60    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
61    #get the 2 significant digits in the esd
62    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
63    #get the number of digits to represent them
64    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
65   
66    mdec = lambda esd: -int(math.log10(abs(esd)))
67    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
68    if esd > 0:
69        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
70        print fmt,ndec(esd),esd*10**(mdec(esd)+1)
71        return fmt%(value,int(esd*10**(mdec(esd)+1)))
72    elif esd < 0:
73         return str(round(value,mdec(esd)))
74    else:
75        text = "%F"%(value)
76        if nTZ:
77            return text.rstrip('0')
78        else:
79            return text
80
81       
82       
83
84def DoPeakFit(peaks,background,limits,inst,data):
85   
86    def backgroundPrint(background,sigback):
87        if background[1]:
88            print 'Background coefficients for',background[0],'function'
89            ptfmt = "%12.5f"
90            ptstr =  'values:'
91            sigstr = 'esds  :'
92            for i,back in enumerate(background[3:]):
93                ptstr += ptfmt % (back)
94                sigstr += ptfmt % (sigback[i+3])
95            print ptstr
96            print sigstr
97        else:
98            print 'Background not refined'
99   
100    def instPrint(instVal,siginst,insLabels):
101        print 'Instrument Parameters:'
102        ptfmt = "%12.6f"
103        ptlbls = 'names :'
104        ptstr =  'values:'
105        sigstr = 'esds  :'
106        for i,value in enumerate(instVal):
107            ptlbls += "%s" % (insLabels[i].center(12))
108            ptstr += ptfmt % (value)
109            if siginst[i]:
110                sigstr += ptfmt % (siginst[i])
111            else:
112                sigstr += 12*' '
113        print ptlbls
114        print ptstr
115        print sigstr
116   
117    def peaksPrint(peaks,sigpeaks):
118        print 'Peak list='
119
120    begin = time.time()
121    Np = len(peaks[0])
122    DataType = inst[1][0]
123    instVal = inst[1][1:]
124    Insref = inst[2][1:]
125    insLabels = inst[3][1:]
126    Ka2 = False
127    Ioff = 3
128    if len(instVal) == 12:
129        lamratio = instVal[1]/instVal[0]
130        Ka2 = True
131        Ioff = 5
132    insref = Insref[len(Insref)-7:-1]               #just U,V,W,X,Y,SH/L
133    for peak in peaks:
134        dip = []
135        dip.append(tand(peak[0]/2.0)**2)
136        dip.append(tand(peak[0]/2.0))
137        dip.append(1.0/cosd(peak[0]/2.0))
138        dip.append(tand(peak[0]/2.0))
139        peak.append(dip)
140    B = background[2]
141    bcof = background[3:3+B]
142    Bv = 0
143    if background[1]:
144        Bv = B
145    x,y,w,yc,yb,yd = data               #these are numpy arrays!
146    V = []
147    A = []
148    swobs = 0.0
149    smin = 0.0
150    first = True
151    GoOn = True
152    Go = True
153    dlg = wx.ProgressDialog("Elapsed time","Fitting peaks to pattern",len(x), \
154        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
155    screenSize = wx.DisplaySize()
156    Size = dlg.GetSize()
157    dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
158    try:
159        i = 0
160        for xi in x:
161            Go = dlg.Update(i)[0]
162            if GoOn:
163                GoOn = Go
164            if limits[0] <= xi <= limits[1]:
165                yb[i] = 0.0
166                dp = []
167                for j in range(B):
168                    t = (xi-limits[0])**j
169                    yb[i] += t*bcof[j]
170                    if background[1]:
171                        dp.append(t)
172                yc[i] = yb[i]
173                Iv = 0
174                for j in range(6):
175                    if insref[j]:
176                        dp.append(0.0)
177                        Iv += 1
178                for peak in peaks:
179                    dip = peak[-1]
180                    f = pyp.pypsvfcj(peak[2],xi-peak[0],peak[0],peak[4],peak[6],peak[8],0.0)
181                    yc[i] += f[0]*peak[2]
182                    if f[0] > 0.0:
183                        j = 0
184                        if insref[0]:              #U
185                            dp[Bv+j] += f[3]*dip[0]
186                            j += 1
187                        if insref[1]:              #V
188                            dp[Bv+j] += f[3]*dip[1]
189                            j += 1
190                        if insref[2]:              #W
191                            dp[Bv+j] += f[3]
192                            j += 1
193                        if insref[3]:              #X
194                            dp[Bv+j] += f[4]*dip[2]
195                            j += 1
196                        if insref[4]:              #Y
197                            dp[Bv+j] += f[4]*dip[3]
198                            j += 1
199                        if insref[5]:              #SH/L
200                            dp[Bv+j] += f[5]
201                    if Ka2:
202                       pos2 = 2.0*asind(lamratio*sind(peak[0]/2.0))
203                       f2 = pyp.pypsvfcj(peak[2],xi-pos2,peak[0],peak[4],peak[6],peak[8],0.0)
204                       yc[i] += f2[0]*peak[2]*instVal[3]
205                       if f[0] > 0.0:
206                           j = 0
207                           if insref[0]:              #U
208                               dp[Bv+j] += f2[3]*dip[0]*instVal[3]
209                               j += 1
210                           if insref[1]:              #V
211                               dp[Bv+j] += f2[3]*dip[1]*instVal[3]
212                               j += 1
213                           if insref[2]:              #W
214                               dp[Bv+j] += f2[3]*instVal[3]
215                               j += 1
216                           if insref[3]:              #X
217                               dp[Bv+j] += f2[4]*dip[2]*instVal[3]
218                               j += 1
219                           if insref[4]:              #Y
220                               dp[Bv+j] += f2[4]*dip[3]*instVal[3]
221                               j += 1
222                           if insref[5]:              #SH/L
223                               dp[Bv+j] += f2[5]*instVal[3]                       
224                    for j in range(0,Np,2):          #14s
225                        if peak[j+1]: dp.append(f[j/2+1])
226                yd[i] = y[i]-yc[i]
227                swobs += w[i]*y[i]**2
228                t2 = w[i]*yd[i]
229                smin += t2*yd[i]                 
230                if first:
231                    first = False
232                    M = len(dp)
233                    A = np.zeros(shape=(M,M))
234                    V = np.zeros(shape=(M))
235                A,V = pyp.buildmv(t2,w[i],M,dp,A,V)
236            i += 1
237    finally:
238        dlg.Destroy()
239    print GoOn
240    Rwp = smin/swobs
241    Rwp = math.sqrt(Rwp)*100.0
242    norm = np.diag(A)
243    for i,elm in enumerate(norm):
244        if elm <= 0.0:
245            print norm
246            return False,0,0,0
247    for i in xrange(len(V)):
248        norm[i] = 1.0/math.sqrt(norm[i])
249        V[i] *= norm[i]
250        a = A[i]
251        for j in xrange(len(V)):
252            a[j] *= norm[i]
253    A = np.transpose(A)
254    for i in xrange(len(V)):
255        a = A[i]
256        for j in xrange(len(V)):
257            a[j] *= norm[i]
258    b = nl.solve(A,V)
259    A = nl.inv(A)
260    sig = np.diag(A)
261    for i in xrange(len(V)):
262        b[i] *= norm[i]
263        sig[i] *= norm[i]*norm[i]
264        sig[i] = math.sqrt(abs(sig[i]))
265    sigback = [0,0,0]
266    for j in range(Bv):
267        background[j+3] += b[j]
268        sigback.append(sig[j])
269    backgroundPrint(background,sigback)
270    k = 0
271    delt = []
272    if Ka2:
273        siginst = [0,0,0,0,0]
274    else:
275        siginst = [0,0,0]
276    for j in range(6):
277        if insref[j]:
278            instVal[j+Ioff] += b[Bv+k]*0.5
279            siginst.append(sig[Bv+k])
280            delt.append(b[Bv+k]*0.5)
281            k += 1
282        else:
283            delt.append(0.0)
284            siginst.append(0.0)
285    delt.append(0.0)                    #dummies for azm
286    siginst.append(0.0)
287    instPrint(instVal,siginst,insLabels)
288    inst[1] = [DataType,]
289    for val in instVal:
290        inst[1].append(val)
291    B = Bv+Iv
292    for peak in peaks:
293        del peak[-1]                        # remove dip from end
294        delsig = delt[0]*tand(peak[0]/2.0)**2+delt[1]*tand(peak[0]/2.0)+delt[2]
295        delgam = delt[3]/cosd(peak[0]/2.0)+delt[4]*tand(peak[0]/2.0)
296        delshl = delt[5]   
297        for j in range(0,len(peak[:-1]),2):
298            if peak[j+1]: 
299                peak[j] += b[B]*0.5
300                B += 1
301        peak[4] += delsig
302        if peak[4] < 0.0:
303            print 'ERROR - negative sigma'
304            return False,0,0,0,False           
305        peak[6] += delgam
306        if peak[6] < 0.0:
307            print 'ERROR - negative gamma'
308            return False,0,0,0,False
309        peak[8] += delshl           
310        if peak[8] < 0.0:
311            print 'ERROR - negative SH/L'
312            return False,0,0,0,False           
313    runtime = time.time()-begin   
314    data = [x,y,w,yc,yb,yd]
315    return True,smin,Rwp,runtime,GoOn
316   
317#some cell utilities
318#for these: H = [h,k,l]; A is as used in calc_rDsq; G - inv metric tensor, g - metric tensor;
319#           cell - a,b,c,alp,bet,gam in A & deg
320   
321def calc_rDsq(H,A):
322    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]
323    return rdsq
324   
325def calc_rDsqZ(H,A,Z,tth,lam):
326    rpd = math.pi/180.
327    rdsq = calc_rDsq(H,A)+Z*math.sin(tth*rpd)*2.0*rpd/(lam*lam)
328    return rdsq
329   
330def calc_rVsq(A):
331    rVsq = A[0]*A[1]*A[2]+0.25*(A[3]*A[4]*A[5]-A[0]*A[5]**2-A[1]*A[4]**2-A[2]*A[3]**2)
332    if rVsq < 0:
333        return 1
334    return rVsq
335   
336def calc_rV(A):
337    return math.sqrt(calc_rVsq(A))
338   
339def calc_V(A):
340    return 1./calc_rV(A)
341   
342def scaleAbyV(A,V):
343    v = calc_V(A)
344    scale = math.exp(math.log(v/V)/3.)**2
345    for i in range(6):
346        A[i] *= scale
347   
348def ranaxis(dmin,dmax):
349    import random as rand
350    return rand.random()*(dmax-dmin)+dmin
351   
352def ran2axis(k,N):
353    import random as rand
354    T = 1.5+0.49*k/N
355#    B = 0.99-0.49*k/N
356#    B = 0.99-0.049*k/N
357    B = 0.99-0.149*k/N
358    R = (T-B)*rand.random()+B
359    return R
360   
361def ranNaxis(k,N):
362    import random as rand
363    T = 1.0+1.0*k/N
364    B = 1.0-1.0*k/N
365    R = (T-B)*rand.random()+B
366    return R
367   
368def ranAbyV(Bravais,dmin,dmax,V):
369    cell = [0,0,0,0,0,0]
370    bad = True
371    while bad:
372        bad = False
373        cell = rancell(Bravais,dmin,dmax)
374        G,g = cell2Gmat(cell)
375        A = Gmat2A(G)
376        if calc_rVsq(A) < 1:
377            scaleAbyV(A,V)
378            cell = A2cell(A)
379            for i in range(3):
380                bad |= cell[i] < dmin
381    return A
382   
383def ranAbyR(Bravais,A,k,N,ranFunc):
384    R = ranFunc(k,N)
385    if Bravais in [0,1,2]:          #cubic - not used
386        A[0] = A[1] = A[2] = A[0]*R
387        A[3] = A[4] = A[5] = 0.
388    elif Bravais in [3,4]:          #hexagonal/trigonal
389        A[0] = A[1] = A[3] = A[0]*R
390        A[2] *= R
391        A[4] = A[5] = 0.       
392    elif Bravais in [5,6]:          #tetragonal
393        A[0] = A[1] = A[0]*R
394        A[2] *= R
395        A[3] = A[4] = A[5] = 0.       
396    elif Bravais in [7,8,9,10]:     #orthorhombic
397        A[0] *= R
398        A[1] *= R
399        A[2] *= R
400        A[3] = A[4] = A[5] = 0.       
401    elif Bravais in [11,12]:        #monoclinic
402        A[0] *= R
403        A[1] *= R
404        A[2] *= R
405        A[4] *= R
406        A[3] = A[5] = 0.       
407    else:                           #triclinic
408        A[0] *= R
409        A[1] *= R
410        A[2] *= R
411        A[3] *= R
412        A[4] *= R
413        A[5] *= R
414    return A
415   
416def rancell(Bravais,dmin,dmax):
417    if Bravais in [0,1,2]:          #cubic
418        a = b = c = ranaxis(dmin,dmax)
419        alp = bet = gam = 90
420    elif Bravais in [3,4]:          #hexagonal/trigonal
421        a = b = ranaxis(dmin,dmax)
422        c = ranaxis(dmin,dmax)
423        alp = bet =  90
424        gam = 120
425    elif Bravais in [5,6]:          #tetragonal
426        a = b = ranaxis(dmin,dmax)
427        c = ranaxis(dmin,dmax)
428        alp = bet = gam = 90
429    elif Bravais in [7,8,9,10]:       #orthorhombic - F,I,C,P - a<b<c convention
430        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
431        abc.sort()
432        a = abc[0]
433        b = abc[1]
434        c = abc[2]
435        alp = bet = gam = 90
436    elif Bravais in [11,12]:        #monoclinic - C,P - a<c convention
437        ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
438        ac.sort()
439        a = ac[0]
440        b = ranaxis(dmin,dmax)
441        c = ac[1]
442        alp = gam = 90
443        bet = ranaxis(90.,130.)
444    else:                           #triclinic - a<b<c convention
445        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
446        abc.sort()
447        a = abc[0]
448        b = abc[1]
449        c = abc[2]
450        r = 0.5*b/c
451        alp = ranaxis(acosd(r),acosd(-r))
452        r = 0.5*a/c
453        bet = ranaxis(acosd(r),acosd(-r))
454        r = 0.5*a/b
455        gam = ranaxis(acosd(r),acosd(-r)) 
456    return [a,b,c,alp,bet,gam]
457   
458def A2Gmat(A):
459    G = np.zeros(shape=(3,3))
460    G = [[A[0],A[3]/2.,A[4]/2.], [A[3]/2.,A[1],A[5]/2.], [A[4]/2.,A[5]/2.,A[2]]]
461    g = nl.inv(G)
462    return G,g
463   
464def fillgmat(cell):
465    a,b,c,alp,bet,gam = cell
466    g = np.array([[a*a,a*b*cosd(gam),a*c*cosd(bet)],[a*b*cosd(gam),b*b,b*c*cosd(alp)],
467        [a*c*cosd(bet),b*c*cosd(alp),c*c]])
468    return g
469           
470def cell2Gmat(cell):
471    #returns reciprocal (G) & real (g) metric tensors
472    g = fillgmat(cell)
473    G = nl.inv(g)       
474    return G,g
475   
476def invcell2Gmat(invcell):
477    G = fillgmat(invcell)
478    g = nl.inv(G)
479    return G,g
480   
481def Gmat2cell(g):
482    #returns lattice parameters from real metric tensor (g)
483    a = math.sqrt(max(0,g[0][0]))
484    b = math.sqrt(max(0,g[1][1]))
485    c = math.sqrt(max(0,g[2][2]))
486    alp = acosd(g[2][1]/(b*c))
487    bet = acosd(g[2][0]/(a*c))
488    gam = acosd(g[0][1]/(a*b))
489    return a,b,c,alp,bet,gam
490   
491def Gmat2A(G):
492    return [G[0][0],G[1][1],G[2][2],2.*G[0][1],2.*G[0][2],2.*G[1][2]]
493   
494def cell2A(cell):
495    G,g = cell2Gmat(cell)
496    return Gmat2A(G)
497   
498def A2cell(A):
499    G,g = A2Gmat(A)
500    return Gmat2cell(g)
501   
502def A2invcell(A):
503    ainv = math.sqrt(max(0.,A[0]))
504    binv = math.sqrt(max(0.,A[1]))
505    cinv = math.sqrt(max(0.,A[2]))
506    gaminv = acosd(max(-0.5,min(0.5,0.5*A[3]/(ainv*binv))))
507    betinv = acosd(max(-0.5,min(0.5,0.5*A[4]/(ainv*cinv))))
508    alpinv = acosd(max(-0.5,min(0.5,0.5*A[5]/(binv*cinv))))
509    return ainv,binv,cinv,alpinv,betinv,gaminv
510   
511def cell2AB(cell):
512    #from real lattice parameters - cell
513    # returns A for Cartesian to crystal transformations A*X = x
514    # and inverse B for crystal to Cartesian transformation B*x = X
515    G,g = cell2Gmat(cell)       #reciprocal & real metric tensors
516    cosAlpStar = G[2][1]/math.sqrt(G[1][1]*G[2][2])
517    sinAlpStar = math.sqrt(1.0-cosAlpStar**2)
518    B = np.eye(3)
519    B *= cell[:3]
520    A = np.zeros(shape=(3,3))
521    A[0][0] = 1.0
522    A[0][1] = cosd(cell[5])
523    A[1][1] = sinAlpStar*sind(cell[5])
524    A[1][2] = -cosAlpStar*sind(cell[5])
525    A[0][2] = cosd(cell[4])
526    A[2][2] = sind(cell[4])
527    B = np.dot(A,B)
528    A = nl.inv(B)
529    return A,B
530   
531def makeMat(Angle,Axis):
532    #Make rotation matrix from Angle in degrees,Axis =0 for rotation about x, =1 for about y, etc.
533    cs = cosd(Angle)
534    ss = sind(Angle)
535    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
536    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
537                   
538def MaxIndex(dmin,A):
539    #finds maximum allowed hkl for given A within dmin
540    Hmax = [0,0,0]
541    try:
542        cell = A2cell(A)
543    except:
544        cell = [1,1,1,90,90,90]
545    for i in range(3):
546        Hmax[i] = int(round(cell[i]/dmin))
547    return Hmax
548   
549def sortHKLd(HKLd,ifreverse,ifdup):
550    #HKLd is a list of [h,k,l,d,...]; ifreverse=True for largest d first
551    #ifdup = True if duplicate d-spacings allowed
552    T = []
553    for i,H in enumerate(HKLd):
554        if ifdup:
555            T.append((H[3],i))
556        else:
557            T.append(H[3])           
558    D = dict(zip(T,HKLd))
559    T.sort()
560    if ifreverse:
561        T.reverse()
562    X = []
563    okey = ''
564    for key in T: 
565        if key != okey: X.append(D[key])    #remove duplicate d-spacings
566        okey = key
567    return X
568   
569def sortM20(cells):
570    #cells is M20,X20,Bravais,a,b,c,alp,bet,gam
571    #sort highest M20 1st
572    T = []
573    for i,M in enumerate(cells):
574        T.append((M[0],i))
575    D = dict(zip(T,cells))
576    T.sort()
577    T.reverse()
578    X = []
579    for key in T:
580        X.append(D[key])
581    return X
582               
583 
584def GenHBravais(dmin,Bravais,A):
585# dmin - minimum d-spacing
586# Bravais in range(14) to indicate Bravais lattice; 0-2 cubic, 3,4 - hexagonal/trigonal,
587# 5,6 - tetragonal, 7-10 - orthorhombic, 11,12 - monoclinic, 13 - triclinic
588# A - as defined in calc_rDsq
589# returns HKL = [h,k,l,d,0] sorted so d largest first
590    import math
591    if Bravais in [9,11]:
592        Cent = 'C'
593    elif Bravais in [1,5,8]:
594        Cent = 'I'
595    elif Bravais in [0,7]:
596        Cent = 'F'
597    elif Bravais in [3]:
598        Cent = 'R'
599    else:
600        Cent = 'P'
601    Hmax = MaxIndex(dmin,A)
602    dminsq = 1./(dmin**2)
603    HKL = []
604    if Bravais == 13:                       #triclinic
605        for l in range(-Hmax[2],Hmax[2]+1):
606            for k in range(-Hmax[1],Hmax[1]+1):
607                hmin = 0
608                if (k < 0): hmin = 1
609                if (k ==0 and l < 0): hmin = 1
610                for h in range(hmin,Hmax[0]+1):
611                    H=[h,k,l]
612                    rdsq = calc_rDsq(H,A)
613                    if 0 < rdsq <= dminsq:
614                        HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
615    elif Bravais in [11,12]:                #monoclinic - b unique
616        Hmax = SwapIndx(2,Hmax)
617        for h in range(Hmax[0]+1):
618            for k in range(-Hmax[1],Hmax[1]+1):
619                lmin = 0
620                if k < 0:lmin = 1
621                for l in range(lmin,Hmax[2]+1):
622                    [h,k,l] = SwapIndx(-2,[h,k,l])
623                    H = []
624                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
625                    if H:
626                        rdsq = calc_rDsq(H,A)
627                        if 0 < rdsq <= dminsq:
628                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
629                    [h,k,l] = SwapIndx(2,[h,k,l])
630    elif Bravais in [7,8,9,10]:            #orthorhombic
631        for h in range(Hmax[0]+1):
632            for k in range(Hmax[1]+1):
633                for l in range(Hmax[2]+1):
634                    H = []
635                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
636                    if H:
637                        rdsq = calc_rDsq(H,A)
638                        if 0 < rdsq <= dminsq:
639                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
640    elif Bravais in [5,6]:                  #tetragonal
641        for l in range(Hmax[2]+1):
642            for k in range(Hmax[1]+1):
643                for h in range(k,Hmax[0]+1):
644                    H = []
645                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
646                    if H:
647                        rdsq = calc_rDsq(H,A)
648                        if 0 < rdsq <= dminsq:
649                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
650    elif Bravais in [3,4]:
651        lmin = 0
652        if Bravais == 3: lmin = -Hmax[2]                  #hexagonal/trigonal
653        for l in range(lmin,Hmax[2]+1):
654            for k in range(Hmax[1]+1):
655                hmin = k
656                if l < 0: hmin += 1
657                for h in range(hmin,Hmax[0]+1):
658                    H = []
659                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
660                    if H:
661                        rdsq = calc_rDsq(H,A)
662                        if 0 < rdsq <= dminsq:
663                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
664
665    else:                                   #cubic
666        for l in range(Hmax[2]+1):
667            for k in range(l,Hmax[1]+1):
668                for h in range(k,Hmax[0]+1):
669                    H = []
670                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
671                    if H:
672                        rdsq = calc_rDsq(H,A)
673                        if 0 < rdsq <= dminsq:
674                            HKL.append([h,k,l,rdsq2d(rdsq,6),-1])
675    return sortHKLd(HKL,True,False)
676   
677def SwapXY(x,y):
678    return [y,x]
679   
680def SwapIndx(Axis,H):
681    if Axis in [1,-1]:
682        return H
683    elif Axis in [2,-3]:
684        return [H[1],H[2],H[0]]
685    else:
686        return [H[2],H[0],H[1]]
687       
688def CentCheck(Cent,H):
689    h,k,l = H
690    if Cent == 'A' and (k+l)%2:
691        return False
692    elif Cent == 'B' and (h+l)%2:
693        return False
694    elif Cent == 'C' and (h+k)%2:
695        return False
696    elif Cent == 'I' and (h+k+l)%2:
697        return False
698    elif Cent == 'F' and ((h+k)%2 or (h+l)%2 or (k+l)%2):
699        return False
700    elif Cent == 'R' and (-h+k+l)%3:
701        return False
702    else:
703        return True
704                                   
705def GenHLaue(dmin,Laue,Cent,Axis,A):
706# dmin - minimum d-spacing
707# Laue - Laue group symbol = '-1','2/m','mmmm','4/m','6/m','4/mmm','6/mmm',
708#                            '3m1', '31m', '3', '3R', '3mR', 'm3', 'm3m'
709# Cent - lattice centering = 'P','A','B','C','I','F'
710# Axis - code for unique monoclinic axis = 'a','b','c'
711# A - 6 terms as defined in calc_rDsq
712# returns - HKL = list of [h,k,l,d] sorted with largest d first and is unique
713# part of reciprocal space ignoring anomalous dispersion
714    import math
715    Hmax = MaxIndex(dmin,A)
716    dminsq = 1./(dmin**2)
717    HKL = []
718    if Laue == '-1':                       #triclinic
719        for l in range(-Hmax[2],Hmax[2]+1):
720            for k in range(-Hmax[1],Hmax[1]+1):
721                hmin = 0
722                if (k < 0) or (k ==0 and l < 0): hmin = 1
723                for h in range(hmin,Hmax[0]+1):
724                    H = []
725                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
726                    rdsq = calc_rDsq(H,A)
727                    if 0 < rdsq <= dminsq:
728                        HKL.append([h,k,l,1/math.sqrt(rdsq)])
729    elif Laue == '2/m':                #monoclinic
730        Hmax = SwapIndx(Axis,Hmax)
731        for h in range(Hmax[0]+1):
732            for k in range(-Hmax[1],Hmax[1]+1):
733                lmin = 0
734                if k < 0:lmin = 1
735                for l in range(lmin,Hmax[2]+1):
736                    [h,k,l] = SwapIndx(-Axis,[h,k,l])
737                    H = []
738                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
739                    if H:
740                        rdsq = calc_rDsq(H,A)
741                        if 0 < rdsq <= dminsq:
742                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
743                    [h,k,l] = SwapIndx(Axis,[h,k,l])
744    elif Laue in ['mmm','4/m','6/m']:            #orthorhombic
745        for l in range(Hmax[2]+1):
746            for h in range(Hmax[0]+1):
747                kmin = 1
748                if Laue == 'mmm' or h ==0: kmin = 0
749                for k in range(kmin,Hmax[2]+1):
750                    H = []
751                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
752                    if H:
753                        rdsq = calc_rDsq(H,A)
754                        if 0 < rdsq <= dminsq:
755                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
756    elif Laue in ['4/mmm','6/mmm']:                  #tetragonal & hexagonal
757        for l in range(Hmax[2]+1):
758            for h in range(Hmax[0]+1):
759                for k in range(h+1):
760                    H = []
761                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
762                    if H:
763                        rdsq = calc_rDsq(H,A)
764                        if 0 < rdsq <= dminsq:
765                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
766    elif Laue in ['3m1','31m','3','3R','3mR']:                  #trigonals
767        for l in range(-Hmax[2],Hmax[2]+1):
768            hmin = 0
769            if l < 0: hmin = 1
770            for h in range(hmin,Hmax[0]+1):
771                if Laue in ['3R','3']:
772                    kmax = h
773                    kmin = -int((h-1)/2.)
774                else:
775                    kmin = 0
776                    kmax = h
777                    if Laue in ['3m1','3mR'] and l < 0: kmax = h-1
778                    if Laue == '31m' and l < 0: kmin = 1
779                for k in range(kmin,kmax+1):
780                    H = []
781                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
782                    if H:
783                        rdsq = calc_rDsq(H,A)
784                        if 0 < rdsq <= dminsq:
785                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
786    else:                                   #cubic
787        for h in range(Hmax[0]+1):
788            for k in range(h+1):
789                lmin = 0
790                lmax = k
791                if Laue =='m3':
792                    lmax = h-1
793                    if h == k: lmax += 1
794                for l in range(lmin,lmax+1):
795                    H = []
796                    if CentCheck(Cent,[h,k,l]): H=[h,k,l]
797                    if H:
798                        rdsq = calc_rDsq(H,A)
799                        if 0 < rdsq <= dminsq:
800                            HKL.append([h,k,l,1/math.sqrt(rdsq)])
801    return sortHKLd(HKL,True,True)
802   
803def calc_M20(peaks,HKL):
804    diff = 0
805    X20 = 0
806    for Nobs20,peak in enumerate(peaks):
807        if peak[3]:
808            Qobs = 1.0/peak[7]**2
809            Qcalc = 1.0/peak[8]**2
810            diff += abs(Qobs-Qcalc)
811        elif peak[2]:
812            X20 += 1
813        if Nobs20 == 19: 
814            d20 = peak[7]
815            break
816    else:
817        d20 = peak[7]
818        Nobs20 = len(peaks)
819    for N20,hkl in enumerate(HKL):
820        if hkl[3] < d20:
821            break               
822    eta = diff/Nobs20
823    Q20 = 1.0/d20**2
824    if diff:
825        M20 = Q20/(2.0*diff)
826    else:
827        M20 = 0
828    M20 /= (1.+X20)
829    return M20,X20
830   
831def IndexPeaks(peaks,HKL):
832    import bisect
833    hklds = [1000.0]                                    #make bounded list of available d-spacings
834    N = len(HKL)
835    if N == 0: return False
836    for hkl in HKL:
837        hklds.append(hkl[3])
838    hklds.append(0.0)
839    hklds.sort()                                        # ascending sort - upper bound at end
840    hklmax = [0,0,0]
841    for ipk,peak in enumerate(peaks):
842        if peak[2]:
843            i = bisect.bisect_right(hklds,peak[7])          # find peak position in hkl list
844            dm = peak[7]-hklds[i-1]                         # peak to neighbor hkls in list
845            dp = hklds[i]-peak[7]
846            pos = N-i                                       # reverse the order
847            if dp > dm: pos += 1                            # closer to upper than lower
848            hkl = HKL[pos]                                 # put in hkl
849            if hkl[4] >= 0:                                 # peak already assigned - test if this one better
850                opeak = peaks[hkl[4]]
851                dold = abs(opeak[7]-hkl[3])
852                dnew = min(dm,dp)
853                if dold > dnew:                             # new better - zero out old
854                    for j in range(3):
855                        opeak[j+4] = 0
856                    opeak[8] = 0.
857                else:                                       # old better - do nothing
858                    continue               
859            hkl[4] = ipk
860            for j in range(3):
861                peak[j+4] = hkl[j]
862            peak[8] = hkl[3]                                # fill in d-calc
863    for peak in peaks:
864        peak[3] = False
865        if peak[2]:
866            if peak[8] > 0.:
867                for j in range(3):
868                    if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
869                peak[3] = True
870    if hklmax[0]*hklmax[1]*hklmax[2] > 0:
871        return True
872    else:
873        return False
874   
875def FitHKL(ibrav,peaks,A,wtP):
876    def ShiftTest(a,b):
877        if b < -0.1*a: 
878            b = -0.0001*a
879        return b
880    smin = 0.
881    first = True
882    for peak in peaks:
883        if peak[2] and peak[3]:
884            h,k,l = H = peak[4:7]
885            Qo = 1./peak[7]**2
886            Qc = calc_rDsq(H,A)
887            try:
888                peak[8] = 1./math.sqrt(Qc)
889            except:
890                print A2invcell(A)
891            delt = Qo-Qc
892            smin += delt**2
893            dp = []
894            if ibrav in [0,1,2]:            #m3m
895                dp.append(h*h+k*k+l*l)
896            elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
897                dp.append(h*h+k*k+h*k)
898                dp.append(l*l)
899            elif ibrav in [5,6]:            #4/mmm
900                dp.append(h*h+k*k)
901                dp.append(l*l)
902            elif ibrav in [7,8,9,10]:       #mmm
903                dp.append(h*h)
904                dp.append(k*k)
905                dp.append(l*l)
906            elif ibrav in [11,12]:          #2/m
907                dp.append(h*h)
908                dp.append(k*k)
909                dp.append(l*l)
910                dp.append(h*l)
911            else:                           #1
912#    derivatives for a0*h^2+a1*k^2+a2*l^2+a3*h*k+a4*h*l+a5*k*l
913                dp.append(h*h)
914                dp.append(k*k)
915                dp.append(l*l)
916                dp.append(h*k)
917                dp.append(h*l)
918                dp.append(k*l)
919            if first:
920                first = False
921                M = len(dp)
922                B = np.zeros(shape=(M,M))
923                V = np.zeros(shape=(M))
924            dwt = peak[7]**wtP
925            B,V = pyp.buildmv(delt*dwt,dwt,M,dp,B,V)
926    if nl.det(B) > 0.0:
927        try:
928            b = nl.solve(B,V)
929            B = nl.inv(B)
930            sig = np.diag(B)
931        except SingularMatrix:
932            return False,0
933        if ibrav in [0,1,2]:            #m3m
934            A[0] += ShiftTest(A[0],b[0])
935            A[1] = A[2] = A[0]
936        elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
937            A[0] += ShiftTest(A[0],b[0])
938            A[1] = A[3] = A[0]
939            A[2] += ShiftTest(A[2],b[1])
940        elif ibrav in [5,6]:            #4/mmm
941            A[0] += ShiftTest(A[0],b[0])
942            A[1] = A[0]
943            A[2] += ShiftTest(A[2],b[1])
944        elif ibrav in [7,8,9,10]:       #mmm
945            for i in range(3):
946                A[i] += ShiftTest(A[i],b[i])
947        elif ibrav in [11,12]:          #2/m
948            for i in range(3):
949                A[i] += ShiftTest(A[i],b[i])
950            A[4] += ShiftTest(A[4],b[3])
951            A[4] = min(1.4*math.sqrt(A[0]*A[2]),A[4])   #min beta star = 45
952        else:                           #1
953            oldV = math.sqrt(1./calc_rVsq(A))
954            oldA = A[:]
955            for i in range(6):
956                A[i] += b[i]*0.2
957            A[3] = min(1.1*math.sqrt(max(0,A[1]*A[2])),A[3])
958            A[4] = min(1.1*math.sqrt(max(0,A[0]*A[2])),A[4])
959            A[5] = min(1.1*math.sqrt(max(0,A[0]*A[1])),A[5])
960            ratio = math.sqrt(1./calc_rVsq(A))/oldV
961            if 0.9 > ratio or ratio > 1.1:
962                A = oldA
963#    else:
964#        print B
965#        print V
966#        for peak in peaks:
967#            print peak
968    return True,smin
969       
970def FitHKLZ(ibrav,peaks,Z,A):
971    return A,Z
972   
973def rotOrthoA(A):
974    return [A[1],A[2],A[0],0,0,0]
975   
976def swapMonoA(A):
977    return [A[2],A[1],A[0],0,A[4],0]
978   
979def oddPeak(indx,peaks):
980    noOdd = True
981    for peak in peaks:
982        H = peak[4:7]
983        if H[indx] % 2:
984            noOdd = False
985    return noOdd
986   
987def halfCell(ibrav,A,peaks):
988    if ibrav in [0,1,2]:
989        if oddPeak(0,peaks):
990            A[0] *= 2
991            A[1] = A[2] = A[0]
992    elif ibrav in [3,4,5,6]:
993        if oddPeak(0,peaks):
994            A[0] *= 2
995            A[1] = A[0]
996        if oddPeak(2,peaks):
997            A[2] *=2
998    else:
999        if oddPeak(0,peaks):
1000            A[0] *=2
1001        if oddPeak(1,peaks):
1002            A[1] *=2
1003        if oddPeak(2,peaks):
1004            A[2] *=2
1005    return A
1006   
1007def getDmin(peaks):
1008    return peaks[-1][7]
1009   
1010def getDmax(peaks):
1011    return peaks[0][7]
1012   
1013def refinePeaks(peaks,ibrav,A):
1014    dmin = getDmin(peaks)
1015    smin = 1.0e10
1016    pwr = 6
1017    maxTries = 10
1018    if ibrav == 13:
1019        pwr = 4
1020        maxTries = 10
1021    OK = False
1022    tries = 0
1023    HKL = GenHBravais(dmin,ibrav,A)
1024    while IndexPeaks(peaks,HKL):
1025        Pwr = pwr - 2*(tries % 2)
1026        HKL = []
1027        tries += 1
1028        osmin = smin
1029        oldA = A
1030        OK,smin = FitHKL(ibrav,peaks,A,Pwr)
1031        for a in A[:3]:
1032            if a < 0:
1033                A = oldA
1034                OK = False
1035                break
1036        if OK:
1037            HKL = GenHBravais(dmin,ibrav,A)
1038        if len(HKL) == 0: break                         #absurd cell obtained!
1039        rat = (osmin-smin)/smin
1040        if abs(rat) < 1.0e-5 or not OK: break
1041        if tries > maxTries: break
1042    if OK:
1043        OK,smin = FitHKL(ibrav,peaks,A,4)
1044    M20,X20 = calc_M20(peaks,HKL)
1045    return len(HKL),M20,X20
1046       
1047def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
1048# dlg & ncMax are used for wx progress bar
1049# A != 0 find the best A near input A,
1050# A = 0 for random cell, volume normalized to V1;
1051# returns number of generated hkls, M20, X20 & A for best found
1052    mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10]
1053    dmin = getDmin(peaks)-0.05
1054    amin = 2.5
1055    amax = 5.*getDmax(peaks)
1056    Asave = []
1057    GoOn = True
1058    if A:
1059        HKL = GenHBravais(dmin,ibrav,A[:])
1060        if len(HKL) > mHKL[ibrav]:
1061            IndexPeaks(peaks,HKL)
1062            Asave.append([calc_M20(peaks,HKL),A[:]])
1063    tries = 0
1064    while tries < Ntries:
1065        if A:
1066            Anew = ranAbyR(ibrav,A[:],tries+1,Ntries,ran2axis)
1067            if ibrav in [11,12,13]:
1068                Anew = ranAbyR(ibrav,A[:],tries/10+1,Ntries,ran2axis)
1069        else:
1070            Anew = ranAbyV(ibrav,amin,amax,V1)
1071        HKL = GenHBravais(dmin,ibrav,Anew)
1072        if len(HKL) > mHKL[ibrav] and IndexPeaks(peaks,HKL):
1073            Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
1074            Asave.append([calc_M20(peaks,HKL),Anew[:]])
1075            if ibrav == 9:                          #C-centered orthorhombic
1076                for i in range(2):
1077                    Anew = rotOrthoA(Anew[:])
1078                    Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
1079                    HKL = GenHBravais(dmin,ibrav,Anew)
1080                    IndexPeaks(peaks,HKL)
1081                    Asave.append([calc_M20(peaks,HKL),Anew[:]])
1082            elif ibrav == 11:                      #C-centered monoclinic
1083                Anew = swapMonoA(Anew[:])
1084                Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
1085                HKL = GenHBravais(dmin,ibrav,Anew)
1086                IndexPeaks(peaks,HKL)
1087                Asave.append([calc_M20(peaks,HKL),Anew[:]])
1088        else:
1089            break
1090        Nc = len(HKL)
1091        if Nc >= ncMax:
1092            GoOn = False
1093        elif dlg:
1094            GoOn = dlg.Update(Nc)[0]
1095            if not GoOn:
1096                break
1097        tries += 1   
1098    X = sortM20(Asave)
1099    if X:
1100        Lhkl,M20,X20 = refinePeaks(peaks,ibrav,X[0][1])
1101        return GoOn,Lhkl,M20,X20,X[0][1]
1102    else:
1103        return GoOn,0,0,0,0
1104       
1105def monoCellReduce(ibrav,A):
1106    a,b,c,alp,bet,gam = A2cell(A)
1107    G,g = A2Gmat(A)
1108    if ibrav in [11]:
1109        u = [0,0,-1]
1110        v = [1,0,2]
1111        anew = math.sqrt(np.dot(np.dot(v,g),v))
1112        if anew < a:
1113            cang = np.dot(np.dot(u,g),v)/(anew*c)
1114            beta = acosd(-abs(cang))
1115            A = cell2A([anew,b,c,90,beta,90])
1116    else:
1117        u = [-1,0,0]
1118        v = [1,0,1]
1119        cnew = math.sqrt(np.dot(np.dot(v,g),v))
1120        if cnew < c:
1121            cang = np.dot(np.dot(u,g),v)/(a*cnew)
1122            beta = acosd(-abs(cang))
1123            A = cell2A([a,b,cnew,90,beta,90])
1124    return A
1125
1126def DoIndexPeaks(peaks,inst,controls,bravais):
1127   
1128    def peakDspace(peaks,A):
1129        for peak in peaks:
1130            if peak[3]:
1131                dsq = calc_rDsq(peak[4:7],A)
1132                if dsq > 0:
1133                    peak[8] = 1./math.sqrt(dsq)
1134        return
1135    delt = 0.05                                     #lowest d-spacing cushion - can be fixed?
1136    amin = 2.5
1137    amax = 5.0*getDmax(peaks)
1138    dmin = getDmin(peaks)-delt
1139    bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
1140        'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C',
1141        'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic']
1142    tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
1143    N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,  50,50, 200]
1144    N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,   2,2,   4]
1145    Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,   2,2,   4]
1146    Nobs = len(peaks)
1147    wave = inst[1]
1148    if len(inst) > 10:
1149        zero = inst[3]
1150    else:
1151        zero = inst[2]
1152    print "%s %8.5f %6.3f" % ('wavelength, zero =',wave,zero)
1153    print "%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax)
1154    ifzero,maxzero,ncno = controls[:3]
1155    ncMax = Nobs*ncno
1156    print "%s %d %s %d %s %d" % ('change zero =',ifzero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs)
1157    cells = []
1158    for ibrav in range(14):
1159        begin = time.time()
1160        if bravais[ibrav]:
1161            print 'cell search for ',bravaisNames[ibrav]
1162            print '      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test'
1163            V1 = controls[3]
1164            bestM20 = 0
1165            topM20 = 0
1166            cycle = 0
1167            while cycle < 5:
1168                dlg = wx.ProgressDialog("Generated reflections",tries[cycle]+" cell search for "+bravaisNames[ibrav],ncMax, \
1169                    style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
1170                screenSize = wx.DisplaySize()
1171                Size = dlg.GetSize()
1172                dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
1173                try:
1174                    GoOn = True
1175                    while GoOn:                                                 #Loop over increment of volume
1176                        N2 = 0
1177                        while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
1178                            if ibrav > 2:
1179                                if not N2:
1180                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
1181                                if A:
1182                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
1183                            else:
1184                                GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
1185                            if Nc >= ncMax:
1186                                GoOn = False
1187                                break
1188                            elif 3*Nc < Nobs:
1189                                N2 = 10
1190                                break
1191                            else:
1192                                if not GoOn:
1193                                    break
1194                                if M20 > 1.0:
1195                                    bestM20 = max(bestM20,M20)
1196                                    A = halfCell(ibrav,A[:],peaks)
1197                                    if ibrav in [12]:
1198                                        A = monoCellReduce(ibrav,A[:])
1199                                    HKL = GenHBravais(dmin,ibrav,A)
1200                                    IndexPeaks(peaks,HKL)
1201                                    a,b,c,alp,bet,gam = A2cell(A)
1202                                    V = calc_V(A)
1203                                    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)
1204                                    if M20 >= 2.0:
1205                                        cells.append([M20,X20,ibrav,a,b,c,alp,bet,gam,V,False])
1206                            if not GoOn:
1207                                break
1208                            N2 += 1
1209                        if ibrav < 11:
1210                            V1 *= 1.1
1211                        elif ibrav in range(11,14):
1212                            V1 *= 1.05
1213                        if not GoOn:
1214                            if bestM20 > topM20:
1215                                topM20 = bestM20
1216                                if cells:
1217                                    V1 = cells[0][9]
1218                                else:
1219                                    V1 = 25
1220                                ncMax += Nobs
1221                                cycle += 1
1222                                print 'Restart search, new Max Nc = ',ncMax
1223                            else:
1224                                cycle = 10
1225                finally:
1226                    dlg.Destroy()
1227            print '%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
1228                ', elapsed time = ',sec2HMS(time.time()-begin))
1229           
1230    if cells:
1231        cells = sortM20(cells)
1232        cells[0][-1] = True
1233        return True,dmin,cells
1234    else:
1235        return False,0,0
1236       
1237def FitRing(ring):
1238    Err,parms = FitCircle(ring)
1239    Err /= len(ring)
1240#    print 'circle error:','%8f'%(Err)
1241    if Err > 20000.:
1242        eparms = FitEllipse(ring)
1243        if eparms:
1244            parms = eparms
1245    return parms
1246       
1247def FitCircle(ring):
1248    import numpy.linalg as nl
1249   
1250    def makeParmsCircle(B):
1251        cent = [-B[0]/2,-B[1]/2]
1252        phi = 0.
1253        sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2-B[2])
1254        return cent,phi,[sr1,sr2]
1255       
1256    ring = np.array(ring)
1257    x = np.asarray(ring.T[0])
1258    y = np.asarray(ring.T[1])
1259   
1260    M = np.array((x,y,np.ones_like(x)))
1261    B = np.array(-(x**2+y**2))
1262    result = nl.lstsq(M.T,B)
1263    return result[1],makeParmsCircle(result[0])
1264       
1265def FitEllipse(ring):
1266    import numpy.linalg as nl
1267           
1268    def makeParmsEllipse(B):
1269        det = 4.*(1.-B[0]**2)-B[1]**2
1270        if det < 0.:
1271            print 'hyperbola!'
1272            return 0
1273        elif det == 0.:
1274            print 'parabola!'
1275            return 0
1276        cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \
1277            (B[1]*B[2]-2.*(1.+B[0])*B[3])/det]
1278        phi = 0.5*atand(0.5*B[1]/B[0])
1279       
1280        a = (1.+B[0])/cosd(2*phi)
1281        b = 2.-a
1282        f = (1.+B[0])*cent[0]**2+(1.-B[0])*cent[1]**2+B[1]*cent[0]*cent[1]-B[4]
1283        if f/a < 0 or f/b < 0:
1284            return 0
1285        sr1 = math.sqrt(f/a)
1286        sr2 = math.sqrt(f/b)
1287        if sr1 > sr2:
1288            sr1,sr2 = SwapXY(sr1,sr2)
1289            phi -= 90.
1290            if phi < -90.:
1291                phi += 180.
1292        return cent,phi,[sr1,sr2]
1293               
1294    ring = np.array(ring)
1295    x = np.asarray(ring.T[0])
1296    y = np.asarray(ring.T[1])
1297    M = np.array((x**2-y**2,x*y,x,y,np.ones_like(x)))
1298    B = np.array(-(x**2+y**2))
1299    result = nl.lstsq(M.T,B)
1300    return makeParmsEllipse(result[0])
1301   
1302def FitDetector(rings,p0,wave):
1303    from scipy.optimize import leastsq
1304    def ellipseCalc(B,xyd,wave):
1305        x = xyd[0]
1306        y = xyd[1]
1307        dsp = xyd[2]
1308        dist,x0,y0,phi,tilt = B
1309        tth = 2.0*npasind(wave/(2.*dsp))
1310        ttth = nptand(tth)
1311        radius = dist*ttth
1312        stth = npsind(tth)
1313        cosb = npcosd(tilt)
1314        R1 = dist*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
1315        R0 = np.sqrt(R1*radius*cosb)
1316        zdis = R1*ttth*nptand(tilt)
1317        X = x-x0+zdis*npsind(phi)
1318        Y = y-y0-zdis*npcosd(phi)
1319        XR = X*npcosd(phi)-Y*npsind(phi)
1320        YR = X*npsind(phi)+Y*npcosd(phi)
1321        return (XR/R0)**2+(YR/R1)**2-1
1322    result = leastsq(ellipseCalc,p0,args=(rings.T,wave))
1323    return result[0]
1324           
1325def ImageLocalMax(image,w,Xpix,Ypix):
1326    w2 = w*2
1327    size = len(image)
1328    if (w < Xpix < size-w) and (w < Ypix < size-w) and image[Ypix,Xpix]:
1329        Z = image[Ypix-w:Ypix+w,Xpix-w:Xpix+w]
1330        Zmax = np.argmax(Z)
1331        Zmin = np.argmin(Z)
1332        Xpix += Zmax%w2-w
1333        Ypix += Zmax/w2-w
1334        return Xpix,Ypix,np.ravel(Z)[Zmax],np.ravel(Z)[Zmin]
1335    else:
1336        return 0,0,0,0
1337   
1338def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
1339    cent,phi,radii = ellipse
1340    cphi = cosd(phi)
1341    sphi = sind(phi)
1342    ring = []
1343    for a in range(0,360,2):
1344        x = radii[0]*cosd(a)
1345        y = radii[1]*sind(a)
1346        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
1347        Y = (sphi*x+cphi*y+cent[1])*scaley
1348        X,Y,I,J = ImageLocalMax(image,pix,X,Y)     
1349        if I and J and I/J > reject:
1350            X /= scalex                         #convert to mm
1351            Y /= scaley
1352            ring.append([X,Y,dsp])
1353    if len(ring) < 45:             #want more than 1/4 of a circle
1354        return []
1355    return ring
1356   
1357def calcDist(radii,tth):
1358    stth = sind(tth)
1359    ctth = cosd(tth)
1360    ttth = tand(tth)
1361    return math.sqrt(radii[0]**4/(ttth**2*((radii[0]*ctth)**2+(radii[1]*stth)**2)))
1362   
1363def calcZdisCosB(radius,tth,radii):
1364    cosB = sinb = radii[0]**2/(radius*radii[1])
1365    if cosB > 1.:
1366        return 0.,1.
1367    else:
1368        cosb = math.sqrt(1.-sinb**2)
1369        ttth = tand(tth)
1370        zdis = radii[1]*ttth*cosb/sinb
1371        return zdis,cosB
1372   
1373def GetEllipse(dsp,data):
1374    dist = data['distance']
1375    cent = data['center']
1376    tilt = data['tilt']
1377    phi = data['rotation']
1378    radii = [0,0]
1379    tth = 2.0*asind(data['wavelength']/(2.*dsp))
1380    ttth = tand(tth)
1381    stth = sind(tth)
1382    ctth = cosd(tth)
1383    cosb = cosd(tilt)
1384    radius = dist*ttth
1385    radii[1] = dist*stth*ctth*cosb/(cosb**2-stth**2)
1386    if radii[1] > 0:
1387        radii[0] = math.sqrt(radii[1]*radius*cosb)
1388        zdis = radii[1]*ttth*tand(tilt)
1389        elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)]
1390        return elcent,phi,radii
1391    else:
1392        return False
1393       
1394def GetDetectorXY(dsp,azm,data):
1395    from scipy.optimize import fsolve
1396    def func(xy,*args):
1397       azm,phi,R0,R1,A,B = args
1398       cp = cosd(phi)
1399       sp = sind(phi)
1400       x,y = xy
1401       out = []
1402       out.append(y-x*tand(azm))
1403       out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
1404       return out
1405    elcent,phi,radii = GetEllipse(dsp,data)
1406    cent = data['center']
1407    tilt = data['tilt']
1408    phi = data['rotation']
1409    wave = data['wavelength']
1410    dist = data['distance']
1411    tth = 2.0*asind(wave/(2.*dsp))
1412    ttth = tand(tth)
1413    radius = dist*ttth
1414    stth = sind(tth)
1415    cosb = cosd(tilt)
1416    R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
1417    R0 = math.sqrt(R1*radius*cosb)
1418    zdis = R1*ttth*tand(tilt)
1419    A = zdis*sind(phi)
1420    B = -zdis*cosd(phi)
1421    xy0 = [radius*cosd(azm),radius*sind(azm)]
1422    xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent
1423    return xy
1424                   
1425def GetTthDspAzm(x,y,data):
1426    wave = data['wavelength']
1427    dist = data['distance']
1428    cent = data['center']
1429    tilt = data['tilt']
1430    phi = data['rotation']
1431    dx = np.array(x-cent[0],dtype=np.float32)
1432    dy = np.array(y-cent[1],dtype=np.float32)
1433    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
1434    X = np.dot(X,makeMat(phi,2))
1435    Z = np.dot(X,makeMat(tilt,0)).T[2]
1436    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
1437    dsp = wave/(2.*npsind(tth/2.))
1438    azm = npatan2d(dy,dx)
1439    return tth,azm,dsp
1440   
1441def GetTth(x,y,data):
1442    return GetTthDspAzm(x,y,data)[0]
1443   
1444def GetTthAzm(x,y,data):
1445    return GetTthDspAzm(x,y,data)[0:2]
1446   
1447def GetDsp(x,y,data):
1448    return GetTthDspAzm(x,y,data)[2]
1449       
1450def ImageCompress(image,scale):
1451    if scale == 1:
1452        return image
1453    else:
1454        return image[::scale,::scale]
1455       
1456def ImageCalibrate(self,data):
1457    import copy
1458    import ImageCalibrants as calFile
1459    print 'image calibrate'
1460    ring = data['ring']
1461    pixelSize = data['pixelSize']
1462    scalex = 1000./pixelSize[0]
1463    scaley = 1000./pixelSize[1]
1464    cutoff = data['cutoff']
1465    if len(ring) < 5:
1466        print 'not enough inner ring points for ellipse'
1467        return False
1468       
1469    #fit start points on inner ring
1470    data['ellipses'] = []
1471    outE = FitRing(ring)
1472    if outE:
1473        print 'start ellipse:',outE
1474        ellipse = outE
1475    else:
1476        return False
1477       
1478    #setup 180 points on that ring for "good" fit
1479    Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ)
1480    if Ring:
1481        ellipse = FitRing(Ring)
1482        Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ)    #do again
1483        ellipse = FitRing(Ring)
1484    else:
1485        print '1st ring not sufficiently complete to proceed'
1486        return False
1487    print 'inner ring:',ellipse
1488    data['center'] = copy.copy(ellipse[0])           #not right!! (but useful for now)
1489    data['ellipses'].append(ellipse[:]+('r',))
1490    G2plt.PlotImage(self)
1491   
1492    #setup for calibration
1493    data['rings'] = []
1494    data['ellipses'] = []
1495    if not data['calibrant']:
1496        print 'no calibration material selected'
1497        return True
1498       
1499    Bravais,cell = calFile.Calibrants[data['calibrant']]
1500    A = cell2A(cell)
1501    wave = data['wavelength']
1502    cent = data['center']
1503    pixLimit = data['pixLimit']
1504    elcent,phi,radii = ellipse
1505    HKL = GenHBravais(0.5,Bravais,A)
1506    dsp = HKL[0][3]
1507    tth = 2.0*asind(wave/(2.*dsp))
1508    ttth = tand(tth)
1509    data['distance'] = dist = calcDist(radii,tth)
1510    radius = dist*tand(tth)
1511    zdis,cosB = calcZdisCosB(radius,tth,radii)
1512    cent1 = []
1513    cent2 = []
1514    xSum = 0
1515    ySum = 0
1516    zxSum = 0
1517    zySum = 0
1518    phiSum = 0
1519    tiltSum = 0
1520    distSum = 0
1521    Zsum = 0
1522    for i,H in enumerate(HKL):
1523        dsp = H[3]
1524        tth = 2.0*asind(0.5*wave/dsp)
1525        stth = sind(tth)
1526        ctth = cosd(tth)
1527        ttth = tand(tth)
1528        radius = dist*ttth
1529        elcent,phi,radii = ellipse
1530        radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2)
1531        radii[0] = math.sqrt(radii[1]*radius*cosB)
1532        zdis,cosB = calcZdisCosB(radius,tth,radii)
1533        zsinp = zdis*sind(phi)
1534        zcosp = zdis*cosd(phi)
1535        cent = data['center']
1536        elcent = [cent[0]+zsinp,cent[1]-zcosp]
1537        ratio = radii[1]/radii[0]
1538        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
1539        if Ring:
1540            numZ = len(Ring)
1541            data['rings'].append(np.array(Ring))
1542            ellipse = FitRing(Ring)
1543            elcent,phi,radii = ellipse               
1544            if abs(phi) > 45. and phi < 0.:
1545                phi += 180.
1546            dist = calcDist(radii,tth)
1547            distR = 1.-dist/data['distance']
1548            if distR > 0.001:
1549                print 'Wavelength too large?'
1550            elif distR < -0.001:
1551                print 'Wavelength too small?'
1552            else:
1553                if abs((radii[1]/radii[0]-ratio)/ratio) > 0.01:
1554                    print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i)
1555                    return False
1556            zdis,cosB = calcZdisCosB(radius,tth,radii)
1557            Tilt = acosd(cosB)          # 0 <= tilt <= 90
1558            zsinp = zdis*sind(ellipse[1])
1559            zcosp = zdis*cosd(ellipse[1])
1560            cent1.append(np.array([elcent[0]+zsinp,elcent[1]-zcosp]))
1561            cent2.append(np.array([elcent[0]-zsinp,elcent[1]+zcosp]))
1562            if i:
1563                d1 = cent1[-1]-cent1[-2]        #get shift of 2 possible center solutions
1564                d2 = cent2[-1]-cent2[-2]
1565                if np.dot(d2,d2) > np.dot(d1,d1):  #right solution is the larger shift
1566                    data['center'] = cent1[-1]
1567                else:
1568                    data['center'] = cent2[-1]
1569                Zsum += numZ
1570                phiSum += numZ*phi
1571                distSum += numZ*dist
1572                xSum += numZ*data['center'][0]
1573                ySum += numZ*data['center'][1]
1574                tiltSum += numZ*abs(Tilt)
1575            cent = data['center']
1576            print 'for ring # %2i dist %.3f rotate %6.2f tilt %6.2f Xcent %.3f Ycent %.3f Npts %d' \
1577                %(i,dist,phi,Tilt,cent[0],cent[1],numZ)
1578            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
1579            G2plt.PlotImage(self)
1580        else:
1581            break
1582    fullSize = len(self.ImageZ)/scalex
1583    if 2*radii[1] < .9*fullSize:
1584        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
1585    if not Zsum:
1586        print 'Only one ring fitted. Check your wavelength.'
1587        return False
1588    cent = data['center'] = [xSum/Zsum,ySum/Zsum]
1589    wave = data['wavelength']
1590    dist = data['distance'] = distSum/Zsum
1591   
1592    #possible error if no. of rings < 3! Might need trap here
1593    d1 = cent1[-1]-cent1[1]             #compare last ring to 2nd ring
1594    d2 = cent2[-1]-cent2[1]
1595    Zsign = 1
1596    len1 = math.sqrt(np.dot(d1,d1))
1597    len2 = math.sqrt(np.dot(d2,d2))
1598    t1 = d1/len1
1599    t2 = d2/len2
1600    if len2 > len1:
1601        if -135. < atan2d(t2[1],t2[0]) < 45.:
1602            Zsign = -1
1603    else:
1604        if -135. < atan2d(t1[1],t1[0]) < 45.:
1605            Zsign = -1
1606   
1607    tilt = data['tilt'] = Zsign*tiltSum/Zsum
1608    phi = data['rotation'] = phiSum/Zsum
1609    rings = np.concatenate((data['rings']),axis=0)
1610    p0 = [dist,cent[0],cent[1],phi,tilt]
1611    result = FitDetector(rings,p0,wave)
1612    data['distance'] = result[0]
1613    data['center'] = result[1:3]
1614    data['rotation'] = np.mod(result[3],360.0)
1615    data['tilt'] = result[4]
1616    N = len(data['ellipses'])
1617    data['ellipses'] = []           #clear away individual ellipse fits
1618    for H in HKL[:N]:
1619        ellipse = GetEllipse(H[3],data)
1620        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
1621    G2plt.PlotImage(self)       
1622    return True
1623   
1624def ImageIntegrate(self,data):
1625    print 'image integrate'
1626    pixelSize = data['pixelSize']
1627    scalex = pixelSize[0]/1000.
1628    scaley = pixelSize[1]/1000.
1629    LUtth = data['IOtth']
1630    if data['fullIntegrate']:
1631        LRazm = [-180,180]
1632    else:
1633        LRazm = data['LRazimuth']
1634    numAzms = data['outAzimuths']
1635    numChans = data['outChannels']
1636    outGrid = np.zeros(shape=(numAzms,numChans))
1637    outNum = np.zeros(shape=(numAzms,numChans))           
1638    imageN = len(self.ImageZ)
1639    t0 = time.time()
1640    print 'Create ',imageN,' X ',imageN,' 2-theta,azimuth map'
1641    tax,tay = np.mgrid[0:imageN,0:imageN]
1642    tax = np.asfarray(tax)
1643    tay = np.asfarray(tay)
1644    tax *= scalex
1645    tay *= scaley
1646    t1 = time.time()
1647    print "Elapsed time:","%8.3f"%(t1-t0), "s"
1648    print 'Fill map with 2-theta/azimuth values'
1649    self.TA = GetTthAzm(tay,tax,data)           #2-theta & azimuth arrays
1650    self.TA = np.reshape(self.TA,(2,imageN,imageN))
1651    self.TA = np.dstack((self.TA[1],self.TA[0],self.ImageZ))    #azimuth, 2-theta, intensity order
1652    t2 = time.time()
1653    print "Elapsed time:","%8.3f"%(t2-t1), "s"
1654    G2plt.PlotTRImage(self,newPlot=True)
1655    print 'Form 1-D histograms for ',numAzms,' azimuthal angles'
1656    print 'Integration limits:',LUtth,LRazm
1657    tax,tay,taz = np.dsplit(self.TA,3)    #azimuth, 2-theta, intensity
1658    NST = np.histogram2d(tax.flatten(),tay.flatten(),normed=False, \
1659        bins=(numAzms,numChans),range=[LRazm,LUtth])
1660    HST = np.histogram2d(tax.flatten(),tay.flatten(),normed=False, \
1661        bins=(numAzms,numChans),weights=taz.flatten(),range=[LRazm,LUtth])
1662    t3 = time.time()
1663    print "Elapsed time:","%8.3f"%(t3-t2), "s"
1664    self.Integrate = [HST[0]/NST[0],HST[1],HST[2]]
1665    G2plt.PlotIntegration(self,newPlot=True)
1666   
1667       
1668def test():
1669    cell = [5,5,5,90,90,90]
1670    A = cell2A(cell)
1671    assert ( calc_V(A) == 125. )
1672   
1673    print 'test passed'
1674
1675print __name__
1676if __name__ == '__main__':
1677    test()
1678   
1679       
Note: See TracBrowser for help on using the repository browser.