source: trunk/GSASIIindex.py @ 2052

Last change on this file since 2052 was 2052, checked in by vondreele, 10 years ago

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