source: trunk/GSASIIindex.py @ 2038

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

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