Changeset 4910 for trunk/GSASIImath.py


Ignore:
Timestamp:
May 20, 2021 11:19:53 AM (7 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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.