source: trunk/GSASIIindex.py @ 3136

Last change on this file since 3136 was 3136, checked in by vondreele, 4 years ago

make GSAS-II python 3.6 compliant & preserve python 2.7 use;changes:
do from future import division, print_function for all GSAS-II py sources
all menu items revised to be py 2.7/3.6 compliant
all wx.OPEN --> wx.FD_OPEN in file dialogs
all integer divides (typically for image pixel math) made explicit with ; ambiguous ones made floats as appropriate
all print "stuff" --> print (stuff)
all print >> pFile,'stuff' --> pFile.writeCIFtemplate('stuff')
all read file opens made explicit 'r' or 'rb'
all cPickle imports made for py2.7 or 3.6 as cPickle or _pickle; test for '2' platform.version_tuple[0] for py 2.7
define cPickleload to select load(fp) or load(fp,encoding='latin-1') for loading gpx files; provides cross compatibility between py 2.7/3.6 gpx files
make dict.keys() as explicit list(dict.keys()) as needed (NB: possible source of remaining py3.6 bugs)
make zip(a,b) as explicit list(zip(a,b)) as needed (NB: possible source of remaining py3.6 bugs)
select unichr/chr according test for '2' platform.version_tuple[0] for py 2.7 (G2pwdGUI * G2plot) for special characters
select wg.EVT_GRID_CELL_CHANGE (classic) or wg.EVT_GRID_CELL_CHANGED (phoenix) in grid Bind
maxint --> maxsize; used in random number stuff
raise Exception,"stuff" --> raise Exception("stuff")
wx 'classic' sizer.DeleteWindows?() or 'phoenix' sizer.Clear(True)
wx 'classic' SetToolTipString?(text) or 'phoenix' SetToolTip?(wx.ToolTip?(text)); define SetToolTipString?(self,text) to handle the choice in plots
status.SetFields? --> status.SetStatusText?
'classic' AddSimpleTool? or 'phoenix' self.AddTool? for plot toolbar; Bind different as well
define GetItemPydata? as it doesn't exist in wx 'phoenix'
allow python versions 2.7 & 3.6 to run GSAS-II
Bind override commented out - no logging capability (NB: remove all logging code?)
all import ContentsValidator? open filename & test if valid then close; filepointer removed from Reader
binary importers (mostly images) test for 'byte' type & convert as needed to satisfy py 3.6 str/byte rules

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 42.5 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: 2017-10-23 16:39:16 +0000 (Mon, 23 Oct 2017) $
6# $Author: vondreele $
7# $Revision: 3136 $
8# $URL: trunk/GSASIIindex.py $
9# $Id: GSASIIindex.py 3136 2017-10-23 16:39:16Z 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: 3136 $")
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]:     #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 [11,12]:        #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]:       #orthorhombic - F,I,C,P - a<b<c convention
143        abc = [ranaxis(dmin,dmax),ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
144        if Bravais in [7,8,10]:
145            abc.sort()
146        a = abc[0]
147        b = abc[1]
148        c = abc[2]
149        alp = bet = gam = 90
150    elif Bravais in [11,12]:        #monoclinic - C,P - a<c convention
151        ac = [ranaxis(dmin,dmax),ranaxis(dmin,dmax)]
152        if Bravais == 12:
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.02       #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[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]:
424        return [values[0],values[1],values[2],0,0,0]
425    elif ibrav in [11,12]:
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]:
437        return [A[0],A[1],A[2]]
438    elif ibrav in [11,12]:
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]:
447        Nskip = 3
448    elif ibrav in [11,12]:
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]:
479            derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
480        elif ibrav in [11,12]:
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]:
511        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
512    elif ibrav in [11,12]:
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    result = so.leastsq(errFitZ,values,Dfun=dervFitZ,full_output=True,ftol=0.0001,
529        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Z,Zref))
530    A = Values2A(ibrav,result[0][:6])
531    if Zref:
532        Z = result[0][-1]
533    chisq = np.sum(errFitZ(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],wave,Z,Zref)**2)
534    return True,chisq,A,Z,result
535   
536def errFitZSS(values,ibrav,d,H,tth,wave,vec,Vref,Z,Zref):
537    Zero = Z
538    if Zref:   
539        Zero = values[-1]
540    A = Values2A(ibrav,values)
541    Vec = Values2Vec(ibrav,vec,Vref,values)
542    Qo = 1./d**2
543    Qc = G2lat.calc_rDsqZSS(H,A,Vec,Zero,tth,wave)
544    return (Qo-Qc)
545   
546def dervFitZSS(values,ibrav,d,H,tth,wave,vec,Vref,Z,Zref):
547    A = Values2A(ibrav,values)
548    Vec = Values2Vec(ibrav,vec,Vref,values)
549    HM = H[:3]+(H[3][:,np.newaxis]*Vec).T
550    if ibrav in [3,4,]:
551        derv = [HM[0]*HM[0]+HM[1]*HM[1]+HM[0]*HM[1],HM[2]*HM[2]]
552    elif ibrav in [5,6]:
553        derv = [HM[0]*HM[0]+HM[1]*HM[1],HM[2]*HM[2]]
554    elif ibrav in [7,8,9,10]:
555        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2]]
556    elif ibrav in [11,12]:
557        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[2]]
558    else:
559        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]]
560    if Vref[0]:
561        derv.append(2.*A[0]*HM[0]*H[3]+A[3]*HM[1]*H[3]+A[4]*HM[2]*H[3])
562    if Vref[1]:
563        derv.append(2.*A[1]*HM[1]*H[3]+A[3]*HM[0]*H[3]+A[5]*HM[2]*H[3])
564    if Vref[2]:
565        derv.append(2.*A[2]*HM[2]*H[3]+A[4]*HM[1]*H[3]+A[5]*HM[0]*H[3])   
566    if Zref:
567        derv.append(npsind(tth)*2.0*rpd/wave**2)
568    derv = -np.array(derv)
569    return derv.T
570   
571def FitHKLZSS(wave,ibrav,peaks,A,V,Vref,Z,Zref):
572    'needs a doc string'
573   
574    Peaks = np.array(peaks).T   
575    values = A2values(ibrav,A)
576    for v,r in zip(V,Vref):
577        if r:
578            values.append(v)
579    if Zref:
580        values.append(Z)
581    result = so.leastsq(errFitZSS,values,Dfun=dervFitZSS,full_output=True,ftol=1.e-6,
582        args=(ibrav,Peaks[8],Peaks[4:8],Peaks[0],wave,V,Vref,Z,Zref))
583    A = Values2A(ibrav,result[0])
584    Vec = Values2Vec(ibrav,V,Vref,result[0])
585    if Zref:
586        Z = result[0][-1]
587    chisq = np.sum(errFitZSS(result[0],ibrav,Peaks[8],Peaks[4:8],Peaks[0],wave,Vec,Vref,Z,Zref)**2) 
588    return True,chisq,A,Vec,Z,result
589   
590def errFitT(values,ibrav,d,H,tof,difC,Z,Zref):
591    Zero = Z
592    if Zref:   
593        Zero = values[-1]
594    A = Values2A(ibrav,values)
595    Qo = 1./d**2
596    Qc = G2lat.calc_rDsqT(H,A,Zero,tof,difC)
597    return (Qo-Qc)
598   
599def dervFitT(values,ibrav,d,H,tof,difC,Z,Zref):
600    if ibrav in [0,1,2]:
601        derv = [H[0]*H[0]+H[1]*H[1]+H[2]*H[2],]
602    elif ibrav in [3,4,]:
603        derv = [H[0]*H[0]+H[1]*H[1]+H[0]*H[1],H[2]*H[2]]
604    elif ibrav in [5,6]:
605        derv = [H[0]*H[0]+H[1]*H[1],H[2]*H[2]]
606    elif ibrav in [7,8,9,10]:
607        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2]]
608    elif ibrav in [11,12]:
609        derv = [H[0]*H[0],H[1]*H[1],H[2]*H[2],H[0]*H[2]]
610    else:
611        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]]
612    if Zref:
613        derv.append(np.ones_like(d)/difC)
614    derv = np.array(derv)
615    return derv.T
616   
617def FitHKLT(difC,ibrav,peaks,A,Z,Zref):
618    'needs a doc string'
619   
620    Peaks = np.array(peaks).T   
621    values = A2values(ibrav,A)
622    if Zref:
623        values.append(Z)
624    result = so.leastsq(errFitT,values,Dfun=dervFitT,full_output=True,ftol=0.0001,
625        args=(ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Z,Zref))
626    A = Values2A(ibrav,result[0])
627    if Zref:
628        Z = result[0][-1]
629    chisq = np.sum(errFitT(result[0],ibrav,Peaks[7],Peaks[4:7],Peaks[0],difC,Z,Zref)**2)
630    return True,chisq,A,Z,result
631   
632def errFitTSS(values,ibrav,d,H,tof,difC,vec,Vref,Z,Zref):
633    Zero = Z
634    if Zref:   
635        Zero = values[-1]
636    A = Values2A(ibrav,values)
637    Vec = Values2Vec(ibrav,vec,Vref,values)
638    Qo = 1./d**2
639    Qc = G2lat.calc_rDsqTSS(H,A,Vec,Zero,tof,difC)
640    return (Qo-Qc)
641   
642def dervFitTSS(values,ibrav,d,H,tof,difC,vec,Vref,Z,Zref):
643    A = Values2A(ibrav,values)
644    Vec = Values2Vec(ibrav,vec,Vref,values)
645    HM = H[:3]+(H[3][:,np.newaxis]*Vec).T
646    if ibrav in [3,4,]:
647        derv = [HM[0]*HM[0]+HM[1]*HM[1]+HM[0]*HM[1],HM[2]*HM[2]]
648    elif ibrav in [5,6]:
649        derv = [HM[0]*HM[0]+HM[1]*HM[1],HM[2]*HM[2]]
650    elif ibrav in [7,8,9,10]:
651        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2]]
652    elif ibrav in [11,12]:
653        derv = [HM[0]*HM[0],HM[1]*HM[1],HM[2]*HM[2],HM[0]*HM[2]]
654    else:
655        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]]
656    if Vref[0]:
657        derv.append(2.*A[0]*HM[0]*H[3]+A[3]*HM[1]*H[3]+A[4]*HM[2]*H[3])
658    if Vref[1]:
659        derv.append(2.*A[1]*HM[1]*H[3]+A[3]*HM[0]*H[3]+A[5]*HM[2]*H[3])
660    if Vref[2]:
661        derv.append(2.*A[2]*HM[2]*H[3]+A[4]*HM[1]*H[3]+A[5]*HM[0]*H[3])   
662    if Zref:
663        derv.append(np.ones_like(d)/difC)
664    derv = np.array(derv)
665    return derv.T
666   
667def FitHKLTSS(difC,ibrav,peaks,A,V,Vref,Z,Zref):
668    'needs a doc string'
669   
670    Peaks = np.array(peaks).T   
671    values = A2values(ibrav,A)
672    for v,r in zip(V,Vref):
673        if r:
674            values.append(v)
675    if Zref:
676        values.append(Z)
677    result = so.leastsq(errFitTSS,values,Dfun=dervFitTSS,full_output=True,ftol=0.0001,
678        args=(ibrav,Peaks[8],Peaks[4:8],Peaks[0],difC,V,Vref,Z,Zref))
679    A = Values2A(ibrav,result[0])
680    Vec = Values2Vec(ibrav,V,Vref,result[0])
681    if Zref:
682        Z = result[0][-1]
683    chisq = np.sum(errFitTSS(result[0],ibrav,Peaks[8],Peaks[4:8],Peaks[0],difC,V,Vref,Z,Zref)**2)
684    return True,chisq,A,Vec,Z,result
685               
686def rotOrthoA(A):
687    'needs a doc string'
688    return [A[1],A[2],A[0],0,0,0]
689   
690def swapMonoA(A):
691    'needs a doc string'
692    return [A[2],A[1],A[0],0,A[4],0]
693   
694def oddPeak(indx,peaks):
695    'needs a doc string'
696    noOdd = True
697    for peak in peaks:
698        H = peak[4:7]
699        if H[indx] % 2:
700            noOdd = False
701    return noOdd
702   
703def halfCell(ibrav,A,peaks):
704    'needs a doc string'
705    if ibrav in [0,1,2]:
706        if oddPeak(0,peaks):
707            A[0] *= 2
708            A[1] = A[2] = A[0]
709    elif ibrav in [3,4,5,6]:
710        if oddPeak(0,peaks):
711            A[0] *= 2
712            A[1] = A[0]
713        if oddPeak(2,peaks):
714            A[2] *=2
715    else:
716        if oddPeak(0,peaks):
717            A[0] *=2
718        if oddPeak(1,peaks):
719            A[1] *=2
720        if oddPeak(2,peaks):
721            A[2] *=2
722    return A
723   
724def getDmin(peaks):
725    'needs a doc string'
726    return peaks[-1][-2]
727   
728def getDmax(peaks):
729    'needs a doc string'
730    return peaks[0][-2]
731   
732def refinePeaksZ(peaks,wave,ibrav,A,Zero,ZeroRef):
733    'needs a doc string'
734    dmin = getDmin(peaks)
735    OK,smin,Aref,Z,result = FitHKLZ(wave,ibrav,peaks,A,Zero,ZeroRef)
736    Peaks = np.array(peaks).T
737    H = Peaks[4:7]
738    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqZ(H,Aref,Z,Peaks[0],wave))
739    peaks = Peaks.T   
740    HKL = G2lat.GenHBravais(dmin,ibrav,A)
741    M20,X20 = calc_M20(peaks,HKL)
742    return len(HKL),M20,X20,Aref,Z
743   
744def refinePeaksT(peaks,difC,ibrav,A,Zero,ZeroRef):
745    'needs a doc string'
746    dmin = getDmin(peaks)
747    OK,smin,Aref,Z,result = FitHKLT(difC,ibrav,peaks,A,Zero,ZeroRef)
748    Peaks = np.array(peaks).T
749    H = Peaks[4:7]
750    Peaks[8] = 1./np.sqrt(G2lat.calc_rDsqT(H,Aref,Z,Peaks[0],difC))
751    peaks = Peaks.T   
752    HKL = G2lat.GenHBravais(dmin,ibrav,A)
753    M20,X20 = calc_M20(peaks,HKL)
754    return len(HKL),M20,X20,Aref,Z
755   
756def refinePeaksZSS(peaks,wave,Inst,SGData,SSGData,maxH,ibrav,A,vec,vecRef,Zero,ZeroRef):
757    'needs a doc string'
758    dmin = getDmin(peaks)
759    OK,smin,Aref,Vref,Z,result = FitHKLZSS(wave,ibrav,peaks,A,vec,vecRef,Zero,ZeroRef)
760    Peaks = np.array(peaks).T
761    H = Peaks[4:8]
762    Peaks[9] = 1./np.sqrt(G2lat.calc_rDsqZSS(H,Aref,Vref,Z,Peaks[0],wave))  #H,A,vec,Z,tth,lam
763    peaks = Peaks.T   
764    HKL =  G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vref,maxH,Aref)
765    M20,X20 = calc_M20SS(peaks,HKL)
766    return len(HKL),M20,X20,Aref,Vref,Z
767   
768def refinePeaksTSS(peaks,difC,Inst,SGData,SSGData,maxH,ibrav,A,vec,vecRef,Zero,ZeroRef):
769    'needs a doc string'
770    dmin = getDmin(peaks)
771    OK,smin,Aref,Vref,Z,result = FitHKLTSS(difC,ibrav,peaks,A,vec,vecRef,Zero,ZeroRef)
772    Peaks = np.array(peaks).T
773    H = Peaks[4:8]
774    Peaks[9] = 1./np.sqrt(G2lat.calc_rDsqTSS(H,Aref,Vref,Z,Peaks[0],difC))
775    peaks = Peaks.T   
776    HKL =  G2pwd.getHKLMpeak(dmin,Inst,SGData,SSGData,Vref,maxH,Aref)
777    HKL = G2lat.GenHBravais(dmin,ibrav,A)
778    M20,X20 = calc_M20SS(peaks,HKL)
779    return len(HKL),M20,X20,Aref,Vref,Z
780   
781def refinePeaks(peaks,ibrav,A,ifX20=True):
782    'needs a doc string'
783    dmin = getDmin(peaks)
784    smin = 1.0e10
785    pwr = 8
786    maxTries = 10
787    OK = False
788    tries = 0
789    HKL = G2lat.GenHBravais(dmin,ibrav,A)
790    while len(HKL) > 2 and IndexPeaks(peaks,HKL)[0]:
791        Pwr = pwr - (tries % 2)
792        HKL = []
793        tries += 1
794        osmin = smin
795        oldA = A[:]
796        Vold = G2lat.calc_V(oldA)
797        OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
798        Vnew = G2lat.calc_V(A)
799        if Vnew > 2.0*Vold or Vnew < 2.:
800            A = ranAbyR(ibrav,oldA,tries+1,maxTries,ran2axis)
801            OK = False
802            continue
803        try:
804            HKL = G2lat.GenHBravais(dmin,ibrav,A)
805        except FloatingPointError:
806            A = oldA
807            OK = False
808            break
809        if len(HKL) == 0: break                         #absurd cell obtained!
810        rat = (osmin-smin)/smin
811        if abs(rat) < 1.0e-5 or not OK: break
812        if tries > maxTries: break
813    if OK:
814        OK,smin,A,result = FitHKL(ibrav,peaks,A,2)
815        Peaks = np.array(peaks).T
816        H = Peaks[4:7]
817        try:
818            Peaks[8] = 1./np.sqrt(G2lat.calc_rDsq(H,A))
819            peaks = Peaks.T
820        except FloatingPointError:
821            A = oldA
822       
823    M20,X20 = calc_M20(peaks,HKL,ifX20)
824    return len(HKL),M20,X20,A
825       
826def findBestCell(dlg,ncMax,A,Ntries,ibrav,peaks,V1,ifX20=True):
827    'needs a doc string'
828# dlg & ncMax are used for wx progress bar
829# A != 0 find the best A near input A,
830# A = 0 for random cell, volume normalized to V1;
831# returns number of generated hkls, M20, X20 & A for best found
832    mHKL = [3,3,3, 5,5, 5,5, 7,7,7,7, 9,9, 10]
833    dmin = getDmin(peaks)-0.05
834    amin = 2.5
835    amax = 5.*getDmax(peaks)
836    Asave = []
837    GoOn = True
838    Skip = False
839    if A:
840        HKL = G2lat.GenHBravais(dmin,ibrav,A[:])
841        if len(HKL) > mHKL[ibrav]:
842            peaks = IndexPeaks(peaks,HKL)[1]
843            Asave.append([calc_M20(peaks,HKL,ifX20),A[:]])
844    tries = 0
845    while tries < Ntries and GoOn:
846        if A:
847            Abeg = ranAbyR(ibrav,A,tries+1,Ntries,ran2axis)
848            if ibrav in [11,12,13]:         #monoclinic & triclinic
849                Abeg = ranAbyR(ibrav,A,tries/10+1,Ntries,ran2axis)
850        else:
851            Abeg = ranAbyV(ibrav,amin,amax,V1)
852        HKL = G2lat.GenHBravais(dmin,ibrav,Abeg)
853        Nc = len(HKL)
854        if Nc >= ncMax:
855            GoOn = False
856        else:
857            if dlg:
858#                GoOn,Skip = dlg.Update(100*Nc/ncMax)   #wx error doesn't work in 32 bit versions!
859                GoOn = dlg.Update(100*Nc/ncMax)[0]
860                if Skip or not GoOn:
861                    GoOn = False
862                    break
863       
864        if IndexPeaks(peaks,HKL)[0] and len(HKL) > mHKL[ibrav]:
865            Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20)
866            Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]])
867            if ibrav == 9:                          #C-centered orthorhombic
868                for i in range(2):
869                    Abeg = rotOrthoA(Abeg[:])
870                    Lhkl,M20,X20,Aref = refinePeaks(peaks,ibrav,Abeg,ifX20)
871                    HKL = G2lat.GenHBravais(dmin,ibrav,Aref)
872                    peaks = IndexPeaks(peaks,HKL)[1]
873                    Asave.append([calc_M20(peaks,HKL,ifX20),Aref[:]])
874            elif ibrav == 11:                      #C-centered monoclinic
875                Abeg = swapMonoA(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        else:
881            break
882        Nc = len(HKL)
883        tries += 1
884    X = sortM20(Asave)
885    if X:
886        Lhkl,M20,X20,A = refinePeaks(peaks,ibrav,X[0][1],ifX20)
887        return GoOn,Skip,Lhkl,M20,X20,A       
888    else:
889        return GoOn,Skip,0,0,0,0
890       
891def monoCellReduce(ibrav,A):
892    'needs a doc string'
893    a,b,c,alp,bet,gam = G2lat.A2cell(A)
894    G,g = G2lat.A2Gmat(A)
895    if ibrav in [11]:
896        u = [0,0,-1]
897        v = [1,0,2]
898        anew = math.sqrt(np.dot(np.dot(v,g),v))
899        if anew < a:
900            cang = np.dot(np.dot(u,g),v)/(anew*c)
901            beta = acosd(-abs(cang))
902            A = G2lat.cell2A([anew,b,c,90,beta,90])
903    else:
904        u = [-1,0,0]
905        v = [1,0,1]
906        cnew = math.sqrt(np.dot(np.dot(v,g),v))
907        if cnew < c:
908            cang = np.dot(np.dot(u,g),v)/(a*cnew)
909            beta = acosd(-abs(cang))
910            A = G2lat.cell2A([a,b,cnew,90,beta,90])
911    return A
912
913def DoIndexPeaks(peaks,controls,bravais,dlg,ifX20=True):
914    'needs a doc string'
915   
916    delt = 0.005                                     #lowest d-spacing cushion - can be fixed?
917    amin = 2.5
918    amax = 5.0*getDmax(peaks)
919    dmin = getDmin(peaks)-delt
920    bravaisNames = ['Cubic-F','Cubic-I','Cubic-P','Trigonal-R','Trigonal/Hexagonal-P',
921        'Tetragonal-I','Tetragonal-P','Orthorhombic-F','Orthorhombic-I','Orthorhombic-C',
922        'Orthorhombic-P','Monoclinic-C','Monoclinic-P','Triclinic']
923    tries = ['1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th']
924    N1s = [1,1,1,   5,5,  5,5, 50,50,50,50,  50,50, 200]
925    N2s = [1,1,1,   2,2,  2,2,     2,2,2,2,   2,2,   4]
926    Nm  = [1,1,1,   1,1,  1,1,     1,1,1,1,   2,2,   4]
927    notUse = 0
928    for peak in peaks:
929        if not peak[2]:
930            notUse += 1
931    Nobs = len(peaks)-notUse
932    zero,ncno = controls[1:3]
933    ncMax = Nobs*ncno
934    print ("%s %8.3f %8.3f" % ('lattice parameter range = ',amin,amax))
935    print ("%s %.4f %s %d %s %d" % ('Zero =',zero,'Nc/No max =',ncno,' Max Nc =',ncno*Nobs))
936    cells = []
937    lastcell = np.zeros(7)
938    for ibrav in range(14):
939        begin = time.time()
940        if bravais[ibrav]:
941            print ('cell search for ',bravaisNames[ibrav])
942            print ('      M20  X20  Nc       a          b          c        alpha       beta      gamma     volume      V-test')
943            V1 = controls[3]
944            bestM20 = 0
945            topM20 = 0
946            cycle = 0
947            while cycle < 5:
948                if dlg:
949                    dlg.Update(0,newmsg=tries[cycle]+" cell search for "+bravaisNames[ibrav])
950                try:
951                    GoOn = True
952                    while GoOn:                                                 #Loop over increment of volume
953                        N2 = 0
954                        while N2 < N2s[ibrav]:                                  #Table 2 step (iii)               
955                            if ibrav > 2:
956                                if not N2:
957                                    A = []
958                                    GoOn,Skip,Nc,M20,X20,A = findBestCell(dlg,ncMax,A,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1,ifX20)
959                                    if Skip:
960                                        break
961                                if A:
962                                    GoOn,Skip,Nc,M20,X20,A = findBestCell(dlg,ncMax,A[:],N1s[ibrav],ibrav,peaks,0,ifX20)
963                            else:
964                                GoOn,Skip,Nc,M20,X20,A = findBestCell(dlg,ncMax,0,Nm[ibrav]*N1s[ibrav],ibrav,peaks,V1,ifX20)
965                            if Skip:
966                                break
967                            elif Nc >= ncMax:
968                                GoOn = False
969                                break
970                            elif 3*Nc < Nobs:
971                                N2 = 10
972                                break
973                            else:
974                                if not GoOn:
975                                    break
976                                if 1.e6 > M20 > 1.0:    #exclude nonsense
977                                    bestM20 = max(bestM20,M20)
978                                    A = halfCell(ibrav,A[:],peaks)
979                                    if ibrav in [12]:
980                                        A = monoCellReduce(ibrav,A[:])
981                                    HKL = G2lat.GenHBravais(dmin,ibrav,A)
982                                    peaks = IndexPeaks(peaks,HKL)[1]
983                                    a,b,c,alp,bet,gam = G2lat.A2cell(A)
984                                    V = G2lat.calc_V(A)
985                                    if M20 >= 2.0:
986                                        cell = [M20,X20,ibrav,a,b,c,alp,bet,gam,V,False,False]
987                                        newcell = np.array(cell[3:10])
988                                        if not np.allclose(newcell,lastcell):
989                                            print ("%10.3f %3d %3d %10.5f %10.5f %10.5f %10.3f %10.3f %10.3f %10.2f %10.2f"  \
990                                                %(M20,X20,Nc,a,b,c,alp,bet,gam,V,V1))
991                                            cells.append(cell)
992                                        lastcell = np.array(cell[3:10])
993                            if not GoOn:
994                                break
995                            N2 += 1
996                        if Skip:
997                            cycle = 10
998                            GoOn = False
999                            break
1000                        if ibrav < 11:
1001                            V1 *= 1.1
1002                        elif ibrav in range(11,14):
1003                            V1 *= 1.05
1004                        if not GoOn:
1005                            if bestM20 > topM20:
1006                                topM20 = bestM20
1007                                if cells:
1008                                    V1 = cells[0][9]
1009                                else:
1010                                    V1 = 25
1011                                ncMax += Nobs
1012                                cycle += 1
1013                                print ('Restart search, new Max Nc = %d'%ncMax)
1014                            else:
1015                                cycle = 10
1016                finally:
1017                    pass
1018#                dlg.Destroy()
1019            print ('%s%s%s%s'%('finished cell search for ',bravaisNames[ibrav], \
1020                ', elapsed time = ',G2lat.sec2HMS(time.time()-begin)))
1021           
1022    if cells:
1023        return True,dmin,cells
1024    else:
1025        return False,0,[]
1026       
1027       
1028NeedTestData = True
1029def TestData():
1030    'needs a doc string'
1031    global NeedTestData
1032    NeedTestData = False
1033    global TestData
1034    TestData = [12, [7.,8.70,10.86,90.,102.95,90.], [7.76006,8.706215,10.865679,90.,102.947,90.],3,
1035        [[2.176562137832974, 761.60902227696033, True, True, 0, 0, 1, 10.591300714328161, 10.589436], 
1036        [3.0477561489789498, 4087.2956049071572, True, True, 1, 0, 0, 7.564238997554908, 7.562777], 
1037        [3.3254921120068524, 1707.0253890991009, True, True, 1, 0, -1, 6.932650301411212, 6.932718], 
1038        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.725106], 
1039        [4.0379791325512118, 1598.4321673135987, True, True, 1, 1, 0, 5.709789097440156, 5.70946], 
1040        [4.2511182350743937, 473.10955149057577, True, True, 1, 1, -1, 5.423637972781876, 5.42333], 
1041        [4.354684330373451, 569.88528280256071, True, True, 0, 0, 2, 5.2947091882172534, 5.294718],
1042        [4.723324574319177, 342.73882372499997, True, True, 1, 0, -2, 4.881681587039431, 4.881592], 
1043        [4.9014773581253994, 5886.3516356615492, True, True, 1, 1, 1, 4.704350709093183, 4.70413], 
1044        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.523829], 
1045        [5.2971997607389518, 1290.0229964239879, True, True, 0, 2, 0, 4.353139557169142, 4.353108], 
1046        [5.4161306205553847, 1472.5726977257755, True, True, 1, 1, -2, 4.257619398422479, 4.257944], 
1047        [5.7277364698554161, 1577.8791668322888, True, True, 0, 2, 1, 4.026169751907777, 4.026193], 
1048        [5.8500213058834163, 420.74210142657131, True, True, 1, 0, 2, 3.9420803081518443, 3.942219],
1049        [6.0986764166731708, 163.02160537058708, True, True, 2, 0, 0, 3.7814965150452537, 3.781389], 
1050        [6.1126665157702753, 943.25461245706833, True, True, 1, 2, 0, 3.772849962062199, 3.772764], 
1051        [6.2559260555056957, 250.55355015505376, True, True, 1, 2, -1, 3.6865353266375283, 3.686602], 
1052        [6.4226243128279892, 5388.5560141098349, True, True, 1, 1, 2, 3.5909481979190283, 3.591214], 
1053        [6.5346132446561134, 1951.6070344509026, True, True, 0, 0, 3, 3.5294722429440584, 3.529812], 
1054        [6.5586952135236443, 259.65938178131034, True, True, 2, 1, -1, 3.516526936765838, 3.516784], 
1055        [6.6509216222783722, 93.265376597376573, True, True, 2, 1, 0, 3.4678179073694952, 3.468369], 
1056        [6.7152737044107722, 289.39386813803162, True, True, 1, 2, 1, 3.4346235125812807, 3.434648], 
1057        [6.8594130457361899, 603.54959764648322, True, True, 0, 2, 2, 3.362534044860622, 3.362553], 
1058        [7.0511627728884454, 126.43246447656593, True, True, 0, 1, 3, 3.2712038721790675, 3.271181], 
1059        [7.077700845503319, 125.49742760019636, True, True, 1, 1, -3, 3.2589538988480626, 3.259037], 
1060        [7.099393757363675, 416.55444885434633, True, True, 1, 2, -2, 3.2490085228959193, 3.248951], 
1061        [7.1623933278642742, 369.27397110921817, True, True, 2, 1, -2, 3.2204673608202383, 3.220487], 
1062        [7.4121734953058924, 482.84120827021826, True, True, 2, 1, 1, 3.1120858221599876, 3.112308]]
1063        ]
1064    global TestData2
1065    TestData2 = [12, [0.15336547830008007, 0.017345499139401827, 0.008122368657493792, 0, 0.02893538955687591, 0], 3,
1066        [[2.176562137832974, 761.6090222769603, True, True, 0, 0, 1, 10.591300714328161, 11.095801], 
1067        [3.0477561489789498, 4087.295604907157, True, True, 0, 1, 0, 7.564238997554908, 7.592881], 
1068        [3.3254921120068524, 1707.025389099101, True, False, 0, 0, 0, 6.932650301411212, 0.0], 
1069        [3.428121546163426, 2777.5082170150563, True, True, 0, 1, 1, 6.725163158013632, 6.266192], 
1070        [4.037979132551212, 1598.4321673135987, True, False, 0, 0, 0, 5.709789097440156, 0.0], 
1071        [4.251118235074394, 473.10955149057577, True, True, 0, 0, 2, 5.423637972781876, 5.5479], 
1072        [4.354684330373451, 569.8852828025607, True, True, 0, 0, 2, 5.2947091882172534, 5.199754], 
1073        [4.723324574319177, 342.738823725, True, False, 0, 0, 0, 4.881681587039431, 0.0], 
1074        [4.901477358125399, 5886.351635661549, True, False, 0, 0, 0, 4.704350709093183, 0.0], 
1075        [5.0970774474587275, 3459.7541692903033, True, True, 0, 1, 2, 4.523933797797693, 4.479534], 
1076        [5.297199760738952, 1290.022996423988, True, True, 0, 1, 0, 4.353139557169142, 4.345087],
1077        [5.416130620555385, 1472.5726977257755, True, False, 0, 0, 0, 4.257619398422479, 0.0], 
1078        [5.727736469855416, 1577.8791668322888, True, False, 0, 0, 0, 4.026169751907777, 0.0], 
1079        [5.850021305883416, 420.7421014265713, True, False, 0, 0, 0, 3.9420803081518443, 0.0], 
1080        [6.098676416673171, 163.02160537058708, True, True, 0, 2, 0, 3.7814965150452537, 3.796441], 
1081        [6.112666515770275, 943.2546124570683, True, False, 0, 0, 0, 3.772849962062199, 0.0], 
1082        [6.255926055505696, 250.55355015505376, True, True, 0, 0, 3, 3.6865353266375283, 3.6986], 
1083        [6.422624312827989, 5388.556014109835, True, True, 0, 2, 1, 3.5909481979190283, 3.592005], 
1084        [6.534613244656113, 191.6070344509026, True, True, 1, 0, -1, 3.5294722429440584, 3.546166], 
1085        [6.558695213523644, 259.65938178131034, True, True, 0, 0, 3, 3.516526936765838, 3.495428], 
1086        [6.650921622278372, 93.26537659737657, True, True, 0, 0, 3, 3.4678179073694952, 3.466503], 
1087        [6.715273704410772, 289.3938681380316, True, False, 0, 0, 0, 3.4346235125812807, 0.0], 
1088        [6.85941304573619, 603.5495976464832, True, True, 0, 1, 3, 3.362534044860622, 3.32509], 
1089        [7.051162772888445, 126.43246447656593, True, True, 0, 1, 2, 3.2712038721790675, 3.352121], 
1090        [7.077700845503319, 125.49742760019636, True, False, 0, 0, 0, 3.2589538988480626, 0.0], 
1091        [7.099393757363675, 416.55444885434633, True, False, 0, 0, 0, 3.2490085228959193, 0.0], 
1092        [7.162393327864274, 369.27397110921817, True, False, 0, 0, 0, 3.2204673608202383, 0.0], 
1093        [7.412173495305892, 482.84120827021826, True, True, 0, 2, 2, 3.112085822159976, 3.133096]]
1094        ]
1095#
1096def test0():
1097    if NeedTestData: TestData()
1098    ibrav,cell,bestcell,Pwr,peaks = TestData
1099    print ('best cell:',bestcell)
1100    print ('old cell:',cell)
1101    Peaks = np.array(peaks)
1102    HKL = Peaks[4:7]
1103    print (calc_M20(peaks,HKL))
1104    A = G2lat.cell2A(cell)
1105    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
1106    print ('new cell:',G2lat.A2cell(A))   
1107    print ('x:',result[0])
1108    print ('cov_x:',result[1])
1109    print ('infodict:')
1110    for item in result[2]:
1111        print (item,result[2][item])
1112    print ('msg:',result[3])
1113    print ('ier:',result[4])
1114    result = refinePeaks(peaks,ibrav,A)
1115    N,M20,X20,A = result
1116    print ('refinePeaks:',N,M20,X20,G2lat.A2cell(A))
1117    print ('compare bestcell:',bestcell)
1118#
1119def test1():
1120    if NeedTestData: TestData()
1121    ibrav,A,Pwr,peaks = TestData2
1122    print ('bad cell:',G2lat.A2cell(A))
1123    print ('FitHKL')
1124    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
1125    result = refinePeaks(peaks,ibrav,A)
1126    N,M20,X20,A = result
1127    print ('refinePeaks:',N,M20,X20,A)
1128#    Peaks = np.array(peaks)
1129#    HKL = Peaks[4:7]
1130#    print calc_M20(peaks,HKL)
1131#    OK,smin,A,result = FitHKL(ibrav,peaks,A,Pwr)
1132#    print 'new cell:',G2lat.A2cell(A)   
1133#    print 'x:',result[0]
1134#    print 'cov_x:',result[1]
1135#    print 'infodict:'
1136#    for item in result[2]:
1137#        print item,result[2][item]
1138#    print 'msg:',result[3]
1139#    print 'ier:',result[4]
1140   
1141#
1142if __name__ == '__main__':
1143    test0()
1144    test1()
1145#    test2()
1146#    test3()
1147#    test4()
1148#    test5()
1149#    test6()
1150#    test7()
1151#    test8()
1152    print ("OK")
Note: See TracBrowser for help on using the repository browser.