Changeset 594


Ignore:
Timestamp:
Nov 4, 2009 2:19:47 PM (13 years ago)
Author:
srkline
Message:

Correcting my merge mistakes in CSParallelepiped

Location:
sans/XOP_Dev/SANSAnalysis/lib
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c

    r507 r594  
    29262926 
    29272927double 
    2928 BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum){ 
    2929  
     2928BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum) 
     2929{ 
    29302930        double dr1,dr2,dr3; 
    29312931        double besarg1,besarg2; 
     
    29592959 
    29602960double 
    2961 CSParallelepiped(double dp[], double q){ 
    2962  
    2963         double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg,inten; 
    2964         double yyy,summ,va,vb,zi; 
    2965         int i; 
    2966         int nord=76; 
    2967          
    2968         scale = dp[0]; 
    2969         aa = dp[1]; 
    2970         bb = dp[2]; 
    2971         cc = dp[3]; 
    2972         ta  = dp[4]; 
    2973         tb  = dp[5]; 
    2974         tc  = dp[6];   // is 0 at the moment   
    2975         rhoA=dp[7];   //rim A SLD 
    2976         rhoB=dp[8];   //rim B SLD 
    2977         rhoC=dp[9];    //rim C SLD 
    2978         rhoP = dp[10];   //Parallelpiped core SLD 
    2979         rhosolv=dp[11];  // Solvent SLD 
    2980         bkg = dp[12]; 
    2981          
    2982         //inten = IntegrateFn20(CSPP_Outer,0,1,w,x) 
    2983         //      inten = IntegrateFn76(PP_Outer,0,1,w,x) 
    2984          
    2985         //Do integral of outer 
    2986         va = 0; 
    2987         vb = 1; 
    2988          
    2989         summ = 0.0; 
    2990          
    2991         for(i=0;i<nord;i++){ 
    2992                 zi = (Gauss76Z[i]*(vb-va)+vb+va)/2.0; 
    2993                 yyy = Gauss76Wt[i]*CSPP_Outer(dp,q,zi); 
    2994                 summ += yyy; 
    2995         } 
    2996         inten = (vb-va)/2.0*summ; 
    2997          
    2998         inten /= (aa*bb*cc+2*ta*bb*cc+2*aa*tb*cc+2*aa*bb*tc);           //divide by outer volume (=Volume of core+edges) 
    2999         inten *= 1e8;           //convert to cm^-1 
    3000         inten *= scale; 
    3001         inten += bkg; 
    3002          
    3003         return (inten);  
    3004          
    3005 } 
    3006  
    3007  
    3008 double 
    3009 CSPP_Outer(double dp[], double q, double dum){ 
    3010         double retVal,mu,aa,bb,cc,mudum,arg ;  
    3011         double summ,yyy,va,vb,zi; 
    3012         double retval; 
    3013         int i; 
    3014         int nord=76; 
    3015         aa = dp[1]; 
    3016         bb = dp[2]; 
    3017         cc = dp[3]; 
    3018          
    3019         mu= bb*q; 
    3020         mudum = mu*sqrt(1-(dum*dum)); 
    3021          
    3022         va = 0; 
    3023         vb = 1; 
    3024          
    3025         summ = 0.0; 
    3026         retval = 0.0; 
    3027  
    3028         //Do Inner integral 
    3029         for(i=0;i<nord;i++){ 
    3030                 zi = (Gauss76Z[i]*(vb-va)+vb+va)/2.0; 
    3031                 yyy = Gauss76Wt[i]*CSPP_Inner(dp,mudum,zi); 
    3032                 summ += yyy; 
    3033         } 
    3034         retval = (vb-va)/2.0*summ; 
    3035          
    3036          
    3037         cc = cc/bb; 
    3038         arg = mu*cc*dum/2; 
    3039         if(arg==0){ 
    3040                 retval *= 1; 
    3041         } else { 
    3042                 retval *= (sin(arg)/arg)*(sin(arg)/arg); 
    3043         } 
    3044                          
    3045         return(retVal); 
    3046  
    3047 } 
    3048  
    3049 double 
    3050 CSPP_Inner(double dp[], double mu, double uu){ 
    3051          
     2961CSPPKernel(double dp[], double mu, double uu) 
     2962{        
    30522963        double aa,bb,cc, ta,tb,tc;  
    30532964        double Vin,Vot,V1,V2; 
    30542965        double rhoA,rhoB,rhoC, rhoP, rhosolv; 
    3055         double dr0, drA,drB, drC,retVal; 
     2966        double dr0, drA,drB, drC; 
    30562967        double arg1,arg2,arg3,arg4,t1,t2, t3, t4; 
    3057         double Pi,retval; 
     2968        double Pi,retVal; 
    30582969 
    30592970        aa = dp[1]; 
     
    30732984        drC=rhoC-rhosolv;  
    30742985        Vin=(aa*bb*cc); 
    3075         Vot=(aa*bb*cc+2*ta*bb*cc+2*aa*tb*cc+2*aa*bb*tc); 
    3076         V1=(2*ta*bb*cc);   //  incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    3077         V2=(2*aa*tb*cc);  // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
     2986        Vot=(aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc); 
     2987        V1=(2.0*ta*bb*cc);   //  incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
     2988        V2=(2.0*aa*tb*cc);  // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    30782989        aa = aa/bb; 
    3079         ta=(aa+2*ta)/bb; 
    3080         tb=(aa+2*tb)/bb; 
     2990        ta=(aa+2.0*ta)/bb; 
     2991        tb=(aa+2.0*tb)/bb; 
    30812992         
    30822993        Pi = 4.0*atan(1.0); 
    30832994         
    3084         arg1 = (mu*aa/2)*sin(Pi*uu/2); 
    3085         arg2 = (mu/2)*cos(Pi*uu/2); 
    3086         arg3=  (mu*ta/2)*sin(Pi*uu/2); 
    3087         arg4=  (mu*tb/2)*cos(Pi*uu/2); 
     2995        arg1 = (mu*aa/2.0)*sin(Pi*uu/2.0); 
     2996        arg2 = (mu/2.0)*cos(Pi*uu/2.0); 
     2997        arg3=  (mu*ta/2.0)*sin(Pi*uu/2.0); 
     2998        arg4=  (mu*tb/2.0)*cos(Pi*uu/2.0); 
    30882999                          
    30893000        if(arg1==0){ 
    3090                 t1 = 1; 
     3001                t1 = 1.0; 
    30913002        } else { 
    30923003                t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)    
    30933004        } 
    30943005        if(arg2==0){ 
    3095                 t2 = 1; 
     3006                t2 = 1.0; 
    30963007        } else { 
    30973008                t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)     
    30983009        }        
    30993010        if(arg3==0){ 
    3100                 t3 = 1; 
     3011                t3 = 1.0; 
    31013012        } else { 
    31023013                t3 = sin(arg3)/arg3; 
    31033014        } 
    31043015        if(arg4==0){ 
    3105                 t4 = 1; 
     3016                t4 = 1.0; 
    31063017        } else { 
    31073018                t4 = sin(arg4)/arg4; 
    31083019        } 
    3109         retval =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors 
     3020        retVal =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors 
    31103021        return(retVal);  
    3111          
    3112  
    3113  
    3114 } 
     3022 
     3023} 
     3024 
     3025/*      CSParallelepiped  :  calculates the form factor of a Parallelepiped with a core-shell structure 
     3026 -- different SLDs can be used for the face and rim 
     3027 
     3028Uses 76 pt Gaussian quadrature for both integrals 
     3029*/ 
     3030double 
     3031CSParallelepiped(double dp[], double q) 
     3032{ 
     3033        int i,j; 
     3034        double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg;         //local variables of coefficient wave 
     3035        int nordi=76;                   //order of integration 
     3036        int nordj=76; 
     3037        double va,vb;           //upper and lower integration limits 
     3038        double summ,yyy,answer;                 //running tally of integration 
     3039        double summj,vaj,vbj;                   //for the inner integration 
     3040        double mu,mudum,arg,sigma,uu,vol; 
     3041         
     3042         
     3043        //      Pi = 4.0*atan(1.0); 
     3044        va = 0.0; 
     3045        vb = 1.0;               //orintational average, outer integral 
     3046        vaj = 0.0; 
     3047        vbj = 1.0;              //endpoints of inner integral 
     3048         
     3049        summ = 0.0;                     //initialize intergral 
     3050         
     3051        scale = dp[0]; 
     3052        aa = dp[1]; 
     3053        bb = dp[2]; 
     3054        cc = dp[3]; 
     3055        ta  = dp[4]; 
     3056        tb  = dp[5]; 
     3057        tc  = dp[6];   // is 0 at the moment   
     3058        rhoA=dp[7];   //rim A SLD 
     3059        rhoB=dp[8];   //rim B SLD 
     3060        rhoC=dp[9];    //rim C SLD 
     3061        rhoP = dp[10];   //Parallelpiped core SLD 
     3062        rhosolv=dp[11];  // Solvent SLD 
     3063        bkg = dp[12]; 
     3064         
     3065        mu = q*bb; 
     3066        vol = aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc;          //calculate volume before rescaling 
     3067         
     3068        // do the rescaling here, not in the kernel 
     3069        // normalize all WRT bb 
     3070        aa = aa/bb; 
     3071        cc = cc/bb; 
     3072         
     3073        for(i=0;i<nordi;i++) { 
     3074                //setup inner integral over the ellipsoidal cross-section 
     3075                summj=0; 
     3076                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy 
     3077                 
     3078                for(j=0;j<nordj;j++) { 
     3079                        //76 gauss points for the inner integral 
     3080                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy 
     3081                        mudum = mu*sqrt(1.0-sigma*sigma); 
     3082                        yyy = Gauss76Wt[j] * CSPPKernel(dp,mudum,uu); 
     3083                        summj += yyy; 
     3084                } 
     3085                //now calculate the value of the inner integral 
     3086                answer = (vbj-vaj)/2.0*summj; 
     3087                 
     3088                //finish the outer integral cc already scaled 
     3089                arg = mu*cc*sigma/2.0; 
     3090                if ( arg == 0 ) { 
     3091                        answer *= 1.0; 
     3092                } else { 
     3093                        answer *= sin(arg)*sin(arg)/arg/arg; 
     3094                } 
     3095                 
     3096                //now sum up the outer integral 
     3097                yyy = Gauss76Wt[i] * answer; 
     3098                summ += yyy; 
     3099        }               //final scaling is done at the end of the function, after the NT_FP64 case 
     3100         
     3101        answer = (vb-va)/2.0*summ; 
     3102 
     3103        //normalize by volume 
     3104        answer /= vol; 
     3105        //convert to [cm-1] 
     3106        answer *= 1.0e8; 
     3107        //Scale 
     3108        answer *= scale; 
     3109        // add in the background 
     3110        answer += bkg; 
     3111         
     3112        return answer; 
     3113} 
     3114 
  • sans/XOP_Dev/SANSAnalysis/lib/libCylinder.h

    r501 r594  
    6464double BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum); 
    6565double BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length); 
    66 double CSPP_Outer(double *, double q, double dum); 
    67 double CSPP_Inner(double *, double mu, double uu); 
     66double CSPPKernel(double dp[], double mu, double uu); 
    6867 
    6968/////////functions for WRC implementation of flexible cylinders 
Note: See TracChangeset for help on using the changeset viewer.