Changeset 1185


Ignore:
Timestamp:
Jan 8, 2014 3:26:10 PM (9 years ago)
Author:
vondreele
Message:

revamp strain fitting - now does rings one at a time giving each one strain coeff.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIimage.py

    r1184 r1185  
    192192       
    193193    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
    194     fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.3f','%12.5f']
     194    fmt = ['%12.3f','%12.3f','%12.3f','%12.3f','%12.3f','%12.3f','%12.5f']
    195195    p0 = [parmDict[key] for key in varyList]
    196196    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
    … …  
    198198    parmDict.update(zip(varyList,result[0]))
    199199    vals = list(result[0])
    200     chi = np.sqrt(np.sum(ellipseCalcD(result[0],rings.T,varyList,parmDict)**2))
    201     sig = list(chi*np.sqrt(np.diag(result[1])))
     200    sig = list(np.sqrt(chisq*np.diag(result[1])))
    202201    sigList = np.zeros(7)
    203202    for i,name in enumerate(names):
    … …  
    207206    if Print:
    208207        CalibPrint(ValSig)
    209     return chi,chisq
     208    return chisq
    210209                   
    211210def ImageLocalMax(image,w,Xpix,Ypix):
    … …  
    636635    varyList = ['dist','det-X','det-Y','tilt','phi']
    637636    if len(Ringp) > 10:
    638         chip,chisqp = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)
     637        chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)
    639638        tiltp = parmDict['tilt']
    640639        phip = parmDict['phi']
    … …  
    647646    if len(Ringm) > 10:
    648647        parmDict['tilt'] *= -1
    649         chim,chisqm = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)
     648        chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)
    650649        tiltm = parmDict['tilt']
    651650        phim = parmDict['phi']
    … …  
    688687            rings = np.concatenate((data['rings']),axis=0)
    689688            if i:
    690                 chi,chisq = FitDetector(rings,varyList,parmDict,False)
     689                chisq = FitDetector(rings,varyList,parmDict,False)
    691690                data['distance'] = parmDict['dist']
    692691                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
    … …  
    859858    scaley = 1000./pixSize[1]
    860859    Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)).T   #returns x,y,dsp for each point in ring
    861     ring['ImxyObs'] = Ring[:2]
    862     TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
    863     TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
    864     ring['ImtaObs'] = TA
    865     ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
    866     Ring[0] = TA[0]
    867     Ring[1] = TA[1]
    868     return Ring,ring
     860    if len(Ring):
     861        ring['ImxyObs'] = copy.copy(Ring[:2])
     862        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
     863        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
     864        ring['ImtaObs'] = TA
     865        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
     866        Ring[0] = TA[0]
     867        Ring[1] = TA[1]
     868        return Ring,ring
     869    else:
     870        return [],[]    #bad ring; no points found
    869871   
    870872def FitStrSta(Image,StrSta,Controls):
    … …  
    875877    wave = Controls['wavelength']
    876878    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
    877     rings = []
    878879
    879880    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
     881        dset = ring['Dset']
    880882        Ring,R = MakeStrStaRing(ring,Image,StaControls)
    881         ring.update(R)
    882         if len(rings):
    883             rings = np.concatenate((rings,Ring),axis=1)
    884         else:
    885             rings = np.array(Ring)     
    886     E = StrSta['strain']
    887     p0 = [E[0][0],E[0][1],E[1][1]]
    888     E = FitStrain(rings,p0,wave,phi)
    889     StrSta['strain'] = E
     883        if len(Ring):
     884            ring.update(R)
     885            p0 = ring['Emat']
     886            ring['Emat'] = FitStrain(Ring,p0,dset,wave,phi)
    890887    CalcStrSta(StrSta,Controls)
    891888   
    … …  
    896893    E = StrSta['strain']
    897894    for ring in StrSta['d-zero']:
     895        Eij = ring['Emat']
     896        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
    898897        th,azm = ring['ImtaObs']
    899898        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
    … …  
    923922    return -Fij*nptand(th)
    924923
    925 def FitStrain(rings,p0,wave,phi):
    926     'Needs a doc string'
    927     def StrainPrint(ValSig):
    928         print 'Strain tensor:'
     924def FitStrain(rings,p0,dset,wave,phi):
     925    'Needs a doc string'
     926    def StrainPrint(ValSig,dset):
     927        print 'Strain tensor for Dset: %.6f'%(dset)
    929928        ptlbls = 'names :'
    930929        ptstr =  'values:'
    … …  
    941940        print sigstr
    942941       
    943     def strainCalc(p,xyd,wave,phi):
     942    def strainCalc(p,xyd,dset,wave,phi):
    944943        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
    945944        dspo,azm,dsp = xyd
    946945        th = npasind(wave/(2.0*dspo))
    947946        V = 1.+np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
    948         dspc = dsp*V
    949         return dspo-dspc
    950        
    951     names = ['e11','e12','e22','e33','e13','e23']
    952     fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.5f']
    953     p1 = [-20000.,0,5000.] 
    954     result = leastsq(strainCalc,p1,args=(rings,wave,phi),full_output=True)
     947        dspc = dset*V
     948        return dspo**2-dspc**2
     949       
     950    names = ['e11','e12','e22']
     951    fmt = ['%12.2f','%12.2f','%12.2f']
     952    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi),full_output=True)
    955953    vals = list(result[0])
    956     chi = np.sqrt(np.sum(strainCalc(result[0],rings,wave,phi)**2))
    957     sig = list(chi*np.sqrt(np.diag(result[1])))
     954    chisq = np.sum(result[2]['fvec']**2)
     955    sig = list(np.sqrt(chisq*np.diag(result[1])))
    958956    ValSig = zip(names,fmt,vals,sig)
    959     StrainPrint(ValSig)
    960     return np.array([[vals[0],vals[1],0],[vals[1],vals[2],0],[0,0,0]])
    961    
    962        
     957    StrainPrint(ValSig,dset)
     958    return vals
     959   
     960       

  • trunk/GSASIIimgGUI.py

    r1184 r1185  
    13221322    def OnAppendDzero(event):
    13231323        data['d-zero'].append({'Dset':1.0,'Dcalc':0.0,'pixLimit':10,'cutoff':10.0,
    1324             'ImxyObs':[[],[]],'ImtaObs':[[],[]],'ImtaCalc':[[],[]]})
     1324            'ImxyObs':[[],[]],'ImtaObs':[[],[]],'ImtaCalc':[[],[]],'Emat':[1.0,1.0,1.0]})
    13251325        UpdateStressStrain(G2frame,data)
    13261326           
    … …  
    13711371                File = open(filename,'w')
    13721372                save = {}
    1373                 keys = ['Type','Sample phi','Sample z','strain']
    1374                 keys2 = ['Dset','Dcalc','pixLimit','cutoff']
     1373                keys = ['Type','Sample phi','Sample z']
     1374                keys2 = ['Dset','Dcalc','pixLimit','cutoff','Emat']
    13751375                File.write('{\n\t')
    13761376                for key in keys:
    … …  
    13961396            G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Image Controls'))
    13971397        G2img.FitStrSta(G2frame.ImageZ,data,Controls)
     1398        print 'Strain fitting finished'
    13981399        UpdateStressStrain(G2frame,data)
    13991400        G2plt.PlotExposedImage(G2frame,event=event)
    … …  
    14481449            Obj = event.GetEventObject()
    14491450            try:
    1450                 value = min(10.0,max(1.0,float(Obj.GetValue())))
     1451                value = min(20.0,max(0.25,float(Obj.GetValue())))
    14511452            except ValueError:
    14521453                value = 1.0
    … …  
    14541455            data['d-zero'][Indx[Obj.GetId()]]['Dset'] = value
    14551456            Ring,R = G2img.MakeStrStaRing(data['d-zero'][Indx[Obj.GetId()]],G2frame.ImageZ,Controls)
    1456             data['d-zero'][Indx[Obj.GetId()]].update(R)
     1457            if len(Ring):
     1458                data['d-zero'][Indx[Obj.GetId()]].update(R)
     1459            else:
     1460                G2frame.ErrorDialog('Strain peak selection','WARNING - No points found for this ring selection')
     1461        #sort them on d-spacing?
    14571462            UpdateStressStrain(G2frame,data)
    14581463            G2plt.PlotExposedImage(G2frame,event=event)
    … …  
    15191524            delIndx.append(dzeroDelete)
    15201525            dzeroSizer.Add(dzeroDelete,0,wx.ALIGN_CENTER_VERTICAL)
     1526           
     1527            dzeroSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,label=(' Strain tensor:')),
     1528                wx.ALIGN_CENTER_VERTICAL)
     1529            names = ['e11','e12','e22']
     1530            for i in range(3):
     1531                dzeroSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,label=names[i]),0,wx.ALIGN_CENTER_VERTICAL)
     1532                tensorElem = wx.TextCtrl(G2frame.dataDisplay,-1,value='%.2f'%(dzero['Emat'][i]),style=wx.TE_READONLY)
     1533                tensorElem.SetBackgroundColour(VERY_LIGHT_GREY)
     1534                dzeroSizer.Add(tensorElem,0,wx.ALIGN_CENTER_VERTICAL)
     1535            dzeroSizer.Add((5,5),0)             
    15211536        return dzeroSizer
    15221537       
    1523     def StrainSizer():
    1524        
    1525         def OnTensor(event):
    1526             Obj = event.GetEventObject()
    1527             try:
    1528                 value = float(Obj.GetValue())
    1529             except ValueError:
    1530                 value = 1.0
    1531             i,j = Indx[Obj.GetId()]
    1532             data['strain'][i][j] = value
    1533             if i != j:
    1534                 data['strain'][j][i] = value
    1535             UpdateStressStrain(G2frame,data)
    1536             G2img.CalcStrSta(data,Controls)
    1537             G2plt.PlotStrain(G2frame,data,newPlot=False)
    1538 
    1539         Indx = {}
    1540         strainSizer = wx.BoxSizer(wx.VERTICAL)
    1541         strainSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,label=(' Strain tensor:')),
    1542             0,wx.ALIGN_CENTER_VERTICAL)
    1543         tensorSizer = wx.FlexGridSizer(3,6,5,5)
    1544         names = [[' e11','e12','e13'],[' e21','e22','e23'],[' e31','e32','e33']]
    1545         for i in range(3):
    1546             for j in range(3):
    1547                 tensorSizer.Add(wx.StaticText(G2frame.dataDisplay,-1,label=names[i][j]),0,wx.ALIGN_CENTER_VERTICAL)
    1548 #                tensorElem = wx.TextCtrl(G2frame.dataDisplay,-1,value='%.2f'%(data['strain'][i][j]),style=wx.TE_READONLY)
    1549 #                tensorElem.SetBackgroundColour(VERY_LIGHT_GREY)
    1550                 tensorElem = wx.TextCtrl(G2frame.dataDisplay,-1,value='%.2f'%(data['strain'][i][j]),style=wx.TE_PROCESS_ENTER)
    1551                 tensorElem.Bind(wx.EVT_TEXT_ENTER,OnTensor)
    1552                 tensorElem.Bind(wx.EVT_KILL_FOCUS,OnTensor)
    1553                 Indx[tensorElem.GetId()] = [i,j]
    1554                 tensorSizer.Add(tensorElem,0,wx.ALIGN_CENTER_VERTICAL)
    1555         strainSizer.Add(tensorSizer)
    1556         return strainSizer
    1557 
    15581538    if G2frame.dataDisplay:
    15591539        G2frame.dataDisplay.Destroy()
    … …  
    15751555    mainSizer.Add((5,10),0)
    15761556    mainSizer.Add(DzeroSizer())
    1577     mainSizer.Add((5,10),0)
    1578     mainSizer.Add(StrainSizer())
    15791557   
    15801558    mainSizer.Layout()   

  • trunk/GSASIIplot.py

    r1184 r1185  
    12641264            Page.canvas.SetCursor(wx.CROSS_CURSOR)
    12651265            try:
    1266                 G2frame.G2plotNB.status.SetStatusText('d-spacing =%9.5f Azimuth =%9.3f'%(xpos,ypos),1)                   
     1266                G2frame.G2plotNB.status.SetStatusText('d-spacing =%9.5f Azimuth =%9.3f'%(ypos,xpos),1)                   
    12671267            except TypeError:
    12681268                G2frame.G2plotNB.status.SetStatusText('Select Strain pattern first',1)
    … …  
    12871287    G2frame.G2plotNB.status.DestroyChildren()
    12881288    Plot.set_title('Strain')
    1289     Plot.set_xlabel(r'd-spacing',fontsize=14)
    1290     Plot.set_ylabel(r'Azimuth',fontsize=14)
     1289    Plot.set_ylabel(r'd-spacing',fontsize=14)
     1290    Plot.set_xlabel(r'Azimuth',fontsize=14)
    12911291    colors=['b','g','r','c','m','k']
    12921292    for N,item in enumerate(data['d-zero']):
    1293         X,Y = np.array(item['ImtaObs'])
     1293        Y,X = np.array(item['ImtaObs'])         #plot azimuth as X & d-spacing as Y
    12941294        Plot.plot(X,Y,colors[N%6]+'+',picker=False)
    1295         X,Y = np.array(item['ImtaCalc'])
     1295        Y,X = np.array(item['ImtaCalc'])
    12961296        Plot.plot(X,Y,colors[N%6],picker=False)
    12971297    if not newPlot:
    … …  
    22922292        except TypeError:
    22932293            return
    2294         if PickName not in ['Image Controls','Masks']:
     2294        if PickName not in ['Image Controls','Masks','Stress/Strain']:
    22952295            return
    22962296        pixelSize = Data['pixelSize']
    … …  
    23652365            G2imG.UpdateMasks(G2frame,Masks)
    23662366            PlotImage(G2frame,newImage=False)
     2367        elif PickName == 'Stress/Strain':
     2368            Xpos,Ypos = [event.xdata,event.ydata]
     2369            if not Xpos or not Ypos or Page.toolbar._active:  #got point out of frame or zoom/pan selected
     2370                return
     2371            dsp = G2img.GetDsp(Xpos,Ypos,Data)
     2372            StrSta['d-zero'].append({'Dset':dsp,'Dcalc':0.0,'pixLimit':10,'cutoff':10.0,
     2373                'ImxyObs':[[],[]],'ImtaObs':[[],[]],'ImtaCalc':[[],[]],'Emat':[1.0,1.0,1.0]})
     2374            R,r = G2img.MakeStrStaRing(StrSta['d-zero'][-1],G2frame.ImageZ,Data)
     2375            if not len(R):
     2376                del StrSta['d-zero'][-1]
     2377                G2frame.ErrorDialog('Strain peak selection','WARNING - No points found for this ring selection')
     2378            G2imG.UpdateStressStrain(G2frame,StrSta)
     2379            PlotStrain(G2frame,StrSta)
     2380            PlotImage(G2frame,newPlot=False)           
    23672381        else:
    23682382            Xpos,Ypos = [event.xdata,event.ydata]
    … …  
    25872601                for ring in Data['rings']:
    25882602                    xring,yring = np.array(ring).T[:2]
    2589                     Plot.plot(xring,yring,'+',color=colors[N%6])
     2603                    Plot.plot(xring,yring,'.',color=colors[N%6])
    25902604                    N += 1           
    25912605            for ellipse in Data['ellipses']:      #what about hyperbola?
    … …  
    25952609                    Plot.text(cent[0],cent[1],'+',color=col,ha='center',va='center')
    25962610        if G2frame.PatternTree.GetItemText(G2frame.PickId) in 'Stress/Strain':
    2597             print 'plot stress/strain stuff'
    25982611            for N,ring in enumerate(StrSta['d-zero']):
    25992612                xring,yring = ring['ImxyObs']
    2600                 Plot.plot(xring,yring,colors[N%6]+'o')
     2613                Plot.plot(xring,yring,colors[N%6]+'.')
    26012614        #masks - mask lines numbered after integration limit lines
    26022615        spots = Masks['Points']
Note: See TracChangeset for help on using the changeset viewer.