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

Last change on this file since 252 was 252, checked in by ajj, 15 years ago

Renaming files with _v40 - see trac ticket #108

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