source: trunk/GSASIIimage.py @ 1405

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

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