source: trunk/GSASIIindex.py @ 230

Last change on this file since 230 was 230, checked in by vondreele, 11 years ago

add svn keywords & add log intensity plotting

  • Property svn:keywords set to Date Author Revision URL Id
File size: 20.8 KB
Line 
1#GSASII cell indexing program: variation on that of A. Coehlo
2#   includes cell refinement from peak positions (not zero as yet)
3########### SVN repository information ###################
4# $Date: 2011-01-07 19:27:24 +0000 (Fri, 07 Jan 2011) $
5# $Author: vondreele $
6# $Revision: 230 $
7# $URL: trunk/GSASIIindex.py $
8# $Id: GSASIIindex.py 230 2011-01-07 19:27:24Z vondreele $
9########### SVN repository information ###################
10import math
11import wx
12import time
13import numpy as np
14import numpy.linalg as nl
15import GSASIIpath
16import pypowder as pyp              #assumes path has been amended to include correctr bin directory
17import GSASIIlattice as G2lat
18
19# trig functions in degrees
20sind = lambda x: math.sin(x*math.pi/180.)
21asind = lambda x: 180.*math.asin(x)/math.pi
22tand = lambda x: math.tan(x*math.pi/180.)
23atand = lambda x: 180.*math.atan(x)/math.pi
24atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
25cosd = lambda x: math.cos(x*math.pi/180.)
26acosd = lambda x: 180.*math.acos(x)/math.pi
27rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
28#numpy versions
29npsind = lambda x: np.sin(x*np.pi/180.)
30npasind = lambda x: 180.*np.arcsin(x)/math.pi
31npcosd = lambda x: np.cos(x*math.pi/180.)
32nptand = lambda x: np.tan(x*math.pi/180.)
33npatand = lambda x: 180.*np.arctan(x)/np.pi
34npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
35   
36def scaleAbyV(A,V):
37    v = G2lat.calc_V(A)
38    scale = math.exp(math.log(v/V)/3.)**2
39    for i in range(6):
40        A[i] *= scale
41   
42def ranaxis(dmin,dmax):
43    import random as rand
44    return rand.random()*(dmax-dmin)+dmin
45   
46def ran2axis(k,N):
47    import random as rand
48    T = 1.5+0.49*k/N
49#    B = 0.99-0.49*k/N
50#    B = 0.99-0.049*k/N
51    B = 0.99-0.149*k/N
52    R = (T-B)*rand.random()+B
53    return R
54   
55#def ranNaxis(k,N):
56#    import random as rand
57#    T = 1.0+1.0*k/N
58#    B = 1.0-1.0*k/N
59#    R = (T-B)*rand.random()+B
60#    return R
61   
62def ranAbyV(Bravais,dmin,dmax,V):
63    cell = [0,0,0,0,0,0]
64    bad = True
65    while bad:
66        bad = False
67        cell = rancell(Bravais,dmin,dmax)
68        G,g = G2lat.cell2Gmat(cell)
69        A = G2lat.Gmat2A(G)
70        if G2lat.calc_rVsq(A) < 1:
71            scaleAbyV(A,V)
72            cell = G2lat.A2cell(A)
73            for i in range(3):
74                bad |= cell[i] < dmin
75    return A
76   
77def ranAbyR(Bravais,A,k,N,ranFunc):
78    R = ranFunc(k,N)
79    if Bravais in [0,1,2]:          #cubic - not used
80        A[0] = A[1] = A[2] = A[0]*R
81        A[3] = A[4] = A[5] = 0.
82    elif Bravais in [3,4]:          #hexagonal/trigonal
83        A[0] = A[1] = A[3] = A[0]*R
84        A[2] *= R
85        A[4] = A[5] = 0.       
86    elif Bravais in [5,6]:          #tetragonal
87        A[0] = A[1] = A[0]*R
88        A[2] *= R
89        A[3] = A[4] = A[5] = 0.       
90    elif Bravais in [7,8,9,10]:     #orthorhombic
91        A[0] *= R
92        A[1] *= R
93        A[2] *= R
94        A[3] = A[4] = A[5] = 0.       
95    elif Bravais in [11,12]:        #monoclinic
96        A[0] *= R
97        A[1] *= R
98        A[2] *= R
99        A[4] *= R
100        A[3] = A[5] = 0.       
101    else:                           #triclinic
102        A[0] *= R
103        A[1] *= R
104        A[2] *= R
105        A[3] *= R
106        A[4] *= R
107        A[5] *= R
108    return A
109   
110def rancell(Bravais,dmin,dmax):
111    if Bravais in [0,1,2]:          #cubic
112        a = b = c = ranaxis(dmin,dmax)
113        alp = bet = gam = 90
114    elif Bravais in [3,4]:          #hexagonal/trigonal
115        a = b = ranaxis(dmin,dmax)
116        c = ranaxis(dmin,dmax)
117        alp = bet =  90
118        gam = 120
119    elif Bravais in [5,6]:          #tetragonal
120        a = b = ranaxis(dmin,dmax)
121        c = ranaxis(dmin,dmax)
122        alp = bet = gam = 90
123    elif Bravais in [7,8,9,10]:       #orthorhombic - F,I,C,P - a<b<c convention
124        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
125        abc.sort()
126        a = abc[0]
127        b = abc[1]
128        c = abc[2]
129        alp = bet = gam = 90
130    elif Bravais in [11,12]:        #monoclinic - C,P - a<c convention
131        ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
132        ac.sort()
133        a = ac[0]
134        b = ranaxis(dmin,dmax)
135        c = ac[1]
136        alp = gam = 90
137        bet = ranaxis(90.,130.)
138    else:                           #triclinic - a<b<c convention
139        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
140        abc.sort()
141        a = abc[0]
142        b = abc[1]
143        c = abc[2]
144        r = 0.5*b/c
145        alp = ranaxis(acosd(r),acosd(-r))
146        r = 0.5*a/c
147        bet = ranaxis(acosd(r),acosd(-r))
148        r = 0.5*a/b
149        gam = ranaxis(acosd(r),acosd(-r)) 
150    return [a,b,c,alp,bet,gam]
151   
152def calc_M20(peaks,HKL):
153    diff = 0
154    X20 = 0
155    for Nobs20,peak in enumerate(peaks):
156        if peak[3]:
157            Qobs = 1.0/peak[7]**2
158            Qcalc = 1.0/peak[8]**2
159            diff += abs(Qobs-Qcalc)
160        elif peak[2]:
161            X20 += 1
162        if Nobs20 == 19: 
163            d20 = peak[7]
164            break
165    else:
166        d20 = peak[7]
167        Nobs20 = len(peaks)
168    for N20,hkl in enumerate(HKL):
169        if hkl[3] < d20:
170            break               
171    eta = diff/Nobs20
172    Q20 = 1.0/d20**2
173    if diff:
174        M20 = Q20/(2.0*diff)
175    else:
176        M20 = 0
177    M20 /= (1.+X20)
178    return M20,X20
179   
180def sortM20(cells):
181    #cells is M20,X20,Bravais,a,b,c,alp,bet,gam
182    #sort highest M20 1st
183    T = []
184    for i,M in enumerate(cells):
185        T.append((M[0],i))
186    D = dict(zip(T,cells))
187    T.sort()
188    T.reverse()
189    X = []
190    for key in T:
191        X.append(D[key])
192    return X
193               
194def IndexPeaks(peaks,HKL):
195    import bisect
196    hklds = [1000.0]                                    #make bounded list of available d-spacings
197    N = len(HKL)
198    if N == 0: return False
199    for hkl in HKL:
200        hklds.append(hkl[3])
201    hklds.append(0.0)
202    hklds.sort()                                        # ascending sort - upper bound at end
203    hklmax = [0,0,0]
204    for ipk,peak in enumerate(peaks):
205        if peak[2]:
206            i = bisect.bisect_right(hklds,peak[7])          # find peak position in hkl list
207            dm = peak[7]-hklds[i-1]                         # peak to neighbor hkls in list
208            dp = hklds[i]-peak[7]
209            pos = N-i                                       # reverse the order
210            if dp > dm: pos += 1                            # closer to upper than lower
211            hkl = HKL[pos]                                 # put in hkl
212            if hkl[4] >= 0:                                 # peak already assigned - test if this one better
213                opeak = peaks[hkl[4]]
214                dold = abs(opeak[7]-hkl[3])
215                dnew = min(dm,dp)
216                if dold > dnew:                             # new better - zero out old
217                    for j in range(3):
218                        opeak[j+4] = 0
219                    opeak[8] = 0.
220                else:                                       # old better - do nothing
221                    continue               
222            hkl[4] = ipk
223            for j in range(3):
224                peak[j+4] = hkl[j]
225            peak[8] = hkl[3]                                # fill in d-calc
226    for peak in peaks:
227        peak[3] = False
228        if peak[2]:
229            if peak[8] > 0.:
230                for j in range(3):
231                    if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
232                peak[3] = True
233    if hklmax[0]*hklmax[1]*hklmax[2] > 0:
234        return True
235    else:
236        return False
237   
238def FitHKL(ibrav,peaks,A,wtP):
239    def ShiftTest(a,b):
240        if b < -0.1*a: 
241            b = -0.0001*a
242        return b
243    smin = 0.
244    first = True
245    for peak in peaks:
246        if peak[2] and peak[3]:
247            h,k,l = H = peak[4:7]
248            Qo = 1./peak[7]**2
249            Qc = G2lat.calc_rDsq(H,A)
250            try:
251                peak[8] = 1./math.sqrt(Qc)
252            except:
253                print G2lat.A2invcell(A)
254            delt = Qo-Qc
255            smin += delt**2
256            dp = []
257            if ibrav in [0,1,2]:            #m3m
258                dp.append(h*h+k*k+l*l)
259            elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
260                dp.append(h*h+k*k+h*k)
261                dp.append(l*l)
262            elif ibrav in [5,6]:            #4/mmm
263                dp.append(h*h+k*k)
264                dp.append(l*l)
265            elif ibrav in [7,8,9,10]:       #mmm
266                dp.append(h*h)
267                dp.append(k*k)
268                dp.append(l*l)
269            elif ibrav in [11,12]:          #2/m
270                dp.append(h*h)
271                dp.append(k*k)
272                dp.append(l*l)
273                dp.append(h*l)
274            else:                           #1
275#    derivatives for a0*h^2+a1*k^2+a2*l^2+a3*h*k+a4*h*l+a5*k*l
276                dp.append(h*h)
277                dp.append(k*k)
278                dp.append(l*l)
279                dp.append(h*k)
280                dp.append(h*l)
281                dp.append(k*l)
282            if first:
283                first = False
284                M = len(dp)
285                B = np.zeros(shape=(M,M))
286                V = np.zeros(shape=(M))
287            dwt = peak[7]**wtP
288            B,V = pyp.buildmv(delt*dwt,dwt,M,dp,B,V)
289    if nl.det(B) > 0.0:
290        try:
291            b = nl.solve(B,V)
292            B = nl.inv(B)
293            sig = np.diag(B)
294        except SingularMatrix:
295            return False,0
296        if ibrav in [0,1,2]:            #m3m
297            A[0] += ShiftTest(A[0],b[0])
298            A[1] = A[2] = A[0]
299        elif ibrav in [3,4]:            #R3H, P3/m & P6/mmm
300            A[0] += ShiftTest(A[0],b[0])
301            A[1] = A[3] = A[0]
302            A[2] += ShiftTest(A[2],b[1])
303        elif ibrav in [5,6]:            #4/mmm
304            A[0] += ShiftTest(A[0],b[0])
305            A[1] = A[0]
306            A[2] += ShiftTest(A[2],b[1])
307        elif ibrav in [7,8,9,10]:       #mmm
308            for i in range(3):
309                A[i] += ShiftTest(A[i],b[i])
310        elif ibrav in [11,12]:          #2/m
311            for i in range(3):
312                A[i] += ShiftTest(A[i],b[i])
313            A[4] += ShiftTest(A[4],b[3])
314            A[4] = min(1.4*math.sqrt(A[0]*A[2]),A[4])   #min beta star = 45
315        else:                           #1
316            oldV = math.sqrt(1./G2lat.calc_rVsq(A))
317            oldA = A[:]
318            for i in range(6):
319                A[i] += b[i]*0.2
320            A[3] = min(1.1*math.sqrt(max(0,A[1]*A[2])),A[3])
321            A[4] = min(1.1*math.sqrt(max(0,A[0]*A[2])),A[4])
322            A[5] = min(1.1*math.sqrt(max(0,A[0]*A[1])),A[5])
323            ratio = math.sqrt(1./G2lat.calc_rVsq(A))/oldV
324            if 0.9 > ratio or ratio > 1.1:
325                A = oldA
326#    else:
327#        print B
328#        print V
329#        for peak in peaks:
330#            print peak
331    return True,smin
332       
333def FitHKLZ(ibrav,peaks,Z,A):
334    return A,Z
335   
336def rotOrthoA(A):
337    return [A[1],A[2],A[0],0,0,0]
338   
339def swapMonoA(A):
340    return [A[2],A[1],A[0],0,A[4],0]
341   
342def oddPeak(indx,peaks):
343    noOdd = True
344    for peak in peaks:
345        H = peak[4:7]
346        if H[indx] % 2:
347            noOdd = False
348    return noOdd
349   
350def halfCell(ibrav,A,peaks):
351    if ibrav in [0,1,2]:
352        if oddPeak(0,peaks):
353            A[0] *= 2
354            A[1] = A[2] = A[0]
355    elif ibrav in [3,4,5,6]:
356        if oddPeak(0,peaks):
357            A[0] *= 2
358            A[1] = A[0]
359        if oddPeak(2,peaks):
360            A[2] *=2
361    else:
362        if oddPeak(0,peaks):
363            A[0] *=2
364        if oddPeak(1,peaks):
365            A[1] *=2
366        if oddPeak(2,peaks):
367            A[2] *=2
368    return A
369   
370def getDmin(peaks):
371    return peaks[-1][7]
372   
373def getDmax(peaks):
374    return peaks[0][7]
375   
376def refinePeaks(peaks,ibrav,A):
377    dmin = getDmin(peaks)
378    smin = 1.0e10
379    pwr = 6
380    maxTries = 10
381    if ibrav == 13:
382        pwr = 4
383        maxTries = 10
384    OK = False
385    tries = 0
386    HKL = G2lat.GenHBravais(dmin,ibrav,A)
387    while IndexPeaks(peaks,HKL):
388        Pwr = pwr - 2*(tries % 2)
389        HKL = []
390        tries += 1
391        osmin = smin
392        oldA = A
393        OK,smin = FitHKL(ibrav,peaks,A,Pwr)
394        for a in A[:3]:
395            if a < 0:
396                A = oldA
397                OK = False
398                break
399        if OK:
400            HKL = G2lat.GenHBravais(dmin,ibrav,A)
401        if len(HKL) == 0: break                         #absurd cell obtained!
402        rat = (osmin-smin)/smin
403        if abs(rat) < 1.0e-5 or not OK: break
404        if tries > maxTries: break
405    if OK:
406        OK,smin = FitHKL(ibrav,peaks,A,4)
407    M20,X20 = calc_M20(peaks,HKL)
408    return len(HKL),M20,X20
409       
410def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1):
411# dlg & ncMax are used for wx progress bar
412# A != 0 find the best A near input A,
413# A = 0 for random cell, volume normalized to V1;
414# returns number of generated hkls, M20, X20 & A for best found
415    mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10]
416    dmin = getDmin(peaks)-0.05
417    amin = 2.5
418    amax = 5.*getDmax(peaks)
419    Asave = []
420    GoOn = True
421    if A:
422        HKL = G2lat.GenHBravais(dmin,ibrav,A[:])
423        if len(HKL) > mHKL[ibrav]:
424            IndexPeaks(peaks,HKL)
425            Asave.append([calc_M20(peaks,HKL),A[:]])
426    tries = 0
427    while tries < Ntries:
428        if A:
429            Anew = ranAbyR(ibrav,A[:],tries+1,Ntries,ran2axis)
430            if ibrav in [11,12,13]:
431                Anew = ranAbyR(ibrav,A[:],tries/10+1,Ntries,ran2axis)
432        else:
433            Anew = ranAbyV(ibrav,amin,amax,V1)
434        HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
435        if len(HKL) > mHKL[ibrav] and IndexPeaks(peaks,HKL):
436            Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
437            Asave.append([calc_M20(peaks,HKL),Anew[:]])
438            if ibrav == 9:                          #C-centered orthorhombic
439                for i in range(2):
440                    Anew = rotOrthoA(Anew[:])
441                    Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
442                    HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
443                    IndexPeaks(peaks,HKL)
444                    Asave.append([calc_M20(peaks,HKL),Anew[:]])
445            elif ibrav == 11:                      #C-centered monoclinic
446                Anew = swapMonoA(Anew[:])
447                Lhkl,M20,X20 = refinePeaks(peaks,ibrav,Anew)
448                HKL = G2lat.GenHBravais(dmin,ibrav,Anew)
449                IndexPeaks(peaks,HKL)
450                Asave.append([calc_M20(peaks,HKL),Anew[:]])
451        else:
452            break
453        Nc = len(HKL)
454        if Nc >= ncMax:
455            GoOn = False
456        elif dlg:
457            GoOn = dlg.Update(Nc)[0]
458            if not GoOn:
459                break
460        tries += 1   
461    X = sortM20(Asave)
462    if X:
463        Lhkl,M20,X20 = refinePeaks(peaks,ibrav,X[0][1])
464        return GoOn,Lhkl,M20,X20,X[0][1]
465    else:
466        return GoOn,0,0,0,0
467       
468def monoCellReduce(ibrav,A):
469    a,b,c,alp,bet,gam = G2lat.A2cell(A)
470    G,g = G2lat.A2Gmat(A)
471    if ibrav in [11]:
472        u = [0,0,-1]
473        v = [1,0,2]
474        anew = math.sqrt(np.dot(np.dot(v,g),v))
475        if anew < a:
476            cang = np.dot(np.dot(u,g),v)/(anew*c)
477            beta = acosd(-abs(cang))
478            A = G2lat.cell2A([anew,b,c,90,beta,90])
479    else:
480        u = [-1,0,0]
481        v = [1,0,1]
482        cnew = math.sqrt(np.dot(np.dot(v,g),v))
483        if cnew < c:
484            cang = np.dot(np.dot(u,g),v)/(a*cnew)
485            beta = acosd(-abs(cang))
486            A = G2lat.cell2A([a,b,cnew,90,beta,90])
487    return A
488
489def DoIndexPeaks(peaks,inst,controls,bravais):
490   
491    def peakDspace(peaks,A):
492        for peak in peaks:
493            if peak[3]:
494                dsq = calc_rDsq(peak[4:7],A)
495                if dsq > 0:
496                    peak[8] = 1./math.sqrt(dsq)
497        return
498    delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
499    amin = 2.5
500    amax = 5.0*getDmax(peaks)
501    dmin = getDmin(peaks)-delt
502    bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
503        'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C',
504        'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic']
505    tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
506    N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,  50,50, 200]
507    N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,   2,2,   4]
508    Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,   2,2,   4]
509    Nobs = len(peaks)
510    wave = inst[1]
511    if len(inst) > 10:
512        zero = inst[3]
513    else:
514        zero = inst[2]
515    print "%s %8.5f %6.3f" % ('wavelength, zero =',wave,zero)
516    print "%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax)
517    ifzero,maxzero,ncno = controls[:3]
518    ncMax = Nobs*ncno
519    print "%s %d %s %d %s %d" % ('change zero =',ifzero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs)
520    cells = []
521    for ibrav in range(14):
522        begin = time.time()
523        if bravais[ibrav]:
524            print 'cell search for ',bravaisNames[ibrav]
525            print '      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test'
526            V1 = controls[3]
527            bestM20 = 0
528            topM20 = 0
529            cycle = 0
530            while cycle < 5:
531                dlg = wx.ProgressDialog("Generated reflections",tries[cycle]+" cell search for "+bravaisNames[ibrav],ncMax, 
532                    style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE|wx.PD_REMAINING_TIME|wx.PD_CAN_ABORT)
533                screenSize = wx.DisplaySize()
534                Size = dlg.GetSize()
535                dlg.SetPosition(wx.Point(screenSize[0]-Size[0]-300,0))
536                try:
537                    GoOn = True
538                    while GoOn:                                                 #Loop over increment of volume
539                        N2 = 0
540                        while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
541                            if ibrav > 2:
542                                if not N2:
543                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
544                                if A:
545                                    GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0)
546                            else:
547                                GoOn,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1)
548                            if Nc >= ncMax:
549                                GoOn = False
550                                break
551                            elif 3*Nc < Nobs:
552                                N2 = 10
553                                break
554                            else:
555                                if not GoOn:
556                                    break
557                                if M20 > 1.0:
558                                    bestM20 = max(bestM20,M20)
559                                    A = halfCell(ibrav,A[:],peaks)
560                                    if ibrav in [12]:
561                                        A = monoCellReduce(ibrav,A[:])
562                                    HKL = G2lat.GenHBravais(dmin,ibrav,A)
563                                    IndexPeaks(peaks,HKL)
564                                    a,b,c,alp,bet,gam = G2lat.A2cell(A)
565                                    V = G2lat.calc_V(A)
566                                    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)
567                                    if M20 >= 2.0:
568                                        cells.append([M20,X20,ibrav,a,b,c,alp,bet,gam,V,False])
569                            if not GoOn:
570                                break
571                            N2 += 1
572                        if ibrav < 11:
573                            V1 *= 1.1
574                        elif ibrav in range(11,14):
575                            V1 *= 1.05
576                        if not GoOn:
577                            if bestM20 > topM20:
578                                topM20 = bestM20
579                                if cells:
580                                    V1 = cells[0][9]
581                                else:
582                                    V1 = 25
583                                ncMax += Nobs
584                                cycle += 1
585                                print 'Restart search, new Max Nc = ',ncMax
586                            else:
587                                cycle = 10
588                finally:
589                    dlg.Destroy()
590            print '%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
591                ', elapsed time = ',G2lat.sec2HMS(time.time()-begin))
592           
593    if cells:
594        cells = sortM20(cells)
595        cells[0][-1] = True
596        return True,dmin,cells
597    else:
598        return False,0,0
599       
Note: See TracBrowser for help on using the repository browser.