Changeset 4910 for trunk


Ignore:
Timestamp:
May 20, 2021 11:19:53 AM (4 months ago)
Author:
vondreele
Message:

Add Wilson statistics to Reflection list menu - computes Wilson plot, <E2-1>, etc. Useful for detecting twins, inversion symmetry, etc.

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIdataGUI.py

    r4890 r4910  
    60836083       
    60846084        # PDR / Reflection Lists
    6085         G2G.Define_wxId('wxID_SELECTPHASE','wxID_SHOWHIDEEXTINCT' ) #some wxIDs defined above in PWDR & SASD
     6085        G2G.Define_wxId('wxID_SELECTPHASE','wxID_SHOWHIDEEXTINCT','wxID_WILSONSTAT' ) #some wxIDs defined above in PWDR & SASD
    60866086        self.ReflMenu = wx.MenuBar()
    60876087        self.PrefillDataMenu(self.ReflMenu)
     
    60926092        self.ReflEdit.Append(G2G.wxID_PWDHKLPLOT,'Plot HKLs','Plot HKLs in 2D')
    60936093        self.ReflEdit.Append(G2G.wxID_PWD3DHKLPLOT,'Plot 3D HKLs','Plot HKLs in 3D')
     6094        self.ReflEdit.Append(G2G.wxID_WILSONSTAT,'Wilson statistics')
    60946095        self.HideShow = self.ReflEdit.Append(G2G.wxID_SHOWHIDEEXTINCT,'Show/hide extinct reflections')
    60956096        self.PostfillDataMenu()
  • trunk/GSASIImath.py

    r4903 r4910  
    42824282    return Atoms
    42834283
    4284 ################################################################################
    4285 ##### Dysnomia setup & return stuff
    4286 ################################################################################
    4287    
    4288 
    4289    
     4284def DoWilsonStat(refList,Super,normEle,Inst):
     4285    ns = 0
     4286    if Super:
     4287        ns=1
     4288    Eldata = G2el.GetElInfo(normEle,Inst)
     4289    RefList = np.array(refList).T
     4290    dsp = RefList[4+ns]
     4291    if  'P' in Inst['Type'][0]:
     4292        Fosq = RefList[8+ns]
     4293    else:
     4294        Fosq = RefList[5+ns]
     4295    SQ = 1./(2.*dsp)
     4296    if 'X' in Inst['Type'][0]:
     4297        Esq = Fosq/G2el.ScatFac(Eldata,SQ)**2
     4298    else:
     4299        Esq = Fosq
     4300    SQ2 = SQ**2
     4301    SQ2min = np.min(SQ2)
     4302    SQ2max = np.max(SQ2)
     4303    SQbins = np.linspace(SQ2min,SQ2max,50,True)
     4304    step = SQbins[1]-SQbins[0]
     4305    SQbins += step/2.
     4306    E2bins = np.zeros_like(SQbins)
     4307    nE2bins = np.zeros_like(SQbins)
     4308    for sq,e in zip(list(SQ2),list(Esq)):
     4309        i = np.searchsorted(SQbins,sq)
     4310        E2bins[i] += e
     4311        nE2bins[i] += 1
     4312    E2bins /= nE2bins
     4313    E2bins = np.nan_to_num(E2bins)
     4314    A = np.vstack([SQbins,np.ones_like(SQbins)]).T
     4315    result = nl.lstsq(A,np.log(E2bins),rcond=None)
     4316    twoB,lnscale = result[0]    #twoB = -2B
     4317    scale = np.exp(lnscale)
     4318    E2calc = lnscale+twoB*SQbins
     4319    normE = np.sqrt(np.where(Esq>0.,Esq,0.)/scale)*np.exp(-0.5*twoB*SQ2)
     4320    return [np.mean(normE),np.mean(normE**2),np.mean(np.abs(-1.+normE**2))],[SQbins,np.log(E2bins),E2calc]
     4321
    42904322################################################################################
    42914323##### single peak fitting profile fxn stuff
  • trunk/GSASIIpwdGUI.py

    r4862 r4910  
    5858if '2' in platform.python_version_tuple()[0]:
    5959    GkDelta = unichr(0x0394)
     60    GkSigma = unichr(0x03a3)
     61    GkTheta = unichr(0x03f4)
     62    Gklambda = unichr(0x03bb)
    6063    Pwr10 = unichr(0x0b9)+unichr(0x2070)
    6164    Pwr20 = unichr(0x0b2)+unichr(0x2070)
     
    6871else:
    6972    GkDelta = chr(0x0394)
     73    GkSigma = chr(0x03a3)
     74    GkTheta = chr(0x03f4)
     75    Gklambda = chr(0x03bb)
    7076    Pwr10 = chr(0x0b9)+chr(0x2070)
    7177    Pwr20 = chr(0x0b2)+chr(0x2070)
     
    7480    Pwrm6 = chr(0x207b)+chr(0x2076)
    7581    Pwrm4 = chr(0x207b)+chr(0x2074)
    76     Angstr = chr(0x00c5)   
     82    Angstr = chr(0x00c5)
    7783    superMinusOne = chr(0xaf)+chr(0xb9)
    7884# trig functions in degrees
     
    48294835        G2plt.Plot3DSngl(G2frame,newPlot=True,Data=controls,hklRef=refList,Title=phaseName)
    48304836       
     4837    def OnWilsonStat(event):
     4838        ''' Show Wilson plot for PWDR and HKLF & return Wilson statistics <<E>, <E^2> & <E^2-1> to console
     4839        '''
     4840        phaseName = G2frame.RefList
     4841        if phaseName not in ['Unknown',]:
     4842            pId = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Phases')
     4843            phaseId =  G2gd.GetGPXtreeItemId(G2frame,pId,phaseName)
     4844            General = G2frame.GPXtree.GetItemPyData(phaseId)['General']
     4845            Super = General.get('Super',0)
     4846        else:
     4847            Super = 0
     4848        if 'list' in str(type(data)):   #single crystal data is 2 dict in list
     4849            refList = data[1]['RefList']
     4850        else:                           #powder data is a dict of dicts; each same structure as SC 2nd dict
     4851            if 'RefList' in data[phaseName]:
     4852                refList = np.array(data[phaseName]['RefList'])
     4853            else:
     4854                wx.MessageBox('No reflection list - do Refine first',caption='Reflection plotting')
     4855                return
     4856       
     4857        PE = G2elemGUI.PickElement(G2frame,ifNone=False)
     4858        if PE.ShowModal() == wx.ID_OK:
     4859            normEle = PE.Elem.strip()
     4860        PE.Destroy()
     4861        Estat,Ehist = G2mth.DoWilsonStat(refList,Super,normEle,Inst)
     4862        print(' Wilson statistics   : <|E|>: %.3f, <|E^2|>: %.3f, <|E^2-1|>: %.3f'%(Estat[0]/Estat[1],1.,Estat[2]/Estat[1]))
     4863        print(' Expected: random P-1: <|E|>: 0.798, <|E^2|>: 1.000, <|E^2-1|>: 0.968')
     4864        print('           random P1 : <|E|>: 0.886, <|E^2|>: 1.000, <|E^2-1|>: 0.736')
     4865        XY = [[Ehist[0],Ehist[1]],[Ehist[0],Ehist[2]]]
     4866        G2plt.PlotXY(G2frame,XY,labelX='sin$^2$%s/%s$^2$'%(GkTheta,Gklambda),
     4867            labelY=r'ln(<|F$_o$|$^2$>/%sf$^2$)'%GkSigma,newPlot=True,Title='Wilson plot')
     4868       
     4869       
    48314870    def MakeReflectionTable(phaseName):
    48324871        '''Returns a wx.grid table (G2G.Table) containing a list of all reflections
     
    50205059        G2frame.Bind(wx.EVT_MENU, OnPlot1DHKL, id=G2G.wxID_1DHKLSTICKPLOT)
    50215060        G2frame.Bind(wx.EVT_MENU, OnPlot3DHKL, id=G2G.wxID_PWD3DHKLPLOT)
     5061        G2frame.Bind(wx.EVT_MENU, OnWilsonStat, id=G2G.wxID_WILSONSTAT)
    50225062        G2frame.Bind(wx.EVT_MENU, OnToggleExt, id=G2G.wxID_SHOWHIDEEXTINCT)
    50235063        G2frame.dataWindow.SelectPhase.Enable(False)
     
    50275067        G2frame.Bind(wx.EVT_MENU, OnPlotHKL, id=G2G.wxID_PWDHKLPLOT)
    50285068        G2frame.Bind(wx.EVT_MENU, OnPlot3DHKL, id=G2G.wxID_PWD3DHKLPLOT)
     5069        G2frame.Bind(wx.EVT_MENU, OnWilsonStat, id=G2G.wxID_WILSONSTAT)
    50295070        G2frame.dataWindow.SelectPhase.Enable(False)
    50305071           
Note: See TracChangeset for help on using the changeset viewer.