Changeset 500 for sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c
 Timestamp:
 May 13, 2009 12:48:45 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c
r453 r500 2840 2840 return(val); 2841 2841 } 2842 2843 double 2844 PolyCoreBicelle(double w[], double q){ 2845 2846 int i; 2847 int nord = 20; 2848 double scale, length, sigma, bkg, radius, radthick, facthick; 2849 double rhoc, rhoh, rhor, rhosolv; 2850 double answer, Vpoly; 2851 double lolim,uplim,summ=0,yyy,rad,AR,Rsqr,Rsqrsumm=0,Rsqryyy; 2852 scale = w[0]; 2853 radius = w[1]; 2854 sigma = w[2]; //sigma is the standard mean deviation 2855 length = w[3]; 2856 radthick = w[4]; 2857 facthick= w[5]; 2858 rhoc = w[6]; 2859 rhoh = w[7]; 2860 rhor=w[8]; 2861 rhosolv = w[9]; 2862 bkg = w[10]; 2863 2864 double Pi = 4.0*atan(1.0); 2865 2866 lolim = exp(log(radius)(4.*sigma)); 2867 if (lolim<0) { 2868 lolim=0; //to avoid numerical error when va<0 (ve r value) 2869 } 2870 uplim = exp(log(radius)+(4.*sigma)); 2871 2872 for(i=0;i<nord;i++) { 2873 rad = ( Gauss20Z[i]*(uplimlolim) + uplim + lolim )/2.0; 2874 AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp((0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma))); 2875 yyy = AR* Gauss20Wt[i] * BicelleIntegration(q,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length); 2876 Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick); //SRK normalize to total dimensions 2877 summ += yyy; 2878 Rsqrsumm += Rsqryyy; 2879 } 2880 2881 answer = (uplimlolim)/2.0*summ; 2882 Rsqr = (uplimlolim)/2.0*Rsqrsumm; 2883 //normalize by average cylinder volume 2884 Vpoly = Pi*Rsqr*(length+2*facthick); 2885 answer /= Vpoly; 2886 //convert to [cm1] 2887 answer *= 1.0e8; 2888 //Scale 2889 answer *= scale; 2890 // add in the background 2891 answer += bkg; 2892 2893 return(answer); 2894 2895 } 2896 2897 double 2898 BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length){ 2899 2900 double answer,halfheight,Pi; 2901 double lolim,uplim,summ,yyy,zi; 2902 int nord,i; 2903 2904 // set up the integration end points 2905 Pi = 4.0*atan(1.0); 2906 nord = 76; 2907 lolim = 0; 2908 uplim = Pi/2; 2909 halfheight = length/2.0; 2910 2911 summ = 0.0; // initialize integral 2912 i=0; 2913 for(i=0;i<nord;i++) { 2914 zi = ( Gauss76Z[i]*(uplimlolim) + uplim + lolim )/2.0; 2915 yyy = Gauss76Wt[i] * BicelleKernel(qq, rad, radthick, facthick, rhoc, rhoh, rhor,rhosolv, halfheight, zi); 2916 summ += yyy; 2917 } 2918 2919 // calculate value of integral to return 2920 answer = (uplimlolim)/2.0*summ; 2921 return(answer); 2922 } 2923 2924 double 2925 BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum){ 2926 2927 double dr1,dr2,dr3; 2928 double besarg1,besarg2; 2929 double vol1,vol2,vol3; 2930 double sinarg1,sinarg2; 2931 double t1,t2,t3; 2932 double retval; 2933 2934 double Pi = 4.0*atan(1.0); 2935 2936 dr1 = rhocrhoh; 2937 dr2 = rhorrhosolv; 2938 dr3= rhohrhor; 2939 vol1 = Pi*rad*rad*(2*length); 2940 vol2 = Pi*(rad+radthick)*(rad+radthick)*(2*length+2*facthick); 2941 vol3= Pi*(rad)*(rad)*(2*length+2*facthick); 2942 besarg1 = qq*rad*sin(dum); 2943 besarg2 = qq*(rad+radthick)*sin(dum); 2944 sinarg1 = qq*length*cos(dum); 2945 sinarg2 = qq*(length+facthick)*cos(dum); 2946 2947 t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 2948 t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 2949 t3 = 2*vol3*dr3*sin(sinarg2)/sinarg2*NR_BessJ1(besarg1)/besarg1; 2950 2951 retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); 2952 return retval; 2953 2954 }
Note: See TracChangeset
for help on using the changeset viewer.