Changeset 1309


Ignore:
Timestamp:
Apr 30, 2014 2:05:57 PM (9 years ago)
Author:
vondreele
Message:

Implement structure factors for non-dilute small angle systems

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIpwdGUI.py

    r1300 r1309  
    26852685            data['Particle']['Levels'].append({
    26862686                'Controls':{'FormFact':'Sphere','DistType':'LogNormal','Material':material,
    2687                     'FFargs':{},'NumPoints':50,'Cutoff':0.01,
     2687                    'FFargs':{},'SFargs':{},'NumPoints':50,'Cutoff':0.01,
    26882688                    'SlitSmear':[0.0,False],'StrFact':'Dilute'},    #last 2 not used - future?
    26892689                'LogNormal':{'Volume':[0.05,False],'Mean':[1000.,False],'StdDev':[0.5,False],'MinSize':[10.,False],},
    … …  
    29992999            'Unified tube':{'Length':[100.,False],'Thickness':[10.,False]},}
    30003000               
     3001        StructureFactors = {'Dilute':{},'Hard sphere':{'VolFr':[0.1,False],'Dist':[100.,False]},
     3002            'Sticky hard sphere':{'VolFr':[0.1,False],'Dist':[100.,False],'epis':[0.05,False],'Sticky':[0.2,False]},
     3003            'Square well':{'VolFr':[0.1,False],'Dist':[100.,False],'Depth':[0.1,False],'Width':[1.,False]},
     3004            'InterPrecipitate':{'VolFr':[0.1,False],'Dist':[100.,False]},}
     3005               
    30013006        ffDistChoices =  ['Sphere','Spheroid','Cylinder','Cylinder diam',
    30023007            'Cylinder AR','Unified sphere','Unified rod','Unified rod AR',
    30033008            'Unified disk','Unified tube',]
    30043009               
    3005         ffMonoChoices = ['Sphere','Spheroid','Cylinder','Cylinder AR',
    3006             ]
     3010        ffMonoChoices = ['Sphere','Spheroid','Cylinder','Cylinder AR',]
     3011       
     3012        sfChoices = ['Dilute','Hard sphere','Sticky hard sphere','Square well','InterPrecipitate',]
     3013           
     3014        slMult = 1000.
    30073015                 
    30083016        def OnValue(event):
    30093017            Obj = event.GetEventObject()
    3010             item,key,id,sldrObj = Indx[Obj.GetId()]
     3018            item,key,sldrObj = Indx[Obj.GetId()]
    30113019            try:
    30123020                value = float(Obj.GetValue())
    … …  
    30143022                    raise ValueError
    30153023            except ValueError:
    3016                 value = item[key][id]
    3017             item[key][id] = value
     3024                value = item[key][0]
     3025            item[key][0] = value
    30183026            Obj.SetValue('%.3g'%(value))
    3019             if key == 'P':
    3020                 sldrObj.SetValue(1000.*value)
     3027            if key in ['P','epis','Sticky','Depth','Width','VolFr','Dist']:
     3028                sldrObj.SetValue(slMult*value)
    30213029            else:
    3022                 sldrObj.SetRange(1000.*(np.log10(value)-2),1000.*(np.log10(value)+2))
    3023                 sldrObj.SetValue(1000.*np.log10(value))
     3030                logv = np.log10(value)
     3031                valMinMax = [logv-1,logv+1]
     3032                sldrObj.SetRange(slMult*valMinMax[0],slMult*valMinMax[1])
     3033                sldrObj.SetValue(slMult*logv)
    30243034            G2sasd.ModelFxn(Profile,ProfDict,Limits,Substances,Sample,data)
    30253035            RefreshPlots()
    … …  
    30323042                if 'FormFact' in key :
    30333043                    item['FFargs'] = FormFactors[Obj.GetValue()]
     3044                elif 'StrFact' in key:
     3045                    item['SFargs'] = StructureFactors[Obj.GetValue()]
    30343046                wx.CallAfter(UpdateModelsGrid,G2frame,data)
    30353047                G2sasd.ModelFxn(Profile,ProfDict,Limits,Substances,Sample,data)
    … …  
    30463058        def OnParmSlider(event):
    30473059            Obj = event.GetEventObject()
    3048             item,key,id,pvObj = Indx[Obj.GetId()]
     3060            item,key,pvObj = Indx[Obj.GetId()]
    30493061            slide = Obj.GetValue()
    3050             if key == 'P':
    3051                 value = float(slide/1000.)
     3062            if key in ['P','epis','Sticky','Depth','Width','VolFr','Dist']:
     3063                value = float(slide/slMult)
    30523064            else:
    3053                 value = 10.**float(slide/1000.)
    3054             item[key][id] = value
    3055             pvObj.SetValue('%.3g'%(item[key][id]))
     3065                value = 10.**float(slide/slMult)
     3066            item[key][0] = value
     3067            pvObj.SetValue('%.3g'%(item[key][0]))
    30563068            G2sasd.ModelFxn(Profile,ProfDict,Limits,Substances,Sample,data)
    30573069            RefreshPlots()
    3058            
    3059         def OnCenterSlider(event):
    3060             Obj = event.GetEventObject()
    3061             Obj.SetValue(False)
    3062             sldrObj,value = Indx[Obj.GetId()]
    3063             sldrObj.SetRange(1000*(value-1),1000*(value+1))
    3064             sldrObj.SetValue(1000*value)
    30653070           
    30663071        def SizeSizer():
    … …  
    30723077            Indx[distChoice.GetId()] = [level['Controls'],'DistType']
    30733078            distChoice.Bind(wx.EVT_COMBOBOX,OnSelect)
    3074             sizeSizer.Add(distChoice,0,WACV)
     3079            sizeSizer.Add(distChoice,0,WACV)    #put structure factor choices here
    30753080            if level['Controls']['DistType'] not in ['Bragg','Unified','Porod',]:
    30763081                sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Form Factor: '),0,WACV)
    3077                 if 'Mono' not in level['Controls']['DistType']:
    3078                     ffChoice = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['FormFact'],choices=ffDistChoices,
     3082                if 'Mono' in level['Controls']['DistType']:
     3083                    ffChoice = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['FormFact'],choices=ffMonoChoices,
    30793084                        style=wx.CB_READONLY|wx.CB_DROPDOWN)
    30803085                else:
    3081                     ffChoice = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['FormFact'],choices=ffMonoChoices,
     3086                    ffChoice = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['FormFact'],choices=ffDistChoices,
    30823087                        style=wx.CB_READONLY|wx.CB_DROPDOWN)
    30833088                Indx[ffChoice.GetId()] = [level['Controls'],'FormFact']
    30843089                ffChoice.Bind(wx.EVT_COMBOBOX,OnSelect)
    30853090                sizeSizer.Add(ffChoice,0,WACV)
     3091               
    30863092                sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Material: '),0,WACV)
    30873093                matSel = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['Material'],
    … …  
    30943100                sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Resonant X-ray contrast: '),0,WACV)
    30953101                sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label='  %.2f 10%scm%s'%(contrast,Pwr20,Pwrm4)),0,WACV)
    3096                 sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Num. radii: '),0,WACV)
    3097                 radii = ['25','50','75','100','200']
    3098                 nRadii = wx.ComboBox(G2frame.dataDisplay,value=str(level['Controls']['NumPoints']),choices=radii,
    3099                     style=wx.CB_READONLY|wx.CB_DROPDOWN)
    3100                 Indx[nRadii.GetId()] = [level['Controls'],'NumPoints']
    3101                 nRadii.Bind(wx.EVT_COMBOBOX,OnSelect)
    3102                 sizeSizer.Add(nRadii,0,WACV)
    3103                 sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' R dist. cutoff: '),0,WACV)
    3104                 rCutoff = G2gd.ValidatedTxtCtrl(G2frame.dataDisplay,level['Controls'],'Cutoff',
    3105                     min=0.001,max=0.1,typeHint=float)
    3106                 sizeSizer.Add(rCutoff,0,WACV)
     3102                if 'Mono' not in level['Controls']['DistType']:
     3103                    sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Num. radii: '),0,WACV)
     3104                    radii = ['25','50','75','100','200']
     3105                    nRadii = wx.ComboBox(G2frame.dataDisplay,value=str(level['Controls']['NumPoints']),choices=radii,
     3106                        style=wx.CB_READONLY|wx.CB_DROPDOWN)
     3107                    Indx[nRadii.GetId()] = [level['Controls'],'NumPoints']
     3108                    nRadii.Bind(wx.EVT_COMBOBOX,OnSelect)
     3109                    sizeSizer.Add(nRadii,0,WACV)
     3110                    sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' R dist. cutoff: '),0,WACV)
     3111                    rCutoff = G2gd.ValidatedTxtCtrl(G2frame.dataDisplay,level['Controls'],'Cutoff',
     3112                        min=0.001,max=0.1,typeHint=float)
     3113                    sizeSizer.Add(rCutoff,0,WACV)
    31073114            elif level['Controls']['DistType']  in ['Unified',]:
    31083115                Parms = level['Unified']
    … …  
    31103117                sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Estimated Dist B: %12.4g'%(Best)),0,WACV)
    31113118            return sizeSizer
     3119           
     3120        def ParmSizer():
     3121            parmSizer = wx.FlexGridSizer(1,3,5,5)
     3122            parmSizer.AddGrowableCol(2,1)
     3123            parmSizer.SetFlexibleDirection(wx.HORIZONTAL)
     3124            Parms = level[level['Controls']['DistType']]
     3125            FFargs = level['Controls']['FFargs']
     3126            SFargs = level['Controls'].get('SFargs',{})
     3127            parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
     3128                'PkInt','PkPos','PkSig','PkGam',]
     3129            for parm in parmOrder:
     3130                if parm in Parms:
     3131                    if parm == 'MinSize':
     3132                        parmSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Dist '+parm),0,wx.ALIGN_CENTER)
     3133                    else:
     3134                        parmVar = wx.CheckBox(G2frame.dataDisplay,label='Refine? Dist '+parm)
     3135                        parmVar.SetValue(Parms[parm][1])
     3136                        parmVar.Bind(wx.EVT_CHECKBOX, OnSelect)
     3137                        parmSizer.Add(parmVar,0,WACV)
     3138                        Indx[parmVar.GetId()] = [Parms[parm],1]
     3139                    parmValue = wx.TextCtrl(G2frame.dataDisplay,value='%.3g'%(Parms[parm][0]),
     3140                        style=wx.TE_PROCESS_ENTER)
     3141                    parmValue.Bind(wx.EVT_TEXT_ENTER,OnValue)       
     3142                    parmValue.Bind(wx.EVT_KILL_FOCUS,OnValue)
     3143                    parmSizer.Add(parmValue,0,WACV)
     3144                    if parm == 'P':
     3145                        value = Parms[parm][0]
     3146                        valMinMax = [0.1,4.2]
     3147                    else:
     3148                        value = np.log10(Parms[parm][0])
     3149                        valMinMax = [value-1,value+1]
     3150                    parmSldr = wx.Slider(G2frame.dataDisplay,minValue=slMult*valMinMax[0],
     3151                        maxValue=slMult*valMinMax[1],value=slMult*value)
     3152                    Indx[parmValue.GetId()] = [Parms,parm,parmSldr]
     3153                    Indx[parmSldr.GetId()] = [Parms,parm,parmValue]
     3154                    parmSldr.Bind(wx.EVT_SLIDER,OnParmSlider)
     3155                    parmSizer.Add(parmSldr,1,wx.EXPAND)
     3156            if level['Controls']['DistType'] not in ['Bragg']:
     3157                parmOrder = ['Aspect ratio','Length','Diameter','Thickness','VolFr','Dist','epis','Sticky','Depth','Width']
     3158                fTypes = ['FF ','SF ']
     3159                for iarg,Args in enumerate([FFargs,SFargs]):
     3160                    for parm in parmOrder:
     3161                        if parm in Args:
     3162                            parmVar = wx.CheckBox(G2frame.dataDisplay,label='Refine? '+fTypes[iarg]+parm)
     3163                            parmVar.SetValue(Args[parm][1])
     3164                            Indx[parmVar.GetId()] = [Args[parm],1]
     3165                            parmVar.Bind(wx.EVT_CHECKBOX, OnSelect)
     3166                            parmSizer.Add(parmVar,0,WACV)
     3167                            parmValue = wx.TextCtrl(G2frame.dataDisplay,value='%.3g'%(Args[parm][0]),
     3168                                style=wx.TE_PROCESS_ENTER)
     3169                            parmValue.Bind(wx.EVT_TEXT_ENTER,OnValue)       
     3170                            parmValue.Bind(wx.EVT_KILL_FOCUS,OnValue)
     3171                            parmSizer.Add(parmValue,0,WACV)
     3172                            value = Args[parm][0]
     3173                            if parm == 'epis':
     3174                                valMinMax = [0,.1]
     3175                            elif parm in ['Sticky','Width',]:
     3176                                valMinMax = [0,1.]
     3177                            elif parm == 'Depth':
     3178                                valMinMax = [-2.,2.]
     3179                            elif parm == 'Dist':
     3180                                valMinMax = [100.,1000.]
     3181                            elif parm == 'VolFr':
     3182                                valMinMax = [1.e-4,1.]
     3183                            else:
     3184                                value = np.log10(Args[parm][0])
     3185                                valMinMax = [value-1,value+1]
     3186                            parmSldr = wx.Slider(G2frame.dataDisplay,minValue=slMult*valMinMax[0],
     3187                                maxValue=slMult*valMinMax[1],value=slMult*value)
     3188                            Indx[parmVar.GetId()] = [Args[parm],1]
     3189                            Indx[parmValue.GetId()] = [Args,parm,parmSldr]
     3190                            Indx[parmSldr.GetId()] = [Args,parm,parmValue]
     3191                            parmSldr.Bind(wx.EVT_SLIDER,OnParmSlider)
     3192                            parmSizer.Add(parmSldr,1,wx.EXPAND)
     3193            return parmSizer               
    31123194           
    31133195        Indx = {}
    … …  
    31423224            partSizer.Add(topLevel,0)
    31433225            partSizer.Add(SizeSizer())
    3144            
    3145             lvlSizer = wx.BoxSizer(wx.HORIZONTAL)
    3146             parmSizer = wx.FlexGridSizer(1,4,5,5)
    3147             parmSizer.AddGrowableCol(2,1)
    3148             parmSizer.SetFlexibleDirection(wx.HORIZONTAL)
    3149             Parms = level[level['Controls']['DistType']]
    3150             FFargs = level['Controls']['FFargs']
    3151             parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
    3152                 'PkInt','PkPos','PkSig','PkGam',]
    3153             for parm in parmOrder:
    3154                 if parm in Parms:
    3155                     if parm == 'MinSize':
    3156                         parmSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Dist '+parm),0,wx.ALIGN_CENTER)
    3157                     else:
    3158                         parmVar = wx.CheckBox(G2frame.dataDisplay,label='Refine? Dist '+parm)
    3159                         parmVar.SetValue(Parms[parm][1])
    3160                         parmVar.Bind(wx.EVT_CHECKBOX, OnSelect)
    3161                         parmSizer.Add(parmVar,0,WACV)
    3162                         Indx[parmVar.GetId()] = [Parms[parm],1]
    3163                     parmValue = wx.TextCtrl(G2frame.dataDisplay,value='%.3g'%(Parms[parm][0]),
    3164                         style=wx.TE_PROCESS_ENTER)
    3165                     parmValue.Bind(wx.EVT_TEXT_ENTER,OnValue)       
    3166                     parmValue.Bind(wx.EVT_KILL_FOCUS,OnValue)
    3167                     parmSizer.Add(parmValue,0,WACV)
    3168                     if parm == 'P':
    3169                         value = Parms[parm][0]
    3170                         valMinMax = [0.1,4.2]
    3171                     else:
    3172                         value = np.log10(Parms[parm][0])
    3173                         valMinMax = [value-1,value+1]
    3174                     parmSldr = wx.Slider(G2frame.dataDisplay,minValue=1000.*valMinMax[0],
    3175                         maxValue=1000.*valMinMax[1],value=1000.*value)
    3176                     Indx[parmValue.GetId()] = [Parms,parm,0,parmSldr]
    3177                     Indx[parmSldr.GetId()] = [Parms,parm,0,parmValue]
    3178                     parmSldr.Bind(wx.EVT_SLIDER,OnParmSlider)
    3179                     parmSizer.Add(parmSldr,1,wx.EXPAND)
    3180                     if parm == 'P':
    3181                         parmSizer.Add((5,5),)
    3182                     else:
    3183                         center = wx.CheckBox(G2frame.dataDisplay,label='Center? ')
    3184                         Indx[center.GetId()] = [parmSldr,value]
    3185                         center.Bind(wx.EVT_CHECKBOX,OnCenterSlider)
    3186                         parmSizer.Add(center,0,WACV)
    3187             if level['Controls']['DistType'] not in ['Bragg']:
    3188                 for parm in list(FFargs):
    3189                     parmVar = wx.CheckBox(G2frame.dataDisplay,label='Refine? FF '+parm)
    3190                     parmVar.SetValue(FFargs[parm][1])
    3191                     Indx[parmVar.GetId()] = [FFargs[parm],1]
    3192                     parmVar.Bind(wx.EVT_CHECKBOX, OnSelect)
    3193                     parmSizer.Add(parmVar,0,WACV)
    3194                     parmValue = wx.TextCtrl(G2frame.dataDisplay,value='%.3g'%(FFargs[parm][0]),
    3195                         style=wx.TE_PROCESS_ENTER)
    3196                     parmValue.Bind(wx.EVT_TEXT_ENTER,OnValue)       
    3197                     parmValue.Bind(wx.EVT_KILL_FOCUS,OnValue)
    3198                     parmSizer.Add(parmValue,0,WACV)
    3199                     value = np.log10(FFargs[parm][0])
    3200                     valMinMax = [value-1,value+1]
    3201                     parmSldr = wx.Slider(G2frame.dataDisplay,minValue=1000.*valMinMax[0],
    3202                         maxValue=1000.*valMinMax[1],value=1000.*value)
    3203                     Indx[parmVar.GetId()] = [FFargs[parm],1]
    3204                     Indx[parmValue.GetId()] = [FFargs[parm],0,parmSldr]
    3205                     Indx[parmSldr.GetId()] = [FFargs,parm,0,parmValue]
    3206                     parmSldr.Bind(wx.EVT_SLIDER,OnParmSlider)
    3207                     parmSizer.Add(parmSldr,1,wx.EXPAND)
    3208                     center = wx.CheckBox(G2frame.dataDisplay,label='Center? ')
    3209                     Indx[center.GetId()] = [parmSldr,value]
    3210                     center.Bind(wx.EVT_CHECKBOX,OnCenterSlider)
    3211                     parmSizer.Add(center,0,WACV)
    3212             lvlSizer.Add(parmSizer,1,wx.EXPAND)
    3213             partSizer.Add(lvlSizer,1,wx.EXPAND)
     3226            if level['Controls']['DistType'] not in ['Bragg','Unified','Porod',]:
     3227                topLevel.Add(wx.StaticText(G2frame.dataDisplay,label=' Structure factor: '),0,WACV)
     3228                strfctr = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['StrFact'],
     3229                    choices=sfChoices,style=wx.CB_READONLY|wx.CB_DROPDOWN)
     3230                Indx[strfctr.GetId()] = [level['Controls'],'StrFact']
     3231                strfctr.Bind(wx.EVT_COMBOBOX,OnSelect)
     3232                topLevel.Add(strfctr,0,WACV)
     3233            partSizer.Add(ParmSizer(),1,wx.EXPAND)
    32143234        return partSizer
    32153235       
    … …  
    33193339    G2frame.dataDisplay.SetSize(Size)
    33203340    G2frame.dataFrame.setSizePosLeft(Size)   
    3321        
    33223341   
    33233342################################################################################

  • trunk/GSASIIsasd.py

    r1302 r1309  
    435435    return []
    436436   
     437################################################################################
     438#### Structure factors for condensed systems
     439################################################################################
     440
     441def DiluteSF(Q,args=[]):
     442    ''' Default: no structure factor correction for dilute system
     443    '''
     444    return np.ones_like(Q)  #or return 1.0
     445
     446def HardSpheresSF(Q,args):
     447    ''' Computes structure factor for not dilute monodisperse hard spheres
     448    Refs.: PERCUS,YEVICK PHYS. REV. 110 1 (1958),THIELE J. CHEM PHYS. 39 474 (1968),
     449    WERTHEIM  PHYS. REV. LETT. 47 1462 (1981)
     450   
     451    param float Q: Q value array (A-1)
     452    param array args: [float R, float VolFrac]: interparticle distance & volume fraction
     453    returns numpy array S(Q)
     454    '''
     455   
     456    R,VolFr = args
     457    denom = (1.-VolFr)**4
     458    num = (1.+2.*VolFr)**2
     459    alp = num/denom
     460    bet = -6.*VolFr*(1.+VolFr/2.)**2/denom
     461    gamm = 0.5*VolFr*num/denom       
     462    A = 2.*Q*R
     463    A2 = A**2
     464    A3 = A**3
     465    A4 = A**4
     466    Rca = np.cos(A)
     467    Rsa = np.sin(A)
     468    calp = alp*(Rsa/A2-Rca/A)
     469    cbet = bet*(2.*Rsa/A2-(A2-2.)*Rca/A3-2./A3)
     470    cgam = gamm*(-Rca/A+(4./A)*((3.*A2-6.)*Rca/A4+(A2-6.)*Rsa/A3+6./A4))
     471    pfac = -24.*VolFr/A
     472    C = pfac*(calp+cbet+cgam)
     473    return 1./(1.-C)
     474       
     475def SquareWellSF(Q,args):
     476    ''' Computes structure factor for not dilute monodisperse hard sphere with a
     477    square well potential interaction.
     478    Refs.: SHARMA,SHARMA, PHYSICA 89A,(1977),212
     479   
     480    param float Q: Q value array (A-1)
     481    param array args: [float R, float VolFrac, float depth, float width]:
     482        interparticle distance, volume fraction (<0.08), well depth (e/kT<1.5kT),
     483        well width
     484    returns numpy array S(Q)
     485    well depth > 0 attractive & values outside above limits nonphysical cf.
     486    Monte Carlo simulations
     487    '''
     488    R,VolFr,Depth,Width = args
     489    eta = VolFr
     490    eta2 = eta*eta
     491    eta3 = eta*eta2
     492    eta4 = eta*eta3       
     493    etam1 = 1. - eta
     494    etam14 = etam1**4
     495    alp = (  ( (1. + 2.*eta)**2 ) + eta3*( eta-4.0 )  )/etam14
     496    bet = -(eta/3.0) * ( 18. + 20.*eta - 12.*eta2 + eta4 )/etam14
     497    gam = 0.5*eta*( (1. + 2.*eta)**2 + eta3*(eta-4.) )/etam14
     498
     499    SK = 2.*Q*R
     500    SK2 = SK*SK
     501    SK3 = SK*SK2
     502    SK4 = SK3*SK
     503    T1 = alp * SK3 * ( np.sin(SK) - SK * np.cos(SK) )
     504    T2 = bet * SK2 * ( 2.*SK*np.sin(SK) - (SK2-2.)*np.cos(SK) - 2.0 )
     505    T3 =   ( 4.0*SK3 - 24.*SK ) * np.sin(SK) 
     506    T3 = T3 - ( SK4 - 12.0*SK2 + 24.0 )*np.cos(SK) + 24.0   
     507    T3 = gam*T3
     508    T4 = -Depth*SK3*(np.sin(Width*SK) - Width*SK*np.cos(Width*SK)+ SK*np.cos(SK) - np.sin(SK) )
     509    CK =  -24.0*eta*( T1 + T2 + T3 + T4 )/SK3/SK3
     510    return 1./(1.-CK)
     511   
     512def StickyHardSpheresSF(Q,args):
     513    ''' Computes structure factor for not dilute monodisperse hard spheres
     514    Refs.: PERCUS,YEVICK PHYS. REV. 110 1 (1958),THIELE J. CHEM PHYS. 39 474 (1968),
     515    WERTHEIM  PHYS. REV. LETT. 47 1462 (1981)
     516   
     517    param float Q: Q value array (A-1)
     518    param array args: [float R, float VolFrac]: sphere radius & volume fraction
     519    returns numpy array S(Q)
     520    '''
     521    R,VolFr,epis,tau = args
     522#//      Input (fitting) variables are:
     523#       //radius = w[0]
     524#       //volume fraction = w[1]
     525#       //epsilon (perurbation param) = w[2]
     526#       //tau (stickiness) = w[3]
     527#       Variable rad,phi,eps,tau,eta
     528#       Variable sig,aa,etam1,qa,qb,qc,radic
     529#       Variable lam,lam2,test,mu,alpha,BetaVar
     530#       Variable qv,kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq
     531#       rad = w[0]
     532#       phi = w[1]
     533#       eps = w[2]
     534#       tau = w[3]
     535#       
     536#       eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps)
     537#       
     538#       sig = 2.0 * rad
     539#       aa = sig/(1.0 - eps)
     540#       etam1 = 1.0 - eta
     541#//C
     542#//C  SOLVE QUADRATIC FOR LAMBDA
     543#//C
     544#       qa = eta/12.0
     545#       qb = -1.0*(tau + eta/(etam1))
     546#       qc = (1.0 + eta/2.0)/(etam1*etam1)
     547#       radic = qb*qb - 4.0*qa*qc
     548#       if(radic<0)
     549#               if(x>0.01 && x<0.015)
     550#                       Print "Lambda unphysical - both roots imaginary"
     551#               endif
     552#               return(-1)
     553#       endif
     554#//C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
     555#       lam = (-1.0*qb-sqrt(radic))/(2.0*qa)
     556#       lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa)
     557#       if(lam2<lam)
     558#               lam = lam2
     559#       endif
     560#       test = 1.0 + 2.0*eta
     561#       mu = lam*eta*etam1
     562#       if(mu>test)
     563#               if(x>0.01 && x<0.015)
     564#                Print "Lambda unphysical mu>test"
     565#               endif
     566#               return(-1)
     567#       endif
     568#       alpha = (1.0 + 2.0*eta - mu)/(etam1*etam1)
     569#       BetaVar = (mu - 3.0*eta)/(2.0*etam1*etam1)
     570#       
     571#//C
     572#//C   CALCULATE THE STRUCTURE FACTOR
     573#//C
     574#
     575#       qv = x
     576#       kk = qv*aa
     577#       k2 = kk*kk
     578#       k3 = kk*k2
     579#       ds = sin(kk)
     580#       dc = cos(kk)
     581#       aq1 = ((ds - kk*dc)*alpha)/k3
     582#       aq2 = (BetaVar*(1.0-dc))/k2
     583#       aq3 = (lam*ds)/(12.0*kk)
     584#       aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
     585#//
     586#       bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
     587#       bq2 = BetaVar*(1.0/kk - ds/k2)
     588#       bq3 = (lam/12.0)*((1.0 - dc)/kk)
     589#       bq = 12.0*eta*(bq1+bq2-bq3)
     590#//
     591#       sq = 1.0/(aq*aq +bq*bq)
     592#
     593#       Return (sq)
     594    pass
     595   
     596#def HayterPenfoldSF(Q,args): #big & ugly function - do later (if ever)
     597#    pass
     598   
     599def InterPrecipitateSF(Q,args):
     600    ''' Computes structure factor for precipitates in a matrix
     601    Refs.: E-Wen Huang, Peter K. Liaw, Lionel Porcar, Yun Liu, Yee-Lang Liu,
     602    Ji-Jung Kai, and Wei-Ren Chen,APPLIED PHYSICS LETTERS 93, 161904 (2008)
     603    R. Giordano, A. Grasso, and J. Teixeira, Phys. Rev. A 43, 6894    1991   
     604    param float Q: Q value array (A-1)
     605    param array args: [float R, float VolFr]: "radius" & volume fraction
     606    returns numpy array S(Q)
     607    '''
     608    R,VolFr = args
     609    QV2 = Q**2*VolFr**2
     610    top = 1 - np.exp(-QV2/4)*np.cos(2.*Q*R)
     611    bot = 1-2*np.exp(-QV2/4)*np.cos(2.*Q*R) + np.exp(-QV2/2)
     612    return 2*(top/bot) - 1
    437613         
    438614################################################################################
    … …  
    9211097        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
    9221098        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
     1099           
     1100    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
     1101            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
     1102               
     1103    parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
     1104        'PkInt','PkPos','PkSig','PkGam',]
     1105       
     1106    FFparmOrder = ['Aspect ratio','Length','Diameter','Thickness']
     1107   
     1108    SFparmOrder = ['Dist','VolFr','epis','Sticky','Depth','Width']
    9231109
    9241110    def GetModelParms():
    925         parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',
    926             'PkInt','PkPos','PkSig','PkGam',]
    9271111        parmDict = {'Scale':Sample['Scale'][0]}
    9281112        varyList = []
    … …  
    9471131                parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0]
    9481132                parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1]
     1133                parmDict[cid+'StrFact'] = sfxns[controls['StrFact']]
    9491134                parmDict[cid+'XAnom density'] = Substances['Substances'][controls['Material']].get('XAnom density',0.0)
    950                 for item in ['Aspect ratio','Length','Thickness','Diameter',]:
     1135                for item in FFparmOrder:
    9511136                    if item in controls['FFargs']:
    9521137                        parmDict[cid+item] = controls['FFargs'][item][0]
    … …  
    9541139                            varyList.append(cid+item)
    9551140                            values.append(controls['FFargs'][item][0])
     1141                for item in SFparmOrder:
     1142                    if item in controls.get('SFargs',{}):
     1143                        parmDict[cid+item] = controls['SFargs'][item][0]
     1144                        if controls['SFargs'][item][1]:
     1145                            varyList.append(cid+item)
     1146                            values.append(controls['SFargs'][item][0])
    9561147            distDict = controls['DistType']
    9571148            for item in parmOrder:
    … …  
    9701161        partData = Model['Particle']
    9711162        for i,level in enumerate(partData['Levels']):
    972             Type = level['Controls']['DistType']
    973             print ' Component %d: Type: %s:'%(i,Type)
    974             cid = str(i)+':'
    9751163            controls = level['Controls']
    9761164            Type = controls['DistType']
    9771165            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
    978                 for item in ['Aspect ratio','Length','Thickness','Diameter',]:
     1166                print ' Component %d: Type: %s: Structure Factor: %s'%(i,Type,controls['StrFact'])               
     1167            else:
     1168                print ' Component %d: Type: %s: '%(i,Type,)
     1169            cid = str(i)+':'
     1170            if Type in ['LogNormal','Gaussian','LSW','Schulz-Zimm','Monodisperse']:
     1171                for item in FFparmOrder:
    9791172                    if cid+item in varyList:
    9801173                        controls['FFargs'][item][0] = parmDict[cid+item]
     1174                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
     1175                for item in SFparmOrder:
     1176                    if cid+item in varyList:
     1177                        controls['SFargs'][item][0] = parmDict[cid+item]
    9811178                        print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item])
    9821179            distDict = controls['DistType']
    … …  
    9991196                FFfxn = parmDict[cid+'FormFact']
    10001197                Volfxn = parmDict[cid+'FFVolume']
     1198                SFfxn = parmDict[cid+'StrFact']
    10011199                FFargs = []
    1002                 for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
     1200                SFargs = []
     1201                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter']:
    10031202                    if item in parmDict:
    10041203                        FFargs.append(parmDict[item])
     1204                for item in [cid+'Dist',cid+'VolFr',cid+'epis',cid+'Sticky',cid+'Depth',cid+'Width']:
     1205                    if item in parmDict:
     1206                        SFargs.append(parmDict[item])
    10051207                distDict = {}
    10061208                for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]:
    … …  
    10121214                Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T
    10131215                dist *= parmDict[cid+'Volume']
    1014                 Ic += np.dot(Gmat,dist)
     1216                Ic += np.dot(Gmat,dist)*SFfxn(Q,args=SFargs)
    10151217            elif 'Unified' in Type:
    10161218                Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \
    … …  
    10271229                FFfxn = parmDict[cid+'FormFact']
    10281230                Volfxn = parmDict[cid+'FFVolume']
     1231                SFfxn = parmDict[cid+'StrFact']
    10291232                FFargs = []
     1233                SFargs = []
    10301234                for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]:
    10311235                    if item in parmDict:
    … …  
    10351239                R = parmDict[cid+'Radius']
    10361240                Gmat = G_matrix(Q,R,contrast,FFfxn,Volfxn,FFargs)             
    1037                 Ic += Gmat[0]*parmDict[cid+'Volume']
     1241                Ic += Gmat[0]*parmDict[cid+'Volume']*SFfxn(Q,args=SFargs)
    10381242            elif 'Bragg' in distFxn:
    10391243                Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'],
    … …  
    10841288   
    10851289def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData):
     1290   
    10861291    shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol],
    10871292        'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol],
    … …  
    10891294        'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol],
    10901295        'Unified tube':[UniTubeFF,UniTubeVol],'Cylinder diam':[CylinderDFF,CylinderDVol]}
     1296    sfxns = {'Dilute':DiluteSF,'Hard sphere':HardSpheresSF,'Square well':SquareWellSF,
     1297            'Sticky hard sphere':StickyHardSpheresSF,'InterPrecipitate':InterPrecipitateSF,}
     1298
    10911299#    pdb.set_trace()
    10921300    partData = sasdData['Particle']
    … …  
    11091317        distFxn = controls['DistType']
    11101318        if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']:
     1319            parmDict = level[controls['DistType']]
    11111320            FFfxn = shapes[controls['FormFact']][0]
    11121321            Volfxn = shapes[controls['FormFact']][1]
     1322            SFfxn = sfxns[controls['StrFact']]
    11131323            FFargs = []
     1324            SFargs = []
     1325            for item in ['Dist','VolFr','epis','Sticky','Depth','Width',]:
     1326                if item in controls.get('SFargs',{}):
     1327                    SFargs.append(controls['SFargs'][item][0])
    11141328            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
    11151329                if item in controls['FFargs']:
    … …  
    11171331            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
    11181332            contrast = rho**2-rhoMat**2
    1119             parmDict = level[controls['DistType']]
    11201333            distDict = {}
    11211334            for item in parmDict:
    … …  
    11241337            Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T
    11251338            dist *= level[distFxn]['Volume'][0]
    1126             Ic[Ibeg:Ifin] += np.dot(Gmat,dist)
     1339            Ic[Ibeg:Ifin] += np.dot(Gmat,dist)*SFfxn(Q[Ibeg:Ifin],args=SFargs)
    11271340            Rbins.append(rBins)
    11281341            Dist.append(dist/(4.*dBins))
    … …  
    11451358            Dist.append([])
    11461359        elif 'Mono' in distFxn:
     1360            parmDict = level[controls['DistType']]
     1361            R = level[controls['DistType']]['Radius'][0]
    11471362            FFfxn = shapes[controls['FormFact']][0]
    11481363            Volfxn = shapes[controls['FormFact']][1]
     1364            SFfxn = sfxns[controls['StrFact']]
    11491365            FFargs = []
     1366            SFargs = [parmDict['Volume'][0],R]
     1367            for item in ['epis','Sticky','Depth','Width',]:
     1368                if item in controls.get('SFargs',{}):
     1369                    SFargs.append(controls['SFargs'][item][0])
    11501370            for item in ['Aspect ratio','Length','Thickness','Diameter',]:
    11511371                if item in controls['FFargs']:
    … …  
    11531373            rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0)
    11541374            contrast = rho**2-rhoMat**2
    1155             R = level[controls['DistType']]['Radius'][0]
    11561375            Gmat = G_matrix(Q[Ibeg:Ifin],R,contrast,FFfxn,Volfxn,FFargs)             
    1157             Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]
     1376            Ic[Ibeg:Ifin] += Gmat[0]*level[distFxn]['Volume'][0]*SFfxn(Q[Ibeg:Ifin],args=SFargs)
    11581377            Rbins.append([])
    11591378            Dist.append([])
Note: See TracChangeset for help on using the changeset viewer.