source: sasmodels/sasmodels/models/binary_hard_sphere.c @ 58c3367

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

cleanup of binary hard sphere model - now passes tests and docs works.
also removed some compilser warnings.

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