source: trunk/GSASIIimage.py @ 1344

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

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