source: trunk/GSASIIindex.py @ 3633

Last change on this file since 3633 was 3633, checked in by vondreele, 5 years ago

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