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

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

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