Changeset 1087


Ignore:
Timestamp:
Oct 7, 2013 12:30:10 PM (8 years ago)
Author:
vondreele
Message:

simulated CW powder data now complete with Poisson errors

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASII.py

    r1086 r1087  
    10791079            else:
    10801080                return False
    1081             x = np.arange(start,end+step,step)
     1081            N = int((end-start)/step)+1
     1082            x = np.linspace(start,end,N,True)
    10821083            N = len(x)
    10831084        rd.powderdata = [
    10841085            np.array(x), # x-axis values
    1085             None, # powder pattern intensities
    1086             None, # 1/sig(intensity)^2 values (weights)
    1087             np.zeros(N), # calc. intensities (zero)
    1088             np.zeros(N), # calc. background (zero)
    1089             np.zeros(N), # obs-calc profiles
     1086            np.zeros_like(x), # powder pattern intensities
     1087            np.ones_like(x), # 1/sig(intensity)^2 values (weights)
     1088            np.zeros_like(x), # calc. intensities (zero)
     1089            np.zeros_like(x), # calc. background (zero)
     1090            np.zeros_like(x), # obs-calc profiles
    10901091            ]
    1091         Tmin = min(rd.powderdata[0])
    1092         Tmax = max(rd.powderdata[0])
     1092        Tmin = rd.powderdata[0][0]
     1093        Tmax = rd.powderdata[0][-1]
    10931094        # data are read, now store them in the tree
    10941095        Id = self.PatternTree.AppendItem(parent=self.root,
  • trunk/GSASIIgrid.py

    r1086 r1087  
    30663066            start,end,step = inp
    30673067        step = abs(step)
    3068         newdata = np.arange(start,end+step,step)
     3068        N = int((end-start)/step)+1
     3069        newdata = np.linspace(start,end,N,True)
    30693070        if len(newdata) < 2: return # too small a range - reject
    30703071        data[1][0] = newdata
    3071         data[1][1] = data[1][2] = None
     3072        data[1][1] = np.zeros_like(newdata)
     3073        data[1][2] = np.ones_like(newdata)
     3074        data[1][3] = np.zeros_like(newdata)
     3075        data[1][4] = np.zeros_like(newdata)
     3076        data[1][5] = np.zeros_like(newdata)
     3077        Tmin = newdata[0]
     3078        Tmax = newdata[-1]
     3079        G2frame.PatternTree.SetItemPyData(GetPatternTreeItemId(G2frame,item,'Limits'),
     3080            [(Tmin,Tmax),[Tmin,Tmax]])
    30723081        UpdatePWHKPlot(G2frame,kind,item) # redisplay data screen
    30733082
  • trunk/GSASIIplot.py

    r1086 r1087  
    777777            X = xye[0]-Parms['Zero'][1]
    778778        if not lenX:
    779             lenX = len(X)           
     779            lenX = len(X)
    780780        Y = xye[1]+offset*N
    781781        if LimitId:
  • trunk/GSASIIstrMath.py

    r1078 r1087  
    1818import numpy.linalg as nl
    1919import scipy.optimize as so
     20import scipy.stats as st
    2021import GSASIIpath
    2122GSASIIpath.SetVersionNumber("$Revision$")
     
    18761877            Limits = calcControls[hfx+'Limits']
    18771878            x,y,w,yc,yb,yd = Histogram['Data']
    1878             W = wtFactor*w
    18791879            yc *= 0.0                           #zero full calcd profiles
    18801880            yb *= 0.0
     
    18821882            xB = np.searchsorted(x,Limits[0])
    18831883            xF = np.searchsorted(x,Limits[1])
     1884            yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
     1885                varylist,Histogram,Phases,calcControls,pawleyLookup)
     1886            yc[xB:xF] += yb[xB:xF]
     1887            if not np.any(y):                   #fill dummy data
     1888                rv = st.poisson(yc[xB:xF])
     1889                y[xB:xF] = rv.rvs()
     1890                w[xB:xF] = np.where(y[xB:xF]>0.,1./y[xB:xF],1.0)
     1891            yd[xB:xF] = y[xB:xF]-yc[xB:xF]
     1892            W = wtFactor*w
     1893            wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
    18841894            Histogram['Residuals']['Nobs'] = ma.count(x[xB:xF])
    18851895            Nobs += Histogram['Residuals']['Nobs']
    18861896            Histogram['Residuals']['sumwYo'] = ma.sum(W[xB:xF]*y[xB:xF]**2)
    18871897            SumwYo += Histogram['Residuals']['sumwYo']
    1888             yc[xB:xF],yb[xB:xF] = getPowderProfile(parmDict,x[xB:xF],
    1889                 varylist,Histogram,Phases,calcControls,pawleyLookup)
    1890             yc[xB:xF] += yb[xB:xF]
    1891             yd[xB:xF] = y[xB:xF]-yc[xB:xF]
    1892             wdy = -ma.sqrt(W[xB:xF])*(yd[xB:xF])
    18931898            Histogram['Residuals']['R'] = min(100.,ma.sum(ma.abs(yd[xB:xF]))/ma.sum(y[xB:xF])*100.)
    18941899            Histogram['Residuals']['wR'] = min(100.,ma.sqrt(ma.sum(wdy**2)/Histogram['Residuals']['sumwYo'])*100.)
Note: See TracChangeset for help on using the changeset viewer.