source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/CoreShellCylinder_v40.ipf @ 510

Last change on this file since 510 was 510, checked in by srkline, 14 years ago

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 7.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////
5// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
6// Adopting these into the experiment will insure that they are always present
7////////////////////////////////////////////////
8// this function is for the form factor of a right circular cylinder with core/shell scattering length density profile
9//
10// the core dimensions are given and a constant shell thickness is added to the radius and to dach and of the length
11// this way, the scattering amplitude is simply the difference between two cylinders of different dimensions
12//
13// 06 NOV 98 SRK
14////////////////////////////////////////////////
15
16Proc PlotCoreShellCylinder(num,qmin,qmax)
17        Variable num=128,qmin=0.001,qmax=0.7
18        Prompt num "Enter number of data points for model: "
19        Prompt qmin "Enter minimum q-value (A^-1) for model: "
20        Prompt qmax "Enter maximum q-value (A^-1) for model: "
21       
22        make/o/d/n=(num) xwave_cscyl,ywave_cscyl
23        xwave_cscyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
24        make/o/d coef_cscyl = {1.,20.,10.,400,1.0e-6,4.0e-6,1.0e-6,0.01}
25        make/o/t parameters_cscyl = {"scale","core radius (A)","shell THICKNESS (A)","CORE length (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
26        Edit parameters_cscyl,coef_cscyl
27        Variable/G root:g_cscyl
28        g_cscyl := CoreShellCylinder(coef_cscyl,ywave_cscyl,xwave_cscyl)
29//      ywave_cscyl := CoreShellCylinder(coef_cscyl,xwave_cscyl)
30        Display ywave_cscyl vs xwave_cscyl
31        ModifyGraph log=1,marker=29,msize=2,mode=4
32        Label bottom "q (A\\S-1\\M)"
33        Label left "Intensity (cm\\S-1\\M)"
34        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
35       
36        AddModelToStrings("CoreShellCylinder","coef_cscyl","parameters_cscyl","cscyl")
37End
38
39///////////////////////////////////////////////////////////
40// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
41Proc PlotSmearedCoreShellCylinder(str)                                                         
42        String str
43        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
44       
45        // if any of the resolution waves are missing => abort
46        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
47                Abort
48        endif
49       
50        SetDataFolder $("root:"+str)
51       
52        // Setup parameter table for model function
53        make/o/d smear_coef_cscyl =  {1.,20.,10.,400,1.0e-6,4.0e-6,1.0e-6,0.01}
54        make/o/t smear_parameters_cscyl = {"scale","core radius (A)","shell radius (A)","length (A)","SLD core (A^-2)","SLD shell (A^-2)","SLD solvent (A^-2)","incoh. bkg (cm^-1)"}
55        Edit smear_parameters_cscyl,smear_coef_cscyl
56       
57        // output smeared intensity wave, dimensions are identical to experimental QSIG values
58        // make extra copy of experimental q-values for easy plotting
59        Duplicate/O $(str+"_q") smeared_cscyl,smeared_qvals
60        SetScale d,0,0,"1/cm",smeared_cscyl     
61                                       
62        Variable/G gs_cscyl=0
63        gs_cscyl := fSmearedCoreShellCylinder(smear_coef_cscyl,smeared_cscyl,smeared_qvals)     //this wrapper fills the STRUCT
64       
65        Display smeared_cscyl vs smeared_qvals
66        ModifyGraph log=1,marker=29,msize=2,mode=4
67        Label bottom "q (A\\S-1\\M)"
68        Label left "Intensity (cm\\S-1\\M)"
69        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
70       
71        SetDataFolder root:
72        AddModelToStrings("SmearedCoreShellCylinder","smear_coef_cscyl","smear_parameters_cscyl","cscyl")
73End
74       
75
76//AAO version
77Function CoreShellCylinder(cw,yw,xw) : FitFunc
78        Wave cw,yw,xw
79
80#if exists("CoreShellCylinderX")
81        yw = CoreShellCylinderX(cw,xw)
82#else
83        yw = fCoreShellCylinder(cw,xw)
84#endif
85        return(0)
86End
87
88///////////////////////////////////////////////////////////////
89// unsmeared model calculation
90///////////////////////////
91Function fCoreShellCylinder(w,x) : FitFunc
92        Wave w
93        Variable x
94
95//The input variables are (and output)
96        //[0] scale
97        //[1] cylinder CORE RADIUS (A)
98        //[2] shell Thickness (A)
99        //[3]  cylinder CORE LENGTH (A)
100        //[4] core SLD (A^-2)
101        //[5] shell SLD (A^-2)
102        //[6] solvent SLD (A^-2)
103        //[7] background (cm^-1)
104        Variable scale,length,delrho,bkg,rcore,thick,rhoc,rhos,rhosolv
105        scale = w[0]
106        rcore = w[1]
107        thick = w[2]
108        length = w[3]
109        rhoc = w[4]
110        rhos = w[5]
111        rhosolv = w[6]
112        bkg = w[7]
113//
114// the OUTPUT form factor is <f^2>/Vcyl [cm-1]
115//
116
117// local variables
118        Variable nord,ii,va,vb,contr,vcyl,nden,summ,yyy,zi,qq,halfheight
119        Variable answer
120        String weightStr,zStr
121       
122        weightStr = "gauss76wt"
123        zStr = "gauss76z"
124
125       
126//      if wt,z waves don't exist, create them
127// 20 Gauss points is not enough for cylinder calculation
128       
129        if (WaveExists($weightStr) == 0) // wave reference is not valid,
130                Make/D/N=76 $weightStr,$zStr
131                Wave w76 = $weightStr
132                Wave z76 = $zStr                // wave references to pass
133                Make76GaussPoints(w76,z76)     
134        //                  printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
135        else
136                if(exists(weightStr) > 1)
137                         Abort "wave name is already in use"    // execute if condition is false
138                endif
139                Wave w76 = $weightStr
140                Wave z76 = $zStr                // Not sure why this has to be "declared" twice
141        //          printf "w[0],z[0] = %g %g\r", w76[0],z76[0]
142        endif
143
144
145// set up the integration
146        // end points and weights
147        nord = 76
148        va = 0
149        vb = Pi/2
150      halfheight = length/2.0
151
152// evaluate at Gauss points
153        // remember to index from 0,size-1
154
155        qq = x          //current x point is the q-value for evaluation
156      summ = 0.0                // initialize integral
157      ii=0
158      do
159                // Using 76 Gauss points
160                zi = ( z76[ii]*(vb-va) + vb + va )/2.0         
161                yyy = w76[ii] * CoreShellcyl(qq, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi)
162                summ += yyy
163
164                ii+=1
165        while (ii<nord)                         // end of loop over quadrature points
166//   
167// calculate value of integral to return
168
169      answer = (vb-va)/2.0*summ
170     
171// contrast is now explicitly included in the core-shell calculation
172
173//normalize by cylinder volume
174//NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
175//calculate TOTAL volume
176// length is the total core length
177        vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2*thick)
178        answer /= vcyl
179//convert to [cm-1]
180        answer *= 1.0e8
181//Scale
182        answer *= scale
183// add in the background
184        answer += bkg
185
186        Return (answer)
187       
188End             //End of function CoreShellCylinderForm()
189
190///////////////////////////////////////////////////////////////
191// F(qq, rcore, thick, rhoc,rhos,rhosolv, length, zi)
192//
193Function CoreShellcyl(qq, rcore, thick, rhoc,rhos,rhosolv, length, dum)
194        Variable qq, rcore, thick, rhoc,rhos,rhosolv, length, dum
195       
196// qq is the q-value for the calculation (1/A)
197// rcore is the core radius of the cylinder (A)
198//thick is the uniform thickness
199// rho(n) are the respective SLD's
200
201// length is the *Half* CORE-LENGTH of the cylinder = L (A)
202
203// dum is the dummy variable for the integration (x in Feigin's notation)
204
205   //Local variables
206        Variable dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval
207       
208        dr1 = rhoc-rhos
209        dr2 = rhos-rhosolv
210        vol1 = Pi*rcore*rcore*(2*length)
211        vol2 = Pi*(rcore+thick)*(rcore+thick)*(2*length+2*thick)
212       
213        besarg1 = qq*rcore*sin(dum)
214        besarg2 = qq*(rcore+thick)*sin(dum)
215        sinarg1 = qq*length*cos(dum)
216        sinarg2 = qq*(length+thick)*cos(dum)
217       
218        t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*bessJ(1,besarg1)/besarg1
219        t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*bessJ(1,besarg2)/besarg2
220       
221        retval = ((t1+t2)^2)*sin(dum)
222       
223    return retval
224   
225End     //Function CoreShellcyl()
226
227// this is all there is to the smeared calculation!
228Function SmearedCoreShellCylinder(s) :FitFunc
229        Struct ResSmearAAOStruct &s
230
231////the name of your unsmeared model is the first argument
232        Smear_Model_20(CoreShellCylinder,s.coefW,s.xW,s.yW,s.resW)
233
234        return(0)
235End
236
237
238//wrapper to calculate the smeared model as an AAO-Struct
239// fills the struct and calls the ususal function with the STRUCT parameter
240//
241// used only for the dependency, not for fitting
242//
243Function fSmearedCoreShellCylinder(coefW,yW,xW)
244        Wave coefW,yW,xW
245       
246        String str = getWavesDataFolder(yW,0)
247        String DF="root:"+str+":"
248       
249        WAVE resW = $(DF+str+"_res")
250       
251        STRUCT ResSmearAAOStruct fs
252        WAVE fs.coefW = coefW   
253        WAVE fs.yW = yW
254        WAVE fs.xW = xW
255        WAVE fs.resW = resW
256       
257        Variable err
258        err = SmearedCoreShellCylinder(fs)
259       
260        return (0)
261End
Note: See TracBrowser for help on using the repository browser.