source: sans/Dev/trunk/NCNR_User_Procedures/SANS/Analysis/Models/SmearedRPA_v40.ipf @ 400

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

SmearedRPA was not properly data folder aware, and would balk if a smeared version was plotted before an unsmeared version. Functions now look in the proper data folders.

File size: 16.9 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)
69       
70        AddModelToStrings("RPAForm","coef_rpa","rpa")
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 containing 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_coef_rpa = {5,5,-.0004,1,0}
107                make/o/t smear_parameters_rpa = {"Seg. Length bc","Seg. Length bd","Chi Param. Kcd","scale","Background"}
108                Edit smear_parameters_rpa,smear_coef_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_coef_rpa = {5,5,5,-.0004,-.0004,-.0004,1,0}
116                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"}
117                Edit smear_parameters_rpa,smear_coef_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_coef_rpa = {5,5,5,5,-.0004,-.0004,-.0004,-.0004,-.0004,-.0004,1,0}
125                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"}
126                Edit smear_parameters_rpa,smear_coef_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:
144        AddModelToStrings("SmearedRPAForm","smear_coef_rpa","rpa")
145End
146
147///////////////////////////////////////////////////////////////
148//AAO version, uses XOP if available
149// simply calls the original single point calculation with
150// a wave assignment (this will behave nicely if given point ranges)
151//
152// there is no RPAFormX due to the extra inputValues wave of information...
153//
154Function RPAForm(cw,yw,xw) : FitFunc
155        Wave cw,yw,xw
156       
157#if exists("RPAFormX")
158        yw = RPAFormX(cw,xw)
159#else
160        yw = fRPAForm(cw,xw)
161#endif
162        return(0)
163End
164
165
166Function fRPAForm(w,x) : FitFunc
167        Wave w
168        Variable x
169       
170//      print (GetWavesDataFolder(w,1)+"inputvalues")
171
172        Wave var=$(GetWavesDataFolder(w,1)+"inputvalues")
173        Nvar lCASE=$(GetWavesDataFolder(w,1)+"gCASE")
174//      print lCASE
175//      Variable lCASE
176//      RANDOM PHASE APPROXIMATION FOR A FOUR-BLOCK COPOLYMER A-B-C-D
177//      THIS FORMALISM APPLIES TO MULTICOMPONENT POLYMER MIXTURES
178//      IN THE HOMOGENEOUS (MIXED) PHASE REGION ONLY.
179
180//      THIS GENERAL CASE INCLUDES A MULTITUDE OF (TEN) SPECIAL CASES.
181
182//      HERE ARE THE VARIOUS CASES COVERED:
183//      CASE 1: C/D BINARY MIXTURE OF HOMOPOLYMERS
184//      CASE 2: C-D DIBLOCK COPOLYMER
185//      CASE 3: B/C/D TERNARY MIXTURE OF HOMOPOLYMERS
186//      CASE 4: B/C-D MIXTURE OF HOMOPOLYMER B AND DIBLOCK COPOLYMER C-D
187//      CASE 5: B-C-D TRIBLOCK COPOLYMER
188//      CASE 6: A/B/C/D QUATERNARY MIXTURE OF HOMOPOLYMERS
189//      CASE 7: A/B/C-D MIXTURE OF TWO HOMOPOLYMERS A/B AND A DIBLOCK C-D
190//      CASE 8: A/B-C-D MIXTURE OF A HOMOPOLYMER A AND A TRIBLOCK B-C-D
191//      CASE 9: A-B/C-D MIXTURE OF TWO DIBLOCK COPOLYMERS A-B AND C-D
192//      CASE 10: A-B-C-D FOUR-BLOCK COPOLYMER
193
194//      B. HAMMOUDA, NIST, JULY 1998
195
196        Variable Na,Nb,Nc,Nd,Nab,Nac,Nad,Nba,Nbc,Nbd,Nca,Ncb,Ncd
197        Variable Phia,Phib,Phic,Phid,Phiab,Phiac,Phiad
198        Variable Phiba,Phibc,Phibd,Phica,Phicb,Phicd,Phida,Phidb,Phidc
199        Variable va,vb,vc,vd,vab,vac,vad,vba,vbc,vbd,vca,vcb,vcd,vda,vdb,vdc
200        Variable m
201        Variable ba,bb,bc,bd
202        Variable Q
203        Variable Xa,Xb,Xc,Xd
204        Variable Paa,S0aa,Pab,S0ab,Pac,S0ac,Pad,S0ad
205        Variable Pba,S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd
206        Variable Pca,S0ca,Pcb,S0cb,Pcc,S0cc,Pcd,S0cd
207        Variable Pda,S0da,Pdb,S0db,Pdc,S0dc,Pdd,S0dd
208        Variable Kaa,Kab,Kac,Kad,Kba,Kbb,Kbc,Kbd
209        Variable Kca,Kcb,Kcc,Kcd,Kda,Kdb,Kdc,Kdd
210        Variable Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc
211        Variable DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33
212        Variable Y1,Y2,Y3,X11,X12,X13,X21,X22,X23,X31,X32,X33
213        Variable ZZ,DenQ1,DenQ2,DenQ3,DenQ,Q11,Q12,Q13,Q21,Q22,Q23,Q31,Q32,Q33
214        Variable N11,N12,N13,N21,N22,N23,N31,N32,N33
215        Variable M11,M12,M13,M21,M22,M23,M31,M32,M33
216        Variable S11,S12,S13,S14,S21,S22,S23,S24
217        Variable S31,S32,S33,S34,S41,S42,S43,S44
218        Variable La,Lb,Lc,Ld,Lad,Lbd,Lcd,Nav,Int
219        Variable scale
220        Variable Background
221        Variable ii=0   //to fool curve fitting dialog to let you pick the coefficient wave
222
223        Na=1000
224        Nb=1000
225        Nc=1000
226        Nd=1000 //DEGREE OF POLYMERIZATION     
227        Phia=0.25
228        Phib=0.25
229        Phic=0.25
230        Phid=0.25       //VOL FRACTION
231        Kab=-.0004
232        Kac=-.0004
233        Kad=-.0004      //CHI PARAM
234        Kbc=-.0004
235        Kbd=-.0004
236        Kcd=-.0004     
237        La=1.E-12
238        Lb=1.E-12
239        Lc=1.E-12
240        Ld=0.E-12       //SCATT. LENGTH
241        va=100
242        vb=100
243        vc=100
244        vd=100          //SPECIFIC VOLUME
245        ba=5
246        bb=5
247        bc=5
248        bd=5            //SEGMENT LENGTH
249
250        IF (lCASE <= 2)
251                Phia=0.0000001
252                Phib=0.0000001
253                Phic=0.5
254                Phid=0.5
255                Nc=var[0]
256                Phic=var[1]
257                vc=var[2]
258                Lc=var[3]
259                bc=w[ii+0]
260                Nd=var[4]
261                Phid=var[5]
262                vd=var[6]
263                Ld=var[7]
264                bd=w[ii+1]
265                Kcd=w[ii+2]
266                scale=w[ii+3]   
267                background=w[ii+4]
268        ENDIF
269
270        IF ((lCASE > 2) %& (lCASE <= 5))
271                Phia=0.0000001
272                Phib=0.333333
273                Phic=0.333333
274                Phid=0.333333
275                Nb=var[0]
276                Phib=var[1]
277                vb=var[2]
278                Lb=var[3]
279                bb=w[ii+0]
280                Nc=var[4]
281                Phic=var[5]
282                vc=var[6]
283                Lc=var[7]
284                bc=w[ii+1]
285                Nd=var[8]
286                Phid=var[9]
287                vd=var[10]
288                Ld=var[11]
289                bd=w[ii+2]
290                Kbc=w[ii+3]
291                Kbd=w[ii+4]
292                Kcd=w[ii+5]
293                scale=w[ii+6]   
294                background=w[ii+7]
295        ENDIF
296
297        IF (lCASE > 5)
298                Phia=0.25
299                Phib=0.25
300                Phic=0.25
301                Phid=0.25
302                Na=var[0]
303                Phia=var[1]
304                va=var[2]
305                La=var[3]
306                ba=w[ii+0]
307                Nb=var[4]
308                Phib=var[5]
309                vb=var[6]
310                Lb=var[7]
311                bb=w[ii+1]
312                Nc=var[8]
313                Phic=var[9]
314                vc=var[10]
315                Lc=var[11]
316                bc=w[ii+2]
317                Nd=var[12]
318                Phid=var[13]
319                vd=var[14]
320                Ld=var[15]
321                bd=w[ii+3]
322                Kab=w[ii+4]
323                Kac=w[ii+5]
324                Kad=w[ii+6]
325                Kbc=w[ii+7]
326                Kbd=w[ii+8]
327                Kcd=w[ii+9]
328                scale=w[ii+10] 
329                background=w[ii+11]
330        ENDIF
331
332        Nab=(Na*Nb)^(0.5)
333        Nac=(Na*Nc)^(0.5)
334        Nad=(Na*Nd)^(0.5)
335        Nbc=(Nb*Nc)^(0.5)
336        Nbd=(Nb*Nd)^(0.5)
337        Ncd=(Nc*Nd)^(0.5)
338
339        vab=(va*vb)^(0.5)
340        vac=(va*vc)^(0.5)
341        vad=(va*vd)^(0.5)
342        vbc=(vb*vc)^(0.5)
343        vbd=(vb*vd)^(0.5)
344        vcd=(vc*vd)^(0.5)
345
346        Phiab=(Phia*Phib)^(0.5)
347        Phiac=(Phia*Phic)^(0.5)
348        Phiad=(Phia*Phid)^(0.5)
349        Phibc=(Phib*Phic)^(0.5)
350        Phibd=(Phib*Phid)^(0.5)
351        Phicd=(Phic*Phid)^(0.5)
352
353        Q=x
354        Xa=Q^2*ba^2*Na/6
355        Xb=Q^2*bb^2*Nb/6
356        Xc=Q^2*bc^2*Nc/6
357        Xd=Q^2*bd^2*Nd/6
358
359        Paa=2*(Exp(-Xa)-1+Xa)/Xa^2
360        S0aa=Na*Phia*va*Paa
361        Pab=((1-Exp(-Xa))/Xa)*((1-Exp(-Xb))/Xb)
362        S0ab=(Phiab*vab*Nab)*Pab
363        Pac=((1-Exp(-Xa))/Xa)*Exp(-Xb)*((1-Exp(-Xc))/Xc)
364        S0ac=(Phiac*vac*Nac)*Pac
365        Pad=((1-Exp(-Xa))/Xa)*Exp(-Xb-Xc)*((1-Exp(-Xd))/Xd)
366        S0ad=(Phiad*vad*Nad)*Pad
367
368        S0ba=S0ab
369        Pbb=2*(Exp(-Xb)-1+Xb)/Xb^2
370        S0bb=Nb*Phib*vb*Pbb
371        Pbc=((1-Exp(-Xb))/Xb)*((1-Exp(-Xc))/Xc)
372        S0bc=(Phibc*vbc*Nbc)*Pbc
373        Pbd=((1-Exp(-Xb))/Xb)*Exp(-Xc)*((1-Exp(-Xd))/Xd)
374        S0bd=(Phibd*vbd*Nbd)*Pbd
375
376        S0ca=S0ac
377        S0cb=S0bc
378        Pcc=2*(Exp(-Xc)-1+Xc)/Xc^2
379        S0cc=Nc*Phic*vc*Pcc
380        Pcd=((1-Exp(-Xc))/Xc)*((1-Exp(-Xd))/Xd)
381        S0cd=(Phicd*vcd*Ncd)*Pcd
382
383        S0da=S0ad
384        S0db=S0bd
385        S0dc=S0cd
386        Pdd=2*(Exp(-Xd)-1+Xd)/Xd^2
387        S0dd=Nd*Phid*vd*Pdd
388
389        IF(lCASE == 1)
390                S0aa=0.000001
391                S0ab=0.000002           
392                S0ac=0.000003           
393                S0ad=0.000004           
394                S0bb=0.000005           
395                S0bc=0.000006           
396                S0bd=0.000007           
397                S0cd=0.000008
398        ENDIF
399
400        IF(lCASE == 2)
401                S0aa=0.000001
402                S0ab=0.000002           
403                S0ac=0.000003           
404                S0ad=0.000004           
405                S0bb=0.000005           
406                S0bc=0.000006           
407                S0bd=0.000007           
408        ENDIF
409       
410        IF(lCASE == 3)
411                S0aa=0.000001           
412                S0ab=0.000002           
413                S0ac=0.000003           
414                S0ad=0.000004           
415                S0bc=0.000005           
416                S0bd=0.000006           
417                S0cd=0.000007           
418        ENDIF
419
420        IF(lCASE == 4)
421                S0aa=0.000001           
422                S0ab=0.000002           
423                S0ac=0.000003           
424                S0ad=0.000004           
425                S0bc=0.000005           
426                S0bd=0.000006           
427        ENDIF
428
429        IF(lCASE == 5)
430                S0aa=0.000001           
431                S0ab=0.000002           
432                S0ac=0.000003           
433                S0ad=0.000004           
434        ENDIF
435
436        IF(lCASE == 6)
437                S0ab=0.000001           
438                S0ac=0.000002           
439                S0ad=0.000003           
440                S0bc=0.000004           
441                S0bd=0.000005           
442                S0cd=0.000006
443        ENDIF
444
445        IF(lCASE == 7)
446                S0ab=0.000001           
447                S0ac=0.000002           
448                S0ad=0.000003           
449                S0bc=0.000004           
450                S0bd=0.000005           
451        ENDIF
452
453        IF(lCASE == 8)
454                S0ab=0.000001           
455                S0ac=0.000002           
456                S0ad=0.000003           
457        ENDIF
458
459        IF(lCASE == 9)
460                S0ac=0.000001           
461                S0ad=0.000002           
462                S0bc=0.000003           
463                S0bd=0.000004
464        ENDIF
465
466        S0ba=S0ab               
467        S0ca=S0ac               
468        S0cb=S0bc               
469        S0da=S0ad               
470        S0db=S0bd               
471        S0dc=S0cd               
472
473        Kaa=0
474        Kbb=0
475        Kcc=0
476        Kdd=0
477
478        Kba=Kab
479        Kca=Kac
480        Kcb=Kbc
481        Kda=Kad
482        Kdb=Kbd
483        Kdc=Kcd
484
485        Zaa=Kaa-Kad-Kad 
486        Zab=Kab-Kad-Kbd 
487        Zac=Kac-Kad-Kcd
488        Zba=Kba-Kbd-Kad 
489        Zbb=Kbb-Kbd-Kbd 
490        Zbc=Kbc-Kbd-Kcd
491        Zca=Kca-Kcd-Kad 
492        Zcb=Kcb-Kcd-Kbd 
493        Zcc=Kcc-Kcd-Kcd
494
495        DenT=(-(S0ac*S0bb*S0ca) + S0ab*S0bc*S0ca + S0ac*S0ba*S0cb - S0aa*S0bc*S0cb - S0ab*S0ba*S0cc + S0aa*S0bb*S0cc)
496
497        T11= (-(S0bc*S0cb) + S0bb*S0cc)/DenT
498        T12= (S0ac*S0cb - S0ab*S0cc)/DenT
499        T13= (-(S0ac*S0bb) + S0ab*S0bc)/DenT
500        T21= (S0bc*S0ca - S0ba*S0cc)/DenT
501        T22= (-(S0ac*S0ca) + S0aa*S0cc)/DenT
502        T23= (S0ac*S0ba - S0aa*S0bc)/DenT
503        T31= (-(S0bb*S0ca) + S0ba*S0cb)/DenT
504        T32= (S0ab*S0ca - S0aa*S0cb)/DenT
505        T33= (-(S0ab*S0ba) + S0aa*S0bb)/DenT
506
507        Y1=T11*S0ad+T12*S0bd+T13*S0cd+1
508        Y2=T21*S0ad+T22*S0bd+T23*S0cd+1
509        Y3=T31*S0ad+T32*S0bd+T33*S0cd+1
510
511        X11=Y1*Y1
512        X12=Y1*Y2
513        X13=Y1*Y3
514        X21=Y2*Y1
515        X22=Y2*Y2
516        X23=Y2*Y3
517        X31=Y3*Y1
518        X32=Y3*Y2
519        X33=Y3*Y3
520
521        ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd)
522
523        m=1/(S0dd-ZZ)
524
525        N11=m*X11+Zaa
526        N12=m*X12+Zab
527        N13=m*X13+Zac
528        N21=m*X21+Zba
529        N22=m*X22+Zbb
530        N23=m*X23+Zbc
531        N31=m*X31+Zca
532        N32=m*X32+Zcb
533        N33=m*X33+Zcc
534
535        M11= N11*S0aa + N12*S0ab + N13*S0ac
536        M12= N11*S0ab + N12*S0bb + N13*S0bc
537        M13= N11*S0ac + N12*S0bc + N13*S0cc
538        M21= N21*S0aa + N22*S0ab + N23*S0ac
539        M22= N21*S0ab + N22*S0bb + N23*S0bc
540        M23= N21*S0ac + N22*S0bc + N23*S0cc
541        M31= N31*S0aa + N32*S0ab + N33*S0ac
542        M32= N31*S0ab + N32*S0bb + N33*S0bc
543        M33= N31*S0ac + N32*S0bc + N33*S0cc
544
545        DenQ1=1+M11-M12*M21+M22+M11*M22-M13*M31-M13*M22*M31
546        DenQ2=  M12*M23*M31+M13*M21*M32-M23*M32-M11*M23*M32+M33+M11*M33
547        DenQ3=  -M12*M21*M33+M22*M33+M11*M22*M33
548        DenQ=DenQ1+DenQ2+DenQ3
549
550        Q11= (1 + M22-M23*M32 + M33 + M22*M33)/DenQ
551        Q12= (-M12 + M13*M32 - M12*M33)/DenQ
552        Q13= (-M13 - M13*M22 + M12*M23)/DenQ
553        Q21= (-M21 + M23*M31 - M21*M33)/DenQ
554        Q22= (1 + M11 - M13*M31 + M33 + M11*M33)/DenQ
555        Q23= (M13*M21 - M23 - M11*M23)/DenQ
556        Q31= (-M31 - M22*M31 + M21*M32)/DenQ
557        Q32= (M12*M31 - M32 - M11*M32)/DenQ
558        Q33= (1 + M11 - M12*M21 + M22 + M11*M22)/DenQ
559
560        S11= Q11*S0aa + Q21*S0ab + Q31*S0ac
561        S12= Q12*S0aa + Q22*S0ab + Q32*S0ac
562        S13= Q13*S0aa + Q23*S0ab + Q33*S0ac
563        S14=-S11-S12-S13
564        S21= Q11*S0ba + Q21*S0bb + Q31*S0bc
565        S22= Q12*S0ba + Q22*S0bb + Q32*S0bc
566        S23= Q13*S0ba + Q23*S0bb + Q33*S0bc
567        S24=-S21-S22-S23
568        S31= Q11*S0ca + Q21*S0cb + Q31*S0cc
569        S32= Q12*S0ca + Q22*S0cb + Q32*S0cc
570        S33= Q13*S0ca + Q23*S0cb + Q33*S0cc
571        S34=-S31-S32-S33
572        S41=S14
573        S42=S24
574        S43=S34
575        S44=S11+S22+S33+2*S12+2*S13+2*S23
576
577        Nav=6.022045E23
578        Lad=(La/va-Ld/vd)*SQRT(Nav)
579        Lbd=(Lb/vb-Ld/vd)*SQRT(Nav)
580        Lcd=(Lc/vc-Ld/vd)*SQRT(Nav)
581
582        Int=Lad^2*S11+Lbd^2*S22+Lcd^2*S33+2*Lad*Lbd*S12+2*Lbd*Lcd*S23+2*Lad*Lcd*S13
583        Int=Int*scale + background
584
585//      print Q,Int,Nd,var[5],lCASE
586
587        return Int
588       
589End
590
591// this is all there is to the smeared calculation!
592Function SmearedRPAForm(s) :FitFunc
593
594        Struct ResSmearAAOStruct &s
595
596////the name of your unsmeared model is the first argument
597        Smear_Model_20(RPAForm,s.coefW,s.xW,s.yW,s.resW)
598
599        return(0)
600End
601
602//wrapper to calculate the smeared model as an AAO-Struct
603// fills the struct and calls the ususal function with the STRUCT parameter
604//
605// used only for the dependency, not for fitting
606//
607Function fSmearedRPAForm(coefW,yW,xW)
608        Wave coefW,yW,xW
609       
610        String str = getWavesDataFolder(yW,0)
611        String DF="root:"+str+":"
612       
613        WAVE resW = $(DF+str+"_res")
614       
615        STRUCT ResSmearAAOStruct fs
616        WAVE fs.coefW = coefW   
617        WAVE fs.yW = yW
618        WAVE fs.xW = xW
619        WAVE fs.resW = resW
620       
621        Variable err
622        err = SmearedRPAForm(fs)
623       
624        return (0)
625End
Note: See TracBrowser for help on using the repository browser.