source: trunk/GSASIIindex.py @ 1592

Last change on this file since 1592 was 1592, checked in by vondreele, 8 years ago

finish modulation vector determination - use so.brute - slow but works
add peak indexing option to not scale M20 by 1+X20 - might be useful

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