source: trunk/GSASIIcomp.py @ 38

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

convert to python 2.6 from Enthought
small fixes to plotting ellipses

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