source: trunk/User form factors for Irena/FlexibleCylinderForIrena.ipf

Last change on this file was 663, checked in by ilavsky, 4 years ago

Add User Form factors for Irena
  • Property svn:eol-style set to native
File size: 15.0 KB
Line 
1#pragma TextEncoding = "UTF-8"
2#pragma rtGlobals=3             // Use modern global access method and strict wave access.
3
4
5//NOTE: this code is copy of code from FexibleCylinder_v40.ipf from NIST Igor package, modified to use in Irena.
6// Below are original notes from the NIST package:
7//CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
8//CCCCCCCC
9//C      SUBROUTINE FOR THE CALCULATION OF THE SCATTERING FUNCTIONS
10//C      OF RODLIKE MICELLES.  METHODLOGY FOLLOWS THAT OF PEDERSEN AND
11//C      SCHURTENBERGER, MACORMOLECULES, VOL 29,PG 7602, 1996.
12//C      WITH EXCULDED VOLUME EFFECTS (METHOD 3)
13//
14// - copied directly from FORTRAN code supplied by Jan Pedersen
15//              SRK - 2002, but shows discontinuity at Qlb = 3.1
16//
17//  Jan 2006 - re-worked FORTRAN correcting typos in paper: now is smooth, but
18// the splicing is actually at Qlb = 2, which is not what the paper
19// says is to be done (and not from earlier models)
20//
21// July 2006 - now is CORRECT with Wei-Ren's changes to the code
22// Matlab code was not too difficult to convert to Igor (only a few hours...)
23//
24// June 2018 modified for Irena use by Jan Ilavsky, all glory for the code goes to Jan Pedersen and Steven Kline (NIST group). I did nothing useful... 
25
26//*************************************************************************************************
27//*************************************************************************************************
28//USE in Irena :
29//In Modeling II select User form factor
30//In panel put in "Name of ForFactor function this string:    IR1T_FlexExclVolCyl
31//In Panel put in Name of volume FF function this string:    IR1T_FlexExclVolCylVol
32//
33// Parameter 1 is the length of the cylinder
34// Parameter 2 is the Kuhn length
35// other parameters are not being used.
36//*************************************************************************************************
37//*************************************************************************************************
38
39
40//****   IR1T_FlexExclVolCyl ***** is simply the form factor, which normalizes to 1 at Q=0
41Function IR1T_FlexExclVolCyl(Qval,radius, par1,par2,par3,par4,par5)
42        variable Qval, radius, par1,par2,par3,par4,par5                                                                                         
43        //create coef wave needed for the parameters
44        Make/O/D/Free coef_fle = {1.,par1,par2,radius,2,1,0}
45        //this is the meaning of those parameters...
46        //make/o/t parameters_fle = {"scale","Contour Length (A)","Kuhn Length, b (A)","Radius (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","bkgd (cm^-1)"}
47        //Irena handles scaling and background on its own, so those need to be set to 1 and 0
48        // this wave also contains the radius, which is parameter Irena changes...
49        //Irena handles contrasts typically also, but in this case it may be easier to stick them here and then use 1 for Contrast in Irena panel.
50        //so par 1 is Contour length
51        //par 2 is Kuhn length
52        //SLD cylinder (A^-2) I took out the slds from calculations, use irena contrast
53        //SLD solvent (A^-2) I took out the slds from calcualtions, use irena contrast
54        //this will now calculate the Form factor for a combination of Q and R 
55        //note, at this moemnt (with removeal of contrasts, volume and other unit corrections in fFlexExclVolCyl below For factor gives 1 for Q=0, asi it shoudl do... 
56#if exists("FlexExclVolCylX")
57        //      flex *= Pi*rad*rad*L    //need to fix this by dividing by (Pi*radius*radius*par1)
58        //      flex *= cont^2                  //cont is set to 1 above, so this does not apply.
59        //      flex *= 1.0e8                   //also, fis this... Convert back to Form factor...
60        return sqrt(FlexExclVolCylX(coef_fle,Qval) / ((Pi*radius*radius*par1)/1e8) )/1e8
61#else
62        //I have cleaned this out of the form factor, no scaling needed...
63        return sqrt(IR1T_fFlexExclVolCyl(coef_fle,Qval))
64#endif
65        return(0)
66end
67//****   IR1T_FlexExclVolCylVol ***** is simply the volume of the cylinder
68Function IR1T_FlexExclVolCylVol(radius, par1,par2,par3,par4,par5)               //returns the sphere volume
69        variable radius, par1,par2,par3,par4,par5
70        // par1 = length, so it may be as simple as this:
71        return (pi*radius*radius*par1)
72end
73//*********************************************************************************************************
74//*********************************************************************************************************
75//*********************************************************************************************************
76//notes and edits is likely Steven Kline, NIST, he developed the code.
77//These are all work functions, nothing really useful here for users...
78static Function IR1T_fFlexExclVolCyl(ww,x)
79        Wave ww
80        Variable x
81
82        //nice names to the input params
83        //ww[0] = scale
84        //ww[1] = L [A]
85        //ww[2] = B [A]
86        //ww[3] = rad [A] cross-sectional radius
87        //ww[4] = sld cylinder
88        //ww[5] = sld solvent [A^-2]
89        //ww[6] = bkg [cm-1]
90        Variable scale,L,B,bkg,rad,qr,cont,sldc,slds
91       
92        scale = ww[0]
93        L = ww[1]
94        B = ww[2]
95        rad = ww[3]
96        sldc = ww[4]
97        slds = ww[5]
98        bkg = ww[6]
99        qr = x*rad              //used for cross section contribution only
100       
101        cont = sldc-slds
102       
103        Variable flex,crossSect
104        flex = IR1T_Sk_WR(x,L,B)
105   
106    if(qr == 0)
107            crossSect = 1
108    else
109        crossSect = (2*bessJ(1,qr)/qr)^2
110    endif
111       
112        //normalize form factor by multiplying by cylinder volume * cont^2
113   // then convert to cm-1 by multiplying by 10^8
114   // then scale = phi
115
116        flex *= crossSect
117        //flex *= Pi*rad*rad*L
118        //flex *= cont^2
119        //flex *= 1.0e-8
120//      print x, flex
121   return (scale*flex + bkg)
122End
123
124
125//////////////////WRC corrected code below
126// main function
127static function IR1T_Sk_WR(q,L,b)
128        Variable q,L,b
129        //
130        Variable p1,p2,p1short,p2short,q0,qconnect
131        Variable c,epsilon,ans,q0short,Sexvmodify
132       
133        p1 = 4.12
134        p2 = 4.42
135        p1short = 5.36
136        p2short = 5.62
137        q0 = 3.1
138        qconnect = q0/b
139        //     
140        q0short = max(1.9/sqrt(IR1T_Rgsquareshort(q,L,b)),3)
141       
142        //
143        if(L/b > 10)
144                C = 3.06/(L/b)^0.44
145                epsilon = 0.176
146        else
147                C = 1
148                epsilon = 0.170
149        endif
150        //
151       
152        if( L > 4*b ) // Longer Chains
153                if (q*b <= 3.1)
154                        //Modified by Yun on Oct. 15,
155                        Sexvmodify = IR1T_Sexvnew(q, L, b)
156                        ans = Sexvmodify + C * (4/15 + 7./(15*IR1T_u_WR(q,L,b)) - (11/15 + 7./(15*IR1T_u_WR(q,L,b)))*exp(-IR1T_u_WR(q,L,b)))*(b/L)
157                else //q(i)*b > 3.1
158                        ans = IR1T_a1long(q, L, b, p1, p2, q0)/((q*b)^p1) + IR1T_a2long(q, L, b, p1, p2, q0)/((q*b)^p2) + pi/(q*L)
159                endif
160        else //L <= 4*b Shorter Chains
161                if (q*b <= max(1.9/sqrt(IR1T_Rgsquareshort(q,L,b)),3) )
162                        if (q*b<=0.01)
163                                ans = 1 - IR1T_Rgsquareshort(q,L,b)*(q^2)/3
164                        else
165                                ans = IR1T_Sdebye1(q,L,b)
166                        endif
167                else    //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
168                        ans = IR1T_a1short(q,L,b,p1short,p2short,q0short)/((q*b)^p1short) + IR1T_a2short(q,L,b,p1short,p2short,q0short)/((q*b)^p2short) + pi/(q*L)
169                endif
170        endif
171       
172        return(ans)
173end
174
175//WR named this w (too generic)
176static Function IR1T_w_WR(x)
177    Variable x
178
179    //C4 = 1.523;
180    //C5 = 0.1477;
181    Variable yy
182    yy = 0.5*(1 + tanh((x - 1.523)/0.1477))
183
184    return (yy)
185end
186
187//
188static function IR1T_u1(q,L,b)
189    Variable q,L,b
190    Variable yy
191
192    yy = IR1T_Rgsquareshort(q,L,b)*q^2
193   
194    return yy
195end
196
197// was named u
198static function IR1T_u_WR(q,L,b)
199    Variable q,L,b
200   
201    Variable yy
202    yy = IR1T_Rgsquare(q,L,b)*(q^2)
203    return yy
204end
205
206
207
208//
209static function IR1T_Rgsquarezero(q,L,b)
210    Variable q,L,b
211   
212    Variable yy
213    yy = (L*b/6) * (1 - 1.5*(b/L) + 1.5*(b/L)^2 - 0.75*(b/L)^3*(1 - exp(-2*(L/b))))
214   
215    return yy
216end
217
218//
219static function IR1T_Rgsquareshort(q,L,b)
220    Variable q,L,b
221   
222    Variable yy
223    yy = IR1T_AlphaSquare(L/b) * IR1T_Rgsquarezero(q,L,b)
224   
225    return yy
226end
227
228//
229static function IR1T_Rgsquare(q,L,b)
230    Variable q,L,b
231   
232    Variable yy
233    yy = IR1T_AlphaSquare(L/b)*L*b/6
234   
235    return yy
236end
237
238//
239static function IR1T_AlphaSquare(x)
240    Variable x
241   
242    Variable yy
243    yy = (1 + (x/3.12)^2 + (x/8.67)^3)^(0.176/3)
244
245    return yy
246end
247
248//
249static function IR1T_miu(x)
250    Variable x
251   
252    Variable yy
253    yy = (1/8)*(9*x - 2 + 2*log(1 + x)/x)*exp(1/2.565*(1/x + (1 - 1/(x^2))*log(1 + x)))
254   
255    return yy
256end
257
258//
259static function IR1T_Sdebye(q,L,b)
260    Variable q,L,b
261   
262    Variable yy
263    yy = 2*(exp(-IR1T_u_WR(q,L,b)) + IR1T_u_WR(q,L,b) -1)/((IR1T_u_WR(q,L,b))^2)
264
265    return yy
266end
267
268//
269static function IR1T_Sdebye1(q,L,b)
270    Variable q,L,b
271   
272    Variable yy
273    yy = 2*(exp(-IR1T_u1(q,L,b)) + IR1T_u1(q,L,b) -1)/((IR1T_u1(q,L,b))^2)
274   
275    return yy
276end
277
278//
279static function IR1T_Sexv(q,L,b)
280    Variable q,L,b
281   
282    Variable yy,C1,C2,C3,miu,Rg2
283    C1=1.22
284    C2=0.4288
285    C3=-1.651
286    miu = 0.585
287
288    Rg2 = IR1T_Rgsquare(q,L,b)
289   
290    yy = (1 - IR1T_w_WR(q*sqrt(Rg2)))*IR1T_Sdebye(q,L,b) + IR1T_w_WR(q*sqrt(Rg2))*(C1*(q*sqrt(Rg2))^(-1/miu) +  C2*(q*sqrt(Rg2))^(-2/miu) +    C3*(q*sqrt(Rg2))^(-3/miu))
291   
292    return yy
293end
294
295// this must be WR modified version
296static function IR1T_Sexvnew(q,L,b)
297    Variable q,L,b
298   
299    Variable yy,C1,C2,C3,miu
300    C1=1.22
301    C2=0.4288
302    C3=-1.651
303    miu = 0.585
304
305    //calculating the derivative to decide on the corection (cutoff) term?
306    // I have modified this from WRs original code
307    Variable del=1.05,C_star2,Rg2
308    if( (IR1T_Sexv(q*del,L,b)-IR1T_Sexv(q,L,b))/(q*del - q) >= 0 )
309        C_star2 = 0
310    else
311        C_star2 = 1
312    endif
313
314    Rg2 = IR1T_Rgsquare(q,L,b)
315   
316    yy = (1 - IR1T_w_WR(q*sqrt(Rg2)))*IR1T_Sdebye(q,L,b) + C_star2*IR1T_w_WR(q*sqrt(Rg2))*(C1*(q*sqrt(Rg2))^(-1/miu) + C2*(q*sqrt(Rg2))^(-2/miu) + C3*(q*sqrt(Rg2))^(-3/miu))
317
318    return yy
319end
320
321
322
323// these are the messy ones
324static function IR1T_a2short(q, L, b, p1short, p2short, q0)
325    Variable q, L, b, p1short, p2short, q0
326   
327    Variable yy,Rg2_sh
328    Rg2_sh = IR1T_Rgsquareshort(q,L,b)
329   
330    Variable t1
331    t1 = ((q0^2*Rg2_sh)/b^2)
332   
333    //E is the number e
334    yy = ((-(1/(L*((p1short - p2short))*Rg2_sh^2)*((b*E^(-t1)*q0^(-4 + p2short)*((8*b^3*L - 8*b^3*E^t1*L - 2*b^3*L*p1short + 2*b^3*E^t1*L*p1short + 4*b*L*q0^2*Rg2_sh + 4*b*E^t1*L*q0^2*Rg2_sh - 2*b*E^t1*L*p1short*q0^2*Rg2_sh - E^t1*pi*q0^3*Rg2_sh^2 + E^t1*p1short*pi*q0^3*Rg2_sh^2)))))))
335         
336    return yy
337end
338
339//
340static function IR1T_a1short(q, L, b, p1short, p2short, q0)
341    Variable q, L, b, p1short, p2short, q0
342   
343    Variable yy,Rg2_sh
344    Rg2_sh = IR1T_Rgsquareshort(q,L,b)
345
346    Variable t1
347    t1 = ((q0^2*Rg2_sh)/b^2)
348   
349    yy = ((1/(L*((p1short - p2short))*Rg2_sh^2)*((b*E^(-t1)*q0^((-4) + p1short)*((8*b^3*L - 8*b^3*E^t1*L - 2*b^3*L*p2short + 2*b^3*E^t1*L*p2short + 4*b*L*q0^2*Rg2_sh + 4*b*E^t1*L*q0^2*Rg2_sh - 2*b*E^t1*L*p2short*q0^2*Rg2_sh - E^t1*pi*q0^3*Rg2_sh^2 + E^t1*p2short*pi*q0^3*Rg2_sh^2))))))
350       
351    return yy
352end
353
354// this one will be lots of trouble
355static function IR1T_a2long(q, L, b, p1, p2, q0)
356    variable q, L, b, p1, p2, q0
357
358    Variable yy,c1,c2,c3,c4,c5,miu,c,Rg2,rRg
359    Variable t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13
360
361    if( L/b > 10)
362        C = 3.06/(L/b)^0.44
363    else
364        C = 1
365    endif
366
367    C1 = 1.22
368    C2 = 0.4288
369    C3 = -1.651
370    C4 = 1.523
371    C5 = 0.1477
372    miu = 0.585
373
374    Rg2 = IR1T_Rgsquare(q,L,b)
375    t1 = (1/(b* p1*q0^((-1) - p1 - p2) - b*p2*q0^((-1) - p1 - p2)))
376    t2 = (b*C*(((-1*((14*b^3)/(15*q0^3*Rg2))) + (14*b^3*E^(-((q0^2*Rg2)/b^2)))/(15*q0^3*Rg2) + (2*E^(-((q0^2*Rg2)/b^2))*q0*((11/15 + (7*b^2)/(15*q0^2*Rg2)))*Rg2)/b)))/L
377    t3 = (sqrt(Rg2)*((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))*IR1T_sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2)/(2*C5)
378    t4 = (b^4*sqrt(Rg2)*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*IR1T_sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2)/(C5*q0^4*Rg2^2)
379    t5 = (2*b^4*(((2*q0*Rg2)/b - (2*E^(-((q0^2*Rg2)/b^2))*q0*Rg2)/b))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q0^4*Rg2^2)
380    t6 = (8*b^5*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q0^5*Rg2^2)
381    t7 = (((-((3*C3*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 3/miu))/miu)) - (2*C2*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 2/miu))/miu - (C1*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 1/miu))/miu))
382    t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))
383    t9 = (b*C*((4/15 - E^(-((q0^2*Rg2)/b^2))*((11/15 + (7*b^2)/(15*q0^2*Rg2))) + (7*b^2)/(15*q0^2*Rg2))))/L
384    t10 = (2*b^4*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q0^4*Rg2^2)
385 
386   
387    yy = ((-1*(t1* (((-q0^(-p1))*(((b^2*pi)/(L*q0^2) + t2 + t3 - t4 + t5 - t6 + 1/2*t7*t8)) - b*p1*q0^((-1) - p1)*(((-((b*pi)/(L*q0))) + t9 + t10 + 1/2*((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))))
388
389
390    return yy
391end
392
393//need to define this on my own
394static Function IR1T_sech_WR(x)
395        variable x
396       
397        return(1/cosh(x))
398end
399//
400static function IR1T_a1long(q, L, b, p1, p2, q0)
401    Variable q, L, b, p1, p2, q0
402   
403    Variable yy,c,c1,c2,c3,c4,c5,miu,Rg2,rRg
404    Variable t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16
405   
406    if( L/b > 10)
407        C = 3.06/(L/b)^0.44
408    else
409        C = 1
410    endif
411
412    C1 = 1.22
413    C2 = 0.4288
414    C3 = -1.651
415    C4 = 1.523
416    C5 = 0.1477
417    miu = 0.585
418
419    Rg2 = IR1T_Rgsquare(q,L,b)
420    t1 = (b*C*((4/15 - E^(-((q0^2*Rg2)/b^2))*((11/15 + (7*b^2)/(15*q0^2*Rg2))) + (7*b^2)/(15*q0^2*Rg2))))
421    t2 = (2*b^4*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))
422    t3 = ((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))
423    t4 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))
424    t5 = (1/(b*p1*q0^((-1) - p1 - p2) - b*p2*q0^((-1) - p1 - p2)))
425    t6 = (b*C*(((-((14*b^3)/(15*q0^3*Rg2))) + (14*b^3*E^(-((q0^2*Rg2)/b^2)))/(15*q0^3*Rg2) + (2*E^(-((q0^2*Rg2)/b^2))*q0*((11/15 + (7*b^2)/(15*q0^2*Rg2)))*Rg2)/b)))
426    t7 = (sqrt(Rg2)*((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))*IR1T_sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2)
427    t8 = (b^4*sqrt(Rg2)*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*IR1T_sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5)^2)
428    t9 = (2*b^4*(((2*q0*Rg2)/b - (2*E^(-((q0^2*Rg2)/b^2))*q0*Rg2)/b))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))
429    t10 = (8*b^5*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))
430    t11 = (((-((3*C3*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 3/miu))/miu)) - (2*C2*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 2/miu))/miu - (C1*sqrt(Rg2)*(((sqrt(Rg2)*q0)/b))^((-1) - 1/miu))/miu))
431    t12 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))
432    t13 = (b*C*((4/15 - E^(-((q0^2*Rg2)/b^2))*((11/15 + (7*b^2)/(15*q0^2* Rg2))) + (7*b^2)/(15*q0^2*Rg2))))
433    t14 = (2*b^4*(((-1) + E^(-((q0^2*Rg2)/b^2)) + (q0^2*Rg2)/b^2))*((1 + 1/2*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))
434    t15 = ((C3*(((sqrt(Rg2)*q0)/b))^((-3)/miu) + C2*(((sqrt(Rg2)*q0)/b))^((-2)/miu) + C1*(((sqrt(Rg2)*q0)/b))^((-1)/miu)))
435
436   
437    yy = (q0^p1*(((-((b*pi)/(L*q0))) +t1/L +t2/(q0^4*Rg2^2) + 1/2*t3*t4)) + (t5*((q0^(p1 - p2)*(((-q0^(-p1))*(((b^2*pi)/(L*q0^2) +t6/L +t7/(2*C5) -t8/(C5*q0^4*Rg2^2) +t9/(q0^4*Rg2^2) -t10/(q0^5*Rg2^2) + 1/2*t11*t12)) - b*p1*q0^((-1) - p1)*(((-((b*pi)/(L*q0))) +t13/L +t14/(q0^4*Rg2^2) + 1/2*t15*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))))
438
439   
440    return yy
441end
Note: See TracBrowser for help on using the repository browser.