source: trunk/GSASIIimage.py @ 1235

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

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