source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/SmearedRPA.ipf @ 153

Last change on this file since 153 was 151, checked in by srkline, 15 years ago

(1) - cursors can now be used to select a subrange of USANS data to fit. This is done by th fit wrapper, assigning a subrange of resW to the struct

(2) all of the smeared model functions are now in the latest form of Smear_Model_N() that is NOT a pointwise calculation anymore, since the USANS matrix smearing in inherently not so.

File size: 16.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4// currently, there is NO XOP version of RPA, since there is an extra input parameter
5// wave that must be carried into the function.
6// this can be done with the STRUCT, if extra space is allocated for such
7//
8//
9Proc PlotRPAForm(num,qmin,qmax,nCASE)
10        Variable num=100,qmin=0.001,qmax=0.5
11        Variable/g gCASE                //Global Variable
12        Variable nCASE
13        Prompt num "Enter number of data points: "
14        Prompt qmin "Enter minimum Q value (A^-1): "
15        Prompt qmax "Enter maximum Q value (A^-1): "
16        Prompt nCASE, "Choose one of the following cases:",popup "CASE 1: C/D BINARY MIXTURE OF HOMOPOLYMERS;"
17        "CASE 2: C-D DIBLOCK COPOLYMER;"
18        "CASE 3: B/C/D TERNARY MIXTURE OF HOMOPOLYMERS;"
19        "CASE 4: B/C-D MIXTURE OF HOMOPOLYMER B AND DIBLOCK COPOLYMER C-D;"
20        "CASE 5: B-C-D TRIBLOCK COPOLYMER;"
21        "CASE 6: A/B/C/D QUATERNARY MIXTURE OF HOMOPOLYMERS;"
22        "CASE 7: A/B/C-D MIXTURE OF TWO HOMOPOLYMERS A/B AND A DIBLOCK C-D;"
23        "CASE 8: A/B-C-D MIXTURE OF A HOMOPOLYMER A AND A TRIBLOCK B-C-D;"
24        "CASE 9: A-B/C-D MIXTURE OF TWO DIBLOCK COPOLYMERS A-B AND C-D;"
25        "CASE 10: A-B-C-D FOUR-BLOCK COPOLYMER"
26
27        Make/O/D/n=(num) xwave_rpa,ywave_rpa
28        //xwave_rpa=qmin+x*(qmax-qmin)/num
29        //switch to log-scaling of the q-values
30        xwave_rpa=alog(log(qmin) + x*((log(qmax)-log(qmin))/num))       
31       
32        gCASE = nCASE
33//      print gCASE
34
35        IF(gCASE <=2)
36                Make/O/D inputvalues={1000,0.5,100,1.e-12,1000,0.5,100,0.e-12}
37                make/o/t inputnames={"Deg. Polym. Nc","Vol. Frac. Phic","Spec. Vol. Vc","Scatt. Length Lc","Deg. Polym. Nd","Vol. Frac. Phid","Spec. Vol. Vd","Scatt. LengthLd"}
38                Edit inputnames,inputvalues
39                Make/O/D coef_rpa = {5,5,-.0004,1,0}
40                make/o/t parameters_rpa = {"Seg. Length bc","Seg. Length bd","Chi Param. Kcd","scale","Background"}
41                Edit parameters_rpa,coef_rpa
42        ENDIF
43
44        IF((gCASE >2) %& (gCASE<=5))
45                Make/O/D inputvalues={1000,0.33,100,1.e-12,1000,0.33,100,1.e-12,1000,0.33,100,0.e-12}
46                make/o/t inputnames={"Deg. Polym. Nb","Vol. Frac. Phib","Spec. Vol. Vb","Scatt. Length Lb","Deg. Polym. Nc","Vol. Frac. Phic","Spec. Vol. Vc","Scatt. Length Lc","Deg. Polym. Nd","Vol. Frac. Phid","Spec. Vol. Vd","Scatt. Length Ld"}
47                Edit inputnames,inputvalues
48                Make/O/D coef_rpa = {5,5,5,-.0004,-.0004,-.0004,1,0}
49                make/o/t parameters_rpa = {"Seg. Length bb","Seg. Length bc","Seg. Length bd","Chi Param. Kbc","Chi Param. Kbd","Chi Param. Kcd","scale","Background"}
50                Edit parameters_rpa,coef_rpa
51        ENDIF
52
53        IF(gCASE >5)
54                Make/O/D inputvalues={1000,0.25,100,1.e-12,1000,0.25,100,1.e-12,1000,0.25,100,1.e-12,1000,0.25,100,0.e-12}
55                make/o/t inputnames={"Deg. Polym. Na","Vol. Frac. Phia","Spec. Vol. Va","Scatt. Length La","Deg. Polym. Nb","Vol. Frac. Phib","Spec. Vol. Vb","Scatt. Length Lb","Deg. Polym. Nc","Vol. Frac. Phic","Spec. Vol. Vc","Scatt. Length Lc","Deg. Polym. Nd","Vol. Frac. Phid","Spec. Vol. Vd","Scatt. Length Ld"}
56                Edit inputnames,inputvalues
57                Make/O/D coef_rpa = {5,5,5,5,-.0004,-.0004,-.0004,-.0004,-.0004,-.0004,1,0}
58                make/o/t parameters_rpa = {"Seg. Length ba","Seg. Length bb","Seg. Length bc","Seg. Length bd","Chi Param. Kab","Chi Param. Kac","Chi Param. Kad","Chi Param. Kbc","Chi Param. Kbd","Chi Param. Kcd","scale","Background"}
59                Edit parameters_rpa,coef_rpa
60        ENDIF
61       
62        Variable/G root:g_rpa
63        g_rpa := RPAForm(coef_rpa,ywave_rpa,xwave_rpa)
64        Display ywave_rpa vs xwave_rpa
65//      ModifyGraph log=1,marker=29,msize=2,mode=4
66        Label bottom "Q (A\\S-1\\M)"
67        Label left "Intensity (cm\\S-1\\M)"
68        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
69End
70
71// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
72Proc PlotSmearedRPAForm(str,nCASE)                                                             
73        String str
74        Variable nCASE
75        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
76        Prompt nCASE, "Choose one of the following cases:",popup "CASE 1: C/D BINARY MIXTURE OF HOMOPOLYMERS;"
77        "CASE 2: C-D DIBLOCK COPOLYMER;"
78        "CASE 3: B/C/D TERNARY MIXTURE OF HOMOPOLYMERS;"
79        "CASE 4: B/C-D MIXTURE OF HOMOPOLYMER B AND DIBLOCK COPOLYMER C-D;"
80        "CASE 5: B-C-D TRIBLOCK COPOLYMER;"
81        "CASE 6: A/B/C/D QUATERNARY MIXTURE OF HOMOPOLYMERS;"
82        "CASE 7: A/B/C-D MIXTURE OF TWO HOMOPOLYMERS A/B AND A DIBLOCK C-D;"
83        "CASE 8: A/B-C-D MIXTURE OF A HOMOPOLYMER A AND A TRIBLOCK B-C-D;"
84        "CASE 9: A-B/C-D MIXTURE OF TWO DIBLOCK COPOLYMERS A-B AND C-D;"
85        "CASE 10: A-B-C-D FOUR-BLOCK COPOLYMER"
86       
87        SetDataFolder $("root:"+str)
88
89        Variable/g gCASE                //Global Variable
90
91        // if any of the resolution waves are missing => abort
92        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
93                Abort
94        endif
95       
96        gCASE = nCASE
97//      print gCASE
98
99       
100        IF(gCASE <=2)
101                Make/O/D inputvalues={1000,0.5,100,1.e-12,1000,0.5,100,0.e-12}
102                make/o/t inputnames={"Deg. Polym. Nc","Vol. Frac. Phic","Spec. Vol. Vc","Scatt. Length Lc","Deg. Polym. Nd","Vol. Frac. Phid","Spec. Vol. Vd","Scatt. LengthLd"}
103                Edit inputnames,inputvalues
104                Make/O/D smear_coef_rpa = {5,5,-.0004,1,0}
105                make/o/t smear_parameters_rpa = {"Seg. Length bc","Seg. Length bd","Chi Param. Kcd","scale","Background"}
106                Edit smear_parameters_rpa,smear_coef_rpa
107        ENDIF
108
109        IF((gCASE >2) %& (gCASE<=5))
110                Make/O/D inputvalues={1000,0.33,100,1.e-12,1000,0.33,100,1.e-12,1000,0.33,100,0.e-12}
111                make/o/t inputnames={"Deg. Polym. Nb","Vol. Frac. Phib","Spec. Vol. Vb","Scatt. Length Lb","Deg. Polym. Nc","Vol. Frac. Phic","Spec. Vol. Vc","Scatt. Length Lc","Deg. Polym. Nd","Vol. Frac. Phid","Spec. Vol. Vd","Scatt. Length Ld"}
112                Edit inputnames,inputvalues
113                Make/O/D smear_coef_rpa = {5,5,5,-.0004,-.0004,-.0004,1,0}
114                make/o/t smear_parameters_rpa = {"Seg. Length bb","Seg. Length bc","Seg. Length bd","Chi Param. Kbc","Chi Param. Kbd","Chi Param. Kcd","scale","Background"}
115                Edit smear_parameters_rpa,smear_coef_rpa
116        ENDIF
117
118        IF(gCASE >5)
119                Make/O/D inputvalues={1000,0.25,100,1.e-12,1000,0.25,100,1.e-12,1000,0.25,100,1.e-12,1000,0.25,100,0.e-12}
120                make/o/t inputnames={"Deg. Polym. Na","Vol. Frac. Phia","Spec. Vol. Va","Scatt. Length La","Deg. Polym. Nb","Vol. Frac. Phib","Spec. Vol. Vb","Scatt. Length Lb","Deg. Polym. Nc","Vol. Frac. Phic","Spec. Vol. Vc","Scatt. Length Lc","Deg. Polym. Nd","Vol. Frac. Phid","Spec. Vol. Vd","Scatt. Length Ld"}
121                Edit inputnames,inputvalues
122                Make/O/D smear_coef_rpa = {5,5,5,5,-.0004,-.0004,-.0004,-.0004,-.0004,-.0004,1,0}
123                make/o/t smear_parameters_rpa = {"Seg. Length ba","Seg. Length bb","Seg. Length bc","Seg. Length bd","Chi Param. Kab","Chi Param. Kac","Chi Param. Kad","Chi Param. Kbc","Chi Param. Kbd","Chi Param. Kcd","scale","Background"}
124                Edit smear_parameters_rpa,smear_coef_rpa
125        ENDIF
126       
127        // output smeared intensity wave, dimensions are identical to experimental QSIG values
128        // make extra copy of experimental q-values for easy plotting
129        Duplicate/O $(str+"_q") smeared_rpa,smeared_qvals       
130        SetScale d,0,0,"1/cm",smeared_rpa                                       
131               
132        Variable/G gs_rpa=0
133        gs_rpa := fSmearedRPAForm(smear_coef_rpa,smeared_rpa,smeared_qvals)     //this wrapper fills the STRUCT
134       
135        Display smeared_rpa vs smeared_qvals
136        ModifyGraph log=1,marker=29,msize=2,mode=4
137        Label bottom "q (A\\S-1\\M)"
138        Label left "Intensity (cm\\S-1\\M)"
139        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
140       
141        SetDataFolder root:
142End
143
144///////////////////////////////////////////////////////////////
145//AAO version, uses XOP if available
146// simply calls the original single point calculation with
147// a wave assignment (this will behave nicely if given point ranges)
148Function RPAForm(cw,yw,xw) : FitFunc
149        Wave cw,yw,xw
150       
151#if exists("RPAFormX")
152        yw = RPAFormX(cw,xw)
153#else
154        yw = fRPAForm(cw,xw)
155#endif
156        return(0)
157End
158
159
160Function fRPAForm(w,x) : FitFunc
161        Wave w
162        Variable x
163       
164        Wave var=$"inputvalues"
165        Nvar lCASE=gCASE
166//      print lCASE
167//      Variable lCASE
168//      RANDOM PHASE APPROXIMATION FOR A FOUR-BLOCK COPOLYMER A-B-C-D
169//      THIS FORMALISM APPLIES TO MULTICOMPONENT POLYMER MIXTURES
170//      IN THE HOMOGENEOUS (MIXED) PHASE REGION ONLY.
171
172//      THIS GENERAL CASE INCLUDES A MULTITUDE OF (TEN) SPECIAL CASES.
173
174//      HERE ARE THE VARIOUS CASES COVERED:
175//      CASE 1: C/D BINARY MIXTURE OF HOMOPOLYMERS
176//      CASE 2: C-D DIBLOCK COPOLYMER
177//      CASE 3: B/C/D TERNARY MIXTURE OF HOMOPOLYMERS
178//      CASE 4: B/C-D MIXTURE OF HOMOPOLYMER B AND DIBLOCK COPOLYMER C-D
179//      CASE 5: B-C-D TRIBLOCK COPOLYMER
180//      CASE 6: A/B/C/D QUATERNARY MIXTURE OF HOMOPOLYMERS
181//      CASE 7: A/B/C-D MIXTURE OF TWO HOMOPOLYMERS A/B AND A DIBLOCK C-D
182//      CASE 8: A/B-C-D MIXTURE OF A HOMOPOLYMER A AND A TRIBLOCK B-C-D
183//      CASE 9: A-B/C-D MIXTURE OF TWO DIBLOCK COPOLYMERS A-B AND C-D
184//      CASE 10: A-B-C-D FOUR-BLOCK COPOLYMER
185
186//      B. HAMMOUDA, NIST, JULY 1998
187
188        Variable Na,Nb,Nc,Nd,Nab,Nac,Nad,Nba,Nbc,Nbd,Nca,Ncb,Ncd
189        Variable Phia,Phib,Phic,Phid,Phiab,Phiac,Phiad
190        Variable Phiba,Phibc,Phibd,Phica,Phicb,Phicd,Phida,Phidb,Phidc
191        Variable va,vb,vc,vd,vab,vac,vad,vba,vbc,vbd,vca,vcb,vcd,vda,vdb,vdc
192        Variable m
193        Variable ba,bb,bc,bd
194        Variable Q
195        Variable Xa,Xb,Xc,Xd
196        Variable Paa,S0aa,Pab,S0ab,Pac,S0ac,Pad,S0ad
197        Variable Pba,S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd
198        Variable Pca,S0ca,Pcb,S0cb,Pcc,S0cc,Pcd,S0cd
199        Variable Pda,S0da,Pdb,S0db,Pdc,S0dc,Pdd,S0dd
200        Variable Kaa,Kab,Kac,Kad,Kba,Kbb,Kbc,Kbd
201        Variable Kca,Kcb,Kcc,Kcd,Kda,Kdb,Kdc,Kdd
202        Variable Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc
203        Variable DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33
204        Variable Y1,Y2,Y3,X11,X12,X13,X21,X22,X23,X31,X32,X33
205        Variable ZZ,DenQ1,DenQ2,DenQ3,DenQ,Q11,Q12,Q13,Q21,Q22,Q23,Q31,Q32,Q33
206        Variable N11,N12,N13,N21,N22,N23,N31,N32,N33
207        Variable M11,M12,M13,M21,M22,M23,M31,M32,M33
208        Variable S11,S12,S13,S14,S21,S22,S23,S24
209        Variable S31,S32,S33,S34,S41,S42,S43,S44
210        Variable La,Lb,Lc,Ld,Lad,Lbd,Lcd,Nav,Int
211        Variable scale
212        Variable Background
213        Variable ii=0   //to fool curve fitting dialog to let you pick the coefficient wave
214
215        Na=1000
216        Nb=1000
217        Nc=1000
218        Nd=1000 //DEGREE OF POLYMERIZATION     
219        Phia=0.25
220        Phib=0.25
221        Phic=0.25
222        Phid=0.25       //VOL FRACTION
223        Kab=-.0004
224        Kac=-.0004
225        Kad=-.0004      //CHI PARAM
226        Kbc=-.0004
227        Kbd=-.0004
228        Kcd=-.0004     
229        La=1.E-12
230        Lb=1.E-12
231        Lc=1.E-12
232        Ld=0.E-12       //SCATT. LENGTH
233        va=100
234        vb=100
235        vc=100
236        vd=100          //SPECIFIC VOLUME
237        ba=5
238        bb=5
239        bc=5
240        bd=5            //SEGMENT LENGTH
241
242        IF (lCASE <= 2)
243                Phia=0.0000001
244                Phib=0.0000001
245                Phic=0.5
246                Phid=0.5
247                Nc=var[0]
248                Phic=var[1]
249                vc=var[2]
250                Lc=var[3]
251                bc=w[ii+0]
252                Nd=var[4]
253                Phid=var[5]
254                vd=var[6]
255                Ld=var[7]
256                bd=w[ii+1]
257                Kcd=w[ii+2]
258                scale=w[ii+3]   
259                background=w[ii+4]
260        ENDIF
261
262        IF ((lCASE > 2) %& (lCASE <= 5))
263                Phia=0.0000001
264                Phib=0.333333
265                Phic=0.333333
266                Phid=0.333333
267                Nb=var[0]
268                Phib=var[1]
269                vb=var[2]
270                Lb=var[3]
271                bb=w[ii+0]
272                Nc=var[4]
273                Phic=var[5]
274                vc=var[6]
275                Lc=var[7]
276                bc=w[ii+1]
277                Nd=var[8]
278                Phid=var[9]
279                vd=var[10]
280                Ld=var[11]
281                bd=w[ii+2]
282                Kbc=w[ii+3]
283                Kbd=w[ii+4]
284                Kcd=w[ii+5]
285                scale=w[ii+6]   
286                background=w[ii+7]
287        ENDIF
288
289        IF (lCASE > 5)
290                Phia=0.25
291                Phib=0.25
292                Phic=0.25
293                Phid=0.25
294                Na=var[0]
295                Phia=var[1]
296                va=var[2]
297                La=var[3]
298                ba=w[ii+0]
299                Nb=var[4]
300                Phib=var[5]
301                vb=var[6]
302                Lb=var[7]
303                bb=w[ii+1]
304                Nc=var[8]
305                Phic=var[9]
306                vc=var[10]
307                Lc=var[11]
308                bc=w[ii+2]
309                Nd=var[12]
310                Phid=var[13]
311                vd=var[14]
312                Ld=var[15]
313                bd=w[ii+3]
314                Kab=w[ii+4]
315                Kac=w[ii+5]
316                Kad=w[ii+6]
317                Kbc=w[ii+7]
318                Kbd=w[ii+8]
319                Kcd=w[ii+9]
320                scale=w[ii+10] 
321                background=w[ii+11]
322        ENDIF
323
324        Nab=(Na*Nb)^(0.5)
325        Nac=(Na*Nc)^(0.5)
326        Nad=(Na*Nd)^(0.5)
327        Nbc=(Nb*Nc)^(0.5)
328        Nbd=(Nb*Nd)^(0.5)
329        Ncd=(Nc*Nd)^(0.5)
330
331        vab=(va*vb)^(0.5)
332        vac=(va*vc)^(0.5)
333        vad=(va*vd)^(0.5)
334        vbc=(vb*vc)^(0.5)
335        vbd=(vb*vd)^(0.5)
336        vcd=(vc*vd)^(0.5)
337
338        Phiab=(Phia*Phib)^(0.5)
339        Phiac=(Phia*Phic)^(0.5)
340        Phiad=(Phia*Phid)^(0.5)
341        Phibc=(Phib*Phic)^(0.5)
342        Phibd=(Phib*Phid)^(0.5)
343        Phicd=(Phic*Phid)^(0.5)
344
345        Q=x
346        Xa=Q^2*ba^2*Na/6
347        Xb=Q^2*bb^2*Nb/6
348        Xc=Q^2*bc^2*Nc/6
349        Xd=Q^2*bd^2*Nd/6
350
351        Paa=2*(Exp(-Xa)-1+Xa)/Xa^2
352        S0aa=Na*Phia*va*Paa
353        Pab=((1-Exp(-Xa))/Xa)*((1-Exp(-Xb))/Xb)
354        S0ab=(Phiab*vab*Nab)*Pab
355        Pac=((1-Exp(-Xa))/Xa)*Exp(-Xb)*((1-Exp(-Xc))/Xc)
356        S0ac=(Phiac*vac*Nac)*Pac
357        Pad=((1-Exp(-Xa))/Xa)*Exp(-Xb-Xc)*((1-Exp(-Xd))/Xd)
358        S0ad=(Phiad*vad*Nad)*Pad
359
360        S0ba=S0ab
361        Pbb=2*(Exp(-Xb)-1+Xb)/Xb^2
362        S0bb=Nb*Phib*vb*Pbb
363        Pbc=((1-Exp(-Xb))/Xb)*((1-Exp(-Xc))/Xc)
364        S0bc=(Phibc*vbc*Nbc)*Pbc
365        Pbd=((1-Exp(-Xb))/Xb)*Exp(-Xc)*((1-Exp(-Xd))/Xd)
366        S0bd=(Phibd*vbd*Nbd)*Pbd
367
368        S0ca=S0ac
369        S0cb=S0bc
370        Pcc=2*(Exp(-Xc)-1+Xc)/Xc^2
371        S0cc=Nc*Phic*vc*Pcc
372        Pcd=((1-Exp(-Xc))/Xc)*((1-Exp(-Xd))/Xd)
373        S0cd=(Phicd*vcd*Ncd)*Pcd
374
375        S0da=S0ad
376        S0db=S0bd
377        S0dc=S0cd
378        Pdd=2*(Exp(-Xd)-1+Xd)/Xd^2
379        S0dd=Nd*Phid*vd*Pdd
380
381        IF(lCASE == 1)
382                S0aa=0.000001
383                S0ab=0.000002           
384                S0ac=0.000003           
385                S0ad=0.000004           
386                S0bb=0.000005           
387                S0bc=0.000006           
388                S0bd=0.000007           
389                S0cd=0.000008
390        ENDIF
391
392        IF(lCASE == 2)
393                S0aa=0.000001
394                S0ab=0.000002           
395                S0ac=0.000003           
396                S0ad=0.000004           
397                S0bb=0.000005           
398                S0bc=0.000006           
399                S0bd=0.000007           
400        ENDIF
401       
402        IF(lCASE == 3)
403                S0aa=0.000001           
404                S0ab=0.000002           
405                S0ac=0.000003           
406                S0ad=0.000004           
407                S0bc=0.000005           
408                S0bd=0.000006           
409                S0cd=0.000007           
410        ENDIF
411
412        IF(lCASE == 4)
413                S0aa=0.000001           
414                S0ab=0.000002           
415                S0ac=0.000003           
416                S0ad=0.000004           
417                S0bc=0.000005           
418                S0bd=0.000006           
419        ENDIF
420
421        IF(lCASE == 5)
422                S0aa=0.000001           
423                S0ab=0.000002           
424                S0ac=0.000003           
425                S0ad=0.000004           
426        ENDIF
427
428        IF(lCASE == 6)
429                S0ab=0.000001           
430                S0ac=0.000002           
431                S0ad=0.000003           
432                S0bc=0.000004           
433                S0bd=0.000005           
434                S0cd=0.000006
435        ENDIF
436
437        IF(lCASE == 7)
438                S0ab=0.000001           
439                S0ac=0.000002           
440                S0ad=0.000003           
441                S0bc=0.000004           
442                S0bd=0.000005           
443        ENDIF
444
445        IF(lCASE == 8)
446                S0ab=0.000001           
447                S0ac=0.000002           
448                S0ad=0.000003           
449        ENDIF
450
451        IF(lCASE == 9)
452                S0ac=0.000001           
453                S0ad=0.000002           
454                S0bc=0.000003           
455                S0bd=0.000004
456        ENDIF
457
458        S0ba=S0ab               
459        S0ca=S0ac               
460        S0cb=S0bc               
461        S0da=S0ad               
462        S0db=S0bd               
463        S0dc=S0cd               
464
465        Kaa=0
466        Kbb=0
467        Kcc=0
468        Kdd=0
469
470        Kba=Kab
471        Kca=Kac
472        Kcb=Kbc
473        Kda=Kad
474        Kdb=Kbd
475        Kdc=Kcd
476
477        Zaa=Kaa-Kad-Kad 
478        Zab=Kab-Kad-Kbd 
479        Zac=Kac-Kad-Kcd
480        Zba=Kba-Kbd-Kad 
481        Zbb=Kbb-Kbd-Kbd 
482        Zbc=Kbc-Kbd-Kcd
483        Zca=Kca-Kcd-Kad 
484        Zcb=Kcb-Kcd-Kbd 
485        Zcc=Kcc-Kcd-Kcd
486
487        DenT=(-(S0ac*S0bb*S0ca) + S0ab*S0bc*S0ca + S0ac*S0ba*S0cb - S0aa*S0bc*S0cb - S0ab*S0ba*S0cc + S0aa*S0bb*S0cc)
488
489        T11= (-(S0bc*S0cb) + S0bb*S0cc)/DenT
490        T12= (S0ac*S0cb - S0ab*S0cc)/DenT
491        T13= (-(S0ac*S0bb) + S0ab*S0bc)/DenT
492        T21= (S0bc*S0ca - S0ba*S0cc)/DenT
493        T22= (-(S0ac*S0ca) + S0aa*S0cc)/DenT
494        T23= (S0ac*S0ba - S0aa*S0bc)/DenT
495        T31= (-(S0bb*S0ca) + S0ba*S0cb)/DenT
496        T32= (S0ab*S0ca - S0aa*S0cb)/DenT
497        T33= (-(S0ab*S0ba) + S0aa*S0bb)/DenT
498
499        Y1=T11*S0ad+T12*S0bd+T13*S0cd+1
500        Y2=T21*S0ad+T22*S0bd+T23*S0cd+1
501        Y3=T31*S0ad+T32*S0bd+T33*S0cd+1
502
503        X11=Y1*Y1
504        X12=Y1*Y2
505        X13=Y1*Y3
506        X21=Y2*Y1
507        X22=Y2*Y2
508        X23=Y2*Y3
509        X31=Y3*Y1
510        X32=Y3*Y2
511        X33=Y3*Y3
512
513        ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd)
514
515        m=1/(S0dd-ZZ)
516
517        N11=m*X11+Zaa
518        N12=m*X12+Zab
519        N13=m*X13+Zac
520        N21=m*X21+Zba
521        N22=m*X22+Zbb
522        N23=m*X23+Zbc
523        N31=m*X31+Zca
524        N32=m*X32+Zcb
525        N33=m*X33+Zcc
526
527        M11= N11*S0aa + N12*S0ab + N13*S0ac
528        M12= N11*S0ab + N12*S0bb + N13*S0bc
529        M13= N11*S0ac + N12*S0bc + N13*S0cc
530        M21= N21*S0aa + N22*S0ab + N23*S0ac
531        M22= N21*S0ab + N22*S0bb + N23*S0bc
532        M23= N21*S0ac + N22*S0bc + N23*S0cc
533        M31= N31*S0aa + N32*S0ab + N33*S0ac
534        M32= N31*S0ab + N32*S0bb + N33*S0bc
535        M33= N31*S0ac + N32*S0bc + N33*S0cc
536
537        DenQ1=1+M11-M12*M21+M22+M11*M22-M13*M31-M13*M22*M31
538        DenQ2=  M12*M23*M31+M13*M21*M32-M23*M32-M11*M23*M32+M33+M11*M33
539        DenQ3=  -M12*M21*M33+M22*M33+M11*M22*M33
540        DenQ=DenQ1+DenQ2+DenQ3
541
542        Q11= (1 + M22-M23*M32 + M33 + M22*M33)/DenQ
543        Q12= (-M12 + M13*M32 - M12*M33)/DenQ
544        Q13= (-M13 - M13*M22 + M12*M23)/DenQ
545        Q21= (-M21 + M23*M31 - M21*M33)/DenQ
546        Q22= (1 + M11 - M13*M31 + M33 + M11*M33)/DenQ
547        Q23= (M13*M21 - M23 - M11*M23)/DenQ
548        Q31= (-M31 - M22*M31 + M21*M32)/DenQ
549        Q32= (M12*M31 - M32 - M11*M32)/DenQ
550        Q33= (1 + M11 - M12*M21 + M22 + M11*M22)/DenQ
551
552        S11= Q11*S0aa + Q21*S0ab + Q31*S0ac
553        S12= Q12*S0aa + Q22*S0ab + Q32*S0ac
554        S13= Q13*S0aa + Q23*S0ab + Q33*S0ac
555        S14=-S11-S12-S13
556        S21= Q11*S0ba + Q21*S0bb + Q31*S0bc
557        S22= Q12*S0ba + Q22*S0bb + Q32*S0bc
558        S23= Q13*S0ba + Q23*S0bb + Q33*S0bc
559        S24=-S21-S22-S23
560        S31= Q11*S0ca + Q21*S0cb + Q31*S0cc
561        S32= Q12*S0ca + Q22*S0cb + Q32*S0cc
562        S33= Q13*S0ca + Q23*S0cb + Q33*S0cc
563        S34=-S31-S32-S33
564        S41=S14
565        S42=S24
566        S43=S34
567        S44=S11+S22+S33+2*S12+2*S13+2*S23
568
569        Nav=6.022045E23
570        Lad=(La/va-Ld/vd)*SQRT(Nav)
571        Lbd=(Lb/vb-Ld/vd)*SQRT(Nav)
572        Lcd=(Lc/vc-Ld/vd)*SQRT(Nav)
573
574        Int=Lad^2*S11+Lbd^2*S22+Lcd^2*S33+2*Lad*Lbd*S12+2*Lbd*Lcd*S23+2*Lad*Lcd*S13
575        Int=Int*scale + background
576
577//      print Q,Int,Nd,var[5],lCASE
578
579        return Int
580       
581End
582
583// this is all there is to the smeared calculation!
584Function SmearedRPAForm(s) :FitFunc
585
586        Struct ResSmearAAOStruct &s
587
588////the name of your unsmeared model is the first argument
589        Smear_Model_20(RPAForm,s.coefW,s.xW,s.yW,s.resW)
590
591        return(0)
592End
593
594//wrapper to calculate the smeared model as an AAO-Struct
595// fills the struct and calls the ususal function with the STRUCT parameter
596//
597// used only for the dependency, not for fitting
598//
599Function fSmearedRPAForm(coefW,yW,xW)
600        Wave coefW,yW,xW
601       
602        String str = getWavesDataFolder(yW,0)
603        String DF="root:"+str+":"
604       
605        WAVE resW = $(DF+str+"_res")
606       
607        STRUCT ResSmearAAOStruct fs
608        WAVE fs.coefW = coefW   
609        WAVE fs.yW = yW
610        WAVE fs.xW = xW
611        WAVE fs.resW = resW
612       
613        Variable err
614        err = SmearedRPAForm(fs)
615       
616        return (0)
617End
Note: See TracBrowser for help on using the repository browser.