source: trunk/GSASIIimage.py @ 1392

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

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