source: trunk/GSASIIimage.py @ 1649

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

setting up image calibration to allow more control over refined parameters
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 40.2 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2015-02-10 19:24:18 +0000 (Tue, 10 Feb 2015) $
5# $Author: vondreele $
6# $Revision: 1644 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 1644 2015-02-10 19:24:18Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIIimage: Image calc module*
12================================
13
14Ellipse fitting & image integration
15
16'''
17
18import math
19import time
20import numpy as np
21import numpy.linalg as nl
22import numpy.ma as ma
23import polymask as pm
24from scipy.optimize import leastsq
25import copy
26import GSASIIpath
27GSASIIpath.SetVersionNumber("$Revision: 1644 $")
28import GSASIIplot as G2plt
29import GSASIIlattice as G2lat
30import GSASIIpwd as G2pwd
31import GSASIIspc as G2spc
32import fellipse as fel
33
34# trig functions in degrees
35sind = lambda x: math.sin(x*math.pi/180.)
36asind = lambda x: 180.*math.asin(x)/math.pi
37tand = lambda x: math.tan(x*math.pi/180.)
38atand = lambda x: 180.*math.atan(x)/math.pi
39atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
40cosd = lambda x: math.cos(x*math.pi/180.)
41acosd = lambda x: 180.*math.acos(x)/math.pi
42rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
43#numpy versions
44npsind = lambda x: np.sin(x*np.pi/180.)
45npasind = lambda x: 180.*np.arcsin(x)/np.pi
46npcosd = lambda x: np.cos(x*np.pi/180.)
47npacosd = lambda x: 180.*np.arccos(x)/np.pi
48nptand = lambda x: np.tan(x*np.pi/180.)
49npatand = lambda x: 180.*np.arctan(x)/np.pi
50npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
51   
52def pointInPolygon(pXY,xy):
53    'Needs a doc string'
54    #pXY - assumed closed 1st & last points are duplicates
55    Inside = False
56    N = len(pXY)
57    p1x,p1y = pXY[0]
58    for i in range(N+1):
59        p2x,p2y = pXY[i%N]
60        if (max(p1y,p2y) >= xy[1] > min(p1y,p2y)) and (xy[0] <= max(p1x,p2x)):
61            if p1y != p2y:
62                xinters = (xy[1]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
63            if p1x == p2x or xy[0] <= xinters:
64                Inside = not Inside
65        p1x,p1y = p2x,p2y
66    return Inside
67   
68def peneCorr(tth,dep,tilt=0.,azm=0.):
69    'Needs a doc string'
70#    return dep*(1.-npcosd(abs(tilt*npsind(azm))-tth*npcosd(azm)))  #something wrong here
71    return dep*(1.-npcosd(tth))         #best one
72#    return dep*npsind(tth)             #not as good as 1-cos2Q
73       
74def makeMat(Angle,Axis):
75    '''Make rotation matrix from Angle and Axis
76
77    :param float Angle: in degrees
78    :param int Axis: 0 for rotation about x, 1 for about y, etc.
79    '''
80    cs = npcosd(Angle)
81    ss = npsind(Angle)
82    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
83    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
84                   
85def FitEllipse(xy):
86   
87    def ellipse_center(p):
88        ''' gives ellipse center coordinates
89        '''
90        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
91        num = b*b-a*c
92        x0=(c*d-b*f)/num
93        y0=(a*f-b*d)/num
94        return np.array([x0,y0])
95   
96    def ellipse_angle_of_rotation( p ):
97        ''' gives rotation of ellipse major axis from x-axis
98        range will be -90 to 90 deg
99        '''
100        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
101        return 0.5*npatand(2*b/(a-c))
102   
103    def ellipse_axis_length( p ):
104        ''' gives ellipse radii in [minor,major] order
105        '''
106        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
107        up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
108        down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
109        down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
110        res1=np.sqrt(up/down1)
111        res2=np.sqrt(up/down2)
112        return np.array([ res2,res1])
113   
114    xy = np.array(xy)
115    x = np.asarray(xy.T[0])[:,np.newaxis]
116    y = np.asarray(xy.T[1])[:,np.newaxis]
117    D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
118    S = np.dot(D.T,D)
119    C = np.zeros([6,6])
120    C[0,2] = C[2,0] = 2; C[1,1] = -1
121    E, V =  nl.eig(np.dot(nl.inv(S), C))
122    n = np.argmax(np.abs(E))
123    a = V[:,n]
124    cent = ellipse_center(a)
125    phi = ellipse_angle_of_rotation(a)
126    radii = ellipse_axis_length(a)
127    phi += 90.
128    if radii[0] > radii[1]:
129        radii = [radii[1],radii[0]]
130        phi -= 90.
131    return cent,phi,radii
132
133def FitDetector(rings,varyList,parmDict,Print=True):
134    'Needs a doc string'
135       
136    def CalibPrint(ValSig,chisq,Npts):
137        print 'Image Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts)
138        ptlbls = 'names :'
139        ptstr =  'values:'
140        sigstr = 'esds  :'
141        for name,fmt,value,sig in ValSig:
142            ptlbls += "%s" % (name.rjust(12))
143            if name == 'phi':
144                ptstr += fmt % (value%360.)
145            else:
146                ptstr += fmt % (value)
147            if sig:
148                sigstr += fmt % (sig)
149            else:
150                sigstr += 12*' '
151        print ptlbls
152        print ptstr
153        print sigstr       
154       
155    def ellipseCalcD(B,xyd,varyList,parmDict):
156       
157        x,y,dsp = xyd
158        varyDict = dict(zip(varyList,B))
159        parms = {}
160        for parm in parmDict:
161            if parm in varyList:
162                parms[parm] = varyDict[parm]
163            else:
164                parms[parm] = parmDict[parm]
165        phi = parms['phi']-90.               #get rotation of major axis from tilt axis
166        tth = 2.0*npasind(parms['wave']/(2.*dsp))
167        phi0 = npatan2d(y-parms['det-Y'],x-parms['det-X'])
168        dxy = peneCorr(tth,parms['dep'],parms['tilt'],phi0)
169        ttth = nptand(tth)
170        stth = npsind(tth)
171        cosb = npcosd(parms['tilt'])
172        tanb = nptand(parms['tilt'])       
173        tbm = nptand((tth-parms['tilt'])/2.)
174        tbp = nptand((tth+parms['tilt'])/2.)
175        sinb = npsind(parms['tilt'])
176        d = parms['dist']+dxy
177        fplus = d*tanb*stth/(cosb+stth)
178        fminus = d*tanb*stth/(cosb-stth)
179        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
180        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
181        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
182        R1 = (vplus+vminus)/2.                                    #major axis
183        zdis = (fplus-fminus)/2.
184        Robs = np.sqrt((x-parms['det-X'])**2+(y-parms['det-Y'])**2)
185        rsqplus = R0**2+R1**2
186        rsqminus = R0**2-R1**2
187        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
188        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
189        P = 2.*R0**2*zdis*npcosd(phi0-phi)
190        Rcalc = (P+Q)/R
191        M = (Robs-Rcalc)*10.        #why 10? does make "chi**2" more reasonable
192        return M
193       
194    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
195    fmt = ['%12.3f','%12.3f','%12.3f','%12.3f','%12.3f','%12.3f','%12.5f']
196    p0 = [parmDict[key] for key in varyList]
197    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
198    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0))   #reduced chi^2 = M/(Nobs-Nvar)
199    parmDict.update(zip(varyList,result[0]))
200    vals = list(result[0])
201    sig = list(np.sqrt(chisq*np.diag(result[1])))
202    sigList = np.zeros(7)
203    for i,name in enumerate(varyList):
204        sigList[i] = sig[varyList.index(name)]
205    ValSig = zip(varyList,fmt,vals,sig)
206    if Print:
207        CalibPrint(ValSig,chisq,rings.shape[0])
208    return chisq
209                   
210def ImageLocalMax(image,w,Xpix,Ypix):
211    'Needs a doc string'
212    w2 = w*2
213    sizey,sizex = image.shape
214    xpix = int(Xpix)            #get reference corner of pixel chosen
215    ypix = int(Ypix)
216    if (w < xpix < sizex-w) and (w < ypix < sizey-w) and image[ypix,xpix]:
217        Z = image[ypix-w:ypix+w,xpix-w:xpix+w]
218        Zmax = np.argmax(Z)
219        Zmin = np.argmin(Z)
220        xpix += Zmax%w2-w
221        ypix += Zmax/w2-w
222        return xpix,ypix,np.ravel(Z)[Zmax],max(0.0001,np.ravel(Z)[Zmin])   #avoid neg/zero minimum
223    else:
224        return 0,0,0,0     
225   
226def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
227    'Needs a doc string'
228    def ellipseC():
229        'compute estimate of ellipse circumference'
230        if radii[0] < 0:        #hyperbola
231            theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
232            print theta
233            return 0
234        apb = radii[1]+radii[0]
235        amb = radii[1]-radii[0]
236        return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
237       
238    cent,phi,radii = ellipse
239    cphi = cosd(phi-90.)        #convert to major axis rotation
240    sphi = sind(phi-90.)
241    ring = []
242    C = int(ellipseC())         #ring circumference
243    for i in range(0,C,1):      #step around ring in 1mm increments
244        a = 360.*i/C
245        x = radii[1]*cosd(a)        #major axis
246        y = radii[0]*sind(a)
247        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
248        Y = (sphi*x+cphi*y+cent[1])*scaley
249        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
250        if I and J and I/J > reject:
251            X += .5                             #set to center of pixel
252            Y += .5
253            X /= scalex                         #convert back to mm
254            Y /= scaley
255            if [X,Y,dsp] not in ring:           #no duplicates!
256                ring.append([X,Y,dsp])
257    if len(ring) < 10:
258        ring = []
259    return ring
260   
261def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
262    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
263    on output
264    radii[0] (b-minor axis) set < 0. for hyperbola
265   
266    '''
267    radii = [0,0]
268    ttth = tand(tth)
269    stth = sind(tth)
270    ctth = cosd(tth)
271    cosb = cosd(tilt)
272    tanb = tand(tilt)
273    tbm = tand((tth-tilt)/2.)
274    tbp = tand((tth+tilt)/2.)
275    sinb = sind(tilt)
276    d = dist+dxy
277    if tth+abs(tilt) < 90.:      #ellipse
278        fplus = d*tanb*stth/(cosb+stth)
279        fminus = d*tanb*stth/(cosb-stth)
280        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
281        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
282        radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
283        radii[1] = (vplus+vminus)/2.                                    #major axis
284        zdis = (fplus-fminus)/2.
285    else:   #hyperbola!
286        f = d*abs(tanb)*stth/(cosb+stth)
287        v = d*(abs(tanb)+tand(tth-abs(tilt)))
288        delt = d*stth*(1.+stth*cosb)/(abs(sinb)*cosb*(stth+cosb))
289        eps = (v-f)/(delt-v)
290        radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.)                     #-minor axis
291        radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
292        if tilt > 0:
293            zdis = f+radii[1]*eps
294        else:
295            zdis = -f
296#NB: zdis is || to major axis & phi is rotation of minor axis
297#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
298    elcent = [cent[0]+zdis*sind(phi),cent[1]-zdis*cosd(phi)]
299    return elcent,phi,radii
300   
301def GetEllipse(dsp,data):
302    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
303    as given in image controls dictionary (data) and a d-spacing (dsp)
304    '''
305    cent = data['center']
306    tilt = data['tilt']
307    phi = data['rotation']
308    dep = data['DetDepth']
309    tth = 2.0*asind(data['wavelength']/(2.*dsp))
310    dxy = peneCorr(tth,dep,tilt)
311    dist = data['distance']
312    return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
313       
314def GetDetectorXY(dsp,azm,data):
315    'Needs a doc string'
316   
317    elcent,phi,radii = GetEllipse(dsp,data)
318    phi = data['rotation']-90.          #to give rotation of major axis
319    tilt = data['tilt']
320    dist = data['distance']
321    cent = data['center']
322    tth = 2.0*asind(data['wavelength']/(2.*dsp))
323    ttth = tand(tth)
324    stth = sind(tth)
325    ctth = cosd(tth)
326    cosb = cosd(tilt)
327    if radii[0] > 0.:
328        sinb = sind(tilt)
329        tanb = tand(tilt)
330        fplus = dist*tanb*stth/(cosb+stth)
331        fminus = dist*tanb*stth/(cosb-stth)
332        zdis = (fplus-fminus)/2.
333        rsqplus = radii[0]**2+radii[1]**2
334        rsqminus = radii[0]**2-radii[1]**2
335        R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
336        Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
337        P = 2.*radii[0]**2*zdis*cosd(azm-phi)
338        radius = (P+Q)/R
339        xy = np.array([radius*cosd(azm),radius*sind(azm)])
340        xy += cent
341    else:   #hyperbola - both branches (one is way off screen!)
342        sinb = abs(sind(tilt))
343        tanb = abs(tand(tilt))
344        f = dist*tanb*stth/(cosb+stth)
345        v = dist*(tanb+tand(tth-abs(tilt)))
346        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
347        ecc = (v-f)/(delt-v)
348        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
349        if tilt > 0.:
350            offset = 2.*radii[1]*ecc+f      #select other branch
351            xy = [-R*cosd(azm)-offset,R*sind(azm)]
352        else:
353            offset = -f
354            xy = [-R*cosd(azm)-offset,-R*sind(azm)]
355        xy = -np.array([xy[0]*cosd(phi)+xy[1]*sind(phi),xy[0]*sind(phi)-xy[1]*cosd(phi)])
356        xy += cent
357    return xy
358   
359def GetDetXYfromThAzm(Th,Azm,data):
360    'Needs a doc string'
361    dsp = data['wavelength']/(2.0*npsind(Th))   
362    return GetDetectorXY(dsp,azm,data)
363                   
364def GetTthAzmDsp(x,y,data): #expensive
365    'Needs a doc string - checked OK for ellipses & hyperbola'
366    wave = data['wavelength']
367    cent = data['center']
368    tilt = data['tilt']
369    dist = data['distance']/cosd(tilt)
370    x0 = data['distance']*tand(tilt)
371    phi = data['rotation']
372    dep = data['DetDepth']
373    LRazim = data['LRazimuth']
374    azmthoff = data['azmthOff']
375    dx = np.array(x-cent[0],dtype=np.float32)
376    dy = np.array(y-cent[1],dtype=np.float32)
377    D = ((dx-x0)**2+dy**2+data['distance']**2)      #sample to pixel distance
378    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
379    X = np.dot(X,makeMat(phi,2))
380    Z = np.dot(X,makeMat(tilt,0)).T[2]
381    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
382    dxy = peneCorr(tth,dep,tilt,npatan2d(dy,dx))
383    DX = dist-Z+dxy
384    DY = np.sqrt(dx**2+dy**2-Z**2)
385    tth = npatan2d(DY,DX) 
386    dsp = wave/(2.*npsind(tth/2.))
387    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
388    G = D/data['distance']**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
389    return np.array([tth,azm,G,dsp])
390   
391def GetTth(x,y,data):
392    'Give 2-theta value for detector x,y position; calibration info in data'
393    return GetTthAzmDsp(x,y,data)[0]
394   
395def GetTthAzm(x,y,data):
396    'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
397    return GetTthAzmDsp(x,y,data)[0:2]
398   
399def GetTthAzmG(x,y,data):
400    '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
401     calibration info in data - only used in integration
402    '''
403    'Needs a doc string - checked OK for ellipses & hyperbola'
404    tilt = data['tilt']
405    dist = data['distance']/npcosd(tilt)
406    x0 = data['distance']*nptand(tilt)
407    MN = -np.inner(makeMat(data['rotation'],2),makeMat(tilt,0))
408    distsq = data['distance']**2
409    dx = x-data['center'][0]
410    dy = y-data['center'][1]
411    G = ((dx-x0)**2+dy**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
412    X = np.dstack([dx.T,dy.T,np.zeros_like(dx.T)])
413    Z = np.dot(X,MN).T[2]
414    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
415    dxy = peneCorr(tth,data['DetDepth'],tilt,npatan2d(dy,dx))
416    tth = npatan2d(np.sqrt(dx**2+dy**2-Z**2),dist-Z+dxy) 
417    azm = (npatan2d(dy,dx)+data['azmthOff']+720.)%360.
418    return tth,azm,G
419
420def GetDsp(x,y,data):
421    'Give d-spacing value for detector x,y position; calibration info in data'
422    return GetTthAzmDsp(x,y,data)[3]
423       
424def GetAzm(x,y,data):
425    'Give azimuth value for detector x,y position; calibration info in data'
426    return GetTthAzmDsp(x,y,data)[1]
427       
428def ImageCompress(image,scale):
429    'Needs a doc string'
430    if scale == 1:
431        return image
432    else:
433        return image[::scale,::scale]
434       
435def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
436    'Needs a doc string'
437    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
438    curr = np.array([dist,x,y])
439    return abs(avg-curr)/avg < .02
440
441def EdgeFinder(image,data):
442    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
443    Not currently used but might be useful in future?
444    '''
445    import numpy.ma as ma
446    Nx,Ny = data['size']
447    pixelSize = data['pixelSize']
448    edgemin = data['edgemin']
449    scalex = pixelSize[0]/1000.
450    scaley = pixelSize[1]/1000.   
451    tay,tax = np.mgrid[0:Nx,0:Ny]
452    tax = np.asfarray(tax*scalex,dtype=np.float32)
453    tay = np.asfarray(tay*scaley,dtype=np.float32)
454    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
455    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
456    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
457    return zip(tax,tay)
458   
459def MakeFrameMask(data,frame):
460    pixelSize = data['pixelSize']
461    scalex = pixelSize[0]/1000.
462    scaley = pixelSize[1]/1000.
463    blkSize = 512
464    Nx,Ny = data['size']
465    nXBlks = (Nx-1)/blkSize+1
466    nYBlks = (Ny-1)/blkSize+1
467    tam = ma.make_mask_none(data['size'])
468    for iBlk in range(nXBlks):
469        iBeg = iBlk*blkSize
470        iFin = min(iBeg+blkSize,Nx)
471        for jBlk in range(nYBlks):
472            jBeg = jBlk*blkSize
473            jFin = min(jBeg+blkSize,Ny)               
474            nI = iFin-iBeg
475            nJ = jFin-jBeg
476            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
477            tax = np.asfarray(tax*scalex,dtype=np.float32)
478            tay = np.asfarray(tay*scaley,dtype=np.float32)
479            tamp = ma.make_mask_none((1024*1024))
480            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
481                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
482            if tamp.shape:
483                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
484                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
485            else:
486                tam[iBeg:iFin,jBeg:jFin] = True
487    return tam.T
488   
489def ImageRecalibrate(self,data,masks):
490    'Needs a doc string'
491    import ImageCalibrants as calFile
492    print 'Image recalibration:'
493    time0 = time.time()
494    pixelSize = data['pixelSize']
495    scalex = 1000./pixelSize[0]
496    scaley = 1000./pixelSize[1]
497    pixLimit = data['pixLimit']
498    cutoff = data['cutoff']
499    data['rings'] = []
500    data['ellipses'] = []
501    if not data['calibrant']:
502        print 'no calibration material selected'
503        return True   
504    skip = data['calibskip']
505    dmin = data['calibdmin']
506    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
507    HKL = []
508    for bravais,sg,cell in zip(Bravais,SGs,Cells):
509        A = G2lat.cell2A(cell)
510        if sg:
511            SGData = G2spc.SpcGroup(sg)[1]
512            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
513            HKL += hkl
514        else:
515            hkl = G2lat.GenHBravais(dmin,bravais,A)
516            HKL += hkl
517    HKL = G2lat.sortHKLd(HKL,True,False)
518    varyList = [item for item in data['varyList'] if data['varyList'][item]]
519    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
520        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
521    Found = False
522    wave = data['wavelength']
523    frame = masks['Frames']
524    tam = ma.make_mask_none(self.ImageZ.shape)
525    if frame:
526        tam = ma.mask_or(tam,MakeFrameMask(data,frame))
527    for iH,H in enumerate(HKL):
528        print H
529        dsp = H[3]
530        tth = 2.0*asind(wave/(2.*dsp))
531        if tth+abs(data['tilt']) > 90.:
532            print 'next line is a hyperbola - search stopped'
533            break
534        ellipse = GetEllipse(dsp,data)
535        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(self.ImageZ,mask=tam))
536        if Ring:
537            if iH >= skip:
538                data['rings'].append(np.array(Ring))
539            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
540            Found = True
541        elif not Found:         #skipping inner rings, keep looking until ring found
542            continue
543        else:                   #no more rings beyond edge of detector
544            data['ellipses'].append([])
545            continue
546#            break
547    rings = np.concatenate((data['rings']),axis=0)
548    chisq = FitDetector(rings,varyList,parmDict)
549    data['distance'] = parmDict['dist']
550    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
551    data['rotation'] = np.mod(parmDict['phi'],360.0)
552    data['tilt'] = parmDict['tilt']
553    data['DetDepth'] = parmDict['dep']
554    data['chisq'] = chisq
555    N = len(data['ellipses'])
556    data['ellipses'] = []           #clear away individual ellipse fits
557    for H in HKL[:N]:
558        ellipse = GetEllipse(H[3],data)
559        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
560   
561    print 'calibration time = ',time.time()-time0
562    G2plt.PlotImage(self,newImage=True)       
563    return True
564           
565def ImageCalibrate(self,data):
566    'Needs a doc string'
567    import copy
568    import ImageCalibrants as calFile
569    print 'Image calibration:'
570    time0 = time.time()
571    ring = data['ring']
572    pixelSize = data['pixelSize']
573    scalex = 1000./pixelSize[0]
574    scaley = 1000./pixelSize[1]
575    pixLimit = data['pixLimit']
576    cutoff = data['cutoff']
577    if len(ring) < 5:
578        print 'not enough inner ring points for ellipse'
579        return False
580       
581    #fit start points on inner ring
582    data['ellipses'] = []
583    data['rings'] = []
584    outE = FitEllipse(ring)
585    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
586    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi**2: %.3f, Np: %d'
587    if outE:
588        print fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1])
589        ellipse = outE
590    else:
591        return False
592       
593    #setup 360 points on that ring for "good" fit
594    data['ellipses'].append(ellipse[:]+('g',))
595    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
596    if Ring:
597        ellipse = FitEllipse(Ring)
598        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
599        ellipse = FitEllipse(Ring)
600    else:
601        print '1st ring not sufficiently complete to proceed'
602        return False
603    print fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],ellipse[2][0],ellipse[2][1],0.,len(Ring))     #cent,phi,radii
604    data['ellipses'].append(ellipse[:]+('r',))
605    data['rings'].append(np.array(Ring))
606    G2plt.PlotImage(self,newImage=True)
607   
608#setup for calibration
609    data['rings'] = []
610    if not data['calibrant']:
611        print 'no calibration material selected'
612        return True
613   
614    skip = data['calibskip']
615    dmin = data['calibdmin']
616#generate reflection set
617    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
618    HKL = []
619    for bravais,sg,cell in zip(Bravais,SGs,Cells):
620        A = G2lat.cell2A(cell)
621        if sg:
622            SGData = G2spc.SpcGroup(sg)[1]
623            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
624            HKL += hkl
625        else:
626            hkl = G2lat.GenHBravais(dmin,bravais,A)
627            HKL += hkl
628    HKL = G2lat.sortHKLd(HKL,True,False)[skip:]
629    wave = data['wavelength']
630#set up 1st ring
631    elcent,phi,radii = ellipse              #from fit of 1st ring
632    dsp = HKL[0][3]
633    print '1st ring: try %.4f'%(dsp)
634    tth = 2.0*asind(wave/(2.*dsp))
635    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,self.ImageZ)
636    ttth = nptand(tth)
637    stth = npsind(tth)
638    ctth = npcosd(tth)
639#1st estimate of tilt; assume ellipse - don't know sign though
640    tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
641    if not tilt:
642        print 'WARNING - selected ring was fitted as a circle'
643        print ' - if detector was tilted we suggest you skip this ring - WARNING'
644#1st estimate of dist: sample to detector normal to plane
645    data['distance'] = dist = radii[0]**2/(ttth*radii[1])
646#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
647    zdisp = radii[1]*ttth*tand(tilt)
648    zdism = radii[1]*ttth*tand(-tilt)
649#cone axis position; 2 choices. Which is right?     
650#NB: zdisp is || to major axis & phi is rotation of minor axis
651#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
652    centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
653    centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
654#check get same ellipse parms either way
655#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
656    fail = True
657    i2 = 1
658    while fail:
659        dsp = HKL[i2][3]
660        print '2nd ring: try %.4f'%(dsp)
661        tth = 2.0*asind(wave/(2.*dsp))
662        ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
663        print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1])
664        Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,self.ImageZ)
665        parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
666            'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
667        varyList = [item for item in data['varyList'] if data['varyList'][item]]
668        if len(Ringp) > 10:
669            chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)
670            tiltp = parmDict['tilt']
671            phip = parmDict['phi']
672            centp = [parmDict['det-X'],parmDict['det-Y']]
673            fail = False
674        else:
675            chip = 1e6
676        ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
677        print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1])
678        Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,self.ImageZ)
679        if len(Ringm) > 10:
680            parmDict['tilt'] *= -1
681            chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)
682            tiltm = parmDict['tilt']
683            phim = parmDict['phi']
684            centm = [parmDict['det-X'],parmDict['det-Y']]
685            fail = False
686        else:
687            chim = 1e6
688        if fail:
689            i2 += 1
690    if chip < chim:
691        data['tilt'] = tiltp
692        data['center'] = centp
693        data['rotation'] = phip
694    else:
695        data['tilt'] = tiltm
696        data['center'] = centm
697        data['rotation'] = phim
698    data['ellipses'].append(ellipsep[:]+('b',))
699    data['rings'].append(np.array(Ringp))
700    data['ellipses'].append(ellipsem[:]+('r',))
701    data['rings'].append(np.array(Ringm))
702    G2plt.PlotImage(self,newImage=True)
703    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
704        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
705    varyList = [item for item in data['varyList'] if data['varyList'][item]]
706    data['rings'] = []
707    data['ellipses'] = []
708    for i,H in enumerate(HKL):
709        dsp = H[3]
710        tth = 2.0*asind(wave/(2.*dsp))
711        if tth+abs(data['tilt']) > 90.:
712            print 'next line is a hyperbola - search stopped'
713            break
714        print 'HKLD:',H[:4],'2-theta: %.4f'%(tth)
715        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
716        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
717        print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
718        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
719        if Ring:
720            data['rings'].append(np.array(Ring))
721            rings = np.concatenate((data['rings']),axis=0)
722            if i:
723                chisq = FitDetector(rings,varyList,parmDict,False)
724                data['distance'] = parmDict['dist']
725                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
726                data['rotation'] = parmDict['phi']
727                data['tilt'] = parmDict['tilt']
728                data['DetDepth'] = parmDict['dep']
729                data['chisq'] = chisq
730                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
731                print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
732            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
733#            G2plt.PlotImage(self,newImage=True)
734        else:
735            print 'insufficient number of points in this ellipse to fit'
736#            break
737    G2plt.PlotImage(self,newImage=True)
738    fullSize = len(self.ImageZ)/scalex
739    if 2*radii[1] < .9*fullSize:
740        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
741    N = len(data['ellipses'])
742    if N > 2:
743        FitDetector(rings,varyList,parmDict)
744        data['distance'] = parmDict['dist']
745        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
746        data['rotation'] = parmDict['phi']
747        data['tilt'] = parmDict['tilt']
748        data['DetDepth'] = parmDict['dep']
749    for H in HKL[:N]:
750        ellipse = GetEllipse(H[3],data)
751        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
752    print 'calibration time = ',time.time()-time0
753    G2plt.PlotImage(self,newImage=True)       
754    return True
755   
756def Make2ThetaAzimuthMap(data,masks,iLim,jLim,times): #most expensive part of integration!
757    'Needs a doc string'
758    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
759    pixelSize = data['pixelSize']
760    scalex = pixelSize[0]/1000.
761    scaley = pixelSize[1]/1000.
762   
763    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
764    tax = np.asfarray(tax*scalex,dtype=np.float32)
765    tay = np.asfarray(tay*scaley,dtype=np.float32)
766    nI = iLim[1]-iLim[0]
767    nJ = jLim[1]-jLim[0]
768    t0 = time.time()
769    #make position masks here
770    frame = masks['Frames']
771    tam = ma.make_mask_none((nI,nJ))
772    if frame:
773        tamp = ma.make_mask_none((1024*1024))
774        tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
775            tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
776        tam = ma.mask_or(tam.flatten(),tamp)
777    polygons = masks['Polygons']
778    for polygon in polygons:
779        if polygon:
780            tamp = ma.make_mask_none((1024*1024))
781            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
782                tay.flatten(),len(polygon),polygon,tamp)[:nI*nJ])
783            tam = ma.mask_or(tam.flatten(),tamp)
784    if tam.shape: tam = np.reshape(tam,(nI,nJ))
785    spots = masks['Points']
786    for X,Y,diam in spots:
787        tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
788        tam = ma.mask_or(tam,tamp)
789    times[0] += time.time()-t0
790    t0 = time.time()
791    TA = np.array(GetTthAzmG(tax,tay,data))     #includes geom. corr. as dist**2/d0**2 - most expensive step
792    times[1] += time.time()-t0
793    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
794    return np.array(TA),tam           #2-theta, azimuth & geom. corr. arrays & position mask
795
796def Fill2ThetaAzimuthMap(masks,TA,tam,image):
797    'Needs a doc string'
798    Zlim = masks['Thresholds'][1]
799    rings = masks['Rings']
800    arcs = masks['Arcs']
801    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
802    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist**2/d0**2
803    for tth,thick in rings:
804        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
805    for tth,azm,thick in arcs:
806        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
807        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
808        tam = ma.mask_or(tam.flatten(),tamt*tama)
809    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
810    tabs = np.ones_like(taz)
811    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
812    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
813    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
814    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
815    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
816    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
817    return tax,tay,taz,tad,tabs
818   
819def ImageIntegrate(image,data,masks,blkSize=128,dlg=None,returnN=False):
820    'Needs a doc string'    #for q, log(q) bins need data['binType']
821    import histogram2d as h2d
822    print 'Begin image integration'
823    LUtth = np.array(data['IOtth'])
824    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
825    numAzms = data['outAzimuths']
826    numChans = data['outChannels']
827    azmOff = data['azmthOff']
828    Dazm = (LRazm[1]-LRazm[0])/numAzms
829    if 'log(q)' in data['binType']:
830        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
831    elif 'q' == data['binType']:
832        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
833    elif '2-theta' in data['binType']:
834        lutth = LUtth               
835    dtth = (lutth[1]-lutth[0])/numChans
836    muT = data['SampleAbs'][0]
837    if 'SASD' in data['type']:
838        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
839    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
840    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
841    imageN = len(image)
842    Nx,Ny = data['size']
843    nXBlks = (Nx-1)/blkSize+1
844    nYBlks = (Ny-1)/blkSize+1
845    Nup = nXBlks*nYBlks*3+3
846    tbeg = time.time()
847    Nup = 0
848    if dlg:
849        dlg.Update(Nup)
850    times = [0,0,0,0,0]
851    for iBlk in range(nYBlks):
852        iBeg = iBlk*blkSize
853        iFin = min(iBeg+blkSize,Ny)
854        for jBlk in range(nXBlks):
855            jBeg = jBlk*blkSize
856            jFin = min(jBeg+blkSize,Nx)
857            # next is most expensive step!
858            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin),times)           #2-theta & azimuth arrays & create position mask
859            Nup += 1
860            if dlg:
861                dlg.Update(Nup)
862            Block = image[iBeg:iFin,jBeg:jFin]
863            t0 = time.time()
864            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
865            times[2] += time.time()-t0
866            Nup += 1
867            if dlg:
868                dlg.Update(Nup)
869            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
870            tax = np.where(tax < LRazm[0],tax+360.,tax)
871            if data['SampleAbs'][1]:
872                if 'PWDR' in data['type']:
873                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))
874                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
875                elif 'SASD' in data['type']:    #assumes flat plate sample normal to beam
876                    tabs = G2pwd.Absorb('Fixed',muT,tay)
877            if 'log(q)' in data['binType']:
878                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
879            elif 'q' == data['binType']:
880                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
881            t0 = time.time()
882            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
883                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz*tad/tabs,
884                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
885            times[3] += time.time()-t0
886            Nup += 1
887            if dlg:
888                dlg.Update(Nup)
889    t0 = time.time()
890    NST = np.array(NST,dtype=np.float)
891    H0 = np.divide(H0,NST)
892    H0 = np.nan_to_num(H0)
893    H2 = np.array([tth for tth in np.linspace(lutth[0],lutth[1],numChans+1)])
894    if 'log(q)' in data['binType']:
895        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
896    elif 'q' == data['binType']:
897        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
898    if Dazm:       
899        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
900    else:
901        H1 = LRazm
902    H0 /= npcosd(H2[:-1])           #**2? I don't think so, **1 is right for powders
903    if 'SASD' in data['type']:
904        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
905    if data['Oblique'][1]:
906        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
907    if 'SASD' in data['type'] and data['PolaVal'][1]:
908        #NB: in G2pwd.Polarization azm is defined from plane of polarization, not image x axis!
909        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm-90.)[0] for azm in (H1[:-1]+np.diff(H1)/2.)])
910    Nup += 1
911    if dlg:
912        dlg.Update(Nup)
913    times[4] += time.time()-t0
914    print 'Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
915        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4])
916    print "Elapsed time:","%8.3fs"%(time.time()-tbeg)
917    print 'Integration complete'
918    if returnN:     #As requested by Steven Weigand
919        return H0,H1,H2,NST
920    else:
921        return H0,H1,H2
922   
923def MakeStrStaRing(ring,Image,Controls):
924    ellipse = GetEllipse(ring['Dset'],Controls)
925    pixSize = Controls['pixelSize']
926    scalex = 1000./pixSize[0]
927    scaley = 1000./pixSize[1]
928    Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)).T   #returns x,y,dsp for each point in ring
929    if len(Ring):
930        ring['ImxyObs'] = copy.copy(Ring[:2])
931        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
932        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
933        ring['ImtaObs'] = TA
934        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
935        Ring[0] = TA[0]
936        Ring[1] = TA[1]
937        return Ring,ring
938    else:
939        ring['ImxyObs'] = [[],[]]
940        ring['ImtaObs'] = [[],[]]
941        ring['ImtaCalc'] = [[],[]]
942        return [],[]    #bad ring; no points found
943   
944def FitStrSta(Image,StrSta,Controls):
945    'Needs a doc string'
946   
947    StaControls = copy.deepcopy(Controls)
948    phi = StrSta['Sample phi']
949    wave = Controls['wavelength']
950    StaType = StrSta['Type']
951    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
952
953    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
954        dset = ring['Dset']
955        Ring,R = MakeStrStaRing(ring,Image,StaControls)
956        if len(Ring):
957            ring.update(R)
958            p0 = ring['Emat']
959            val,esd = FitStrain(Ring,p0,dset,wave,phi,StaType)
960            ring['Emat'] = val
961            ring['Esig'] = esd
962    CalcStrSta(StrSta,Controls)
963   
964def CalcStrSta(StrSta,Controls):
965
966    wave = Controls['wavelength']
967    phi = StrSta['Sample phi']
968    StaType = StrSta['Type']
969    for ring in StrSta['d-zero']:
970        Eij = ring['Emat']
971        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
972        th,azm = ring['ImtaObs']
973        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
974        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
975        if StaType == 'True':
976            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
977        else:
978            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
979        ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
980
981def calcFij(omg,phi,azm,th):
982    '''Does something...
983
984    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
985
986    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
987        90 when perp. to sample surface
988    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
989    :param azm: his chi = azimuth around incident beam
990    :param th:  his theta = theta
991    '''
992    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
993    b = -npcosd(azm)*npcosd(th)
994    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
995    d = a*npsind(phi)+b*npcosd(phi)
996    e = a*npcosd(phi)-b*npsind(phi)
997    Fij = np.array([
998        [d**2,d*e,c*d],
999        [d*e,e**2,c*e],
1000        [c*d,c*e,c**2]])
1001    return -Fij
1002
1003def FitStrain(rings,p0,dset,wave,phi,StaType):
1004    'Needs a doc string'
1005    def StrainPrint(ValSig,dset):
1006        print 'Strain tensor for Dset: %.6f'%(dset)
1007        ptlbls = 'names :'
1008        ptstr =  'values:'
1009        sigstr = 'esds  :'
1010        for name,fmt,value,sig in ValSig:
1011            ptlbls += "%s" % (name.rjust(12))
1012            ptstr += fmt % (value)
1013            if sig:
1014                sigstr += fmt % (sig)
1015            else:
1016                sigstr += 12*' '
1017        print ptlbls
1018        print ptstr
1019        print sigstr
1020       
1021    def strainCalc(p,xyd,dset,wave,phi,StaType):
1022        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1023        dspo,azm,dsp = xyd
1024        th = npasind(wave/(2.0*dspo))
1025        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1026        if StaType == 'True':
1027            dspc = dset*np.exp(V)
1028        else:
1029            dspc = dset*(V+1.)
1030        return dspo-dspc
1031       
1032    names = ['e11','e12','e22']
1033    fmt = ['%12.2f','%12.2f','%12.2f']
1034    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1035    vals = list(result[0])
1036    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1037    sig = list(np.sqrt(chisq*np.diag(result[1])))
1038    ValSig = zip(names,fmt,vals,sig)
1039    StrainPrint(ValSig,dset)
1040    return vals,sig
1041   
1042       
Note: See TracBrowser for help on using the repository browser.