Changeset 501 for sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c
 Timestamp:
 May 13, 2009 9:18:52 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c
r500 r501 679 679 double summ,yyy,answer,Vpoly; //running tally of integration 680 680 double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr; 681 681 682 682 Pi = 4.0*atan(1.0); 683 683 … … 2841 2841 } 2842 2842 2843 double 2844 PolyCoreBicelle(double w[], double q){ 2845 2843 double PolyCoreBicelle(double dp[], double q) 2844 { 2846 2845 int i; 2847 2846 int nord = 20; … … 2849 2848 double rhoc, rhoh, rhor, rhosolv; 2850 2849 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); 2850 double Pi,lolim,uplim,summ,yyy,rad,AR,Rsqr,Rsqrsumm,Rsqryyy; 2851 2852 scale = dp[0]; 2853 radius = dp[1]; 2854 sigma = dp[2]; //sigma is the standard mean deviation 2855 length = dp[3]; 2856 radthick = dp[4]; 2857 facthick= dp[5]; 2858 rhoc = dp[6]; 2859 rhoh = dp[7]; 2860 rhor=dp[8]; 2861 rhosolv = dp[9]; 2862 bkg = dp[10]; 2863 2864 char temp[256]; 2865 2866 /*sprintf(temp, "%f; %f; %f; %f; %f; %f; %f; %f; %f; %f; %f;\r", scale, radius, sigma,length,radthick,facthick,rhoc,rhoh,rhor,rhosolv,bkg); 2867 XOPNotice(temp);*/ 2868 2869 Pi = 4.0*atan(1.0); 2865 2870 2866 2871 lolim = exp(log(radius)(4.*sigma)); … … 2869 2874 } 2870 2875 uplim = exp(log(radius)+(4.*sigma)); 2876 2877 /*sprintf(temp, "%f; %f;\r", lolim,uplim); 2878 XOPNotice(temp);*/ 2879 2880 summ = 0.0; 2881 Rsqrsumm = 0.0; 2871 2882 2872 2883 for(i=0;i<nord;i++) { … … 2891 2902 answer += bkg; 2892 2903 2893 return(answer); 2904 /*sprintf(temp, "Q = %f : summ = %f : Rsqrsumm = %f : Ans = %f\r", q, summ, Rsqrsumm,answer); 2905 XOPNotice(temp);*/ 2906 2907 return answer; 2894 2908 2895 2909 } … … 2950 2964 2951 2965 retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); 2952 return retval; 2953 2954 } 2966 return(retval); 2967 2968 } 2969 2970 2971 double 2972 CSParallelpiped(double dp[], double q){ 2973 2974 double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg,inten; 2975 double yyy,summ,va,vb,zi; 2976 int i; 2977 int nord=76; 2978 2979 scale = dp[0]; 2980 aa = dp[1]; 2981 bb = dp[2]; 2982 cc = dp[3]; 2983 ta = dp[4]; 2984 tb = dp[5]; 2985 tc = dp[6]; // is 0 at the moment 2986 rhoA=dp[7]; //rim A SLD 2987 rhoB=dp[8]; //rim B SLD 2988 rhoC=dp[9]; //rim C SLD 2989 rhoP = dp[10]; //Parallelpiped core SLD 2990 rhosolv=dp[11]; // Solvent SLD 2991 bkg = dp[12]; 2992 2993 //inten = IntegrateFn20(CSPP_Outer,0,1,w,x) 2994 // inten = IntegrateFn76(PP_Outer,0,1,w,x) 2995 2996 //Do integral of outer 2997 va = 0; 2998 vb = 1; 2999 3000 summ = 0.0; 3001 3002 for(i=0;i<nord;i++){ 3003 zi = (Gauss76Z[i]*(vbva)+vb+va)/2.0; 3004 yyy = Gauss76Wt[i]*CSPP_Outer(dp,q,zi); 3005 summ += yyy; 3006 } 3007 inten = (vbva)/2.0*summ; 3008 3009 inten /= (aa*bb*cc+2*ta*bb*cc+2*aa*tb*cc+2*aa*bb*tc); //divide by outer volume (=Volume of core+edges) 3010 inten *= 1e8; //convert to cm^1 3011 inten *= scale; 3012 inten += bkg; 3013 3014 return (inten); 3015 3016 } 3017 3018 3019 double 3020 CSPP_Outer(double dp[], double q, double dum){ 3021 double retVal,mu,aa,bb,cc,mudum,arg ; 3022 double summ,yyy,va,vb,zi; 3023 double retval; 3024 int i; 3025 int nord=76; 3026 aa = dp[1]; 3027 bb = dp[2]; 3028 cc = dp[3]; 3029 3030 mu= bb*q; 3031 mudum = mu*sqrt(1(dum*dum)); 3032 3033 va = 0; 3034 vb = 1; 3035 3036 //Do Inner integral 3037 for(i=0;i<nord;i++){ 3038 zi = (Gauss76Z[i]*(vbva)+vb+va)/2.0; 3039 yyy = Gauss76Wt[i]*CSPP_Inner(dp,mudum,zi); 3040 summ += yyy; 3041 } 3042 retval = (vbva)/2.0*summ; 3043 3044 3045 cc = cc/bb; 3046 arg = mu*cc*dum/2; 3047 if(arg==0){ 3048 retval *= 1; 3049 } else { 3050 retval *= (sin(arg)/arg)*(sin(arg)/arg); 3051 } 3052 3053 return(retVal); 3054 3055 } 3056 3057 double 3058 CSPP_Inner(double dp[], double mu, double uu){ 3059 3060 double aa,bb,cc, ta,tb,tc; 3061 double Vin,Vot,V1,V2; 3062 double rhoA,rhoB,rhoC, rhoP, rhosolv; 3063 double dr0, drA,drB, drC,retVal; 3064 double arg1,arg2,arg3,arg4,t1,t2, t3, t4; 3065 double Pi,retval; 3066 3067 aa = dp[1]; 3068 bb = dp[2]; 3069 cc = dp[3]; 3070 ta = dp[4]; 3071 tb = dp[5]; 3072 tc = dp[6]; 3073 rhoA=dp[7]; 3074 rhoB=dp[8]; 3075 rhoC=dp[9]; 3076 rhoP=dp[10]; 3077 rhosolv=dp[11]; 3078 dr0=rhoPrhosolv; 3079 drA=rhoArhosolv; 3080 drB=rhoBrhosolv; 3081 drC=rhoCrhosolv; 3082 Vin=(aa*bb*cc); 3083 Vot=(aa*bb*cc+2*ta*bb*cc+2*aa*tb*cc+2*aa*bb*tc); 3084 V1=(2*ta*bb*cc); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 3085 V2=(2*aa*tb*cc); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 3086 aa = aa/bb; 3087 ta=(aa+2*ta)/bb; 3088 tb=(aa+2*tb)/bb; 3089 3090 Pi = 4.0*atan(1.0); 3091 3092 arg1 = (mu*aa/2)*sin(Pi*uu/2); 3093 arg2 = (mu/2)*cos(Pi*uu/2); 3094 arg3= (mu*ta/2)*sin(Pi*uu/2); 3095 arg4= (mu*tb/2)*cos(Pi*uu/2); 3096 3097 if(arg1==0){ 3098 t1 = 1; 3099 } else { 3100 t1 = (sin(arg1)/arg1); //defn for CSPP model sin(arg1)/arg1 test: (sin(arg1)/arg1)*(sin(arg1)/arg1) 3101 } 3102 if(arg2==0){ 3103 t2 = 1; 3104 } else { 3105 t2 = (sin(arg2)/arg2); //defn for CSPP model sin(arg2)/arg2 test: (sin(arg2)/arg2)*(sin(arg2)/arg2) 3106 } 3107 if(arg3==0){ 3108 t3 = 1; 3109 } else { 3110 t3 = sin(arg3)/arg3; 3111 } 3112 if(arg4==0){ 3113 t4 = 1; 3114 } else { 3115 t4 = sin(arg4)/arg4; 3116 } 3117 retval =( dr0*t1*t2*Vin + drA*(t3t1)*t2*V1+ drB*t1*(t4t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3t1)*t2*V1+ drB*t1*(t4t2)*V2 ); // correct FF : square of sum of phase factors 3118 return(retVal); 3119 3120 3121 3122 }
Note: See TracChangeset
for help on using the changeset viewer.