[115d0f0] | 1 | double form_volume(void); |
---|
| 2 | |
---|
| 3 | double 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 | |
---|
| 9 | double 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 | |
---|
[e481a39] | 15 | void calculate_psfs(double qval, |
---|
[115d0f0] | 16 | double r2, double nf2, |
---|
| 17 | double aa, double phi, |
---|
| 18 | double *s11, double *s22, double *s12 |
---|
| 19 | ); |
---|
| 20 | |
---|
| 21 | double form_volume(void) |
---|
| 22 | { |
---|
| 23 | return 1.0; |
---|
| 24 | } |
---|
| 25 | |
---|
| 26 | double 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 |
---|
[e481a39] | 53 | calculate_psfs(q,r2,nf2,aa,phi,&psf11,&psf22,&psf12); |
---|
[115d0f0] | 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 | |
---|
| 91 | double 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 | |
---|
[e481a39] | 104 | void calculate_psfs(double qval, |
---|
[115d0f0] | 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 |
---|
[e481a39] | 113 | double s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4; |
---|
[115d0f0] | 114 | double a1,a2i,a2,b1,b2,b12,gm1,gm12; |
---|
[e481a39] | 115 | double yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3; |
---|
[115d0f0] | 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; |
---|
[e481a39] | 120 | // s1 = aa*s2; why is this never used? check original paper? |
---|
[115d0f0] | 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 | |
---|
[e481a39] | 190 | return; |
---|
[115d0f0] | 191 | } |
---|
| 192 | |
---|