Ignore:
Timestamp:
May 13, 2009 12:48:45 PM (14 years ago)
Author:
ajj
Message:
  • Adding Bicelle model for Andrea Hamill
  • Un "fixing" previous bug fix on SANSModelPicker
File:
1 edited

Legend:

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

    r453 r500  
    28402840        return(val); 
    28412841} 
     2842 
     2843double 
     2844PolyCoreBicelle(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]*(uplim-lolim) + 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 = (uplim-lolim)/2.0*summ; 
     2882        Rsqr = (uplim-lolim)/2.0*Rsqrsumm; 
     2883        //normalize by average cylinder volume 
     2884        Vpoly = Pi*Rsqr*(length+2*facthick); 
     2885        answer /= Vpoly; 
     2886        //convert to [cm-1] 
     2887        answer *= 1.0e8; 
     2888        //Scale 
     2889        answer *= scale; 
     2890        // add in the background 
     2891        answer += bkg; 
     2892         
     2893        return(answer); 
     2894         
     2895} 
     2896 
     2897double 
     2898BicelleIntegration(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]*(uplim-lolim) + 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 = (uplim-lolim)/2.0*summ; 
     2921        return(answer);  
     2922} 
     2923 
     2924double 
     2925BicelleKernel(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 = rhoc-rhoh; 
     2937        dr2 = rhor-rhosolv; 
     2938        dr3=  rhoh-rhor; 
     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.