Changeset 345 for trunk/GSASIIpwd.py


Ignore:
Timestamp:
Aug 18, 2011 8:54:16 PM (11 years ago)
Author:
vondreele
Message:

put new banner in GSASII.BAT
fix for vertical tilt & offset images
speedup of peak calcs
allow change from lam to lam1 & lam2
Pawley refinement fixes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIpwd.py

    r343 r345  
    410410np.seterr(divide='ignore')
    411411
    412 # Voigt function = convolution (norm,cauchy)
    413 
    414 class voigt_gen(st.rv_continuous):
    415     '''
    416     Voigt function
    417     Parameters
    418     -----------------------------------------
    419     '''
    420     def _pdf(self,x,s,g):
    421         z = np.empty([2,len(x)])
    422         z[0] = st.norm.pdf(x,scale=s)
    423         z[1] = st.cauchy.pdf(x,scale=g)
    424         Z = fft(z)
    425         return fftshift(ifft(Z.prod(axis=0))).real
    426 #    def _cdf(self, x):
    427 #    def _ppf(self, q):
    428 #    def _sf(self, x):
    429 #    def _isf(self, q):
    430 #    def _stats(self):
    431 #    def _entropy(self):
    432 
    433 voigt = voigt_gen(name='voigt',shapes='ts,g')
    434        
    435 # Finger-Cox_Jephcoat D(2phi,2th) function for S/L = H/L
     412# Normal distribution
     413
     414# loc = mu, scale = std
     415_norm_pdf_C = 1./math.sqrt(2*math.pi)
     416class norm_gen(st.rv_continuous):
     417       
     418    def pdf(self,x,*args,**kwds):
     419        loc,scale=kwds['loc'],kwds['scale']
     420        x = (x-loc)/scale
     421        return np.exp(-x**2/2.0) * _norm_pdf_C / scale
     422       
     423norm = norm_gen(name='norm',longname='A normal',extradoc="""
     424
     425Normal distribution
     426
     427The location (loc) keyword specifies the mean.
     428The scale (scale) keyword specifies the standard deviation.
     429
     430normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi)
     431""")
     432
     433## Cauchy
     434
     435# median = loc
     436
     437class cauchy_gen(st.rv_continuous):
     438
     439    def pdf(self,x,*args,**kwds):
     440        loc,scale=kwds['loc'],kwds['scale']
     441        x = (x-loc)/scale
     442        return 1.0/np.pi/(1.0+x*x) / scale
     443       
     444cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc="""
     445
     446Cauchy distribution
     447
     448cauchy.pdf(x) = 1/(pi*(1+x**2))
     449
     450This is the t distribution with one degree of freedom.
     451""")
     452   
     453   
     454#GSASII peak fitting routine: Finger, Cox & Jephcoat model       
     455
    436456
    437457class fcjde_gen(st.rv_continuous):
     
    457477    def _pdf(self,x,t,s,dx):
    458478        T = dx*x+t
    459         ax = npcosd(T)**2
     479        ax2 = abs(npcosd(T))
     480        ax = ax2**2
    460481        bx = npcosd(t)**2
    461482        bx = np.where(ax>bx,bx,ax)
    462         fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/np.sqrt(ax),0.0)
     483        fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0)
    463484        fx = np.where(fx > 0.,fx,0.0)
    464485        return fx
    465 #    def _cdf(self, x):
    466 #    def _ppf(self, q):
    467 #    def _sf(self, x):
    468 #    def _isf(self, q):
    469 #    def _stats(self):
    470 #    def _entropy(self):
     486             
     487    def pdf(self,x,*args,**kwds):
     488        loc=kwds['loc']
     489        return self._pdf(x-loc,*args)
    471490       
    472491fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx')
     
    485504    return yb
    486505
    487 def calcPeakFFT(x,fxns,widths,pos,args):
    488     dx = x[1]-x[0]
    489     z = np.empty([len(fxns),len(x)])
    490     for i,(fxn,width,arg) in enumerate(zip(fxns,widths,args)):
    491         z[i] = fxn.pdf(x,loc=pos,scale=width,*arg)*dx
    492     Z = fft(z)
    493     if len(fxns)%2:
    494         return ifft(Z.prod(axis=0)).real
    495     else:
    496         return fftshift(ifft(Z.prod(axis=0))).real
    497                    
    498 def getFCJVoigt(pos,intens,sig,gam,shl,xdata):
    499    
     506def getWidths(pos,sig,gam,shl):
     507    widths = [np.sqrt(sig)/100.,gam/200.]
     508    fwhm = 2.355*widths[0]+2.*widths[1]
     509    fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
     510    fmax = 15.0*fwhm
     511    if pos > 90:
     512        fmin,fmax = [fmax,fmin]         
     513    return widths,fmin,fmax
     514               
     515def getFCJVoigt(pos,intens,sig,gam,shl,xdata):   
    500516    DX = xdata[1]-xdata[0]
    501517    widths,fmin,fmax = getWidths(pos,sig,gam,shl)
    502     x = np.linspace(pos-fmin,pos+fmin,1024)
    503     if pos > 90:
    504         fmin,fmax = [fmax,fmin]         
     518    x = np.linspace(pos-fmin,pos+fmin,256)
    505519    dx = x[1]-x[0]
    506     arg = [pos,shl/10.,dx,]
    507     Df = fcjde.pdf(x,*arg,loc=pos,scale=1.0)
    508     if len(np.nonzero(Df)[0])>5:
    509         fxns = [st.norm,st.cauchy,fcjde]
    510         args = [[],[],arg]
     520    Norm = norm.pdf(x,loc=pos,scale=widths[0])
     521    Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1])
     522    arg = [pos,shl/57.2958,dx,]
     523    FCJ = fcjde.pdf(x,*arg,loc=pos)
     524    if len(np.nonzero(FCJ)[0])>5:
     525        z = np.column_stack([Norm,Cauchy,FCJ]).T
     526        Z = fft(z)
     527        Df = ifft(Z.prod(axis=0)).real
    511528    else:
    512         fxns = [st.norm,st.cauchy]
    513         args = [[],[]]
    514     Df = calcPeakFFT(x,fxns,widths,pos,args)
     529        z = np.column_stack([Norm,Cauchy]).T
     530        Z = fft(z)
     531        Df = fftshift(ifft(Z.prod(axis=0))).real
    515532    Df /= np.sum(Df)
    516533    Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0)
    517534    return intens*Df(xdata)*DX/dx
    518                                    
     535
    519536def getPeakProfile(parmDict,xdata,varyList,bakType):
    520537   
    521538    yb = getBackground('',parmDict,bakType,xdata)
    522539    yc = np.zeros_like(yb)
     540    dx = xdata[1]-xdata[0]
    523541    U = parmDict['U']
    524542    V = parmDict['V']
     
    526544    X = parmDict['X']
    527545    Y = parmDict['Y']
    528     shl = max(parmDict['SH/L'],0.001)
     546    shl = max(parmDict['SH/L'],0.0005)
    529547    Ka2 = False
    530548    if 'Lam1' in parmDict.keys():
    531549        Ka2 = True
    532         lamRatio = parmDict['Lam2']/parmDict['Lam1']
     550        lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1'])
    533551        kRatio = parmDict['I(L2)/I(L1)']
    534552    iPeak = 0
     
    548566            else:
    549567                gam = X/cosd(pos/2.0)+Y*tand(pos/2.0)
    550             gam = max(gam,0.01)             #avoid neg gamma
     568            gam = max(gam,0.001)             #avoid neg gamma
    551569            Wd,fmin,fmax = getWidths(pos,sig,gam,shl)
    552             if pos > 90:
    553                 fmin,fmax = [fmax,fmin]         
    554570            iBeg = np.searchsorted(xdata,pos-fmin)
    555             iFin = np.searchsorted(xdata,pos+fmax)
     571            lenX = len(xdata)
     572            if not iBeg:
     573                iFin = np.searchsorted(xdata,pos+fmin)
     574            elif iBeg == lenX:
     575                iFin = iBeg
     576            else:
     577                iFin = min(lenX,iBeg+int((fmin+fmax)/dx))
    556578            if not iBeg+iFin:       #peak below low limit
    557579                iPeak += 1
     
    561583            yc[iBeg:iFin] += getFCJVoigt(pos,intens,sig,gam,shl,xdata[iBeg:iFin])
    562584            if Ka2:
    563                 pos2 = 2.0*asind(lamRatio*sind(pos/2.0))
    564                 Wd,fmin,fmax = getWidths(pos2,sig,gam,shl)
    565                 if pos > 90:
    566                     fmin,fmax = [fmax,fmin]         
    567                 iBeg = np.searchsorted(xdata,pos2-fmin)
    568                 iFin = np.searchsorted(xdata,pos2+fmax)
    569                 yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
     585                pos2 = pos+lamRatio*tand(pos/2.0)       # + 360/pi * Dlam/lam * tan(th)
     586                kdelt = int((pos2-pos)/dx)               
     587                iBeg = min(lenX,iBeg+kdelt)
     588                iFin = min(lenX,iFin+kdelt)
     589                if iBeg-iFin:
     590                    yc[iBeg:iFin] += getFCJVoigt(pos2,intens*kRatio,sig,gam,shl,xdata[iBeg:iFin])
    570591            iPeak += 1
    571592        except KeyError:        #no more peaks to process
    572593            return yb+yc
    573        
    574 def getWidths(pos,sig,gam,shl):
    575     if shl:
    576         widths = [np.sqrt(sig)/100.,gam/200.,.1]
    577     else:
    578         widths = [np.sqrt(sig)/100.,gam/200.]
    579     fwhm = 2.355*widths[0]+2.*widths[1]
    580     fmin = 10.*(fwhm+shl*abs(npcosd(pos)))
    581     fmax = 15.0*fwhm
    582     return widths,fmin,fmax
    583                
     594
    584595def Dict2Values(parmdict, varylist):
    585596    '''Use before call to leastsq to setup list of values for the parameters
     
    635646                insVary.append(insNames[i])
    636647        instDict = dict(zip(insNames,insVals))
    637         instDict['X'] = max(instDict['X'],0.1)
    638         instDict['Y'] = max(instDict['Y'],0.1)
    639         instDict['SH/L'] = max(instDict['SH/L'],0.001)
     648        instDict['X'] = max(instDict['X'],0.01)
     649        instDict['Y'] = max(instDict['Y'],0.01)
     650        instDict['SH/L'] = max(instDict['SH/L'],0.0001)
    640651        return dataType,instDict,insVary
    641652       
     
    735746        return M
    736747       
    737     Ftol = 0.01
     748    Ftol = 0.0001
    738749    if oneCycle:
    739750        Ftol = 1.0
     
    759770            dlg.SetPosition(wx.Point(screenSize[2]-Size[0]-305,screenSize[1]+5))
    760771            try:
    761                 result = so.leastsq(errPeakProfile,values,full_output=True,ftol=Ftol,
     772                result = so.leastsq(errPeakProfile,values,full_output=True,             #ftol=Ftol,
    762773                    args=(x[xBeg:xFin],y[xBeg:xFin],w[xBeg:xFin],parmDict,varyList,bakType,dlg))
    763774            finally:
Note: See TracChangeset for help on using the changeset viewer.