Ignore:
Timestamp:
May 13, 2009 9:18:52 PM (14 years ago)
Author:
ajj
Message:

Bicelle and Core-Shell Parallelepiped models for Andrea Hamill

File:
1 edited

Legend:

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

    r500 r501  
    679679        double summ,yyy,answer,Vpoly;                   //running tally of integration 
    680680        double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr; 
    681          
     681                 
    682682        Pi = 4.0*atan(1.0); 
    683683         
     
    28412841} 
    28422842 
    2843 double 
    2844 PolyCoreBicelle(double w[], double q){ 
    2845  
     2843double PolyCoreBicelle(double dp[], double q) 
     2844{ 
    28462845        int i; 
    28472846        int nord = 20; 
     
    28492848        double rhoc, rhoh, rhor, rhosolv; 
    28502849        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); 
    28652870         
    28662871        lolim = exp(log(radius)-(4.*sigma)); 
     
    28692874        } 
    28702875        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; 
    28712882         
    28722883        for(i=0;i<nord;i++) { 
     
    28912902        answer += bkg; 
    28922903         
    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; 
    28942908         
    28952909} 
     
    29502964         
    29512965        retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); 
    2952         return retval; 
    2953          
    2954 } 
     2966        return(retval); 
     2967         
     2968} 
     2969 
     2970 
     2971double 
     2972CSParallelpiped(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]*(vb-va)+vb+va)/2.0; 
     3004                yyy = Gauss76Wt[i]*CSPP_Outer(dp,q,zi); 
     3005                summ += yyy; 
     3006        } 
     3007        inten = (vb-va)/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 
     3019double 
     3020CSPP_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]*(vb-va)+vb+va)/2.0; 
     3039                yyy = Gauss76Wt[i]*CSPP_Inner(dp,mudum,zi); 
     3040                summ += yyy; 
     3041        } 
     3042        retval = (vb-va)/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 
     3057double 
     3058CSPP_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=rhoP-rhosolv; 
     3079        drA=rhoA-rhosolv; 
     3080        drB=rhoB-rhosolv; 
     3081        drC=rhoC-rhosolv;  
     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*(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 
     3118        return(retVal);  
     3119         
     3120 
     3121 
     3122} 
Note: See TracChangeset for help on using the changeset viewer.