source: trunk/GSASIIindex.py @ 5112

Last change on this file since 5112 was 5112, checked in by vondreele, 3 years ago

implement energy dispersive x-ray data as a new type 'PXE'; assume Gaussian peak shape only; no sample broadening considered. Peak resolution is only ~0.5% not good enough to see sample broadening
import, plotting, peak fitting, indexing & cell refinement from peak positions "complete"
data is old (from ~20 yrs ago) so no idea about modern ED data.
Start on use in Pawley/LeBail? refinement - Rietveld excluded; currently not working
Some work on ISODISTORT implementation
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 47.0 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII cell indexing program: variation on that of A. Coehlo
3#   includes cell refinement from peak positions (not zero as yet)
4########### SVN repository information ###################
5# $Date: 2021-12-15 22:56:41 +0000 (Wed, 15 Dec 2021) $
6# $Author: vondreele $
7# $Revision: 5112 $
8# $URL: trunk/GSASIIindex.py $
9# $Id: GSASIIindex.py 5112 2021-12-15 22:56:41Z vondreele $
10########### SVN repository information ###################
11'''
12*GSASIIindex: Cell Indexing Module*
13===================================
14
15Cell indexing program: variation on that of A. Coehlo
16includes cell refinement from peak positions
17'''
18
19from __future__ import division, print_function
20import math
21import time
22import numpy as np
23import GSASIIpath
24GSASIIpath.SetVersionNumber("$Revision: 5112 $")
25import GSASIIlattice as G2lat
26import GSASIIpwd as G2pwd
27import GSASIIspc as G2spc
28import GSASIImath as G2mth
29import scipy.optimize as so
30
31# trig functions in degrees
32sind = lambda x: math.sin(x*math.pi/180.)
33asind = lambda x: 180.*math.asin(x)/math.pi
34tand = lambda x: math.tan(x*math.pi/180.)
35atand = lambda x: 180.*math.atan(x)/math.pi
36atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
37cosd = lambda x: math.cos(x*math.pi/180.)
38acosd = lambda x: 180.*math.acos(x)/math.pi
39rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
40#numpy versions
41npsind = lambda x: np.sin(x*np.pi/180.)
42npasind = lambda x: 180.*np.arcsin(x)/math.pi
43npcosd = lambda x: np.cos(x*math.pi/180.)
44nptand = lambda x: np.tan(x*math.pi/180.)
45npatand = lambda x: 180.*np.arctan(x)/np.pi
46npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
47try:  # fails on doc build
48    rpd = np.pi/180.
49except TypeError:
50    pass
51   
52def scaleAbyV(A,V):
53    'needs a doc string'
54    v = G2lat.calc_V(A)
55    scale = math.exp(math.log(v/V)/3.)**2
56    for i in range(6):
57        A[i] *= scale
58   
59def ranaxis(dmin,dmax):
60    'needs a doc string'
61    import random as rand
62    return rand.random()*(dmax-dmin)+dmin
63   
64def ran2axis(k,N):
65    'needs a doc string'
66    import random as rand
67    T = 1.5+0.49*k/N
68#    B = 0.99-0.49*k/N
69#    B = 0.99-0.049*k/N
70    B = 0.99-0.149*k/N
71    R = (T-B)*rand.random()+B
72    return R
73   
74#def ranNaxis(k,N):
75#    import random as rand
76#    T = 1.0+1.0*k/N
77#    B = 1.0-1.0*k/N
78#    R = (T-B)*rand.random()+B
79#    return R
80   
81def ranAbyV(Bravais,dmin,dmax,V):
82    'needs a doc string'
83    cell = [0,0,0,0,0,0]
84    bad = True
85    while bad:
86        bad = False
87        cell = rancell(Bravais,dmin,dmax)
88        G,g = G2lat.cell2Gmat(cell)
89        A = G2lat.Gmat2A(G)
90        if G2lat.calc_rVsq(A) < 1:
91            scaleAbyV(A,V)
92            cell = G2lat.A2cell(A)
93            for i in range(3):
94                bad |= cell[i] < dmin
95    return A
96   
97def ranAbyR(Bravais,A,k,N,ranFunc):
98    'needs a doc string'
99    R = ranFunc(k,N)
100    if Bravais in [0,1,2]:          #cubic - not used
101        A[0] = A[1] = A[2] = A[0]*R
102        A[3] = A[4] = A[5] = 0.
103    elif Bravais in [3,4]:          #hexagonal/trigonal
104        A[0] = A[1] = A[3] = A[0]*R
105        A[2] *= R
106        A[4] = A[5] = 0.       
107    elif Bravais in [5,6]:          #tetragonal
108        A[0] = A[1] = A[0]*R
109        A[2] *= R
110        A[3] = A[4] = A[5] = 0.       
111    elif Bravais in [7,8,9,10,11,12]:     #orthorhombic
112        A[0] *= R
113        A[1] *= R
114        A[2] *= R
115        A[3] = A[4] = A[5] = 0.       
116    elif Bravais in [13,14,15,16]:        #monoclinic
117        A[0] *= R
118        A[1] *= R
119        A[2] *= R
120        A[4] *= R
121        A[3] = A[5] = 0.       
122    else:                           #triclinic
123        A[0] *= R
124        A[1] *= R
125        A[2] *= R
126        A[3] *= R
127        A[4] *= R
128        A[5] *= R
129    return A
130   
131def rancell(Bravais,dmin,dmax):
132    'needs a doc string'
133    if Bravais in [0,1,2]:          #cubic
134        a = b = c = ranaxis(dmin,dmax)
135        alp = bet = gam = 90
136    elif Bravais in [3,4]:          #hexagonal/trigonal
137        a = b = ranaxis(dmin,dmax)
138        c = ranaxis(dmin,dmax)
139        alp = bet =  90
140        gam = 120
141    elif Bravais in [5,6]:          #tetragonal
142        a = b = ranaxis(dmin,dmax)
143        c = ranaxis(dmin,dmax)
144        alp = bet = gam = 90
145    elif Bravais in [7,8,9,10,11,12]:       #orthorhombic - F,I,P - a<b<c convention
146        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
147        if Bravais in [7,8,12]:
148            abc.sort()
149        a = abc[0]
150        b = abc[1]
151        c = abc[2]
152        alp = bet = gam = 90
153    elif Bravais in [13,14,15,16]:        #monoclinic - C,P - a<c convention
154        ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
155        if Bravais in [13,14]:
156            ac.sort()
157        a = ac[0]
158        b = ranaxis(dmin,dmax)
159        c = ac[1]
160        alp = gam = 90
161        bet = ranaxis(90.,140.)
162    else:                           #triclinic - a<b<c convention
163        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
164        abc.sort()
165        a = abc[0]
166        b = abc[1]
167        c = abc[2]
168        r = 0.5*b/c
169        alp = ranaxis(acosd(r),acosd(-r))
170        r = 0.5*a/c
171        bet = ranaxis(acosd(r),acosd(-r))
172        r = 0.5*a/b
173        gam = ranaxis(acosd(r),acosd(-r)) 
174    return [a,b,c,alp,bet,gam]
175   
176def calc_M20(peaks,HKL,ifX20=True):
177    'needs a doc string'
178    diff = 0
179    X20 = 0
180    for Nobs20,peak in enumerate(peaks):
181        if peak[3]:
182            Qobs = 1.0/peak[7]**2
183            Qcalc = 1.0/peak[8]**2
184            diff += abs(Qobs-Qcalc)
185        elif peak[2]:
186            X20 += 1
187        if Nobs20 == 19: 
188            d20 = peak[7]
189            break
190    else:
191        d20 = peak[7]
192        Nobs20 = len(peaks)
193    for N20,hkl in enumerate(HKL):
194        if hkl[3] < d20:
195            break               
196    Q20 = 1.0/d20**2
197    if diff:
198        M20 = Q20/(2.0*diff)
199    else:
200        M20 = 0
201    if ifX20:
202        M20 /= (1.+X20)
203    return M20,X20
204   
205def calc_M20SS(peaks,HKL):
206    'needs a doc string'
207    diff = 0
208    X20 = 0
209    for Nobs20,peak in enumerate(peaks):
210        if peak[3]:
211            Qobs = 1.0/peak[8]**2
212            Qcalc = 1.0/peak[9]**2
213            diff += abs(Qobs-Qcalc)
214        elif peak[2]:
215            X20 += 1
216        if Nobs20 == 19: 
217            d20 = peak[8]
218            break
219    else:
220        d20 = peak[8]
221        Nobs20 = len(peaks)
222    for N20,hkl in enumerate(HKL):
223        if hkl[4] < d20:
224            break               
225    Q20 = 1.0/d20**2
226    if diff:
227        M20 = Q20/(2.0*diff)
228    else:
229        M20 = 0
230    M20 /= (1.+X20)
231    return M20,X20
232   
233def sortM20(cells):
234    'needs a doc string'
235    #cells is M20,X20,Bravais,a,b,c,alp,bet,gam
236    #sort highest M20 1st
237    T = []
238    for i,M in enumerate(cells):
239        T.append((M[0],i))
240    D = dict(zip(T,cells))
241    T.sort()
242    T.reverse()
243    X = []
244    for key in T:
245        X.append(D[key])
246    return X
247               
248def sortCells(cells,col):
249    #cells is M20,X20,Bravais,a,b,c,alp,bet,gam,volume
250    #sort smallest a,b,c,alpha,beta,gamma or volume 1st
251    T = []
252    for i,M in enumerate(cells):
253        T.append((M[col],i))
254    D = dict(zip(T,cells))
255    T.sort()
256    X = []
257    for key in T:
258        X.append(D[key])
259    return X
260   
261def findMV(peaks,controls,ssopt,Inst,dlg):
262       
263    def Val2Vec(vec,Vref,values):
264        Vec = []
265        i = 0
266        for j,r in enumerate(Vref):
267            if r:
268                if values.size > 1:
269                    Vec.append(max(-1.,min(1.0,values[i])))
270                else:
271                    Vec.append(max(0.0,min(1.0,values)))                   
272                i += 1
273            else:
274                Vec.append(vec[j])
275        return np.array(Vec,dtype=float) 
276     
277    def ZSSfunc(values,peaks,dmin,Inst,SGData,SSGData,vec,Vref,maxH,A,wave,Z,dlg=None):
278        Vec = Val2Vec(vec,Vref,values)
279        HKL =  G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A)
280        Peaks = np.array(IndexSSPeaks(peaks,HKL)[1]).T
281        Qo = 1./Peaks[-2]**2
282        Qc = G2lat.calc_rDsqZSS(Peaks[4:8],A,Vec,Z,Peaks[0],wave)
283        chi = np.sum((Qo-Qc)**2)
284        if dlg:
285            dlg.Pulse()   
286        return chi
287       
288    def TSSfunc(values,peaks,dmin,Inst,SGData,SSGData,vec,Vref,maxH,A,difC,Z,dlg=None):
289        Vec = Val2Vec(vec,Vref,values)
290        HKL =  G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A)
291        Peaks = np.array(IndexSSPeaks(peaks,HKL)[1]).T
292        Qo = 1./Peaks[-2]**2
293        Qc = G2lat.calc_rDsqTSS(Peaks[4:8],A,Vec,Z,Peaks[0],difC)
294        chi = np.sum((Qo-Qc)**2)
295        if dlg:
296            dlg.Pulse()   
297        return chi
298
299    if 'T' in Inst['Type'][0]:
300        difC = Inst['difC'][1]
301    else:
302        wave = G2mth.getWave(Inst)
303    SGData = G2spc.SpcGroup(controls[13])[1]
304    SSGData = G2spc.SSpcGroup(SGData,ssopt['ssSymb'])[1]
305    A = G2lat.cell2A(controls[6:12])
306    Z = controls[1]
307    Vref = [True if x in ssopt['ssSymb'] else False for x in ['a','b','g']]
308    values = []
309    ranges = []
310    dT = 0.01       #seems to be a good choice
311    for v,r in zip(ssopt['ModVec'],Vref):
312        if r:
313            ranges += [slice(dT,1.-dT,dT),] #NB: unique part for (00g) & (a0g); (abg)?
314            values += [v,]
315    dmin = getDmin(peaks)-0.005
316    if 'T' in Inst['Type'][0]:   
317        result = so.brute(TSSfunc,ranges,finish=so.fmin_cg,full_output=True,
318            args=(peaks,dmin,Inst,SGData,SSGData,ssopt['ModVec'],Vref,1,A,difC,Z,dlg))
319    else:
320        result = so.brute(ZSSfunc,ranges,finish=so.fmin_cg,full_output=True,
321            args=(peaks,dmin,Inst,SGData,SSGData,ssopt['ModVec'],Vref,1,A,wave,Z,dlg))
322    return Val2Vec(ssopt['ModVec'],Vref,result[0]),result
323               
324def IndexPeaks(peaks,HKL):
325    'needs a doc string'
326    import bisect
327    N = len(HKL)
328    if N == 0: return False,peaks
329    hklds = list(np.array(HKL).T[3])+[1000.0,0.0,]
330    hklds.sort()                                        # ascending sort - upper bound at end
331    hklmax = [0,0,0]
332    for ipk,peak in enumerate(peaks):
333        peak[4:7] = [0,0,0]                           #clear old indexing
334        peak[8] = 0.
335        if peak[2]:
336            i = bisect.bisect_right(hklds,peak[7])          # find peak position in hkl list
337            dm = peak[-2]-hklds[i-1]                         # peak to neighbor hkls in list
338            dp = hklds[i]-peak[-2]
339            pos = N-i                                       # reverse the order
340            if dp > dm: pos += 1                            # closer to upper than lower
341            if pos >= N:
342                break
343            hkl = HKL[pos]                                 # put in hkl
344            if hkl[-1] >= 0:                                 # peak already assigned - test if this one better
345                opeak = peaks[int(hkl[-1])]                 #hkl[-1] needs to be int here
346                dold = abs(opeak[-2]-hkl[3])
347                dnew = min(dm,dp)
348                if dold > dnew:                             # new better - zero out old
349                    opeak[4:7] = [0,0,0]
350                    opeak[8] = 0.
351                else:                                       # old better - do nothing
352                    continue               
353            hkl[-1] = ipk
354            peak[4:7] = hkl[:3]
355            peak[8] = hkl[3]                                # fill in d-calc
356    for peak in peaks:
357        peak[3] = False
358        if peak[2]:
359            if peak[-1] > 0.:
360                for j in range(3):
361                    if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
362                peak[3] = True
363    if hklmax[0]*hklmax[1]*hklmax[2] > 0:
364        return True,peaks
365    else:
366        return False,peaks  #nothing indexed!
367       
368def IndexSSPeaks(peaks,HKL):
369    'needs a doc string'
370    import bisect
371    N = len(HKL)
372    Peaks = np.copy(peaks)
373    if N == 0: return False,Peaks
374    if len(peaks[0]) == 9:      #add m column if missing
375        Peaks = np.insert(Peaks,7,np.zeros_like(Peaks.T[0]),axis=1)
376#        for peak in Peaks:
377#            peak.insert(7,0)
378    hklds = list(np.array(HKL).T[4])+[1000.0,0.0,]
379    hklds.sort()                                        # ascending sort - upper bound at end
380    hklmax = [0,0,0,0]
381    for ipk,peak in enumerate(Peaks):
382        peak[4:8] = [0,0,0,0]                           #clear old indexing
383        peak[9] = 0.
384        if peak[2]: #Use
385            i = bisect.bisect_right(hklds,peak[8])          # find peak position in hkl list
386            dm = peak[8]-hklds[i-1]                         # peak to neighbor hkls in list
387            dp = hklds[i]-peak[8]
388            pos = N-i                                       # reverse the order
389            if dp > dm: pos += 1                            # closer to upper than lower
390            if pos >= N:
391#                print ('%.4f %d'%(pos,N))
392                break
393            hkl = HKL[pos]                                 # put in hkl
394            if hkl[-1] >= 0:                                 # peak already assigned - test if this one better
395                opeak = Peaks[int(hkl[-1])]
396                dold = abs(opeak[-2]-hkl[4])
397                dnew = min(dm,dp)
398                if dold > dnew:                             # new better - zero out old
399                    opeak[4:8] = [0,0,0,0]
400                    opeak[9] = 0.
401                else:                                       # old better - do nothing
402                    continue
403            hkl[-1] = ipk
404            peak[4:8] = hkl[:4]
405            peak[9] = hkl[4]                                # fill in d-calc
406    for peak in Peaks:
407        peak[3] = False
408        if peak[2]:
409            if peak[-1] > 0.:
410                for j in range(4):
411                    if abs(peak[j+4]) > hklmax[j]: hklmax[j] = abs(peak[j+4])
412                peak[3] = True
413    if hklmax[0]*hklmax[1]*hklmax[2]*hklmax[3] > 0:
414        return True,Peaks
415    else:
416        return False,Peaks  #nothing indexed!
417       
418def Values2A(ibrav,values):
419    'needs a doc string'
420    if ibrav in [0,1,2]:
421        return [values[0],values[0],values[0],0,0,0]
422    elif ibrav in [3,4]:
423        return [values[0],values[0],values[1],values[0],0,0]
424    elif ibrav in [5,6]:
425        return [values[0],values[0],values[1],0,0,0]
426    elif ibrav in [7,8,9,10,11,12]:
427        return [values[0],values[1],values[2],0,0,0]
428    elif ibrav in [13,14,15,16]:
429        return [values[0],values[1],values[2],0,values[3],0]
430    else:
431        return list(values[:6])
432       
433def A2values(ibrav,A):
434    'needs a doc string'
435    if ibrav in [0,1,2]:
436        return [A[0],]
437    elif ibrav in [3,4,5,6]:
438        return [A[0],A[2]]
439    elif ibrav in [7,8,9,10,11,12]:
440        return [A[0],A[1],A[2]]
441    elif ibrav in [13,14,15,16]:
442        return [A[0],A[1],A[2],A[4]]
443    else:
444        return A
445       
446def Values2Vec(ibrav,vec,Vref,val):
447    if ibrav in [3,4,5,6]:
448        Nskip = 2
449    elif ibrav in [7,8,9,10,11,12]:
450        Nskip = 3
451    elif ibrav in [13,14,15,16]:
452        Nskip = 4
453    else:
454        Nskip = 6
455    Vec = []
456    i = 0
457    for j,r in enumerate(Vref):
458        if r:
459            Vec.append(val[i+Nskip])
460            i += 1
461        else:
462            Vec.append(vec[j])
463    return np.array(Vec) 
464
465def FitHKL(ibrav,peaks,A,Pwr):
466    'needs a doc string'
467               
468    def errFit(values,ibrav,d,H,Pwr):
469        A = Values2A(ibrav,values)
470        Qo = 1./d**2
471        Qc = G2lat.calc_rDsq(H,A)
472        return (Qo-Qc)*d**Pwr
473       
474    def dervFit(values,ibrav,d,H,Pwr):
475        if ibrav in [0,1,2]:
476            derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],]
477        elif ibrav in [3,4,]:
478            derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]]
479        elif ibrav in [5,6]:
480            derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]]
481        elif ibrav in [7,8,9,10,11,12]:
482            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
483        elif ibrav in [13,14,15,16]:
484            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]]
485        else:
486            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]]
487        derv = -np.array(derv)
488        return (derv*d**Pwr).T
489   
490    Peaks = np.array(peaks).T
491    values = A2values(ibrav,A)
492    result = so.leastsq(errFit,values,Dfun=dervFit,full_output=True,ftol=0.000001,
493        args=(ibrav,Peaks[7],Peaks[4:7],Pwr))
494    A = Values2A(ibrav,result[0])
495    return True,np.sum(errFit(result[0],ibrav,Peaks[7],Peaks[4:7],Pwr)**2),A,result
496           
497def errFitZ(values,ibrav,d,H,tth,wave,Z,Zref):
498    Zero = Z
499    if Zref:   
500        Zero = values[-1]
501    A = Values2A(ibrav,values)
502    Qo = 1./d**2
503    Qc = G2lat.calc_rDsqZ(H,A,Zero,tth,wave)
504    return (Qo-Qc)
505   
506def dervFitZ(values,ibrav,d,H,tth,wave,Z,Zref):
507    if ibrav in [0,1,2]:
508        derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],]
509    elif ibrav in [3,4,]:
510        derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]]
511    elif ibrav in [5,6]:
512        derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]]
513    elif ibrav in [7,8,9,10,11,12]:
514        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
515    elif ibrav in [13,14,15,16]:
516        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]]
517    else:
518        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]]
519    if Zref:
520        derv.append(npsind(tth)*2.0*rpd/wave**2)
521    derv = -np.array(derv)
522    return derv.T
523   
524def FitHKLZ(wave,ibrav,peaks,A,Z,Zref):
525    'needs a doc string'
526   
527    Peaks = np.array(peaks).T   
528    values = A2values(ibrav,A)
529    if Zref:
530        values.append(Z)
531#TODO: try Hessian refinement here - might improve things
532#    result = G2mth.HessianLSQ(errFitZ,values,Hess=peakHess,        #would need "peakHess" routine; returns vec & Hess
533#        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Z,Zref]),ftol=.01,maxcyc=10)
534    result = so.leastsq(errFitZ,values,Dfun=dervFitZ,full_output=True,ftol=0.0001,
535        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Z,Zref))
536    A = Values2A(ibrav,result[0][:6])
537    if Zref:
538        Z = result[0][-1]
539    chisq = np.sum(errFitZ(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Z,Zref)**2)
540    return True,chisq,A,Z,result
541   
542def errFitZSS(values,ibrav,d,H,tth,wave,vec,Vref,Z,Zref):
543    Zero = Z
544    if Zref:   
545        Zero = values[-1]
546    A = Values2A(ibrav,values)
547    Vec = Values2Vec(ibrav,vec,Vref,values)
548    Qo = 1./d**2
549    Qc = G2lat.calc_rDsqZSS(H,A,Vec,Zero,tth,wave)
550    return (Qo-Qc)
551   
552def dervFitZSS(values,ibrav,d,H,tth,wave,vec,Vref,Z,Zref):
553    A = Values2A(ibrav,values)
554    Vec = Values2Vec(ibrav,vec,Vref,values)
555    HM = H[:3]+(H[3][:,np.newaxis]*Vec).T
556    if ibrav in [3,4,]:
557        derv = [HM[0]*HM[0]+HM[1]*HM[1]+HM[0]*HM[1],HM[2]*HM[2]]
558    elif ibrav in [5,6]:
559        derv = [HM[0]*HM[0]+HM[1]*HM[1],HM[2]*HM[2]]
560    elif ibrav in [7,8,9,10,11,12]:
561        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2]]
562    elif ibrav in [13,14,15,16]:
563        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[2]]
564    else:
565        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[1],HM[0]*HM[2],HM[1]*HM[2]]
566    if Vref[0]:
567        derv.append(2.*A[0]*HM[0]*H[3]+A[3]*HM[1]*H[3]+A[4]*HM[2]*H[3])
568    if Vref[1]:
569        derv.append(2.*A[1]*HM[1]*H[3]+A[3]*HM[0]*H[3]+A[5]*HM[2]*H[3])
570    if Vref[2]:
571        derv.append(2.*A[2]*HM[2]*H[3]+A[4]*HM[1]*H[3]+A[5]*HM[0]*H[3])   
572    if Zref:
573        derv.append(npsind(tth)*2.0*rpd/wave**2)
574    derv = -np.array(derv)
575    return derv.T
576   
577def FitHKLZSS(wave,ibrav,peaks,A,V,Vref,Z,Zref):
578    'needs a doc string'
579   
580    Peaks = np.array(peaks).T   
581    values = A2values(ibrav,A)
582    for v,r in zip(V,Vref):
583        if r:
584            values.append(v)
585    if Zref:
586        values.append(Z)
587    result = so.leastsq(errFitZSS,values,Dfun=dervFitZSS,full_output=True,ftol=1.e-6,
588        args=(ibrav,Peaks[8],Peaks[4:8],Peaks[0],wave,V,Vref,Z,Zref))
589    A = Values2A(ibrav,result[0])
590    Vec = Values2Vec(ibrav,V,Vref,result[0])
591    if Zref:
592        Z = result[0][-1]
593    chisq = np.sum(errFitZSS(result[0],ibrav,Peaks[8],Peaks[4:8],Peaks[0],wave,Vec,Vref,Z,Zref)**2) 
594    return True,chisq,A,Vec,Z,result
595   
596def errFitT(values,ibrav,d,H,tof,difC,Z,Zref):
597    Zero = Z
598    if Zref:   
599        Zero = values[-1]
600    A = Values2A(ibrav,values)
601    Qo = 1./d**2
602    Qc = G2lat.calc_rDsqT(H,A,Zero,tof,difC)
603    return (Qo-Qc)
604   
605def dervFitT(values,ibrav,d,H,tof,difC,Z,Zref):
606    if ibrav in [0,1,2]:
607        derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],]
608    elif ibrav in [3,4,]:
609        derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]]
610    elif ibrav in [5,6]:
611        derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]]
612    elif ibrav in [7,8,9,10,11,12]:
613        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
614    elif ibrav in [13,14,15,16]:
615        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]]
616    else:
617        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]]
618    if Zref:
619        derv.append(np.ones_like(d)/difC)
620    derv = np.array(derv)
621    return derv.T
622   
623def FitHKLT(difC,ibrav,peaks,A,Z,Zref):
624    'needs a doc string'
625   
626    Peaks = np.array(peaks).T   
627    values = A2values(ibrav,A)
628    if Zref:
629        values.append(Z)
630    result = so.leastsq(errFitT,values,Dfun=dervFitT,full_output=True,ftol=0.0001,
631        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Z,Zref))
632    A = Values2A(ibrav,result[0])
633    if Zref:
634        Z = result[0][-1]
635    chisq = np.sum(errFitT(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Z,Zref)**2)
636    return True,chisq,A,Z,result
637   
638def FitHKLE(tth,ibrav,peaks,A):
639    'needs a doc string'
640   
641    def errFitE(values,ibrav,d,H,keV,tth):
642        A = Values2A(ibrav,values)
643        Qo = 1./d**2
644        Qc = G2lat.calc_rDsq(H,A)
645        return (Qo-Qc)
646       
647    def dervFitE(values,ibrav,d,H,keV,tth):
648        if ibrav in [0,1,2]:
649            derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],]
650        elif ibrav in [3,4,]:
651            derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]]
652        elif ibrav in [5,6]:
653            derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]]
654        elif ibrav in [7,8,9,10,11,12]:
655            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
656        elif ibrav in [13,14,15,16]:
657            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]]
658        else:
659            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[1],H[0]*H[2],H[1]*H[2]]
660        derv = np.array(derv)
661        return derv.T
662     
663    Peaks = np.array(peaks).T   
664    values = A2values(ibrav,A)
665    result = so.leastsq(errFitE,values,Dfun=dervFitE,full_output=True,ftol=0.0001,
666        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],tth))
667    A = Values2A(ibrav,result[0])
668    chisq = np.sum(errFitE(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],tth)**2)
669    return True,chisq,A,result
670   
671def errFitTSS(values,ibrav,d,H,tof,difC,vec,Vref,Z,Zref):
672    Zero = Z
673    if Zref:   
674        Zero = values[-1]
675    A = Values2A(ibrav,values)
676    Vec = Values2Vec(ibrav,vec,Vref,values)
677    Qo = 1./d**2
678    Qc = G2lat.calc_rDsqTSS(H,A,Vec,Zero,tof,difC)
679    return (Qo-Qc)
680   
681def dervFitTSS(values,ibrav,d,H,tof,difC,vec,Vref,Z,Zref):
682    A = Values2A(ibrav,values)
683    Vec = Values2Vec(ibrav,vec,Vref,values)
684    HM = H[:3]+(H[3][:,np.newaxis]*Vec).T
685    if ibrav in [3,4,]:
686        derv = [HM[0]*HM[0]+HM[1]*HM[1]+HM[0]*HM[1],HM[2]*HM[2]]
687    elif ibrav in [5,6]:
688        derv = [HM[0]*HM[0]+HM[1]*HM[1],HM[2]*HM[2]]
689    elif ibrav in [7,8,9,10,11,12]:
690        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2]]
691    elif ibrav in [13,14,15,16]:
692        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[2]]
693    else:
694        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[1],HM[0]*HM[2],HM[1]*HM[2]]
695    if Vref[0]:
696        derv.append(2.*A[0]*HM[0]*H[3]+A[3]*HM[1]*H[3]+A[4]*HM[2]*H[3])
697    if Vref[1]:
698        derv.append(2.*A[1]*HM[1]*H[3]+A[3]*HM[0]*H[3]+A[5]*HM[2]*H[3])
699    if Vref[2]:
700        derv.append(2.*A[2]*HM[2]*H[3]+A[4]*HM[1]*H[3]+A[5]*HM[0]*H[3])   
701    if Zref:
702        derv.append(np.ones_like(d)/difC)
703    derv = np.array(derv)
704    return derv.T
705   
706def FitHKLTSS(difC,ibrav,peaks,A,V,Vref,Z,Zref):
707    'needs a doc string'
708   
709    Peaks = np.array(peaks).T   
710    values = A2values(ibrav,A)
711    for v,r in zip(V,Vref):
712        if r:
713            values.append(v)
714    if Zref:
715        values.append(Z)
716    print(values)
717    result = so.leastsq(errFitTSS,values,Dfun=dervFitTSS,full_output=True,ftol=0.0001,
718        args=(ibrav,Peaks[8],Peaks[4:8],Peaks[0],difC,V,Vref,Z,Zref))
719    print(result)
720    A = Values2A(ibrav,result[0])
721    Vec = Values2Vec(ibrav,V,Vref,result[0])
722    if Zref:
723        Z = result[0][-1]
724    chisq = np.sum(errFitTSS(result[0],ibrav,Peaks[8],Peaks[4:8],Peaks[0],difC,V,Vref,Z,Zref)**2)
725    return True,chisq,A,Vec,Z,result
726               
727def rotOrthoA(A):
728    'needs a doc string'
729    return [A[1],A[2],A[0],0,0,0]
730   
731def swapMonoA(A):
732    'needs a doc string'
733    return [A[2],A[1],A[0],0,A[4],0]
734   
735def oddPeak(indx,peaks):
736    'needs a doc string'
737    noOdd = True
738    for peak in peaks:
739        H = peak[4:7]
740        if H[indx] % 2:
741            noOdd = False
742    return noOdd
743   
744def halfCell(ibrav,A,peaks):
745    'needs a doc string'
746    if ibrav in [0,1,2]:
747        if oddPeak(0,peaks):
748            A[0] *= 2
749            A[1] = A[2] = A[0]
750    elif ibrav in [3,4,5,6]:
751        if oddPeak(0,peaks):
752            A[0] *= 2
753            A[1] = A[0]
754        if oddPeak(2,peaks):
755            A[2] *=2
756    else:
757        if oddPeak(0,peaks):
758            A[0] *=2
759        if oddPeak(1,peaks):
760            A[1] *=2
761        if oddPeak(2,peaks):
762            A[2] *=2
763    return A
764   
765def getDmin(peaks):
766    'needs a doc string'
767    return peaks[-1][-2]
768   
769def getDmax(peaks):
770    'needs a doc string'
771    return peaks[0][-2]
772   
773def refinePeaksZ(peaks,wave,ibrav,A,Zero,ZeroRef):
774    'needs a doc string'
775    dmin = getDmin(peaks)
776    OK,smin,Aref,Z,result = FitHKLZ(wave,ibrav,peaks,A,Zero,ZeroRef)
777    Peaks = np.array(peaks).T
778    H = Peaks[4:7]
779    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqZ(H,Aref,Z,Peaks[0],wave))
780    peaks = Peaks.T   
781    HKL = G2lat.GenHBravais(dmin,ibrav,A)
782    M20,X20 = calc_M20(peaks,HKL)
783    return len(HKL),M20,X20,Aref,Z
784   
785def refinePeaksT(peaks,difC,ibrav,A,Zero,ZeroRef):
786    'needs a doc string'
787    dmin = getDmin(peaks)
788    OK,smin,Aref,Z,result = FitHKLT(difC,ibrav,peaks,A,Zero,ZeroRef)
789    Peaks = np.array(peaks).T
790    H = Peaks[4:7]
791    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqT(H,Aref,Z,Peaks[0],difC))
792    peaks = Peaks.T   
793    HKL = G2lat.GenHBravais(dmin,ibrav,A)
794    M20,X20 = calc_M20(peaks,HKL)
795    return len(HKL),M20,X20,Aref,Z
796   
797def refinePeaksE(peaks,tth,ibrav,A):
798    'needs a doc string'
799    dmin = getDmin(peaks)
800    OK,smin,Aref,result = FitHKLE(tth,ibrav,peaks,A)
801    Peaks = np.array(peaks).T
802    H = Peaks[4:7]
803    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,Aref))
804    peaks = Peaks.T   
805    HKL = G2lat.GenHBravais(dmin,ibrav,A)
806    M20,X20 = calc_M20(peaks,HKL)
807    return len(HKL),M20,X20,Aref
808   
809def refinePeaksZSS(peaks,wave,Inst,SGData,SSGData,maxH,ibrav,A,vec,vecRef,Zero,ZeroRef):
810    'needs a doc string'
811    dmin = getDmin(peaks)
812    OK,smin,Aref,Vref,Z,result = FitHKLZSS(wave,ibrav,peaks,A,vec,vecRef,Zero,ZeroRef)
813    Peaks = np.array(peaks).T
814    H = Peaks[4:8]
815    Peaks[9] = 1./np.sqrt(G2lat.calc_rDsqZSS(H,Aref,Vref,Z,Peaks[0],wave))  #H,A,vec,Z,tth,lam
816    peaks = Peaks.T   
817    HKL =  G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vref,maxH,Aref)
818    M20,X20 = calc_M20SS(peaks,HKL)
819    return len(HKL),M20,X20,Aref,Vref,Z
820   
821def refinePeaksTSS(peaks,difC,Inst,SGData,SSGData,maxH,ibrav,A,vec,vecRef,Zero,ZeroRef):
822    'needs a doc string'
823    dmin = getDmin(peaks)
824    OK,smin,Aref,Vref,Z,result = FitHKLTSS(difC,ibrav,peaks,A,vec,vecRef,Zero,ZeroRef)
825    Peaks = np.array(peaks).T
826    H = Peaks[4:8]
827    Peaks[9] = 1./np.sqrt(G2lat.calc_rDsqTSS(H,Aref,Vref,Z,Peaks[0],difC))
828    peaks = Peaks.T   
829    HKL =  G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vref,maxH,Aref)
830    HKL = G2lat.GenHBravais(dmin,ibrav,A)
831    M20,X20 = calc_M20SS(peaks,HKL)
832    return len(HKL),M20,X20,Aref,Vref,Z
833   
834def refinePeaks(peaks,ibrav,A,ifX20=True,cctbx_args=None):
835    'needs a doc string'
836    dmin = getDmin(peaks)
837    smin = 1.0e10
838    pwr = 8
839    maxTries = 10
840    OK = False
841    tries = 0
842    HKL = G2lat.GenHBravais(dmin,ibrav,A,cctbx_args)
843    while len(HKL) > 2 and IndexPeaks(peaks,HKL)[0]:
844        Pwr = pwr - (tries % 2)
845        HKL = []
846        tries += 1
847        osmin = smin
848        oldA = A[:]
849        Vold = G2lat.calc_V(oldA)
850        OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
851        Vnew = G2lat.calc_V(A)
852        if Vnew > 2.0*Vold or Vnew < 2.:
853            A = ranAbyR(ibrav,oldA,tries+1,maxTries,ran2axis)
854            OK = False
855            continue
856        try:
857            HKL = G2lat.GenHBravais(dmin,ibrav,A,cctbx_args)
858        except FloatingPointError:
859            A = oldA
860            OK = False
861            break
862        if len(HKL) == 0: break                         #absurd cell obtained!
863        rat = (osmin-smin)/smin
864        if abs(rat) < 1.0e-5 or not OK: break
865        if tries > maxTries: break
866    if OK:
867        OK,smin,A,result = FitHKL(ibrav,peaks,A,2)
868        Peaks = np.array(peaks).T
869        H = Peaks[4:7]
870        try:
871            Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
872            peaks = Peaks.T
873        except FloatingPointError:
874            A = oldA
875       
876    M20,X20 = calc_M20(peaks,HKL,ifX20)
877    return len(HKL),M20,X20,A
878       
879def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1,ifX20=True,cctbx_args=None):
880    'needs a doc string'
881# dlg & ncMax are used for wx progress bar
882# A != 0 find the best A near input A,
883# A = 0 for random cell, volume normalized to V1;
884# returns number of generated hkls, M20, X20 & A for best found
885    mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7,7,7, 9,9,9,9, 10]
886    dmin = getDmin(peaks)-0.05
887    amin = 2.5
888    amax = 5.*getDmax(peaks)
889    Asave = []
890    GoOn = True
891    Skip = False
892    if A:
893        HKL = G2lat.GenHBravais(dmin,ibrav,A[:])
894        if len(HKL) > mHKL[ibrav]:
895            peaks = IndexPeaks(peaks,HKL)[1]
896            Asave.append([calc_M20(peaks,HKL,ifX20),A[:]])
897    tries = 0
898    while tries < Ntries and GoOn:
899        if A:
900            Abeg = ranAbyR(ibrav,A,tries+1,Ntries,ran2axis)
901            if ibrav > 12:         #monoclinic & triclinic
902                Abeg = ranAbyR(ibrav,A,tries/10+1,Ntries,ran2axis)
903        else:
904            Abeg = ranAbyV(ibrav,amin,amax,V1)
905        HKL = G2lat.GenHBravais(dmin,ibrav,Abeg)
906        Nc = len(HKL)
907        if Nc >= ncMax:
908            GoOn = False
909        else:
910            if dlg:
911                dlg.Raise()
912                GoOn = dlg.Update(100*Nc/ncMax)[0]
913                if Skip or not GoOn:
914                    GoOn = False
915                    break
916       
917        if IndexPeaks(peaks,HKL)[0] and len(HKL) > mHKL[ibrav]:
918            Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20,cctbx_args=cctbx_args)
919            Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]])
920            if ibrav in [9,10,11]:                          #A,B,or C-centered orthorhombic
921                for i in range(2):
922                    Abeg = rotOrthoA(Abeg[:])
923                    Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20,cctbx_args=cctbx_args)
924                    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
925                    peaks = IndexPeaks(peaks,HKL)[1]
926                    Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]])
927            # elif ibrav == 15:                      #C-centered monoclinic
928            #     Abeg = swapMonoA(Abeg[:])
929            #     Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20,cctbx_args=cctbx_args)
930            #     HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
931            #     peaks = IndexPeaks(peaks,HKL)[1]
932            #     Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]])
933        else:
934            break
935        Nc = len(HKL)
936        tries += 1
937    X = sortM20(Asave)
938    if X:
939        Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1],ifX20,cctbx_args=cctbx_args)
940        return GoOn,Skip,Lhkl,M20,X20,A       
941    else:
942        return GoOn,Skip,0,0,0,0
943       
944def monoCellReduce(ibrav,A):
945    'needs a doc string'
946    a,b,c,alp,bet,gam = G2lat.A2cell(A)
947    if bet < 90.: 
948        bet = 180.-bet     #always want obtuse beta
949        A = G2lat.cell2A([a,b,c,90.,bet,90.])
950    G,g = G2lat.A2Gmat(A)
951    if ibrav in [14,]:      #A-monoclinic - OK?
952        u = [0,0,-1]
953        v = [1,0,2]
954        anew = math.sqrt(np.dot(np.dot(v,g),v))
955        if anew < a:
956            cang = np.dot(np.dot(u,g),v)/(anew*c)
957            beta = acosd(-abs(cang))
958            if beta < 90.: beta = 180.-beta     #always want obtuse beta
959            A = G2lat.cell2A([anew,b,c,90,beta,90])
960    elif ibrav in [15,]:    #C-monoclinic - OK?
961        u = [-1,0,0]
962        v = [1,0,1]
963        cnew = math.sqrt(np.dot(np.dot(v,g),v))
964        if cnew < c:
965            cang = np.dot(np.dot(u,g),v)/(a*cnew)
966            beta = acosd(-abs(cang))
967            if beta < 90.: beta = 180.-beta     #always want obtuse beta
968            A = G2lat.cell2A([a,b,cnew,90,beta,90])
969    elif ibrav in [13,]:    #I-monoclinic - checked OK
970        uc = [0,0,-1]
971        vc = [1,0,2]
972        ua = [-1,0,0]
973        va = [2,0,1]
974        anew = math.sqrt(np.dot(np.dot(va,g),va))
975        cnew = math.sqrt(np.dot(np.dot(vc,g),vc))
976        if cnew < c:
977            cang = np.dot(np.dot(uc,g),v)/(a*cnew)
978            beta = acosd(-abs(cang))
979            if beta < 90.: beta = 180.-beta     #always want obtuse beta
980            A = G2lat.cell2A([a,b,cnew,90,beta,90])
981        if anew < a:
982            cang = np.dot(np.dot(ua,g),v)/(c*anew)
983            beta = acosd(-abs(cang))
984            if beta < 90.: beta = 180.-beta     #always want obtuse beta
985            A = G2lat.cell2A([anew,b,c,90,beta,90])
986    else:   #P
987        ua = [0,0,-1]
988        uc = [-1,0,0]
989        va = [1,0,1]
990        vc = [1,0,1]
991        anew = math.sqrt(np.dot(np.dot(va,g),va))
992        cnew = math.sqrt(np.dot(np.dot(vc,g),vc))
993        if anew < a:
994            cang = np.dot(np.dot(ua,g),va)/(anew*c)
995            beta = acosd(-abs(cang))
996            if beta < 90.: beta = 180.-beta     #always want obtuse beta
997            A = G2lat.cell2A([anew,b,c,90,beta,90])
998        if cnew < c:
999            cang = np.dot(np.dot(uc,g),vc)/(a*cnew)
1000            beta = acosd(-abs(cang))
1001            if beta < 90.: beta = 180.-beta     #always want obtuse beta
1002            A = G2lat.cell2A([a,b,cnew,90,beta,90])
1003    return A
1004
1005def DoIndexPeaks(peaks,controls,bravais,dlg,ifX20=True,
1006            timeout=None,M20_min=2.0,X20_max=None,return_Nc=False,
1007            cctbx_args=None):
1008    'needs a doc string'
1009   
1010    delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
1011    amin = 2.5
1012    amax = 5.0*getDmax(peaks)
1013    dmin = getDmin(peaks)-delt
1014    bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
1015        'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-A',
1016        'Orthorhombic-B','Orthorhombic-C',
1017        'Orthorhombic-P','Monoclinic-I','Monoclinic-A','Monoclinic-C','Monoclinic-P','Triclinic']
1018    tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
1019    N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,50,50,  100,100,100,100, 200]
1020    N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,2,2,   2,2,2,2,   4]
1021    Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,1,1,   2,2,2,2,   4]
1022    notUse = 0
1023    for peak in peaks:
1024        if not peak[2]:
1025            notUse += 1
1026    Nobs = len(peaks)-notUse
1027    zero,ncno = controls[1:3]
1028    ncMax = Nobs*ncno
1029    print ("%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax))
1030    print ("%s %.4f %s %d %s %d" % ('Zero =',zero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs))
1031    cells = []
1032    lastcell = np.zeros(7)
1033    for ibrav in range(len(bravaisNames)):
1034        begin = time.time()
1035        if bravais[ibrav]:
1036            print ('cell search for ',bravaisNames[ibrav])
1037            print ('      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test')
1038            V1 = controls[3]
1039            bestM20 = 0
1040            topM20 = 0
1041            cycle = 0
1042            while cycle < 5:
1043                if dlg:
1044                    dlg.Raise()
1045                    dlg.Update(0,newmsg=tries[cycle]+" cell search for "+bravaisNames[ibrav])
1046                try:
1047                    GoOn = True
1048                    while GoOn:                                                 #Loop over increment of volume
1049                        N2 = 0
1050                        while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
1051                            if timeout and time.time() - begin > timeout:
1052                                GoOn = False
1053                                break
1054
1055                            if ibrav > 2:
1056                                if not N2:
1057                                    A = []
1058                                    GoOn,Skip,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1,ifX20,cctbx_args=cctbx_args)
1059                                    if Skip:
1060                                        break
1061                                if A:
1062                                    GoOn,Skip,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0,ifX20,cctbx_args=cctbx_args)
1063                            else:
1064                                GoOn,Skip,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1,ifX20,cctbx_args=cctbx_args)
1065                            if Skip:
1066                                break
1067                            elif Nc >= ncMax:
1068                                GoOn = False
1069                                break
1070                            elif 3*Nc < Nobs:
1071                                N2 = 10
1072                                break
1073                            else:
1074                                if not GoOn:
1075                                    break
1076                                if 1.e6 > M20 > 1.0:    #exclude nonsense
1077                                    bestM20 = max(bestM20,M20)
1078                                    A = halfCell(ibrav,A[:],peaks)
1079                                    if ibrav in [13,14,15,16]:
1080                                        A = monoCellReduce(ibrav,A[:])
1081                                    HKL = G2lat.GenHBravais(dmin,ibrav,A)
1082                                    peaks = IndexPeaks(peaks,HKL)[1]
1083                                    a,b,c,alp,bet,gam = G2lat.A2cell(A)
1084                                    V = G2lat.calc_V(A)
1085                                    if (
1086                                        (M20 >= M20_min) and
1087                                        (X20_max is None or X20 <= X20_max)
1088                                    ):
1089                                        cell = [M20,X20,ibrav,a,b,c,alp,bet,gam,V,False,False]
1090                                        if return_Nc: cell.append(Nc)
1091                                        newcell = np.array(cell[3:10])
1092                                        if not np.allclose(newcell,lastcell):
1093                                            print ("%10.3f %3d %3d %10.5f %10.5f %10.5f %10.3f %10.3f %10.3f %10.2f %10.2f %s"
1094                                                %(M20,X20,Nc,a,b,c,alp,bet,gam,V,V1,bravaisNames[ibrav]))
1095                                            cells.append(cell)
1096                                        lastcell = np.array(cell[3:10])
1097                            if not GoOn:
1098                                break
1099                            N2 += 1
1100                        if Skip:
1101                            cycle = 10
1102                            GoOn = False
1103                            break
1104                        if ibrav < 13:
1105                            V1 *= 1.1
1106                        elif ibrav in range(13,18):
1107                            V1 *= 1.025
1108                        if not GoOn:
1109                            if bestM20 > topM20:
1110                                topM20 = bestM20
1111                                if cells:
1112                                    V1 = cells[0][9]
1113                                else:
1114                                    V1 = controls[3]
1115                                ncMax += Nobs
1116                                cycle += 1
1117                                print ('Restart search, new Max Nc = %d'%ncMax)
1118                            else:
1119                                cycle = 10
1120                finally:
1121                    pass
1122#                dlg.Destroy()
1123            print ('%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
1124                ', elapsed time = ',G2lat.sec2HMS(time.time()-begin)))
1125           
1126    if cells:
1127        return True,dmin,cells
1128    else:
1129        return False,0,[]
1130       
1131       
1132NeedTestData = True
1133def TestData():
1134    'needs a doc string'
1135    global NeedTestData
1136    NeedTestData = False
1137    global TestData
1138    TestData = [14, [7.,8.70,10.86,90.,102.95,90.], [7.76006,8.706215,10.865679,90.,102.947,90.],3,
1139        [[2.176562137832974, 761.60902227696033, True, True, 0, 0, 1, 10.591300714328161, 10.589436], 
1140        [3.0477561489789498, 4087.2956049071572, True, True, 1, 0, 0, 7.564238997554908, 7.562777], 
1141        [3.3254921120068524, 1707.0253890991009, True, True, 1, 0, -1, 6.932650301411212, 6.932718], 
1142        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.725106], 
1143        [4.0379791325512118, 1598.4321673135987, True, True, 1, 1, 0, 5.709789097440156, 5.70946], 
1144        [4.2511182350743937, 473.10955149057577, True, True, 1, 1, -1, 5.423637972781876, 5.42333], 
1145        [4.354684330373451, 569.88528280256071, True, True, 0, 0, 2, 5.2947091882172534, 5.294718],
1146        [4.723324574319177, 342.73882372499997, True, True, 1, 0, -2, 4.881681587039431, 4.881592], 
1147        [4.9014773581253994, 5886.3516356615492, True, True, 1, 1, 1, 4.704350709093183, 4.70413], 
1148        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.523829], 
1149        [5.2971997607389518, 1290.0229964239879, True, True, 0, 2, 0, 4.353139557169142, 4.353108], 
1150        [5.4161306205553847, 1472.5726977257755, True, True, 1, 1, -2, 4.257619398422479, 4.257944], 
1151        [5.7277364698554161, 1577.8791668322888, True, True, 0, 2, 1, 4.026169751907777, 4.026193], 
1152        [5.8500213058834163, 420.74210142657131, True, True, 1, 0, 2, 3.9420803081518443, 3.942219],
1153        [6.0986764166731708, 163.02160537058708, True, True, 2, 0, 0, 3.7814965150452537, 3.781389], 
1154        [6.1126665157702753, 943.25461245706833, True, True, 1, 2, 0, 3.772849962062199, 3.772764], 
1155        [6.2559260555056957, 250.55355015505376, True, True, 1, 2, -1, 3.6865353266375283, 3.686602], 
1156        [6.4226243128279892, 5388.5560141098349, True, True, 1, 1, 2, 3.5909481979190283, 3.591214], 
1157        [6.5346132446561134, 1951.6070344509026, True, True, 0, 0, 3, 3.5294722429440584, 3.529812], 
1158        [6.5586952135236443, 259.65938178131034, True, True, 2, 1, -1, 3.516526936765838, 3.516784], 
1159        [6.6509216222783722, 93.265376597376573, True, True, 2, 1, 0, 3.4678179073694952, 3.468369], 
1160        [6.7152737044107722, 289.39386813803162, True, True, 1, 2, 1, 3.4346235125812807, 3.434648], 
1161        [6.8594130457361899, 603.54959764648322, True, True, 0, 2, 2, 3.362534044860622, 3.362553], 
1162        [7.0511627728884454, 126.43246447656593, True, True, 0, 1, 3, 3.2712038721790675, 3.271181], 
1163        [7.077700845503319, 125.49742760019636, True, True, 1, 1, -3, 3.2589538988480626, 3.259037], 
1164        [7.099393757363675, 416.55444885434633, True, True, 1, 2, -2, 3.2490085228959193, 3.248951], 
1165        [7.1623933278642742, 369.27397110921817, True, True, 2, 1, -2, 3.2204673608202383, 3.220487], 
1166        [7.4121734953058924, 482.84120827021826, True, True, 2, 1, 1, 3.1120858221599876, 3.112308]]
1167        ]
1168    global TestData2
1169    TestData2 = [14, [0.15336547830008007, 0.017345499139401827, 0.008122368657493792, 0, 0.02893538955687591, 0], 3,
1170        [[2.176562137832974, 761.6090222769603, True, True, 0, 0, 1, 10.591300714328161, 11.095801], 
1171        [3.0477561489789498, 4087.295604907157, True, True, 0, 1, 0, 7.564238997554908, 7.592881], 
1172        [3.3254921120068524, 1707.025389099101, True, False, 0, 0, 0, 6.932650301411212, 0.0], 
1173        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.266192], 
1174        [4.037979132551212, 1598.4321673135987, True, False, 0, 0, 0, 5.709789097440156, 0.0], 
1175        [4.251118235074394, 473.10955149057577, True, True, 0, 0, 2, 5.423637972781876, 5.5479], 
1176        [4.354684330373451, 569.8852828025607, True, True, 0, 0, 2, 5.2947091882172534, 5.199754], 
1177        [4.723324574319177, 342.738823725, True, False, 0, 0, 0, 4.881681587039431, 0.0], 
1178        [4.901477358125399, 5886.351635661549, True, False, 0, 0, 0, 4.704350709093183, 0.0], 
1179        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.479534], 
1180        [5.297199760738952, 1290.022996423988, True, True, 0, 1, 0, 4.353139557169142, 4.345087],
1181        [5.416130620555385, 1472.5726977257755, True, False, 0, 0, 0, 4.257619398422479, 0.0], 
1182        [5.727736469855416, 1577.8791668322888, True, False, 0, 0, 0, 4.026169751907777, 0.0], 
1183        [5.850021305883416, 420.7421014265713, True, False, 0, 0, 0, 3.9420803081518443, 0.0], 
1184        [6.098676416673171, 163.02160537058708, True, True, 0, 2, 0, 3.7814965150452537, 3.796441], 
1185        [6.112666515770275, 943.2546124570683, True, False, 0, 0, 0, 3.772849962062199, 0.0], 
1186        [6.255926055505696, 250.55355015505376, True, True, 0, 0, 3, 3.6865353266375283, 3.6986], 
1187        [6.422624312827989, 5388.556014109835, True, True, 0, 2, 1, 3.5909481979190283, 3.592005], 
1188        [6.534613244656113, 191.6070344509026, True, True, 1, 0, -1, 3.5294722429440584, 3.546166], 
1189        [6.558695213523644, 259.65938178131034, True, True, 0, 0, 3, 3.516526936765838, 3.495428], 
1190        [6.650921622278372, 93.26537659737657, True, True, 0, 0, 3, 3.4678179073694952, 3.466503], 
1191        [6.715273704410772, 289.3938681380316, True, False, 0, 0, 0, 3.4346235125812807, 0.0], 
1192        [6.85941304573619, 603.5495976464832, True, True, 0, 1, 3, 3.362534044860622, 3.32509], 
1193        [7.051162772888445, 126.43246447656593, True, True, 0, 1, 2, 3.2712038721790675, 3.352121], 
1194        [7.077700845503319, 125.49742760019636, True, False, 0, 0, 0, 3.2589538988480626, 0.0], 
1195        [7.099393757363675, 416.55444885434633, True, False, 0, 0, 0, 3.2490085228959193, 0.0], 
1196        [7.162393327864274, 369.27397110921817, True, False, 0, 0, 0, 3.2204673608202383, 0.0], 
1197        [7.412173495305892, 482.84120827021826, True, True, 0, 2, 2, 3.112085822159976, 3.133096]]
1198        ]
1199#
1200def test0():
1201    if NeedTestData: TestData()
1202    ibrav,cell,bestcell,Pwr,peaks = TestData
1203    print ('best cell:',bestcell)
1204    print ('old cell:',cell)
1205    Peaks = np.array(peaks)
1206    HKL = Peaks[4:7]
1207    print (calc_M20(peaks,HKL))
1208    A = G2lat.cell2A(cell)
1209    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
1210    print ('new cell:',G2lat.A2cell(A))   
1211    print ('x:',result[0])
1212    print ('cov_x:',result[1])
1213    print ('infodict:')
1214    for item in result[2]:
1215        print (item,result[2][item])
1216    print ('msg:',result[3])
1217    print ('ier:',result[4])
1218    result = refinePeaks(peaks,ibrav,A)
1219    N,M20,X20,A = result
1220    print ('refinePeaks:',N,M20,X20,G2lat.A2cell(A))
1221    print ('compare bestcell:',bestcell)
1222#
1223def test1():
1224    if NeedTestData: TestData()
1225    ibrav,A,Pwr,peaks = TestData2
1226    print ('bad cell:',G2lat.A2cell(A))
1227    print ('FitHKL')
1228    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
1229    result = refinePeaks(peaks,ibrav,A)
1230    N,M20,X20,A = result
1231    print ('refinePeaks:',N,M20,X20,A)
1232#    Peaks = np.array(peaks)
1233#    HKL = Peaks[4:7]
1234#    print calc_M20(peaks,HKL)
1235#    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
1236#    print 'new cell:',G2lat.A2cell(A)   
1237#    print 'x:',result[0]
1238#    print 'cov_x:',result[1]
1239#    print 'infodict:'
1240#    for item in result[2]:
1241#        print item,result[2][item]
1242#    print 'msg:',result[3]
1243#    print 'ier:',result[4]
1244   
1245#
1246if __name__ == '__main__':
1247    test0()
1248    test1()
1249#    test2()
1250#    test3()
1251#    test4()
1252#    test5()
1253#    test6()
1254#    test7()
1255#    test8()
1256    print ("OK")
Note: See TracBrowser for help on using the repository browser.