# Changeset 632 for sans/XOP_Dev/SANSAnalysis/lib/libSphere.c

Ignore:
Timestamp:
Mar 1, 2010 1:57:12 PM (12 years ago)
Message:

incorporating Jae-Hie's changes to the library functions to explicitly handle the cases at qr==0. Also a lot of changes to explicitly declare constant values as floating point, rather than chance the interpretation as integer. (probably not necessary)

File:
1 edited

Unmodified
Added
Removed
• ## sans/XOP_Dev/SANSAnalysis/lib/libSphere.c

 r593 { double scale,radius,delrho,bkg,sldSph,sldSolv;          //my local names double bes,f,vol,f2,pi; double bes,f,vol,f2,pi,qr; pi = 4.0*atan(1.0); scale = dp[0]; delrho = sldSph - sldSolv; //handle q==0 separately if(q==0){ f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*scale*1.0e8 + bkg; return(f); //handle qr==0 separately qr = q*radius; if(qr == 0.0){ bes = 1.0; }else{ bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); } bes = 3.0*(sin(q*radius)-q*radius*cos(q*radius))/(q*q*q)/(radius*radius*radius); vol = 4.0*pi/3.0*radius*radius*radius; f = vol*bes*delrho;             // [=]  f = vol*bes*delrho;             // [=] A-1 // normalize to single particle volume, convert to 1/cm f2 = f * f / vol * 1.0e8;               // [=] 1/cm qr=x*rcore; contr = rhocore-rhoshel; bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); if(qr == 0.0){ bes = 1.0; }else{ bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); } vol = 4.0*pi/3.0*rcore*rcore*rcore; f = vol*bes*contr; qr=x*(rcore+thick); contr = rhoshel-rhosolv; bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); if(qr == 0.0){ bes = 1.0; }else{ bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); } vol = 4.0*pi/3.0*pow((rcore+thick),3); f += vol*bes*contr; // normalize to particle volume and rescale from [-1] to [cm-1] // normalize to particle volume and rescale from [A-1] to [cm-1] f2 = f*f/vol*1.0e8; qr=x*rcore; contr = rhocore-rhoshel; bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); if(qr == 0){ bes = 1.0; }else{ bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); } vol = 4.0*pi/3.0*rcore*rcore*rcore; f = vol*bes*contr; qr=x*(rcore+thick); contr = rhoshel-rhosolv; bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); if(qr == 0.0){ bes = 1.0; }else{ bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); } vol = 4.0*pi/3.0*pow((rcore+thick),3); f += vol*bes*contr; // normalize to the particle volume and rescale from [-1] to [cm-1] // normalize to the particle volume and rescale from [A-1] to [cm-1] //note that for the vesicle model, the volume is ONLY the shell volume vol = 4.0*pi/3.0*(pow((rcore+thick),3)-pow(rcore,3)); //small QR limit - use Guinier approx zp8 = zz+8; zp8 = zz+8.0; if(x*ravg < 0.1) { i_zero = scale*delrho*delrho*1.e8*4.*pi/3.*pow(ravg,3); // aa = (zz+1)/x/ravg; aa = (zz+1.0)/x/ravg; at1 = atan(1.0/aa); pq = log10(g1) - 6.0*log10(zp1) + 6.0*log10(ravg); pq = pow(10,pq)*8*pi*pi*delrho*delrho; pq = pow(10,pq)*8.0*pi*pi*delrho*delrho; // // so just use the limiting Guiner value if(qr<0.1) { h1 = scale*cont*cont*1.e8*4.*pi/3*pow(rad,3); h1 = scale*cont*cont*1.e8*4.*pi/3.0*pow(rad,3); h1 *=   (1. + 15.*pow(pd,2) + 27.*pow(pd,4) +27./7.*pow(pd,6) );                                //6th moment h1 /= (1.+3.*pd*pd);    //3rd moment // normal calculation h1 = -0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0; h1 -= 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw); h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw); h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw); h1 += qw*qr*sin(2*qr)*cos(2*qw); h1 -= 5.0/2.0*cos(2.0*qr)*sin(qw)*cos(qw); h1 += 0.5*qr*qr*cos(2.0*qr)*sin(2.0*qw); h1 += 0.5*qw*qw*cos(2.0*qr)*sin(2.0*qw); h1 += qw*qr*sin(2.0*qr)*cos(2.0*qw); h1 += 3.0*qw*(cos(qr)*cos(qw))*(cos(qr)*cos(qw)); h1+= 3.0*qw*(sin(qr)*sin(qw))*(sin(qr)*sin(qw)); va = -4.0*sig + rad; if (va<0) { va=0;           //to avoid numerical error when  va<0 (-ve q-value) if (va<0.0) { va=0.0;         //to avoid numerical error when  va<0 (-ve q-value) } vb = 4.0*sig +rad; va = -3.5*sig + mu; va = exp(va); if (va<0) { va=0;           //to avoid numerical error when  va<0 (-ve q-value) if (va<0.0) { va=0.0;         //to avoid numerical error when  va<0 (-ve q-value) } vb = 3.5*sig*(1.0+sig) +mu; pi = 4.0*atan(1.0); retval = (1/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); return(retval); } double psf11,psf12,psf22; double phi1,phi2,phr,a3; double v1,v2,n1,n2,qr1,qr2,b1,b2; double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2; int err; qr1 = r1*x; qr2 = r2*x; b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*(sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1; b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*(sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2; if (qr1 == 0){ sc1 = 1.0/3.0; }else{ sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1; } if (qr2 == 0){ sc2 = 1.0/3.0; }else{ sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2; } b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*sc1; b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*sc2; inten = n1*b1*b1*psf11; inten += sqrt(n1*n2)*2.0*b1*b2*psf12; double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4; double a1,a2i,a2,b1,b2,b12,gm1,gm12; double err=0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3; double err=0.0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3; double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13; double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1; ii=0; fval=0; fval=0.0; do { ri = rcore + (double)ii*ts + (double)ii*tw; voli = 4*pi/3*ri*ri*ri; voli = 4.0*pi/3.0*ri*ri*ri; sldi = rhocore-rhoshel; fval += voli*sldi*F_func(ri*x); ri += ts; voli = 4*pi/3*ri*ri*ri; voli = 4.0*pi/3.0*ri*ri*ri; sldi = rhoshel-rhocore; fval += voli*sldi*F_func(ri*x); fval *= fval;           //square it fval /= voli;           //normalize by the overall volume fval *= scale*1e8; fval *= scale*1.0e8; fval += bkg; ii=minPairs; fval=0; d1 = 0; d2 = 0; r1 = 0; r2 = 0; distr = 0; first = 1; fval=0.0; d1 = 0.0; d2 = 0.0; r1 = 0.0; r2 = 0.0; distr = 0.0; first = 1.0; do { //make the current values old } while(ii<=maxPairs); fval /= 4*pi/3*avg*avg*avg;             //normalize by the overall volume fval /= 4.0*pi/3.0*avg*avg*avg;         //normalize by the overall volume fval /= distr; fval *= scale; pi = 4.0*atan(1.0); ii=0; fval=0; fval=0.0; do { ri = rcore + (double)ii*ts + (double)ii*tw; voli = 4*pi/3*ri*ri*ri; voli = 4.0*pi/3.0*ri*ri*ri; sldi = rhocore-rhoshel; fval += voli*sldi*F_func(ri*x); ri += ts; voli = 4*pi/3*ri*ri*ri; voli = 4.0*pi/3.0*ri*ri*ri; sldi = rhoshel-rhocore; fval += voli*sldi*F_func(ri*x); fval *= fval; fval /= voli; fval *= 1e8; fval *= 1.0e8; return(fval);       // this result still needs to be multiplied by scale and have background added double dr; dr = zz*log(x) - gammln(zz+1)+(zz+1)*log((zz+1)/avg)-(x/avg*(zz+1)); dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0)); return (exp(dr)); } double F_func(double qr) { return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); double sc; if (qr == 0.0){ sc = 1.0; }else{ sc=(3.0*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); } return sc; } contr = rhocore-rhoshel1; bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); vol = 4.0*pi/3*rcore*rcore*rcore; vol = 4.0*pi/3.0*rcore*rcore*rcore; f = vol*bes*contr; //now the shell (1) contr = rhocore-rhoshel1; bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); vol = 4.0*pi/3*rcore*rcore*rcore; vol = 4.0*pi/3.0*rcore*rcore*rcore; f = vol*bes*contr; //now the shell (1) contr = rhocore-rhoshel1; bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); vol = 4.0*pi/3*rcore*rcore*rcore; vol = 4.0*pi/3.0*rcore*rcore*rcore; f = vol*bes*contr; //now the shell (1) range = 8.0;                    //std deviations for the integration va = rcore*(1.0-range*pd); if (va<0) { va=0;           //otherwise numerical error when pd >= 0.3, making a<0 if (va<0.0) { va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0 } if (pd>0.3) { range = 8.0;                    //std deviations for the integration va = rcore*(1.0-range*pd); if (va<0) { va=0;           //otherwise numerical error when pd >= 0.3, making a<0 if (va<0.0) { va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0 } if (pd>0.3) {
Note: See TracChangeset for help on using the changeset viewer.