source: trunk/GSASIIindex.py @ 1637

Last change on this file since 1637 was 1637, checked in by vondreele, 9 years ago

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