source: trunk/GSASIIimage.py @ 1665

Last change on this file since 1665 was 1665, checked in by vondreele, 7 years ago

fix image calibrate so dist can be fixed & wave unknown

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