Changeset 1309
- Timestamp:
- Apr 30, 2014 2:05:57 PM (9 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIpwdGUI.py
r1300 r1309 2685 2685 data['Particle']['Levels'].append({ 2686 2686 'Controls':{'FormFact':'Sphere','DistType':'LogNormal','Material':material, 2687 'FFargs':{},' NumPoints':50,'Cutoff':0.01,2687 'FFargs':{},'SFargs':{},'NumPoints':50,'Cutoff':0.01, 2688 2688 'SlitSmear':[0.0,False],'StrFact':'Dilute'}, #last 2 not used - future? 2689 2689 'LogNormal':{'Volume':[0.05,False],'Mean':[1000.,False],'StdDev':[0.5,False],'MinSize':[10.,False],}, … … 2999 2999 'Unified tube':{'Length':[100.,False],'Thickness':[10.,False]},} 3000 3000 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 3001 3006 ffDistChoices = ['Sphere','Spheroid','Cylinder','Cylinder diam', 3002 3007 'Cylinder AR','Unified sphere','Unified rod','Unified rod AR', 3003 3008 'Unified disk','Unified tube',] 3004 3009 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. 3007 3015 3008 3016 def OnValue(event): 3009 3017 Obj = event.GetEventObject() 3010 item,key, id,sldrObj = Indx[Obj.GetId()]3018 item,key,sldrObj = Indx[Obj.GetId()] 3011 3019 try: 3012 3020 value = float(Obj.GetValue()) … … 3014 3022 raise ValueError 3015 3023 except ValueError: 3016 value = item[key][ id]3017 item[key][ id] = value3024 value = item[key][0] 3025 item[key][0] = value 3018 3026 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) 3021 3029 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) 3024 3034 G2sasd.ModelFxn(Profile,ProfDict,Limits,Substances,Sample,data) 3025 3035 RefreshPlots() … … 3032 3042 if 'FormFact' in key : 3033 3043 item['FFargs'] = FormFactors[Obj.GetValue()] 3044 elif 'StrFact' in key: 3045 item['SFargs'] = StructureFactors[Obj.GetValue()] 3034 3046 wx.CallAfter(UpdateModelsGrid,G2frame,data) 3035 3047 G2sasd.ModelFxn(Profile,ProfDict,Limits,Substances,Sample,data) … … 3046 3058 def OnParmSlider(event): 3047 3059 Obj = event.GetEventObject() 3048 item,key, id,pvObj = Indx[Obj.GetId()]3060 item,key,pvObj = Indx[Obj.GetId()] 3049 3061 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) 3052 3064 else: 3053 value = 10.**float(slide/ 1000.)3054 item[key][ id] = value3055 pvObj.SetValue('%.3g'%(item[key][ id]))3065 value = 10.**float(slide/slMult) 3066 item[key][0] = value 3067 pvObj.SetValue('%.3g'%(item[key][0])) 3056 3068 G2sasd.ModelFxn(Profile,ProfDict,Limits,Substances,Sample,data) 3057 3069 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)3065 3070 3066 3071 def SizeSizer(): … … 3072 3077 Indx[distChoice.GetId()] = [level['Controls'],'DistType'] 3073 3078 distChoice.Bind(wx.EVT_COMBOBOX,OnSelect) 3074 sizeSizer.Add(distChoice,0,WACV) 3079 sizeSizer.Add(distChoice,0,WACV) #put structure factor choices here 3075 3080 if level['Controls']['DistType'] not in ['Bragg','Unified','Porod',]: 3076 3081 sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Form Factor: '),0,WACV) 3077 if 'Mono' notin level['Controls']['DistType']:3078 ffChoice = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['FormFact'],choices=ff DistChoices,3082 if 'Mono' in level['Controls']['DistType']: 3083 ffChoice = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['FormFact'],choices=ffMonoChoices, 3079 3084 style=wx.CB_READONLY|wx.CB_DROPDOWN) 3080 3085 else: 3081 ffChoice = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['FormFact'],choices=ff MonoChoices,3086 ffChoice = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['FormFact'],choices=ffDistChoices, 3082 3087 style=wx.CB_READONLY|wx.CB_DROPDOWN) 3083 3088 Indx[ffChoice.GetId()] = [level['Controls'],'FormFact'] 3084 3089 ffChoice.Bind(wx.EVT_COMBOBOX,OnSelect) 3085 3090 sizeSizer.Add(ffChoice,0,WACV) 3091 3086 3092 sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Material: '),0,WACV) 3087 3093 matSel = wx.ComboBox(G2frame.dataDisplay,value=level['Controls']['Material'], … … 3094 3100 sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Resonant X-ray contrast: '),0,WACV) 3095 3101 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) 3107 3114 elif level['Controls']['DistType'] in ['Unified',]: 3108 3115 Parms = level['Unified'] … … 3110 3117 sizeSizer.Add(wx.StaticText(G2frame.dataDisplay,label=' Estimated Dist B: %12.4g'%(Best)),0,WACV) 3111 3118 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 3112 3194 3113 3195 Indx = {} … … 3142 3224 partSizer.Add(topLevel,0) 3143 3225 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) 3214 3234 return partSizer 3215 3235 … … 3319 3339 G2frame.dataDisplay.SetSize(Size) 3320 3340 G2frame.dataFrame.setSizePosLeft(Size) 3321 3322 3341 3323 3342 ################################################################################ -
trunk/GSASIIsasd.py
r1302 r1309 435 435 return [] 436 436 437 ################################################################################ 438 #### Structure factors for condensed systems 439 ################################################################################ 440 441 def DiluteSF(Q,args=[]): 442 ''' Default: no structure factor correction for dilute system 443 ''' 444 return np.ones_like(Q) #or return 1.0 445 446 def 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 475 def 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 512 def 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 599 def 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 437 613 438 614 ################################################################################ … … 921 1097 'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol], 922 1098 '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'] 923 1109 924 1110 def GetModelParms(): 925 parmOrder = ['Volume','Radius','Mean','StdDev','MinSize','G','Rg','B','P','Cutoff',926 'PkInt','PkPos','PkSig','PkGam',]927 1111 parmDict = {'Scale':Sample['Scale'][0]} 928 1112 varyList = [] … … 947 1131 parmDict[cid+'FormFact'] = shapes[controls['FormFact']][0] 948 1132 parmDict[cid+'FFVolume'] = shapes[controls['FormFact']][1] 1133 parmDict[cid+'StrFact'] = sfxns[controls['StrFact']] 949 1134 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: 951 1136 if item in controls['FFargs']: 952 1137 parmDict[cid+item] = controls['FFargs'][item][0] … … 954 1139 varyList.append(cid+item) 955 1140 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]) 956 1147 distDict = controls['DistType'] 957 1148 for item in parmOrder: … … 970 1161 partData = Model['Particle'] 971 1162 for i,level in enumerate(partData['Levels']): 972 Type = level['Controls']['DistType']973 print ' Component %d: Type: %s:'%(i,Type)974 cid = str(i)+':'975 1163 controls = level['Controls'] 976 1164 Type = controls['DistType'] 977 1165 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: 979 1172 if cid+item in varyList: 980 1173 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] 981 1178 print ' %15s: %15.4g esd: %15.4g'%(cid+item,parmDict[cid+item],sigDict[cid+item]) 982 1179 distDict = controls['DistType'] … … 999 1196 FFfxn = parmDict[cid+'FormFact'] 1000 1197 Volfxn = parmDict[cid+'FFVolume'] 1198 SFfxn = parmDict[cid+'StrFact'] 1001 1199 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']: 1003 1202 if item in parmDict: 1004 1203 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]) 1005 1207 distDict = {} 1006 1208 for item in [cid+'Volume',cid+'Mean',cid+'StdDev',cid+'MinSize',]: … … 1012 1214 Gmat = G_matrix(Q,rBins,contrast,FFfxn,Volfxn,FFargs).T 1013 1215 dist *= parmDict[cid+'Volume'] 1014 Ic += np.dot(Gmat,dist) 1216 Ic += np.dot(Gmat,dist)*SFfxn(Q,args=SFargs) 1015 1217 elif 'Unified' in Type: 1016 1218 Rg,G,B,P,Rgco = parmDict[cid+'Rg'],parmDict[cid+'G'],parmDict[cid+'B'], \ … … 1027 1229 FFfxn = parmDict[cid+'FormFact'] 1028 1230 Volfxn = parmDict[cid+'FFVolume'] 1231 SFfxn = parmDict[cid+'StrFact'] 1029 1232 FFargs = [] 1233 SFargs = [] 1030 1234 for item in [cid+'Aspect ratio',cid+'Length',cid+'Thickness',cid+'Diameter',]: 1031 1235 if item in parmDict: … … 1035 1239 R = parmDict[cid+'Radius'] 1036 1240 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) 1038 1242 elif 'Bragg' in distFxn: 1039 1243 Ic += parmDict[cid+'PkInt']*G2pwd.getPsVoigt(parmDict[cid+'PkPos'], … … 1084 1288 1085 1289 def ModelFxn(Profile,ProfDict,Limits,Substances,Sample,sasdData): 1290 1086 1291 shapes = {'Spheroid':[SpheroidFF,SpheroidVol],'Cylinder':[CylinderDFF,CylinderDVol], 1087 1292 'Cylinder AR':[CylinderARFF,CylinderARVol],'Unified sphere':[UniSphereFF,UniSphereVol], … … 1089 1294 'Unified disk':[UniDiskFF,UniDiskVol],'Sphere':[SphereFF,SphereVol], 1090 1295 '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 1091 1299 # pdb.set_trace() 1092 1300 partData = sasdData['Particle'] … … 1109 1317 distFxn = controls['DistType'] 1110 1318 if distFxn in ['LogNormal','Gaussian','LSW','Schulz-Zimm']: 1319 parmDict = level[controls['DistType']] 1111 1320 FFfxn = shapes[controls['FormFact']][0] 1112 1321 Volfxn = shapes[controls['FormFact']][1] 1322 SFfxn = sfxns[controls['StrFact']] 1113 1323 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]) 1114 1328 for item in ['Aspect ratio','Length','Thickness','Diameter',]: 1115 1329 if item in controls['FFargs']: … … 1117 1331 rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0) 1118 1332 contrast = rho**2-rhoMat**2 1119 parmDict = level[controls['DistType']]1120 1333 distDict = {} 1121 1334 for item in parmDict: … … 1124 1337 Gmat = G_matrix(Q[Ibeg:Ifin],rBins,contrast,FFfxn,Volfxn,FFargs).T 1125 1338 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) 1127 1340 Rbins.append(rBins) 1128 1341 Dist.append(dist/(4.*dBins)) … … 1145 1358 Dist.append([]) 1146 1359 elif 'Mono' in distFxn: 1360 parmDict = level[controls['DistType']] 1361 R = level[controls['DistType']]['Radius'][0] 1147 1362 FFfxn = shapes[controls['FormFact']][0] 1148 1363 Volfxn = shapes[controls['FormFact']][1] 1364 SFfxn = sfxns[controls['StrFact']] 1149 1365 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]) 1150 1370 for item in ['Aspect ratio','Length','Thickness','Diameter',]: 1151 1371 if item in controls['FFargs']: … … 1153 1373 rho = Substances['Substances'][level['Controls']['Material']].get('XAnom density',0.0) 1154 1374 contrast = rho**2-rhoMat**2 1155 R = level[controls['DistType']]['Radius'][0]1156 1375 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) 1158 1377 Rbins.append([]) 1159 1378 Dist.append([])
Note: See TracChangeset
for help on using the changeset viewer.