source: sasmodels/sasmodels/models/binary_hard_sphere.c @ 115d0f0

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 115d0f0 was 115d0f0, checked in by butler, 8 years ago

conversion of binary hard sphere model. Note that the model is badly
behaved at low (and sometimes high) qr. It may be possible to improve
on this by working on the calculations of the partial structure factors
sf11, sf12, sf22 which are a mass of arithmatic at present. However the
instability is outside the normal useful range. SasView? sasmodels and
samsmodels package do not agree well (delta of 1 or 2) outside the
"good" range.

  • Property mode set to 100644
File size: 5.7 KB
Line 
1double form_volume(void);
2
3double Iq(double q, 
4    double lg_radius, double sm_radius,
5    double lg_vol_frac, double sm_vol_frac,
6    double lg_sld, double sm_sld, double solvent_sld
7    );
8   
9double Iqxy(double qx, double qy,
10    double lg_radius, double sm_radius,
11    double lg_vol_frac, double sm_vol_frac,
12    double lg_sld, double sm_sld, double solvent_sld
13    );
14
15int calculate_psfs(double qval,
16    double r2, double nf2,
17    double aa, double phi,
18    double *s11, double *s22, double *s12
19    );
20
21double form_volume(void)
22{
23    return 1.0;
24}
25
26double Iq(double q,
27    double lg_radius, double sm_radius,
28    double lg_vol_frac, double sm_vol_frac,
29    double lg_sld, double sm_sld, double solvent_sld)
30   
31{
32    double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten;       //my local names
33    double psf11,psf12,psf22;
34    double phi1,phi2,phr,a3;
35    double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2;
36    int err;
37   
38    r2 = lg_radius;
39    r1 = sm_radius;
40    phi2 = lg_vol_frac;
41    phi1 = sm_vol_frac;
42    rho2 = lg_sld;
43    rho1 = sm_sld;
44    rhos = solvent_sld;
45   
46   
47    phi = phi1 + phi2;
48    aa = r1/r2;
49    //calculate the number fraction of larger spheres (eqn 2 in reference)
50    a3=aa*aa*aa;
51    phr=phi2/phi;
52    nf2 = phr*a3/(1.0-phr+phr*a3);
53    // calculate the PSF's here
54    err = calculate_psfs(q,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
55   
56    // /* do form factor calculations  */
57   
58    v1 = 4.0*M_PI/3.0*r1*r1*r1;
59    v2 = 4.0*M_PI/3.0*r2*r2*r2;
60   
61    n1 = phi1/v1;
62    n2 = phi2/v2;
63   
64    qr1 = r1*q;
65    qr2 = r2*q;
66
67    //if (qr1 == 0){
68        //sc1 = 1.0/3.0;
69    //}else{
70        //sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
71    //}
72    //if (qr2 == 0){
73        //sc2 = 1.0/3.0;
74    //}else{
75        //sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
76    //}
77    sc1 = sph_j1c(qr1);
78    sc2 = sph_j1c(qr2);
79    b1 = r1*r1*r1*(rho1-rhos)*4.0/3.0*M_PI*sc1;
80    b2 = r2*r2*r2*(rho2-rhos)*4.0/3.0*M_PI*sc2;
81    inten = n1*b1*b1*psf11;
82    inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
83    inten += n2*b2*b2*psf22;
84    ///* convert I(1/A) to (1/cm)  */
85    inten *= 1.0e8;
86    ///*convert rho^2 in 10^-6A to A*/
87    inten *= 1.0e-12;   
88    return(inten);
89}
90
91
92double Iqxy(double qx, double qy,
93    double lg_radius, double sm_radius,
94    double lg_vol_frac, double sm_vol_frac,
95    double lg_sld, double sm_sld, double solvent_sld)
96   
97{
98    double q = sqrt(qx*qx + qy*qy);
99    return Iq(q,
100        lg_radius, sm_radius,
101        lg_vol_frac, sm_vol_frac,
102        lg_sld, sm_sld, solvent_sld);
103}
104
105int calculate_psfs(double qval,
106    double r2, double nf2,
107    double aa, double phi,
108    double *s11, double *s22, double *s12)
109
110{
111    //  variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
112   
113    //   calculate constant terms
114    double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
115    double a1,a2i,a2,b1,b2,b12,gm1,gm12;
116    double err=0.0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
117    double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
118    double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
119   
120    s2 = 2.0*r2;
121    s1 = aa*s2;
122    v = phi;
123    a3 = aa*aa*aa;
124    v1=((1.-nf2)*a3/(nf2+(1.-nf2)*a3))*v;
125    v2=(nf2/(nf2+(1.-nf2)*a3))*v;
126    g11=((1.+.5*v)+1.5*v2*(aa-1.))/(1.-v)/(1.-v);
127    g22=((1.+.5*v)+1.5*v1*(1./aa-1.))/(1.-v)/(1.-v);
128    g12=((1.+.5*v)+1.5*(1.-aa)*(v1-v2)/(1.+aa))/(1.-v)/(1.-v);
129    wmv = 1/(1.-v);
130    wmv3 = wmv*wmv*wmv;
131    wmv4 = wmv*wmv3;
132    a1=3.*wmv4*((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2))) + ((v1+a3*v2)*(1.+2.*v)+(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)-3.*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
133    a2i=((v1+a3*v2)*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*3*wmv4 + ((v1+a3*v2)*(1.+2.*v)+a3*(1.+v+v*v)-3.*v1*v2*(1.-aa)*(1.-aa)*aa-3.*v1*(1.-aa)*(1.-aa)*(1.+v1+aa*(1.+v2)))*wmv3;
134    a2=a2i/a3;
135    b1=-6.*(v1*g11*g11+.25*v2*(1.+aa)*(1.+aa)*aa*g12*g12);
136    b2=-6.*(v2*g22*g22+.25*v1/a3*(1.+aa)*(1.+aa)*g12*g12);
137    b12=-3.*aa*(1.+aa)*(v1*g11/aa/aa+v2*g22)*g12;
138    gm1=(v1*a1+a3*v2*a2)*.5;
139    gm12=2.*gm1*(1.-aa)/aa;
140    //c 
141    //c   calculate the direct correlation functions and print results
142    //c
143    //  do 20 j=1,npts
144   
145    yy=qval*s2;
146    //c   calculate direct correlation functions
147    //c   ----c11
148    ay=aa*yy;
149    ay2 = ay*ay;
150    ay3 = ay*ay*ay;
151    t1=a1*(sin(ay)-ay*cos(ay));
152    t2=b1*(2.*ay*sin(ay)-(ay2-2.)*cos(ay)-2.)/ay;
153    t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
154    f11=24.*v1*(t1+t2+t3)/ay3;
155   
156    //c ------c22
157    y2=yy*yy;
158    y3=yy*y2;
159    tt1=a2*(sin(yy)-yy*cos(yy));
160    tt2=b2*(2.*yy*sin(yy)-(y2-2.)*cos(yy)-2.)/yy;
161    tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
162    f22=24.*v2*(tt1+tt2+tt3)/y3;
163   
164    //c   -----c12
165    yl=.5*yy*(1.-aa);
166    yl3=yl*yl*yl;
167    wma3 = (1.-aa)*(1.-aa)*(1.-aa);
168    y1=aa*yy;
169    y13 = y1*y1*y1;
170    ttt1=3.*wma3*v*sqrt(nf2)*sqrt(1.-nf2)*a1*(sin(yl)-yl*cos(yl))/((nf2+(1.-nf2)*a3)*yl3);
171    t21=b12*(2.*y1*cos(y1)+(y1*y1-2.)*sin(y1));
172    t22=gm12*((3.*y1*y1-6.)*cos(y1)+(y1*y1*y1-6.*y1)*sin(y1)+6.)/y1;
173    t23=gm1*((4.*y13-24.*y1)*cos(y1)+(y13*y1-12.*y1*y1+24.)*sin(y1))/(y1*y1);
174    t31=b12*(2.*y1*sin(y1)-(y1*y1-2.)*cos(y1)-2.);
175    t32=gm12*((3.*y1*y1-6.)*sin(y1)-(y1*y1*y1-6.*y1)*cos(y1))/y1;
176    t33=gm1*((4.*y13-24.*y1)*sin(y1)-(y13*y1-12.*y1*y1+24.)*cos(y1)+24.)/(y1*y1);
177    t41=cos(yl)*((sin(y1)-y1*cos(y1))/(y1*y1) + (1.-aa)/(2.*aa)*(1.-cos(y1))/y1);
178    t42=sin(yl)*((cos(y1)+y1*sin(y1)-1.)/(y1*y1) + (1.-aa)/(2.*aa)*sin(y1)/y1);
179    ttt2=sin(yl)*(t21+t22+t23)/(y13*y1);
180    ttt3=cos(yl)*(t31+t32+t33)/(y13*y1);
181    ttt4=a1*(t41+t42)/y1;
182    f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
183   
184    c11=f11;
185    c22=f22;
186    c12=f12;
187    *s11=1./(1.+c11-(c12)*c12/(1.+c22));
188    *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
189    *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));   
190   
191    return(err);
192}
193
Note: See TracBrowser for help on using the repository browser.