source: trunk/GSASIIimage.py @ 1653

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

allow all image controls to be varied in calibration and recalibration
stop debug printing in image calibration/recalibration

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 40.9 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2015-02-13 15:25:22 +0000 (Fri, 13 Feb 2015) $
5# $Author: vondreele $
6# $Revision: 1653 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 1653 2015-02-13 15:25:22Z 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: 1653 $")
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    wave = data['wavelength']
639#set up 1st ring
640    elcent,phi,radii = ellipse              #from fit of 1st ring
641    dsp = HKL[0][3]
642    print '1st ring: try %.4f'%(dsp)
643    tth = 2.0*asind(wave/(2.*dsp))
644    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,self.ImageZ)
645    ttth = nptand(tth)
646    stth = npsind(tth)
647    ctth = npcosd(tth)
648#1st estimate of tilt; assume ellipse - don't know sign though
649    if varyDict['tilt']:
650        tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
651        if not tilt:
652            print 'WARNING - selected ring was fitted as a circle'
653            print ' - if detector was tilted we suggest you skip this ring - WARNING'
654    else:
655        tilt = data['tilt']
656#1st estimate of dist: sample to detector normal to plane
657    if varyDict['dist']:
658        data['distance'] = dist = radii[0]**2/(ttth*radii[1])
659    else:
660        dist = data['distance']
661    if varyDict['tilt']:
662#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
663        zdisp = radii[1]*ttth*tand(tilt)
664        zdism = radii[1]*ttth*tand(-tilt)
665#cone axis position; 2 choices. Which is right?     
666#NB: zdisp is || to major axis & phi is rotation of minor axis
667#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
668        centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
669        centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
670#check get same ellipse parms either way
671#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
672        fail = True
673        i2 = 1
674        while fail:
675            dsp = HKL[i2][3]
676            print '2nd ring: try %.4f'%(dsp)
677            tth = 2.0*asind(wave/(2.*dsp))
678            ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
679            print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1])
680            Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,self.ImageZ)
681            parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
682                'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
683            varyList = [item for item in varyDict if varyDict[item]]
684            if len(Ringp) > 10:
685                chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)
686                tiltp = parmDict['tilt']
687                phip = parmDict['phi']
688                centp = [parmDict['det-X'],parmDict['det-Y']]
689                fail = False
690            else:
691                chip = 1e6
692            ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
693            print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1])
694            Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,self.ImageZ)
695            if len(Ringm) > 10:
696                parmDict['tilt'] *= -1
697                chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)
698                tiltm = parmDict['tilt']
699                phim = parmDict['phi']
700                centm = [parmDict['det-X'],parmDict['det-Y']]
701                fail = False
702            else:
703                chim = 1e6
704            if fail:
705                i2 += 1
706        if chip < chim:
707            data['tilt'] = tiltp
708            data['center'] = centp
709            data['rotation'] = phip
710        else:
711            data['tilt'] = tiltm
712            data['center'] = centm
713            data['rotation'] = phim
714        data['ellipses'].append(ellipsep[:]+('b',))
715        data['rings'].append(np.array(Ringp))
716        data['ellipses'].append(ellipsem[:]+('r',))
717        data['rings'].append(np.array(Ringm))
718        G2plt.PlotImage(self,newImage=True)
719    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
720        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
721    varyList = [item for item in varyDict if varyDict[item]]
722    data['rings'] = []
723    data['ellipses'] = []
724    for i,H in enumerate(HKL):
725        dsp = H[3]
726        tth = 2.0*asind(wave/(2.*dsp))
727        if tth+abs(data['tilt']) > 90.:
728            print 'next line is a hyperbola - search stopped'
729            break
730        if debug:   print 'HKLD:',H[:4],'2-theta: %.4f'%(tth)
731        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
732        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
733        if debug:   print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
734        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
735        if Ring:
736            data['rings'].append(np.array(Ring))
737            rings = np.concatenate((data['rings']),axis=0)
738            if i:
739                chisq = FitDetector(rings,varyList,parmDict,False)
740                data['distance'] = parmDict['dist']
741                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
742                data['rotation'] = parmDict['phi']
743                data['tilt'] = parmDict['tilt']
744                data['DetDepth'] = parmDict['dep']
745                data['chisq'] = chisq
746                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
747                if debug:   print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
748            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
749#            G2plt.PlotImage(self,newImage=True)
750        else:
751            if debug:   print 'insufficient number of points in this ellipse to fit'
752#            break
753    G2plt.PlotImage(self,newImage=True)
754    fullSize = len(self.ImageZ)/scalex
755    if 2*radii[1] < .9*fullSize:
756        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
757    N = len(data['ellipses'])
758    if N > 2:
759        FitDetector(rings,varyList,parmDict)
760        data['wavelength'] = parmDict['wave']
761        data['distance'] = parmDict['dist']
762        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
763        data['rotation'] = parmDict['phi']
764        data['tilt'] = parmDict['tilt']
765        data['DetDepth'] = parmDict['dep']
766    for H in HKL[:N]:
767        ellipse = GetEllipse(H[3],data)
768        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
769    print 'calibration time = ',time.time()-time0
770    G2plt.PlotImage(self,newImage=True)       
771    return True
772   
773def Make2ThetaAzimuthMap(data,masks,iLim,jLim,times): #most expensive part of integration!
774    'Needs a doc string'
775    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
776    pixelSize = data['pixelSize']
777    scalex = pixelSize[0]/1000.
778    scaley = pixelSize[1]/1000.
779   
780    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
781    tax = np.asfarray(tax*scalex,dtype=np.float32)
782    tay = np.asfarray(tay*scaley,dtype=np.float32)
783    nI = iLim[1]-iLim[0]
784    nJ = jLim[1]-jLim[0]
785    t0 = time.time()
786    #make position masks here
787    frame = masks['Frames']
788    tam = ma.make_mask_none((nI,nJ))
789    if frame:
790        tamp = ma.make_mask_none((1024*1024))
791        tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
792            tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
793        tam = ma.mask_or(tam.flatten(),tamp)
794    polygons = masks['Polygons']
795    for polygon in polygons:
796        if polygon:
797            tamp = ma.make_mask_none((1024*1024))
798            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
799                tay.flatten(),len(polygon),polygon,tamp)[:nI*nJ])
800            tam = ma.mask_or(tam.flatten(),tamp)
801    if tam.shape: tam = np.reshape(tam,(nI,nJ))
802    spots = masks['Points']
803    for X,Y,diam in spots:
804        tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
805        tam = ma.mask_or(tam,tamp)
806    times[0] += time.time()-t0
807    t0 = time.time()
808    TA = np.array(GetTthAzmG(tax,tay,data))     #includes geom. corr. as dist**2/d0**2 - most expensive step
809    times[1] += time.time()-t0
810    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
811    return np.array(TA),tam           #2-theta, azimuth & geom. corr. arrays & position mask
812
813def Fill2ThetaAzimuthMap(masks,TA,tam,image):
814    'Needs a doc string'
815    Zlim = masks['Thresholds'][1]
816    rings = masks['Rings']
817    arcs = masks['Arcs']
818    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
819    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist**2/d0**2
820    for tth,thick in rings:
821        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
822    for tth,azm,thick in arcs:
823        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
824        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
825        tam = ma.mask_or(tam.flatten(),tamt*tama)
826    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
827    tabs = np.ones_like(taz)
828    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
829    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
830    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
831    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
832    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
833    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
834    return tax,tay,taz,tad,tabs
835   
836def ImageIntegrate(image,data,masks,blkSize=128,dlg=None,returnN=False):
837    'Needs a doc string'    #for q, log(q) bins need data['binType']
838    import histogram2d as h2d
839    print 'Begin image integration'
840    LUtth = np.array(data['IOtth'])
841    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
842    numAzms = data['outAzimuths']
843    numChans = data['outChannels']
844    azmOff = data['azmthOff']
845    Dazm = (LRazm[1]-LRazm[0])/numAzms
846    if 'log(q)' in data['binType']:
847        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
848    elif 'q' == data['binType']:
849        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
850    elif '2-theta' in data['binType']:
851        lutth = LUtth               
852    dtth = (lutth[1]-lutth[0])/numChans
853    muT = data['SampleAbs'][0]
854    if 'SASD' in data['type']:
855        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
856    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
857    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
858    imageN = len(image)
859    Nx,Ny = data['size']
860    nXBlks = (Nx-1)/blkSize+1
861    nYBlks = (Ny-1)/blkSize+1
862    Nup = nXBlks*nYBlks*3+3
863    tbeg = time.time()
864    Nup = 0
865    if dlg:
866        dlg.Update(Nup)
867    times = [0,0,0,0,0]
868    for iBlk in range(nYBlks):
869        iBeg = iBlk*blkSize
870        iFin = min(iBeg+blkSize,Ny)
871        for jBlk in range(nXBlks):
872            jBeg = jBlk*blkSize
873            jFin = min(jBeg+blkSize,Nx)
874            # next is most expensive step!
875            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin),times)           #2-theta & azimuth arrays & create position mask
876            Nup += 1
877            if dlg:
878                dlg.Update(Nup)
879            Block = image[iBeg:iFin,jBeg:jFin]
880            t0 = time.time()
881            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
882            times[2] += time.time()-t0
883            Nup += 1
884            if dlg:
885                dlg.Update(Nup)
886            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
887            tax = np.where(tax < LRazm[0],tax+360.,tax)
888            if data['SampleAbs'][1]:
889                if 'PWDR' in data['type']:
890                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))
891                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
892                elif 'SASD' in data['type']:    #assumes flat plate sample normal to beam
893                    tabs = G2pwd.Absorb('Fixed',muT,tay)
894            if 'log(q)' in data['binType']:
895                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
896            elif 'q' == data['binType']:
897                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
898            t0 = time.time()
899            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
900                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz*tad/tabs,
901                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
902            times[3] += time.time()-t0
903            Nup += 1
904            if dlg:
905                dlg.Update(Nup)
906    t0 = time.time()
907    NST = np.array(NST,dtype=np.float)
908    H0 = np.divide(H0,NST)
909    H0 = np.nan_to_num(H0)
910    H2 = np.array([tth for tth in np.linspace(lutth[0],lutth[1],numChans+1)])
911    if 'log(q)' in data['binType']:
912        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
913    elif 'q' == data['binType']:
914        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
915    if Dazm:       
916        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
917    else:
918        H1 = LRazm
919    H0 /= npcosd(H2[:-1])           #**2? I don't think so, **1 is right for powders
920    if 'SASD' in data['type']:
921        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
922    if data['Oblique'][1]:
923        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
924    if 'SASD' in data['type'] and data['PolaVal'][1]:
925        #NB: in G2pwd.Polarization azm is defined from plane of polarization, not image x axis!
926        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm-90.)[0] for azm in (H1[:-1]+np.diff(H1)/2.)])
927    Nup += 1
928    if dlg:
929        dlg.Update(Nup)
930    times[4] += time.time()-t0
931    print 'Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
932        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4])
933    print "Elapsed time:","%8.3fs"%(time.time()-tbeg)
934    print 'Integration complete'
935    if returnN:     #As requested by Steven Weigand
936        return H0,H1,H2,NST
937    else:
938        return H0,H1,H2
939   
940def MakeStrStaRing(ring,Image,Controls):
941    ellipse = GetEllipse(ring['Dset'],Controls)
942    pixSize = Controls['pixelSize']
943    scalex = 1000./pixSize[0]
944    scaley = 1000./pixSize[1]
945    Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)).T   #returns x,y,dsp for each point in ring
946    if len(Ring):
947        ring['ImxyObs'] = copy.copy(Ring[:2])
948        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
949        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
950        ring['ImtaObs'] = TA
951        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
952        Ring[0] = TA[0]
953        Ring[1] = TA[1]
954        return Ring,ring
955    else:
956        ring['ImxyObs'] = [[],[]]
957        ring['ImtaObs'] = [[],[]]
958        ring['ImtaCalc'] = [[],[]]
959        return [],[]    #bad ring; no points found
960   
961def FitStrSta(Image,StrSta,Controls):
962    'Needs a doc string'
963   
964    StaControls = copy.deepcopy(Controls)
965    phi = StrSta['Sample phi']
966    wave = Controls['wavelength']
967    StaType = StrSta['Type']
968    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
969
970    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
971        dset = ring['Dset']
972        Ring,R = MakeStrStaRing(ring,Image,StaControls)
973        if len(Ring):
974            ring.update(R)
975            p0 = ring['Emat']
976            val,esd = FitStrain(Ring,p0,dset,wave,phi,StaType)
977            ring['Emat'] = val
978            ring['Esig'] = esd
979    CalcStrSta(StrSta,Controls)
980   
981def CalcStrSta(StrSta,Controls):
982
983    wave = Controls['wavelength']
984    phi = StrSta['Sample phi']
985    StaType = StrSta['Type']
986    for ring in StrSta['d-zero']:
987        Eij = ring['Emat']
988        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
989        th,azm = ring['ImtaObs']
990        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
991        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
992        if StaType == 'True':
993            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
994        else:
995            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
996        ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
997
998def calcFij(omg,phi,azm,th):
999    '''Does something...
1000
1001    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1002
1003    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1004        90 when perp. to sample surface
1005    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1006    :param azm: his chi = azimuth around incident beam
1007    :param th:  his theta = theta
1008    '''
1009    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1010    b = -npcosd(azm)*npcosd(th)
1011    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1012    d = a*npsind(phi)+b*npcosd(phi)
1013    e = a*npcosd(phi)-b*npsind(phi)
1014    Fij = np.array([
1015        [d**2,d*e,c*d],
1016        [d*e,e**2,c*e],
1017        [c*d,c*e,c**2]])
1018    return -Fij
1019
1020def FitStrain(rings,p0,dset,wave,phi,StaType):
1021    'Needs a doc string'
1022    def StrainPrint(ValSig,dset):
1023        print 'Strain tensor for Dset: %.6f'%(dset)
1024        ptlbls = 'names :'
1025        ptstr =  'values:'
1026        sigstr = 'esds  :'
1027        for name,fmt,value,sig in ValSig:
1028            ptlbls += "%s" % (name.rjust(12))
1029            ptstr += fmt % (value)
1030            if sig:
1031                sigstr += fmt % (sig)
1032            else:
1033                sigstr += 12*' '
1034        print ptlbls
1035        print ptstr
1036        print sigstr
1037       
1038    def strainCalc(p,xyd,dset,wave,phi,StaType):
1039        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1040        dspo,azm,dsp = xyd
1041        th = npasind(wave/(2.0*dspo))
1042        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1043        if StaType == 'True':
1044            dspc = dset*np.exp(V)
1045        else:
1046            dspc = dset*(V+1.)
1047        return dspo-dspc
1048       
1049    names = ['e11','e12','e22']
1050    fmt = ['%12.2f','%12.2f','%12.2f']
1051    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1052    vals = list(result[0])
1053    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1054    sig = list(np.sqrt(chisq*np.diag(result[1])))
1055    ValSig = zip(names,fmt,vals,sig)
1056    StrainPrint(ValSig,dset)
1057    return vals,sig
1058   
1059       
Note: See TracBrowser for help on using the repository browser.