Changeset 1269 for trunk/GSASIIsasd.py


Ignore:
Timestamp:
Apr 4, 2014 12:38:43 PM (8 years ago)
Author:
vondreele
Message:

more model fitting in SASD
add a possible for plotting hkls from powder data - commented out fro now

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIsasd.py

    r1268 r1269  
    9898    NP = 200
    9999    alp = np.linspace(0,np.pi/2.,NP)
    100     LBessArg = 0.5*L*(Q[:,np.newaxis]*np.cos(alp))
     100    QL = Q[:,np.newaxis]*L
     101    LBessArg = 0.5*QL[:,:,np.newaxis]**np.cos(alp)
    101102    LBess = np.where(LBessArg<1.e-6,1.,np.sin(LBessArg)/LBessArg)
    102103    QR = Q[:,np.newaxis]*R
    103104    SBessArg = QR[:,:,np.newaxis]*np.sin(alp)
    104105    SBess = np.where(SBessArg<1.e-6,0.5,scsp.jv(1,SBessArg)/SBessArg)
    105     LSBess = LBess[:,np.newaxis,:]*SBess
     106    LSBess = LBess*SBess
    106107    return np.sqrt(2.*np.pi*np.sum(np.sin(alp)*LSBess**2,axis=2)/NP)
    107108   
     
    114115    '''
    115116    R = args[0]
    116     return CylinderFF(Q,R,2.*L)   
     117    return CylinderFF(Q,R,args=[2.*L])   
    117118   
    118119def CylinderARFF(Q,R,args):
     
    124125    '''
    125126    AR = args[0]
    126     return CylinderFF(Q,R,2.*R*AR)   
     127    return CylinderFF(Q,R,args=[2.*R*AR])   
    127128   
    128129def UniSphereFF(Q,R,args=0):
     
    165166    '''
    166167    AR = args[0]
    167     return UniRodFF(Q,R,[AR*R,])
     168    return UniRodFF(Q,R,args=[AR*R,])
    168169   
    169170def UniDiskFF(Q,R,args):
     
    204205    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*Rg3/np.sqrt(6)))**3
    205206    FF = np.exp(-Q[:,np.newaxis]**2*Rg3**2/3.)+(B3/QstV)*np.exp(-Q[:,np.newaxis]**2*R**2/3.)
    206     QstV = Q/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
     207    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*R/np.sqrt(6)))**3
    207208    FF += (B2/QstV**2)*np.exp(-Q[:,np.newaxis]**2*T**2/3.)
    208209    QstV = Q[:,np.newaxis]/(scsp.erf(Q[:,np.newaxis]*T/np.sqrt(6)))**3
     
    316317################################################################################
    317318
    318 def LogNormalDist(x,pos,scale,shape):
     319def LogNormalDist(x,pos,args):
    319320    ''' Standard LogNormal distribution - numpy friendly on x axis
    320321    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
     
    325326    returns float: LogNormal distribution
    326327    '''
     328    scale,shape = args
    327329    return np.exp(-np.log((x-pos)/scale)**2/(2.*shape**2))/(np.sqrt(2.*np.pi)*(x-pos)*shape)
    328330   
    329 def GaussDist(x,pos,scale,shape):
     331def GaussDist(x,pos,args):
    330332    ''' Standard Normal distribution - numpy friendly on x axis
    331333    param float x: independent axis (can be numpy array)
     
    335337    returns float: Normal distribution
    336338    '''
    337     return (1./scale*np.sqrt(2.*np.pi))*np.exp(-(x-pos)**2/(2.*scale**2))
    338    
    339 def LSWDist(x,pos,scale,shape):
     339    scale = args[0]
     340    return (1./(scale*np.sqrt(2.*np.pi)))*np.exp(-(x-pos)**2/(2.*scale**2))
     341   
     342def LSWDist(x,pos,args=[]):
    340343    ''' Lifshitz-Slyozov-Wagner Ostwald ripening distribution - numpy friendly on x axis
    341344    ref:
     
    348351    redX = x/pos
    349352    result = (81.*2**(-5/3.))*(redX**2*np.exp(-redX/(1.5-redX)))/((1.5-redX)**(11/3.)*(3.-redX)**(7/3.))
    350     return result/pos
    351    
    352 def SchulzZimmDist(x,pos,scale,shape):
     353   
     354    return np.nan_to_num(result/pos)
     355   
     356def SchulzZimmDist(x,pos,args):
    353357    ''' Schulz-Zimm macromolecule distribution - numpy friendly on x axis
    354358    ref: http://goldbook.iupac.org/S05502.html
     
    359363    returns float: Schulz-Zimm distribution
    360364    '''
     365    scale = args[0]
    361366    b = (2.*pos/scale)**2
    362367    a = b/pos
     
    364369        return (a**(b+1.))*x**b*np.exp(-a*x)/scsp.gamma(b+1.)
    365370    else:
    366         return np.exp((b+1.)*np.log(a)-spsc.gammaln(b+1.)+b*np.log(x)-(a*x))
     371        return np.exp((b+1.)*np.log(a)-scsp.gammaln(b+1.)+b*np.log(x)-(a*x))
    367372           
    368 def LogNormalCume(x,pos,scale,shape):
     373def LogNormalCume(x,pos,args):
    369374    ''' Standard LogNormal cumulative distribution - numpy friendly on x axis
    370375    ref: http://www.itl.nist.gov/div898/handbook/index.htm 1.3.6.6.9
     
    375380    returns float: LogNormal cumulative distribution
    376381    '''
    377     return scsp.erf(np.log((x-pos)/shape)/(np.sqrt(2.)*scale)+1.)/2.
    378    
    379 def GaussCume(x,pos,scale,shape):
     382    scale,shape = args
     383    lredX = np.log((x-pos)/scale)
     384    return (scsp.erf((lredX/shape)/np.sqrt(2.))+1.)/2.
     385   
     386def GaussCume(x,pos,args):
    380387    ''' Standard Normal cumulative distribution - numpy friendly on x axis
    381388    param float x: independent axis (can be numpy array)
     
    385392    returns float: Normal cumulative distribution
    386393    '''
    387     return scsp.erf((x-pos)/(np.sqrt(2.)*scale)+1.)/2.
    388    
    389 def LSWCume(x,pos,scale,shape):
     394    scale = args[0]
     395    redX = (x-pos)/scale
     396    return (scsp.erf(redX/np.sqrt(2.))+1.)/2.
     397   
     398def LSWCume(x,pos,args=[]):
    390399    ''' Lifshitz-Slyozov-Wagner Ostwald ripening cumulative distribution - numpy friendly on x axis
    391400    param float x: independent axis (can be numpy array)
     
    395404    returns float: LSW cumulative distribution
    396405    '''
    397     nP = int(np.ceil(x/30+30))
    398     return []
    399    
    400 def SchulzZimmCume(x,pos,scale,shape):
     406    nP = 500
     407    xMin,xMax = [0.,2.*pos]
     408    X = np.linspace(xMin,xMax,nP,True)
     409    fxn = LSWDist(X,pos)
     410    mat = np.outer(np.ones(nP),fxn)
     411    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
     412    return np.interp(x,X,cume,0,1)
     413   
     414def SchulzZimmCume(x,pos,args):
    401415    ''' Schulz-Zimm cumulative distribution - numpy friendly on x axis
    402416    param float x: independent axis (can be numpy array)
     
    406420    returns float: Normal distribution
    407421    '''
     422    scale = args[0]
     423    nP = 500
     424    xMin = np.max([0.,pos-4.*scale])
     425    xMax = np.min([pos+4.*scale,1.e5])
     426    X = np.linspace(xMin,xMax,nP,True)
     427    fxn = LSWDist(X,pos)
     428    mat = np.outer(np.ones(nP),fxn)
     429    cume = np.sum(np.tril(mat),axis=1)/np.sum(fxn)
     430    return np.interp(x,X,cume,0,1)
     431   
    408432    return []
    409433   
     
    902926    print 'do model fit'
    903927   
    904 def ModelFxn(Q,G,parmDict,SQfxn,args=[]):
    905     return SQfxn(Q,args)
    906    
    907 def MakeDiamDist(DistName,nPoints,prec,pos,scale,shape):
    908    
    909     DistFuncs = {'Gauss':[GaussCume,0.02*scale,pos],
    910         'LogNormal':[LogNormalCume,0.02*scale,pos],
    911         'LSW':[LSWCume,0.3*pos,pos],
    912         'Schulz-Zimm':[SchulzZimmCume,
    913         4*np.sqrt(np.exp(shape**2)*(np.exp(shape**2)-1),pos+scale/np.exp(shape**2))]}
    914    
    915     step,mode = DistFuncs[DistName][1:]       
    916     minXP = 1.  # min cutoff at least 1A
    917     temp = mode
    918     while True:     #stepwise search for min cutoff
    919         temp -= step
    920         temp = min(temp,minXP)
    921         result = eval(DistFuncs[DistName][0](temp,pos,scale,shape))
    922         if (result <= prec) or (temp<=minXP):
    923             break
    924     startX = temp
    925     StartCumTgt = prec
    926     if startX == minXP:
    927         StartCumTgt = eval(DistFuncs[DistName][0](startX,pos,scale,shape))
    928     maxXP = 1.e15   #very high max cutoff (could be 10^5)
    929     temp = mode
    930     while True:     #stepwise search for max cutoff
    931         temp += step
    932         temp = max(temp,maxXP)
    933         result = eval(DistFuncs[DistName][0](temp,pos,scale,shape))
    934         if (result >= 1.-prec) or (temp>=maxXP):
    935             break
    936     endX = temp
    937    
    938     #special one for "Power Law" (future) - see p 27 of IRL2_NLSQCalc.ipf for code
    939    
    940     Diam = np.array([startX+i*(endX-startX/(3*nPoints-1)) for i in range(3*nPoints)])
    941     TCW = eval(DistFuncs[DistName][1](Diam,pos,scale,shape))
    942     CumeTgts = np.array([StartCumTgt+i*(1-prec-StartCumTgt)/(nPoints-1) for i in range(nPoints)])
    943     return np.interp(CumeTgts,Diam,TCW,0,0)
     928def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
     929    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
     930        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
     931        'Unified rod':[UniRodFF,UniRodVol],'Unified rod AR':[UniRodARFF,UniRodARVol],
     932        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol]}
     933    partData = sasdData['Particle']
     934    rhoMat = Substances['Substances'][partData['Matrix']['Name']].get('XAnom density',0.0)
     935    matFrac = partData['Matrix']['VolFrac'] #[value,flag]       
     936    Scale = Sample['Scale']     #[value,flag]
     937    Back = sasdData['Back']
     938    Q,Io,wt,Ic,Ib = Profile[:5]
     939    Qmin = Limits[1][0]
     940    Qmax = Limits[1][1]
     941    wtFactor = ProfDict['wtFactor']
     942    Ibeg = np.searchsorted(Q,Qmin)
     943    Ifin = np.searchsorted(Q,Qmax)
     944    Ib[:] = Back[0]
     945    Ic[:] = 0
     946    Ic[Ibeg:Ifin] += Back[0]
     947    Rbins = []
     948    Dist = []
     949    for level in partData['Levels']:
     950        controls = level['Controls']
     951        FFfxn = shapes[controls['FormFact']][0]
     952        Volfxn = shapes[controls['FormFact']][1]
     953        FFargs = []
     954        for item in ['Aspect ratio','Length','Thickness','Diameter',]:
     955            if item in controls['FFargs']:
     956                FFargs.append(controls['FFargs'][item][0])
     957        distFxn = controls['DistType']
     958        rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
     959        contrast = rho**2-rhoMat**2
     960        parmDict = level[controls['DistType']]
     961        rBins,dist = MakeDiamDist(controls['DistType'],controls['NumPoints'],controls['Cutoff'],parmDict)
     962        Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs)
     963        dist *= level[distFxn]['Volume'][0]
     964        Ic[Ibeg:Ifin] += np.dot(Gmat.T,dist)*100.
     965        Rbins.append(rBins)
     966        Dist.append(dist)
     967       
     968    return Rbins,Dist
     969             
     970   
     971def MakeDiamDist(DistName,nPoints,cutoff,parmDict):
     972   
     973    if 'LogNormal' in DistName:
     974        distFxn = 'LogNormalDist'
     975        cumeFxn = 'LogNormalCume'
     976        pos = parmDict['MinSize'][0]
     977        args = [parmDict['Mean'][0],parmDict['StdDev'][0]]
     978        step = 4.*np.sqrt(np.exp(parmDict['StdDev'][0]**2)*(np.exp(parmDict['StdDev'][0]**2)-1.))
     979        mode = parmDict['MinSize'][0]+parmDict['Mean'][0]/np.exp(parmDict['StdDev'][0]**2)
     980        minX = 1. #pos
     981        maxX = 1.e4 #np.min([mode+nPoints*step*40,1.e5])
     982    elif 'Gauss' in DistName:
     983        distFxn = 'GaussDist'
     984        cumeFxn = 'GaussCume'
     985        pos = parmDict['Mean'][0]
     986        args = [parmDict['StdDev'][0]]
     987        step = 0.02*parmDict['StdDev'][0]
     988        mode = parmDict['Mean'][0]
     989        minX = np.max([mode-4.*parmDict['StdDev'][0],1.])
     990        maxX = np.min([mode+4.*parmDict['StdDev'][0],1.e5])
     991    elif 'LSW' in DistName:
     992        distFxn = 'LSWDist'
     993        cumeFxn = 'LSWCume'
     994        pos = parmDict['Mean'][0]
     995        args = []
     996        minX,maxX = [0.,2.*pos]
     997    elif 'Schulz' in DistName:
     998        distFxn = 'SchulzZimmDist'
     999        cumeFxn = 'SchulzZimmCume'
     1000        pos = parmDict['Mean'][0]
     1001        args = [parmDict['StdDev'][0]]
     1002        minX = np.max([1.,pos-4.*parmDict['StdDev'][0]])
     1003        maxX = np.min([pos+4.*parmDict['StdDev'][0],1.e5])
     1004    nP = 500
     1005    X = np.linspace(minX,maxX,nP,True)
     1006    cume = eval(cumeFxn+'(X,pos,args)')
     1007    Imin = np.searchsorted(cume,cutoff)
     1008    Imax = np.min([np.searchsorted(cume,1-cutoff),nP-1])
     1009    startX =  X[Imin]
     1010    endX = X[Imax]
     1011#    Diam = np.linspace(startX,endX,3*nPoints,True)
     1012    Diam = np.logspace(0.,5.,500,True)
     1013    TCW = eval(cumeFxn+'(Diam,pos,args)')
     1014    CumeTgts = np.linspace(cutoff,(1.-cutoff),nPoints,True)
     1015    rBins = np.interp(CumeTgts,TCW,Diam,0,0)
     1016    return rBins,eval(distFxn+'(rBins,pos,args)')
    9441017
    9451018################################################################################
Note: See TracChangeset for help on using the changeset viewer.