source: trunk/GSASIIimage.py @ 1387

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

fixes to seq display for strain refinements
get sign right on strain calcs
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 38.7 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2014-06-18 17:08:07 +0000 (Wed, 18 Jun 2014) $
5# $Author: vondreele $
6# $Revision: 1387 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 1387 2014-06-18 17:08:07Z 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: 1387 $")
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,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):
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    chisq = 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    data['chisq'] = chisq
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, Np: %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                data['chisq'] = chisq
720                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
721                print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
722            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
723#            G2plt.PlotImage(self,newImage=True)
724        else:
725            print 'insufficient number of points in this ellipse to fit'
726#            break
727    G2plt.PlotImage(self,newImage=True)
728    fullSize = len(self.ImageZ)/scalex
729    if 2*radii[1] < .9*fullSize:
730        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
731    N = len(data['ellipses'])
732    if N > 2:
733        FitDetector(rings,varyList,parmDict)
734        data['distance'] = parmDict['dist']
735        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
736        data['rotation'] = parmDict['phi']
737        data['tilt'] = parmDict['tilt']
738        data['DetDepth'] = parmDict['dep']
739    for H in HKL[:N]:
740        ellipse = GetEllipse(H[3],data)
741        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
742    print 'calibration time = ',time.time()-time0
743    G2plt.PlotImage(self,newImage=True)       
744    return True
745   
746def Make2ThetaAzimuthMap(data,masks,iLim,jLim):
747    'Needs a doc string'
748    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
749    pixelSize = data['pixelSize']
750    scalex = pixelSize[0]/1000.
751    scaley = pixelSize[1]/1000.
752   
753    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
754    tax = np.asfarray(tax*scalex,dtype=np.float32)
755    tay = np.asfarray(tay*scaley,dtype=np.float32)
756    nI = iLim[1]-iLim[0]
757    nJ = jLim[1]-jLim[0]
758    #make position masks here
759    frame = masks['Frames']
760    tam = ma.make_mask_none((nI,nJ))
761    if frame:
762        tamp = ma.make_mask_none((1024*1024))
763        tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
764            tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
765        tam = ma.mask_or(tam.flatten(),tamp)
766    polygons = masks['Polygons']
767    for polygon in polygons:
768        if polygon:
769            tamp = ma.make_mask_none((1024*1024))
770            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
771                tay.flatten(),len(polygon),polygon,tamp)[:nI*nJ])
772            tam = ma.mask_or(tam.flatten(),tamp)
773    if tam.shape: tam = np.reshape(tam,(nI,nJ))
774    spots = masks['Points']
775    for X,Y,diam in spots:
776        tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
777        tam = ma.mask_or(tam,tamp)
778    TA = np.array(GetTthAzmD(tax,tay,data))     #includes geom. corr. as dist**2/d0**2
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    t0 = time.time()
833    Nup = 0
834    if dlg:
835        dlg.Update(Nup)
836    for iBlk in range(nYBlks):
837        iBeg = iBlk*blkSize
838        iFin = min(iBeg+blkSize,Ny)
839        for jBlk in range(nXBlks):
840            jBeg = jBlk*blkSize
841            jFin = min(jBeg+blkSize,Nx)               
842            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
843           
844            Nup += 1
845            if dlg:
846                dlg.Update(Nup)
847            Block = image[iBeg:iFin,jBeg:jFin]
848            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
849            Nup += 1
850            if dlg:
851                dlg.Update(Nup)
852            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
853            tax = np.where(tax < LRazm[0],tax+360.,tax)
854            if data['SampleAbs'][1]:
855                if 'PWDR' in data['type']:
856                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))
857                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
858                elif 'SASD' in data['type']:    #assumes flat plate sample normal to beam
859                    tabs = G2pwd.Absorb('Fixed',muT,tay)
860            if 'log(q)' in data['binType']:
861                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
862            elif 'q' == data['binType']:
863                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
864            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
865                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz*tad/tabs,
866                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
867            Nup += 1
868            if dlg:
869                dlg.Update(Nup)
870    NST = np.array(NST,dtype=np.float)
871    H0 = np.divide(H0,NST)
872    H0 = np.nan_to_num(H0)
873    del NST
874    H2 = np.array([tth for tth in np.linspace(lutth[0],lutth[1],numChans+1)])
875    if 'log(q)' in data['binType']:
876        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
877    elif 'q' == data['binType']:
878        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
879    if Dazm:       
880        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
881    else:
882        H1 = LRazm
883    H0 /= npcosd(H2[:-1])           #**2? I don't think so, **1 is right for powders
884    if 'SASD' in data['type']:
885        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
886    if data['Oblique'][1]:
887        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
888    if 'SASD' in data['type'] and data['PolaVal'][1]:
889        #NB: in G2pwd.Polarization azm is defined from plane of polarization, not image x axis!
890        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm-90.)[0] for azm in (H1[:-1]+np.diff(H1)/2.)])
891    Nup += 1
892    if dlg:
893        dlg.Update(Nup)
894    t1 = time.time()
895    print 'Integration complete'
896    print "Elapsed time:","%8.3f"%(t1-t0), "s"
897    return H0,H1,H2
898   
899def MakeStrStaRing(ring,Image,Controls):
900    ellipse = GetEllipse(ring['Dset'],Controls)
901    pixSize = Controls['pixelSize']
902    scalex = 1000./pixSize[0]
903    scaley = 1000./pixSize[1]
904    Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)).T   #returns x,y,dsp for each point in ring
905    if len(Ring):
906        ring['ImxyObs'] = copy.copy(Ring[:2])
907        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
908        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
909        ring['ImtaObs'] = TA
910        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
911        Ring[0] = TA[0]
912        Ring[1] = TA[1]
913        return Ring,ring
914    else:
915        ring['ImxyObs'] = [[],[]]
916        ring['ImtaObs'] = [[],[]]
917        ring['ImtaCalc'] = [[],[]]
918        return [],[]    #bad ring; no points found
919   
920def FitStrSta(Image,StrSta,Controls):
921    'Needs a doc string'
922   
923    StaControls = copy.deepcopy(Controls)
924    phi = StrSta['Sample phi']
925    wave = Controls['wavelength']
926    StaType = StrSta['Type']
927    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
928
929    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
930        dset = ring['Dset']
931        Ring,R = MakeStrStaRing(ring,Image,StaControls)
932        if len(Ring):
933            ring.update(R)
934            p0 = ring['Emat']
935            val,esd = FitStrain(Ring,p0,dset,wave,phi,StaType)
936            ring['Emat'] = val
937            ring['Esig'] = esd
938    CalcStrSta(StrSta,Controls)
939   
940def CalcStrSta(StrSta,Controls):
941
942    wave = Controls['wavelength']
943    phi = StrSta['Sample phi']
944    StaType = StrSta['Type']
945    for ring in StrSta['d-zero']:
946        Eij = ring['Emat']
947        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
948        th,azm = ring['ImtaObs']
949        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
950        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
951        if StaType == 'True':
952            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
953        else:
954            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
955        ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
956
957def calcFij(omg,phi,azm,th):
958    '''Does something...
959
960    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
961
962    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
963        90 when perp. to sample surface
964    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
965    :param azm: his chi = azimuth around incident beam
966    :param th:  his theta = theta
967    '''
968    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
969    b = -npcosd(azm)*npcosd(th)
970    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
971    d = a*npcosd(phi)-b*npsind(phi)
972    e = a*npsind(phi)+b*npcosd(phi)
973    Fij = np.array([
974        [d**2,d*e,c*d],
975        [d*e,e**2,c*e],
976        [c*d,c*e,c**2]])
977    return -Fij*nptand(th)
978
979def FitStrain(rings,p0,dset,wave,phi,StaType):
980    'Needs a doc string'
981    def StrainPrint(ValSig,dset):
982        print 'Strain tensor for Dset: %.6f'%(dset)
983        ptlbls = 'names :'
984        ptstr =  'values:'
985        sigstr = 'esds  :'
986        for name,fmt,value,sig in ValSig:
987            ptlbls += "%s" % (name.rjust(12))
988            ptstr += fmt % (value)
989            if sig:
990                sigstr += fmt % (sig)
991            else:
992                sigstr += 12*' '
993        print ptlbls
994        print ptstr
995        print sigstr
996       
997    def strainCalc(p,xyd,dset,wave,phi,StaType):
998        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
999        dspo,azm,dsp = xyd
1000        th = npasind(wave/(2.0*dspo))
1001        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1002        if StaType == 'True':
1003            dspc = dset*np.exp(V)
1004        else:
1005            dspc = dset*(V+1.)
1006        return dspo-dspc
1007       
1008    names = ['e11','e12','e22']
1009    fmt = ['%12.2f','%12.2f','%12.2f']
1010    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1011    vals = list(result[0])
1012    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1013    sig = list(np.sqrt(chisq*np.diag(result[1])))
1014    ValSig = zip(names,fmt,vals,sig)
1015    StrainPrint(ValSig,dset)
1016    return vals,sig
1017   
1018       
Note: See TracBrowser for help on using the repository browser.