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

Last change on this file since 127 was 127, checked in by srkline, 16 years ago

LOTS of changes to the analysis ipf files:

-- see sphere.ipf for the simplest example of the changes --

  • #pragma Igor 6
  • #if directive to look for XOP
  • AAO unsmeared functions
  • STRUCT functions for smearing (also AAO)
  • new wrappers for dependencies to struct functions

(2006 models have NOT been completed yet, only the old models)

  • loading data files into data folders (PlotUtils?) + some streamlining of the loaders
  • Smear_Model_N is now AAO + some streamlining of the quadrature code

-- SHS and SW structure factor XOPs are crashing (need DP wave, I may have old XOP)
-- this breaks fitting of the smeared models until wrappers can be devised
-- all packages will be broken due to the new data folder structure
-- multiple instances of functions will now cause problems (MSA)
-- RPA model is a problem with its odd functional form (extra wave)

-- lots of other carnage to follow as the bugs and typos are shaken out

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