Changeset 71b751d in sasmodels for sasmodels/models


Ignore:
Timestamp:
Aug 14, 2018 12:09:34 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
86aa992
Parents:
2f8cbb9
Message:

update remaining form factors to use Fq interface

Location:
sasmodels/models
Files:
65 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/barbell.c

    r108e70e r71b751d  
    6262} 
    6363 
    64 static double 
    65 Iq(double q, double sld, double solvent_sld, 
     64static void 
     65Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 
    6666    double radius_bell, double radius, double length) 
    6767{ 
    … …  
    7272    const double zm = M_PI_4; 
    7373    const double zb = M_PI_4; 
    74     double total = 0.0; 
     74    double total_F1 = 0.0; 
     75    double total_F2 = 0.0; 
    7576    for (int i = 0; i < GAUSS_N; i++){ 
    7677        const double alpha= GAUSS_Z[i]*zm + zb; 
    … …  
    7879        SINCOS(alpha, sin_alpha, cos_alpha); 
    7980        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    80         total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
     81        total_F1 += GAUSS_W[i] * Aq * sin_alpha; 
     82        total_F2 += GAUSS_W[i] * Aq * Aq * sin_alpha; 
    8183    } 
    8284    // translate dx in [-1,1] to dx in [lower,upper] 
    83     const double form = total*zm; 
     85    const double form_avg = total_F1*zm; 
     86    const double form_squared_avg = total_F2*zm; 
    8487 
    8588    //Contrast 
    8689    const double s = (sld - solvent_sld); 
    87     return 1.0e-4 * s * s * form; 
     90    *F1 = 1.0e-2 * s * form_avg; 
     91    *F2 = 1.0e-4 * s * s * form_squared_avg; 
    8892} 
    89  
    9093 
    9194static double 

  • sasmodels/models/barbell.py

    r2d81cfe r71b751d  
    116116 
    117117source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 
     118have_Fq = True 
    118119 
    119120def random(): 

  • sasmodels/models/binary_hard_sphere.c

    r925ad6e r71b751d  
    11double form_volume(void); 
    22 
    3 double Iq(double q,  
     3double Iq(double q, 
    44    double lg_radius, double sm_radius, 
    55    double lg_vol_frac, double sm_vol_frac, 
    66    double lg_sld, double sm_sld, double solvent_sld 
    77    ); 
    8      
     8 
    99void calculate_psfs(double qval, 
    1010    double r2, double nf2, 
    … …  
    2222    double lg_vol_frac, double sm_vol_frac, 
    2323    double lg_sld, double sm_sld, double solvent_sld) 
    24      
    2524{ 
    2625    double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten;       //my local names 
    … …  
    2827    double phi1,phi2,phr,a3; 
    2928    double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2; 
    30      
     29 
    3130    r2 = lg_radius; 
    3231    r1 = sm_radius; 
    … …  
    3635    rho1 = sm_sld; 
    3736    rhos = solvent_sld; 
    38      
    39      
     37 
     38 
    4039    phi = phi1 + phi2; 
    4140    aa = r1/r2; 
    … …  
    4645    // calculate the PSF's here 
    4746    calculate_psfs(q,r2,nf2,aa,phi,&psf11,&psf22,&psf12); 
    48      
     47 
    4948    // /* do form factor calculations  */ 
    50      
     49 
    5150    v1 = M_4PI_3*r1*r1*r1; 
    5251    v2 = M_4PI_3*r2*r2*r2; 
    53      
     52 
    5453    n1 = phi1/v1; 
    5554    n2 = phi2/v2; 
    56      
     55 
    5756    qr1 = r1*q; 
    5857    qr2 = r2*q; 
    … …  
    6867    inten *= 1.0e8; 
    6968    ///*convert rho^2 in 10^-6A to A*/ 
    70     inten *= 1.0e-12;     
     69    inten *= 1.0e-12; 
    7170    return(inten); 
    7271} 
    … …  
    7776    double aa, double phi, 
    7877    double *s11, double *s22, double *s12) 
    79  
    8078{ 
    8179    //  variable qval,r2,nf2,aa,phi,&s11,&s22,&s12 
    82      
     80 
    8381    //   calculate constant terms 
    8482    double s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4; 
    … …  
    8785    double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13; 
    8886    double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1; 
    89      
     87 
    9088    s2 = 2.0*r2; 
    9189//    s1 = aa*s2;  why is this never used?  check original paper? 
    … …  
    108106    gm1=(v1*a1+a3*v2*a2)*.5; 
    109107    gm12=2.*gm1*(1.-aa)/aa; 
    110     //c   
     108    //c 
    111109    //c   calculate the direct correlation functions and print results 
    112110    //c 
    113111    //  do 20 j=1,npts 
    114      
     112 
    115113    yy=qval*s2; 
    116114    //c   calculate direct correlation functions 
    … …  
    123121    t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3; 
    124122    f11=24.*v1*(t1+t2+t3)/ay3; 
    125      
     123 
    126124    //c ------c22 
    127125    y2=yy*yy; 
    … …  
    131129    tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3; 
    132130    f22=24.*v2*(tt1+tt2+tt3)/y3; 
    133      
     131 
    134132    //c   -----c12 
    135133    yl=.5*yy*(1.-aa); 
    … …  
    151149    ttt4=a1*(t41+t42)/y1; 
    152150    f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3); 
    153      
     151 
    154152    c11=f11; 
    155153    c22=f22; 
    156154    c12=f12; 
    157155    *s11=1./(1.+c11-(c12)*c12/(1.+c22)); 
    158     *s22=1./(1.+c22-(c12)*c12/(1.+c11));  
    159     *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));    
    160      
     156    *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
     157    *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12)); 
     158 
    161159    return; 
    162160} 
    163  

  • sasmodels/models/capped_cylinder.c

    r108e70e r71b751d  
    8484} 
    8585 
    86 static double 
    87 Iq(double q, double sld, double solvent_sld, 
     86static void 
     87Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 
    8888    double radius, double radius_cap, double length) 
    8989{ 
    … …  
    9494    const double zm = M_PI_4; 
    9595    const double zb = M_PI_4; 
    96     double total = 0.0; 
     96    double total_F1 = 0.0; 
     97    double total_F2 = 0.0; 
    9798    for (int i=0; i<GAUSS_N ;i++) { 
    9899        const double theta = GAUSS_Z[i]*zm + zb; 
    … …  
    103104        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
    104105        // scale by sin_theta for spherical coord integration 
    105         total += GAUSS_W[i] * Aq * Aq * sin_theta; 
     106        total_F1 += GAUSS_W[i] * Aq * sin_theta; 
     107        total_F2 += GAUSS_W[i] * Aq * Aq * sin_theta; 
    106108    } 
    107109    // translate dx in [-1,1] to dx in [lower,upper] 
    108     const double form = total * zm; 
     110    const double form_avg = total_F1 * zm; 
     111    const double form_squared_avg = total_F2 * zm; 
    109112 
    110113    // Contrast 
    111114    const double s = (sld - solvent_sld); 
    112     return 1.0e-4 * s * s * form; 
     115    *F1 = 1.0e-2 * s * form_avg; 
     116    *F2 = 1.0e-4 * s * s * form_squared_avg; 
    113117} 
    114118 

  • sasmodels/models/capped_cylinder.py

    r2d81cfe r71b751d  
    136136 
    137137source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
     138have_Fq = True 
    138139 
    139140def random(): 

  • sasmodels/models/core_multi_shell.c

    rc3ccaec r71b751d  
    1919} 
    2020 
    21 static double 
    22 Iq(double q, double core_sld, double core_radius, 
     21static void 
     22Fq(double q, double *F1, double *F2, double core_sld, double core_radius, 
    2323   double solvent_sld, double fp_n, double sld[], double thickness[]) 
    2424{ 
    … …  
    3434  } 
    3535  f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sas_3j1x_x(q*r); 
    36   return f * f * 1.0e-4; 
     36  *F1 = 1e-2 * f; 
     37  *F2 = 1e-4 * f * f; 
    3738} 

  • sasmodels/models/core_multi_shell.py

    r2d81cfe r71b751d  
    9999 
    100100source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 
     101have_Fq = True 
    101102 
    102103def random(): 

  • sasmodels/models/core_shell_bicelle.c

    r108e70e r71b751d  
    3636} 
    3737 
    38 static double 
    39 Iq(double q, 
     38static void 
     39Fq(double q, 
     40    double *F1, 
     41    double *F2, 
    4042    double radius, 
    4143    double thick_radius, 
    … …  
    5153    const double halflength = 0.5*length; 
    5254 
    53     double total = 0.0; 
     55    double total_F1 = 0.0; 
     56    double total_F2 = 0.0; 
    5457    for(int i=0;i<GAUSS_N;i++) { 
    5558        double theta = (GAUSS_Z[i] + 1.0)*uplim; 
    5659        double sin_theta, cos_theta; // slots to hold sincos function output 
    5760        SINCOS(theta, sin_theta, cos_theta); 
    58         double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
     61        double form = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
    5962                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
    60         total += GAUSS_W[i]*fq*fq*sin_theta; 
     63        total_F1 += GAUSS_W[i]*form*sin_theta; 
     64        total_F2 += GAUSS_W[i]*form*form*sin_theta; 
    6165    } 
     66    // Correct for integration range 
     67    total_F1 *= uplim; 
     68    total_F2 *= uplim; 
    6269 
    63     // calculate value of integral to return 
    64     double answer = total*uplim; 
    65     return 1.0e-4*answer; 
     70    *F1 = 1.0e-2*total_F1; 
     71    *F2 = 1.0e-4*total_F2; 
    6672} 
    6773 

  • sasmodels/models/core_shell_bicelle.py

    r2d81cfe r71b751d  
    154154source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    155155          "core_shell_bicelle.c"] 
     156have_Fq = True 
    156157 
    157158def random(): 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r108e70e r71b751d  
    1010} 
    1111 
    12 static double 
    13 Iq(double q, 
     12static void 
     13Fq(double q, 
     14    double *F1, 
     15    double *F2, 
    1416    double r_minor, 
    1517    double x_core, 
    … …  
    3638 
    3739    //initialize integral 
    38     double outer_sum = 0.0; 
     40    double outer_total_F1 = 0.0; 
     41    double outer_total_F2 = 0.0; 
    3942    for(int i=0;i<GAUSS_N;i++) { 
    4043        //setup inner integral over the ellipsoidal cross-section 
    41         //const double va = 0.0; 
    42         //const double vb = 1.0; 
    4344        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    4445        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
    … …  
    4849        const double si1 = sas_sinx_x(halfheight*qc); 
    4950        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
    50         double inner_sum=0.0; 
     51        double inner_total_F1 = 0; 
     52        double inner_total_F2 = 0; 
    5153        for(int j=0;j<GAUSS_N;j++) { 
    5254            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    53             // inner integral limits 
    54             //const double vaj=0.0; 
    55             //const double vbj=M_PI; 
    56             //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    57             const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    58             const double rr = sqrt(r2A - r2B*cos(phi)); 
     55            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     56            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
     57            const double rr = sqrt(r2A - r2B*cos(beta)); 
    5958            const double be1 = sas_2J1x_x(rr*qab); 
    6059            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
    61             const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
     60            const double f = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
    6261 
    63             inner_sum += GAUSS_W[j] * fq * fq; 
     62            inner_total_F1 += GAUSS_W[j] * f; 
     63            inner_total_F2 += GAUSS_W[j] * f * f; 
    6464        } 
    6565        //now calculate outer integral 
    66         outer_sum += GAUSS_W[i] * inner_sum; 
     66        outer_total_F1 += GAUSS_W[i] * inner_total_F1; 
     67        outer_total_F2 += GAUSS_W[i] * inner_total_F2; 
    6768    } 
     69    // now complete change of integration variables (1-0)/(1-(-1))= 0.5 
     70    outer_total_F1 *= 0.25; 
     71    outer_total_F2 *= 0.25; 
    6872 
    69     return outer_sum*2.5e-05; 
     73    //convert from [1e-12 A-1] to [cm-1] 
     74    *F1 = 1e-2*outer_total_F1; 
     75    *F2 = 1e-4*outer_total_F2; 
    7076} 
    7177 

  • sasmodels/models/core_shell_bicelle_elliptical.py

    rfc3ae1b r71b751d  
    146146source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    147147          "core_shell_bicelle_elliptical.c"] 
     148have_Fq = True 
    148149 
    149150def random(): 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r108e70e r71b751d  
    1111} 
    1212 
    13 static double 
    14 Iq(double q, 
     13static void 
     14Fq(double q, 
     15        double *F1, 
     16        double *F2, 
    1517        double r_minor, 
    1618        double x_core, 
    … …  
    2426        double sigma) 
    2527{ 
    26     double si1,si2,be1,be2; 
    2728     // core_shell_bicelle_elliptical_belt, RKH 5th Oct 2017, core_shell_bicelle_elliptical 
    2829     // tested briefly against limiting cases of cylinder, hollow cylinder & elliptical cylinder models 
    … …  
    3839    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    3940    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
     41    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
     42    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*halfheight; 
     43    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    4044    // dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face, 
    41     const double dr1 = (-rhor - rhoh + rhoc + rhosolv) *M_PI*r_minor*r_major* 
    42           2.0*halfheight; 
    43     const double dr2 = (rhor-rhosolv) *M_PI*(r_minor+thick_rim)*( 
    44          r_major+thick_rim)* 2.0*halfheight; 
    45     const double dr3 = (rhoh-rhosolv) *M_PI*r_minor*r_major* 
    46          2.0*(halfheight+thick_face); 
     45    const double dr1 = vol1*(-rhor - rhoh + rhoc + rhosolv); 
     46    const double dr2 = vol2*(rhor-rhosolv); 
     47    const double dr3 = vol3*(rhoh-rhosolv); 
     48 
    4749    //initialize integral 
    48     double outer_sum = 0.0; 
     50    double outer_total_F1 = 0.0; 
     51    double outer_total_F2 = 0.0; 
    4952    for(int i=0;i<GAUSS_N;i++) { 
    5053        //setup inner integral over the ellipsoidal cross-section 
    5154        // since we generate these lots of times, why not store them somewhere? 
    52         //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    53         const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0; 
    54         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    55         double inner_sum=0; 
    56         double sinarg1 = q*halfheight*cos_alpha; 
    57         double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
    58         si1 = sas_sinx_x(sinarg1); 
    59         si2 = sas_sinx_x(sinarg2); 
     55        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     56        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
     57        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     58        const double qab = q*sin_theta; 
     59        const double qc = q*cos_theta; 
     60        const double si1 = sas_sinx_x(halfheight*qc); 
     61        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
     62        double inner_total_F1 = 0; 
     63        double inner_total_F2 = 0; 
    6064        for(int j=0;j<GAUSS_N;j++) { 
    6165            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    … …  
    6367            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    6468            const double rr = sqrt(r2A - r2B*cos(beta)); 
    65             double besarg1 = q*rr*sin_alpha; 
    66             double besarg2 = q*(rr+thick_rim)*sin_alpha; 
    67             be1 = sas_2J1x_x(besarg1); 
    68             be2 = sas_2J1x_x(besarg2); 
    69             inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 
    70                                               dr2*si1*be2 + 
    71                                               dr3*si2*be1); 
     69            const double be1 = sas_2J1x_x(rr*qab); 
     70            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
     71            const double f = dr1*si1*be1 + dr2*si1*be2 + dr3*si2*be1; 
     72 
     73            inner_total_F1 += GAUSS_W[j] * f; 
     74            inner_total_F2 += GAUSS_W[j] * f * f; 
    7275        } 
    7376        //now calculate outer integral 
    74         outer_sum += GAUSS_W[i] * inner_sum; 
     77        outer_total_F1 += GAUSS_W[i] * inner_total_F1; 
     78        outer_total_F2 += GAUSS_W[i] * inner_total_F2; 
    7579    } 
     80    // now complete change of integration variables (1-0)/(1-(-1))= 0.5 
     81    outer_total_F1 *= 0.25; 
     82    outer_total_F2 *= 0.25; 
    7683 
    77     return outer_sum*2.5e-05*exp(-0.5*square(q*sigma)); 
     84    //convert from [1e-12 A-1] to [cm-1] 
     85    *F1 = 1e-2*outer_total_F1*exp(-0.25*square(q*sigma)); 
     86    *F2 = 1e-4*outer_total_F2*exp(-0.5*square(q*sigma)); 
    7887} 
    7988 
    … …  
    111120    const double si1 = sas_sinx_x( halfheight*qc ); 
    112121    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 
    113     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1); 
    114     return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 
     122    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1; 
     123    const double atten = exp(-0.5*(qa*qa + qb*qb + qc*qc)*(sigma*sigma)); 
     124    return 1.0e-4 * fq*fq * atten; 
    115125} 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    rfc3ae1b r71b751d  
    159159source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    160160          "core_shell_bicelle_elliptical_belt_rough.c"] 
     161have_Fq = True 
    161162 
    162163demo = dict(scale=1, background=0, 

  • sasmodels/models/core_shell_cylinder.c

    rd86f0fc r71b751d  
    1313} 
    1414 
    15 static double 
    16 Iq(double q, 
     15static void 
     16Fq(double q, 
     17    double *F1, 
     18    double *F2, 
    1719    double core_sld, 
    1820    double shell_sld, 
    … …  
    2931    const double shell_h = (0.5*length + thickness); 
    3032    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    31     double total = 0.0; 
     33    double total_F1 = 0.0; 
     34    double total_F2 = 0.0; 
    3235    for (int i=0; i<GAUSS_N ;i++) { 
    3336        // translate a point in [-1,1] to a point in [0, pi/2] 
    … …  
    4043        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
    4144            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    42         total += GAUSS_W[i] * fq * fq * sin_theta; 
     45        total_F1 += GAUSS_W[i] * fq * sin_theta; 
     46        total_F2 += GAUSS_W[i] * fq * fq * sin_theta; 
    4347    } 
    4448    // translate dx in [-1,1] to dx in [lower,upper] 
    4549    //const double form = (upper-lower)/2.0*total; 
    46     return 1.0e-4 * total * M_PI_4; 
     50    *F1 = 1.0e-2 * total_F1 * M_PI_4; 
     51    *F2 = 1.0e-4 * total_F2 * M_PI_4; 
    4752} 
    48  
    4953 
    5054static double 

  • sasmodels/models/core_shell_cylinder.py

    r2d81cfe r71b751d  
    125125 
    126126source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 
     127have_Fq = True 
    127128 
    128129def ER(radius, thickness, length): 

  • sasmodels/models/core_shell_ellipsoid.c

    r108e70e r71b751d  
    3838} 
    3939 
    40 static double 
    41 Iq(double q, 
     40static void 
     41Fq(double q, 
     42    double *F1, 
     43    double *F2, 
    4244    double radius_equat_core, 
    4345    double x_core, 
    … …  
    5860    const double m = 0.5; 
    5961    const double b = 0.5; 
    60     double total = 0.0;     //initialize intergral 
     62    double total_F1 = 0.0;     //initialize intergral 
     63    double total_F2 = 0.0;     //initialize intergral 
    6164    for(int i=0;i<GAUSS_N;i++) { 
    6265        const double cos_theta = GAUSS_Z[i]*m + b; 
    … …  
    6669            equat_shell, polar_shell, 
    6770            sld_core_shell, sld_shell_solvent); 
    68         total += GAUSS_W[i] * fq * fq; 
     71        total_F1 += GAUSS_W[i] * fq; 
     72        total_F2 += GAUSS_W[i] * fq * fq; 
    6973    } 
    70     total *= m; 
     74    total_F1 *= m; 
     75    total_F2 *= m; 
    7176 
    7277    // convert to [cm-1] 
    73     return 1.0e-4 * total; 
     78    *F1 = 1.0e-2 * total_F1; 
     79    *F2 = 1.0e-4 * total_F2; 
    7480} 
     81 
    7582 
    7683static double 

  • sasmodels/models/core_shell_ellipsoid.py

    r2d81cfe r71b751d  
    145145 
    146146source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
     147have_Fq = True 
    147148 
    148149def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 

  • sasmodels/models/core_shell_parallelepiped.c

    rdbf1a60 r71b751d  
    2727} 
    2828 
    29 static double 
    30 Iq(double q, 
     29static void 
     30Fq(double q, 
     31    double *F1, 
     32    double *F2, 
    3133    double core_sld, 
    3234    double arim_sld, 
    … …  
    6062    // outer integral (with gauss points), integration limits = 0, 1 
    6163    // substitute d_cos_alpha for sin_alpha d_alpha 
    62     double outer_sum = 0; //initialize integral 
     64    double outer_sum_F1 = 0; //initialize integral 
     65    double outer_sum_F2 = 0; //initialize integral 
    6366    for( int i=0; i<GAUSS_N; i++) { 
    6467        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    … …  
    6972        // inner integral (with gauss points), integration limits = 0, 1 
    7073        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta) 
    71         double inner_sum = 0.0; 
     74        double inner_sum_F1 = 0.0; 
     75        double inner_sum_F2 = 0.0; 
    7276        for(int j=0; j<GAUSS_N; j++) { 
    7377            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    … …  
    9195#endif 
    9296 
    93             inner_sum += GAUSS_W[j] * f * f; 
     97            inner_sum_F1 += GAUSS_W[j] * f; 
     98            inner_sum_F2 += GAUSS_W[j] * f * f; 
    9499        } 
    95100        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    96         inner_sum *= 0.5; 
    97         // now sum up the outer integral 
    98         outer_sum += GAUSS_W[i] * inner_sum; 
     101        // and sum up the outer integral 
     102        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * 0.5; 
     103        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * 0.5; 
    99104    } 
    100105    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    101     outer_sum *= 0.5; 
     106    outer_sum_F1 *= 0.5; 
     107    outer_sum_F2 *= 0.5; 
    102108 
    103109    //convert from [1e-12 A-1] to [cm-1] 
    104     return 1.0e-4 * outer_sum; 
     110    *F1 = 1.0e-2 * outer_sum_F1; 
     111    *F2 = 1.0e-4 * outer_sum_F2; 
    105112} 
    106113 

  • sasmodels/models/core_shell_parallelepiped.py

    rf89ec96 r71b751d  
    226226 
    227227source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] 
     228have_Fq = True 
    228229 
    229230 

  • sasmodels/models/core_shell_sphere.c

    r3a48772 r71b751d  
    1 double form_volume(double radius, double thickness); 
    2 double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld); 
     1static double 
     2form_volume(double radius, double thickness) 
     3{ 
     4    return M_4PI_3 * cube(radius + thickness); 
     5} 
    36 
    4 double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) { 
    5  
    6  
    7     double intensity = core_shell_kernel(q, 
     7static void 
     8Fq(double q, double *F1, double *F2, double radius, 
     9   double thickness, double core_sld, double shell_sld, double solvent_sld) { 
     10    double form = core_shell_fq(q, 
    811                              radius, 
    912                              thickness, 
    … …  
    1114                              shell_sld, 
    1215                              solvent_sld); 
    13     return intensity; 
     16    *F1 = 1.0e-2*form; 
     17    *F2 = 1.0e-4*form*form; 
    1418} 
    15  
    16 double form_volume(double radius, double thickness) 
    17 { 
    18     return M_4PI_3 * cube(radius + thickness); 
    19 } 

  • sasmodels/models/core_shell_sphere.py

    rc036ddb r71b751d  
    7777 
    7878source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 
     79have_Fq = True 
    7980 
    8081demo = dict(scale=1, background=0, radius=60, thickness=10, 

  • sasmodels/models/cylinder.c

    r108e70e r71b751d  
    88 
    99static double 
    10 fq(double qab, double qc, double radius, double length) 
     10_fq(double qab, double qc, double radius, double length) 
    1111{ 
    1212    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
    1313} 
    1414 
    15 static double 
    16 orient_avg_1D(double q, double radius, double length) 
     15static void 
     16Fq(double q, 
     17    double *F1, 
     18    double *F2, 
     19    double sld, 
     20    double solvent_sld, 
     21    double radius, 
     22    double length) 
    1723{ 
    1824    // translate a point in [-1,1] to a point in [0, pi/2] 
    … …  
    2026    const double zb = M_PI_4; 
    2127 
    22     double total = 0.0; 
     28    double total_F1 = 0.0; 
     29    double total_F2 = 0.0; 
    2330    for (int i=0; i<GAUSS_N ;i++) { 
    2431        const double theta = GAUSS_Z[i]*zm + zb; 
    … …  
    2633        // theta (theta,phi) the projection of the cylinder on the detector plane 
    2734        SINCOS(theta , sin_theta, cos_theta); 
    28         const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
    29         total += GAUSS_W[i] * form * form * sin_theta; 
     35        const double form = _fq(q*sin_theta, q*cos_theta, radius, length); 
     36        total_F1 += GAUSS_W[i] * form * sin_theta; 
     37        total_F2 += GAUSS_W[i] * form * form * sin_theta; 
    3038    } 
    3139    // translate dx in [-1,1] to dx in [lower,upper] 
    32     return total*zm; 
     40    total_F1 *= zm; 
     41    total_F2 *= zm; 
     42    const double s = (sld - solvent_sld) * form_volume(radius, length); 
     43    *F1 = 1e-2 * s * total_F1; 
     44    *F2 = 1e-4 * s * s * total_F2; 
    3345} 
    3446 
    35 static double 
    36 Iq(double q, 
    37     double sld, 
    38     double solvent_sld, 
    39     double radius, 
    40     double length) 
    41 { 
    42     const double s = (sld - solvent_sld) * form_volume(radius, length); 
    43     return 1.0e-4 * s * s * orient_avg_1D(q, radius, length); 
    44 } 
     47 
    4548 
    4649static double 
    … …  
    5154    double length) 
    5255{ 
     56    const double form = _fq(qab, qc, radius, length); 
    5357    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    54     const double form = fq(qab, qc, radius, length); 
    5558    return 1.0e-4 * square(s * form); 
    5659} 
     60 

  • sasmodels/models/cylinder.py

    r2d81cfe r71b751d  
    138138 
    139139source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
     140have_Fq = True 
    140141 
    141142def ER(radius, length): 

  • sasmodels/models/ellipsoid.c

    rc036ddb r71b751d  
    44    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    55} 
    6  
    7 /* Fq overrides Iq 
    8 static  double 
    9 Iq(double q, 
    10     double sld, 
    11     double sld_solvent, 
    12     double radius_polar, 
    13     double radius_equatorial) 
    14 { 
    15     // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955) 
    16     //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT 
    17     //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT 
    18     //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT 
    19     // u-substitution of 
    20     //     u = sin, du = cos dT 
    21     //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
    22     const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
    23  
    24     // translate a point in [-1,1] to a point in [0, 1] 
    25     // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    26     const double zm = 0.5; 
    27     const double zb = 0.5; 
    28     double total = 0.0; 
    29     for (int i=0;i<GAUSS_N;i++) { 
    30         const double u = GAUSS_Z[i]*zm + zb; 
    31         const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
    32         const double f = sas_3j1x_x(q*r); 
    33         total += GAUSS_W[i] * f * f; 
    34     } 
    35     // translate dx in [-1,1] to dx in [lower,upper] 
    36     const double form = total*zm; 
    37     const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    38     return 1.0e-4 * s * s * form; 
    39 } 
    40 */ 
    416 
    427static void 
    … …  
    7136    } 
    7237    // translate dx in [-1,1] to dx in [lower,upper] 
    73     const double form_squared_avg = total_F2*zm; 
    74     const double form_avg = total_F1*zm; 
     38    total_F1 *= zm; 
     39    total_F2 *= zm; 
    7540    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    76     *F1 = 1e-2 * s * form_avg; 
    77     *F2 = 1e-4 * s * s * form_squared_avg; 
     41    *F1 = 1e-2 * s * total_F1; 
     42    *F2 = 1e-4 * s * s * total_F2; 
    7843} 
    7944 

  • sasmodels/models/ellipsoid.py

    rc036ddb r71b751d  
    163163 
    164164source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
    165  
    166165have_Fq = True 
    167166 

  • sasmodels/models/elliptical_cylinder.c

    r108e70e r71b751d  
    55} 
    66 
    7 static double 
    8 Iq(double q, double radius_minor, double r_ratio, double length, 
     7static void 
     8Fq(double q, double *F1, double *F2, double radius_minor, double r_ratio, double length, 
    99   double sld, double solvent_sld) 
    1010{ 
    … …  
    2121 
    2222    //initialize integral 
    23     double outer_sum = 0.0; 
     23    double outer_sum_F1 = 0.0; 
     24    double outer_sum_F2 = 0.0; 
    2425    for(int i=0;i<GAUSS_N;i++) { 
    2526        //setup inner integral over the ellipsoidal cross-section 
    … …  
    2728        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    2829        //const double arg = radius_minor*sin_val; 
    29         double inner_sum=0; 
     30        double inner_sum_F1 = 0.0; 
     31        double inner_sum_F2 = 0.0; 
    3032        for(int j=0;j<GAUSS_N;j++) { 
    3133            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3234            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3335            const double be = sas_2J1x_x(q*r); 
    34             inner_sum += GAUSS_W[j] * be * be; 
     36            inner_sum_F1 += GAUSS_W[j] * be; 
     37            inner_sum_F2 += GAUSS_W[j] * be * be; 
    3538        } 
    3639        //now calculate the value of the inner integral 
    37         inner_sum *= 0.5*(vbj-vaj); 
     40        inner_sum_F1 *= 0.5*(vbj-vaj); 
     41        inner_sum_F2 *= 0.5*(vbj-vaj); 
    3842 
    3943        //now calculate outer integral 
    4044        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    41         outer_sum += GAUSS_W[i] * inner_sum * si * si; 
     45        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * si; 
     46        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * si * si; 
    4247    } 
    43     outer_sum *= 0.5*(vb-va); 
    44  
    45     //divide integral by Pi 
    46     const double form = outer_sum/M_PI; 
     48    // correct limits and divide integral by pi 
     49    outer_sum_F1 *= 0.5*(vb-va)/M_PI; 
     50    outer_sum_F2 *= 0.5*(vb-va)/M_PI; 
    4751 
    4852    // scale by contrast and volume, and convert to to 1/cm units 
    49     const double vol = form_volume(radius_minor, r_ratio, length); 
    50     const double delrho = sld - solvent_sld; 
    51     return 1.0e-4*square(delrho*vol)*form; 
     53    const double volume = form_volume(radius_minor, r_ratio, length); 
     54    const double contrast = sld - solvent_sld; 
     55    const double s = contrast*volume; 
     56    *F1 = 1.0e-2*s*outer_sum_F1; 
     57    *F2 = 1.0e-4*s*s*outer_sum_F2; 
    5258} 
    5359 
    … …  
    6369    const double be = sas_2J1x_x(qr); 
    6470    const double si = sas_sinx_x(qc*0.5*length); 
    65     const double Aq = be * si; 
    66     const double delrho = sld - solvent_sld; 
    67     const double vol = form_volume(radius_minor, r_ratio, length); 
    68     return 1.0e-4 * square(delrho * vol * Aq); 
     71    const double fq = be * si; 
     72    const double contrast = sld - solvent_sld; 
     73    const double volume = form_volume(radius_minor, r_ratio, length); 
     74    return 1.0e-4 * square(contrast * volume * fq); 
    6975} 

  • sasmodels/models/elliptical_cylinder.py

    ra261a83 r71b751d  
    122122 
    123123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
     124have_Fq = True 
    124125 
    125126demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 

  • sasmodels/models/fcc_paracrystal.c

    r108e70e r71b751d  
    7979} 
    8080 
    81  
    8281static double Iqabc(double qa, double qb, double qc, 
    8382    double dnn, double d_factor, double radius, 

  • sasmodels/models/flexible_cylinder_elliptical.c

    r74768cb r71b751d  
    7272    return result; 
    7373} 
    74  

  • sasmodels/models/fractal.c

    r4788822 r71b751d  
    1919 
    2020    //calculate P(q) for the spherical subunits 
    21     const double pq = square(form_volume(radius) * (sld_block-sld_solvent) 
    22                       *sas_3j1x_x(q*radius)); 
     21    const double fq = form_volume(radius) * (sld_block-sld_solvent) 
     22                      *sas_3j1x_x(q*radius); 
    2323 
    2424    // scale to units cm-1 sr-1 (assuming data on absolute scale) 
    … …  
    2727    //    combined: 1e-4 * I(q) 
    2828 
    29     return 1.e-4 * volfraction * sq * pq; 
     29    return 1.e-4 * volfraction * sq * fq * fq; 
    3030} 
    31  

  • sasmodels/models/fractal_core_shell.c

    rbdd08df r71b751d  
    1616   double cor_length) 
    1717{ 
    18     //The radius for the building block of the core shell particle that is  
     18    //The radius for the building block of the core shell particle that is 
    1919    //needed by the Teixeira fractal S(q) is the radius of the whole particle. 
    2020    const double cs_radius = radius + thickness; 
    2121    const double sq = fractal_sq(q, cs_radius, fractal_dim, cor_length); 
    22     const double pq = core_shell_kernel(q, radius, thickness, 
    23                                         core_sld, shell_sld, solvent_sld); 
     22    const double fq = core_shell_fq(q, radius, thickness, 
     23                                    core_sld, shell_sld, solvent_sld); 
    2424 
    2525    // Note: core_shell_kernel already performs the 1e-4 unit conversion 
    26     return volfraction * sq * pq; 
     26    return 1.0e-4 * volfraction * sq * fq * fq; 
    2727} 
    28  

  • sasmodels/models/fuzzy_sphere.py

    r2d81cfe r71b751d  
    8383 
    8484source = ["lib/sas_3j1x_x.c"] 
     85have_Fq = True 
    8586 
    86 # No volume normalization despite having a volume parameter 
    87 # This should perhaps be volume normalized? 
    88 form_volume = """ 
     87c_code = """ 
     88static double form_volume(double radius) 
     89{ 
    8990    return M_4PI_3*cube(radius); 
    90     """ 
     91} 
    9192 
    92 Iq = """ 
     93static void Fq(double q, double *F1, double *F2, double sld, double sld_solvent, 
     94               double radius, double fuzziness) 
     95{ 
    9396    const double qr = q*radius; 
    9497    const double bes = sas_3j1x_x(qr); 
    95     const double qf = q*fuzziness; 
    96     const double fq = bes * (sld - sld_solvent) * form_volume(radius) * exp(-0.5*qf*qf); 
    97     return 1.0e-4*fq*fq; 
    98     """ 
     98    const double qf = exp(-0.5*square(q*fuzziness)); 
     99    const double contrast = (sld - sld_solvent); 
     100    const double form = contrast * form_volume(radius) * bes * qf; 
     101    *F1 = 1.0e-2*form; 
     102    *F2 = 1.0e-4*form*form; 
     103} 
     104""" 
    99105 
    100106def ER(radius): 

  • sasmodels/models/hardsphere.py

    r2d81cfe r71b751d  
    6262 
    6363#             ["name", "units", default, [lower, upper], "type","description"], 
    64 parameters = [["radius_effective", "Ang", 50.0, [0, inf], "volume", 
     64parameters = [["radius_effective", "Ang", 50.0, [0, inf], "", 
    6565               "effective radius of hard sphere"], 
    6666              ["volfraction", "", 0.2, [0, 0.74], "", 
    6767               "volume fraction of hard spheres"], 
    6868             ] 
    69  
    70 # No volume normalization despite having a volume parameter 
    71 # This should perhaps be volume normalized? 
    72 form_volume = """ 
    73     return 1.0; 
    74     """ 
    7569 
    7670Iq = r""" 

  • sasmodels/models/hollow_cylinder.c

    r108e70e r71b751d  
    11//#define INVALID(v) (v.radius_core >= v.radius) 
    2  
    3 // From Igor library 
    4 static double 
    5 _hollow_cylinder_scaling(double integrand, double delrho, double volume) 
    6 { 
    7     return 1.0e-4 * square(volume * delrho) * integrand; 
    8 } 
    92 
    103static double 
    … …  
    3023 
    3124 
    32 static double 
    33 Iq(double q, double radius, double thickness, double length, 
     25static void 
     26Fq(double q, double *F1, double *F2, double radius, double thickness, double length, 
    3427    double sld, double solvent_sld) 
    3528{ 
    … …  
    3730    const double upper = 1.0;        //limits of numerical integral 
    3831 
    39     double summ = 0.0;            //initialize intergral 
     32    double total_F1 = 0.0;            //initialize intergral 
     33    double total_F2 = 0.0; 
    4034    for (int i=0;i<GAUSS_N;i++) { 
    4135        const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper ); 
    … …  
    4337        const double form = _fq(q*sin_theta, q*cos_theta, 
    4438                                radius, thickness, length); 
    45         summ += GAUSS_W[i] * form * form; 
     39        total_F1 += GAUSS_W[i] * form; 
     40        total_F2 += GAUSS_W[i] * form * form; 
    4641    } 
     42    total_F1 *= 0.5*(upper-lower); 
     43    total_F2 *= 0.5*(upper-lower); 
     44    const double s = (sld - solvent_sld) * form_volume(radius, thickness, length); 
     45    *F1 = 1e-2 * s * total_F1; 
     46    *F2 = 1e-4 * s*s * total_F2; 
     47} 
    4748 
    48     const double Aq = 0.5*summ*(upper-lower); 
    49     const double volume = form_volume(radius, thickness, length); 
    50     return _hollow_cylinder_scaling(Aq, solvent_sld - sld, volume); 
    51 } 
    5249 
    5350static double 
    … …  
    5754{ 
    5855    const double form = _fq(qab, qc, radius, thickness, length); 
    59  
    60     const double vol = form_volume(radius, thickness, length); 
    61     return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
     56    const double s = (sld - solvent_sld) * form_volume(radius, thickness, length); 
     57    return 1.0e-4*square(s * form); 
    6258} 

  • sasmodels/models/hollow_cylinder.py

    r2d81cfe r71b751d  
    8989 
    9090source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
     91have_Fq = True 
    9192 
    9293# pylint: disable=W0613 

  • sasmodels/models/hollow_rectangular_prism.c

    rd86f0fc r71b751d  
    1313} 
    1414 
    15 static double 
    16 Iq(double q, 
     15static void 
     16Fq(double q, 
     17    double *F1, 
     18    double *F2, 
    1719    double sld, 
    1820    double solvent_sld, 
    … …  
    3638    const double v2b = M_PI_2;  //phi integration limits 
    3739 
    38     double outer_sum = 0.0; 
     40    double outer_sum_F1 = 0.0; 
     41    double outer_sum_F2 = 0.0; 
    3942    for(int i=0; i<GAUSS_N; i++) { 
    4043 
    … …  
    4649        const double termC2 = sas_sinx_x(q * (c_half-thickness)*cos(theta)); 
    4750 
    48         double inner_sum = 0.0; 
     51        double inner_sum_F1 = 0.0; 
     52        double inner_sum_F2 = 0.0; 
    4953        for(int j=0; j<GAUSS_N; j++) { 
    5054 
    … …  
    6468            const double AP2 = vol_core * termA2 * termB2 * termC2; 
    6569 
    66             inner_sum += GAUSS_W[j] * square(AP1-AP2); 
     70            inner_sum_F1 += GAUSS_W[j] * (AP1-AP2); 
     71            inner_sum_F2 += GAUSS_W[j] * square(AP1-AP2); 
    6772        } 
    68         inner_sum *= 0.5 * (v2b-v2a); 
     73        inner_sum_F1 *= 0.5 * (v2b-v2a); 
     74        inner_sum_F2 *= 0.5 * (v2b-v2a); 
    6975 
    70         outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
     76        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * sin(theta); 
     77        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * sin(theta); 
    7178    } 
    72     outer_sum *= 0.5*(v1b-v1a); 
     79    outer_sum_F1 *= 0.5*(v1b-v1a); 
     80    outer_sum_F2 *= 0.5*(v1b-v1a); 
    7381 
    7482    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 
    7583    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 
    76     const double form = outer_sum/M_PI_2; 
     84    const double form_avg = outer_sum_F1/M_PI_2; 
     85    const double form_squared_avg = outer_sum_F2/M_PI_2; 
    7786 
    7887    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
    79     const double delrho = sld - solvent_sld; 
     88    const double contrast = sld - solvent_sld; 
    8089 
    8190    // Convert from [1e-12 A-1] to [cm-1] 
    82     return 1.0e-4 * delrho * delrho * form; 
     91    *F1 = 1.0e-2 * contrast * form_avg; 
     92    *F2 = 1.0e-4 * contrast * contrast * form_squared_avg; 
    8393} 
    8494 

  • sasmodels/models/hollow_rectangular_prism.py

    r0e55afe r71b751d  
    142142 
    143143source = ["lib/gauss76.c", "hollow_rectangular_prism.c"] 
     144have_Fq = True 
    144145 
    145146def ER(length_a, b2a_ratio, c2a_ratio, thickness): 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rd86f0fc r71b751d  
    88} 
    99 
    10 static double 
    11 Iq(double q, 
     10static void 
     11Fq(double q, 
     12    double *F1, 
     13    double *F2, 
    1214    double sld, 
    1315    double solvent_sld, 
    … …  
    2830    const double v2b = M_PI_2;  //phi integration limits 
    2931 
    30     double outer_sum = 0.0; 
     32    double outer_sum_F1 = 0.0; 
     33    double outer_sum_F2 = 0.0; 
    3134    for(int i=0; i<GAUSS_N; i++) { 
    3235        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
    … …  
    4144        const double termAT_theta = 8.0 * sin_c / (q*q*sin_theta*cos_theta); 
    4245 
    43         double inner_sum = 0.0; 
     46        double inner_sum_F1 = 0.0; 
     47        double inner_sum_F2 = 0.0; 
    4448        for(int j=0; j<GAUSS_N; j++) { 
    4549            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
    … …  
    6064                * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 
    6165 
    62             inner_sum += GAUSS_W[j] * square(AL+AT); 
     66            inner_sum_F1 += GAUSS_W[j] * (AL+AT); 
     67            inner_sum_F2 += GAUSS_W[j] * square(AL+AT); 
    6368        } 
    6469 
    65         inner_sum *= 0.5 * (v2b-v2a); 
    66         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     70        inner_sum_F1 *= 0.5 * (v2b-v2a); 
     71        inner_sum_F2 *= 0.5 * (v2b-v2a); 
     72        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * sin_theta; 
     73        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * sin_theta; 
    6774    } 
    6875 
    69     outer_sum *= 0.5*(v1b-v1a); 
     76    outer_sum_F1 *= 0.5*(v1b-v1a); 
     77    outer_sum_F2 *= 0.5*(v1b-v1a); 
    7078 
    7179    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 
    7280    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 
    73     double answer = outer_sum/M_PI_2; 
     81    const double form_avg = outer_sum_F1/M_PI_2; 
     82    const double form_squared_avg = outer_sum_F2/M_PI_2; 
    7483 
    7584    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
    76     answer *= square(sld-solvent_sld); 
     85    const double contrast = sld - solvent_sld; 
    7786 
    7887    // Convert from [1e-12 A-1] to [cm-1] 
    79     answer *= 1.0e-4; 
    80  
    81     return answer; 
     88    *F1 = 1e-2 * contrast * form_avg; 
     89    *F2 = 1e-4 * contrast * contrast * form_squared_avg; 
    8290} 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    r2d81cfe r71b751d  
    102102 
    103103source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"] 
     104have_Fq = True 
    104105 
    105106def ER(length_a, b2a_ratio, c2a_ratio): 

  • sasmodels/models/lamellar_hg_stack_caille.c

    r1c6286d r71b751d  
    5656  return inten; 
    5757} 
    58  

  • sasmodels/models/lamellar_stack_caille.c

    r1c6286d r71b751d  
    5252  return(inten); 
    5353} 
    54  

  • sasmodels/models/lamellar_stack_paracrystal.c

    r5467cd8 r71b751d  
    1818        int n2 = n1 + 1; 
    1919        const double xn = (double)n2 - fp_Nlayers;      //fractional contribution of n1 
    20          
     20 
    2121        const double ww = exp(-0.5*square(qval*pd*davg)); 
    2222 
    … …  
    2727 
    2828        Znq = xn*Snq; 
    29          
     29 
    3030        //calculate the n2 contribution 
    3131        an = paraCryst_an(ww,qval,davg,n2); 
    … …  
    3333 
    3434        Znq += (1.0-xn)*Snq; 
    35          
     35 
    3636        //and the independent contribution 
    3737        Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg)); 
    38          
     38 
    3939        //the limit when Nlayers approaches infinity 
    4040//      Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg)) 
    41          
     41 
    4242        const double xi = th/2.0;               //use 1/2 the bilayer thickness 
    4343        const double Pbil = square(sas_sinx_x(qval*xi)); 
    44          
     44 
    4545        const double contr = sld - solvent_sld; 
    4646        const double inten = 2.0*M_PI*contr*contr*Pbil*Znq/(qval*qval); 
    … …  
    5252double 
    5353paraCryst_sn(double ww, double qval, double davg, int Nlayers, double an) { 
    54          
     54 
    5555        double Snq; 
    5656 
    5757        Snq = an/( (double)Nlayers*square(1.0+ww*ww-2.0*ww*cos(qval*davg)) ); 
    58          
     58 
    5959        return Snq; 
    6060} 
    … …  
    6363paraCryst_an(double ww, double qval, double davg, int Nlayers) { 
    6464        double an; 
    65          
     65 
    6666        an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg); 
    6767        an -= 4.0*pow(ww,(Nlayers+2))*cos((double)Nlayers*qval*davg); 
    6868        an += 2.0*pow(ww,(Nlayers+3))*cos((double)(Nlayers-1)*qval*davg); 
    6969        an += 2.0*pow(ww,(Nlayers+1))*cos((double)(Nlayers+1)*qval*davg); 
    70          
     70 
    7171        return an; 
    7272} 

  • sasmodels/models/lib/core_shell.c

    r925ad6e r71b751d  
    77********************************************************************/ 
    88static 
    9 double core_shell_kernel(double q, 
     9double core_shell_fq(double q, 
    1010                         double radius, 
    1111                         double thickness, 
    1212                         double core_sld, 
    1313                         double shell_sld, 
    14                          double solvent_sld) { 
     14                         double solvent_sld) 
     15{ 
    1516    // Core first, then add in shell 
    1617    const double core_qr = q * radius; 
    … …  
    2627    const double shell_volume = M_4PI_3 * cube(radius + thickness); 
    2728    f += shell_volume * shell_bes * shell_contrast; 
    28     return f * f * 1.0e-4; 
     29    return f; 
    2930} 
     31 
     32// Deprecated function: use core_shell_fq instead. 
     33static 
     34double core_shell_kernel(double q, 
     35                         double radius, 
     36                         double thickness, 
     37                         double core_sld, 
     38                         double shell_sld, 
     39                         double solvent_sld) 
     40{ 
     41    const double fq = core_shell_fq(q, radius, thickness, core_sld, shell_sld, solvent_sld); 
     42    return 1.0e-4 * fq*fq; 
     43} 

  • sasmodels/models/linear_pearls.c

    rc3ccaec r71b751d  
    4646    double psi = sas_3j1x_x(q * radius); 
    4747 
    48     // N pearls interaction terms  
     48    // N pearls interaction terms 
    4949    double structure_factor = (double)num_pearls; 
    5050    for(int num=1; num<num_pearls; num++) { 

  • sasmodels/models/mass_surface_fractal.py

    r7994359 r71b751d  
    3838 
    3939    The surface ( $D_s$ ) and mass ( $D_m$ ) fractal dimensions are only 
    40     valid if $0 < surface\_dim < 6$ , $0 < mass\_dim < 6$ , and 
    41     $(surface\_dim + mass\_dim ) < 6$ .  
     40    valid if $0 < surface\_dim < 6$, $0 < mass\_dim < 6$, and 
     41    $(surface\_dim + mass\_dim ) < 6$. 
    4242    Older versions of sasview may have the default primary particle radius 
    43     larger than the cluster radius, this was an error, also present in the  
    44     Schmidt review paper below. The primary particle should be the smaller  
    45     as described in the original Hurd et.al. who also point out that  
    46     polydispersity in the primary particle sizes may affect their  
     43    larger than the cluster radius, this was an error, also present in the 
     44    Schmidt review paper below. The primary particle should be the smaller 
     45    as described in the original Hurd, et al., who also point out that 
     46    polydispersity in the primary particle sizes may affect their 
    4747    apparent surface fractal dimension. 
    48      
     48 
    4949 
    5050References 

  • sasmodels/models/multilayer_vesicle.c

    rec1d4bc r71b751d  
    1212static double 
    1313multilayer_vesicle_kernel(double q, 
    14           double volfraction, 
    1514          double radius, 
    1615          double thick_shell, 
    … …  
    4544                               //unilamellar vesicles (C. Glinka, 11/24/03) 
    4645 
    47     return 1.0e-4*volfraction*fval*fval;  // Volume normalization happens in caller 
     46    return fval;  // Volume normalization happens in caller 
    4847} 
    4948 
    50 static double 
    51 Iq(double q, 
     49static void 
     50Fq(double q, 
     51          double *F1, 
     52          double *F2, 
    5253          double volfraction, 
    5354          double radius, 
    … …  
    5960{ 
    6061    int n_shells = (int)(fp_n_shells + 0.5); 
    61     return multilayer_vesicle_kernel(q, 
    62            volfraction, 
     62    const double fq = multilayer_vesicle_kernel(q, 
    6363           radius, 
    6464           thick_shell, 
    … …  
    6767           sld, 
    6868           n_shells); 
     69    // See comment in vesicle.c regarding volfraction normalization. 
     70    *F1 = 1.0e-2 * sqrt(volfraction)*fq; 
     71    *F2 = 1.0e-4 * volfraction*fq*fq; 
    6972} 
    7073 

  • sasmodels/models/multilayer_vesicle.py

    ref07e95 r71b751d  
    145145 
    146146source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
     147have_Fq = True 
    147148 
    148149def ER(radius, thick_shell, thick_solvent, n_shells): 

  • sasmodels/models/onion.c

    rc3ccaec r71b751d  
    4040} 
    4141 
    42 static double 
    43 Iq(double q, double sld_core, double radius_core, double sld_solvent, 
     42static void 
     43Fq(double q, double *F1, double *F2, double sld_core, double radius_core, double sld_solvent, 
    4444    double n_shells, double sld_in[], double sld_out[], double thickness[], 
    4545    double A[]) 
    … …  
    5555  } 
    5656  f -= f_exp(q, r_out, sld_solvent, 0.0, 0.0, 0.0, 0.0); 
    57   const double f2 = f * f * 1.0e-4; 
    5857 
    59   return f2; 
     58  *F1 = 1e-2 * f; 
     59  *F2 = 1e-4 * f * f; 
    6060} 

  • sasmodels/models/onion.py

    r2d81cfe r71b751d  
    315315source = ["lib/sas_3j1x_x.c", "onion.c"] 
    316316single = False 
     317have_Fq = True 
    317318 
    318319profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 

  • sasmodels/models/parallelepiped.c

    rdbf1a60 r71b751d  
    66 
    77 
    8 static double 
    9 Iq(double q, 
     8static void 
     9Fq(double q, 
     10    double *F1, 
     11    double *F2, 
    1012    double sld, 
    1113    double solvent_sld, 
    … …  
    2123 
    2224    // outer integral (with gauss points), integration limits = 0, 1 
    23     double outer_total = 0; //initialize integral 
    24  
     25    double outer_total_F1 = 0.0; //initialize integral 
     26    double outer_total_F2 = 0.0; //initialize integral 
    2527    for( int i=0; i<GAUSS_N; i++) { 
    2628        const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    … …  
    2931        // inner integral (with gauss points), integration limits = 0, 1 
    3032        // corresponding to angles from 0 to pi/2. 
    31         double inner_total = 0.0; 
     33        double inner_total_F1 = 0.0; 
     34        double inner_total_F2 = 0.0; 
    3235        for(int j=0; j<GAUSS_N; j++) { 
    3336            const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    … …  
    3639            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    3740            const double si2 = sas_sinx_x(mu_proj * cos_uu); 
    38             inner_total += GAUSS_W[j] * square(si1 * si2); 
     41            const double fq = si1 * si2; 
     42            inner_total_F1 += GAUSS_W[j] * fq; 
     43            inner_total_F2 += GAUSS_W[j] * fq * fq; 
    3944        } 
    4045        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    41         inner_total *= 0.5; 
     46        inner_total_F1 *= 0.5; 
     47        inner_total_F2 *= 0.5; 
    4248 
    4349        const double si = sas_sinx_x(mu * c_scaled * sigma); 
    44         outer_total += GAUSS_W[i] * inner_total * si * si; 
     50        outer_total_F1 += GAUSS_W[i] * inner_total_F1 * si; 
     51        outer_total_F2 += GAUSS_W[i] * inner_total_F2 * si * si; 
    4552    } 
    4653    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    47     outer_total *= 0.5; 
     54    outer_total_F1 *= 0.5; 
     55    outer_total_F2 *= 0.5; 
    4856 
    4957    // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1] 
    5058    const double V = form_volume(length_a, length_b, length_c); 
    51     const double drho = (sld-solvent_sld); 
    52     return 1.0e-4 * square(drho * V) * outer_total; 
     59    const double contrast = (sld-solvent_sld); 
     60    const double s = contrast * V; 
     61    *F1 = 1.0e-2 * s * outer_total_F1; 
     62    *F2 = 1.0e-4 * s * s * outer_total_F2; 
    5363} 
    54  
    5564 
    5665static double 

  • sasmodels/models/parallelepiped.py

    rf89ec96 r71b751d  
    230230 
    231231source = ["lib/gauss76.c", "parallelepiped.c"] 
     232have_Fq = True 
    232233 
    233234def ER(length_a, length_b, length_c): 

  • sasmodels/models/raspberry.c

    r925ad6e r71b751d  
    11double form_volume(double radius_lg); 
    22 
    3 double Iq(double q,  
     3double Iq(double q, 
    44          double sld_lg, double sld_sm, double sld_solvent, 
    55          double volfraction_lg, double volfraction_sm, double surf_fraction, 
    … …  
    2727    double sfLS,sfSS; 
    2828    double slT; 
    29   
     29 
    3030    vfL = volfraction_lg; 
    3131    rL = radius_lg; 
    … …  
    3737    deltaS = penetration; 
    3838    sldSolv = sld_solvent; 
    39      
     39 
    4040    delrhoL = fabs(sldL - sldSolv); 
    41     delrhoS = fabs(sldS - sldSolv);  
    42       
     41    delrhoS = fabs(sldS - sldSolv); 
     42 
    4343    VL = M_4PI_3*rL*rL*rL; 
    4444    VS = M_4PI_3*rS*rS*rS; 
    … …  
    7676    //I(q) for free small particles 
    7777    f2+= vfS*(1.0-fSs)*delrhoS*delrhoS*VS*psiS*psiS; 
    78      
     78 
    7979    // normalize to single particle volume and convert to 1/cm 
    8080    f2 *= 1.0e8;        // [=] 1/cm 
    8181    f2 *= 1.0e-12;      // convert for (1/A^-6)^2 to (1/A)^2 
    82      
     82 
    8383    return f2; 
    8484} 

  • sasmodels/models/rectangular_prism.c

    rd86f0fc r71b751d  
    6868} 
    6969 
     70static void 
     71Fq(double q, 
     72    double *F1, 
     73    double *F2, 
     74    double sld, 
     75    double solvent_sld, 
     76    double length_a, 
     77    double b2a_ratio, 
     78    double c2a_ratio) 
     79{ 
     80    const double length_b = length_a * b2a_ratio; 
     81    const double length_c = length_a * c2a_ratio; 
     82    const double a_half = 0.5 * length_a; 
     83    const double b_half = 0.5 * length_b; 
     84    const double c_half = 0.5 * length_c; 
     85 
     86   //Integration limits to use in Gaussian quadrature 
     87    const double v1a = 0.0; 
     88    const double v1b = M_PI_2;  //theta integration limits 
     89    const double v2a = 0.0; 
     90    const double v2b = M_PI_2;  //phi integration limits 
     91 
     92    double outer_sum_F1 = 0.0; 
     93    double outer_sum_F2 = 0.0; 
     94    for(int i=0; i<GAUSS_N; i++) { 
     95        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     96        double sin_theta, cos_theta; 
     97        SINCOS(theta, sin_theta, cos_theta); 
     98 
     99        const double termC = sas_sinx_x(q * c_half * cos_theta); 
     100 
     101        double inner_sum_F1 = 0.0; 
     102        double inner_sum_F2 = 0.0; 
     103        for(int j=0; j<GAUSS_N; j++) { 
     104            double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     105            double sin_phi, cos_phi; 
     106            SINCOS(phi, sin_phi, cos_phi); 
     107 
     108            // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 
     109            const double termA = sas_sinx_x(q * a_half * sin_theta * sin_phi); 
     110            const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi); 
     111            const double AP = termA * termB * termC; 
     112            inner_sum_F1 += GAUSS_W[j] * AP; 
     113            inner_sum_F2 += GAUSS_W[j] * AP * AP; 
     114        } 
     115        inner_sum_F1 = 0.5 * (v2b-v2a) * inner_sum_F1; 
     116        inner_sum_F2 = 0.5 * (v2b-v2a) * inner_sum_F2; 
     117        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * sin_theta; 
     118        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * sin_theta; 
     119    } 
     120 
     121    outer_sum_F1 *= 0.5*(v1b-v1a); 
     122    outer_sum_F2 *= 0.5*(v1b-v1a); 
     123 
     124    // Normalize by Pi (Eqn. 16). 
     125    // The term (ABC)^2 does not appear because it was introduced before on 
     126    // the definitions of termA, termB, termC. 
     127    // The factor 2 appears because the theta integral has been defined between 
     128    // 0 and pi/2, instead of 0 to pi. 
     129    outer_sum_F1 /= M_PI_2; 
     130    outer_sum_F2 /= M_PI_2; 
     131 
     132    // Multiply by contrast and volume 
     133    const double s = (sld-solvent_sld) * (length_a * length_b * length_c); 
     134 
     135    // Convert from [1e-12 A-1] to [cm-1] 
     136    *F1 = 1e-2 * s * outer_sum_F1; 
     137    *F2 = 1e-4 * s * s * outer_sum_F2; 
     138} 
     139 
    70140 
    71141static double 
    … …  
    82152    const double b_half = 0.5 * length_b; 
    83153    const double c_half = 0.5 * length_c; 
    84     const double volume = length_a * length_b * length_c; 
    85154 
    86155    // Amplitude AP from eqn. (13) 
    87  
    88156    const double termA = sas_sinx_x(qa * a_half); 
    89157    const double termB = sas_sinx_x(qb * b_half); 
    90158    const double termC = sas_sinx_x(qc * c_half); 
    91  
    92159    const double AP = termA * termB * termC; 
    93160 
    94     // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
    95     const double delrho = sld - solvent_sld; 
     161    // Multiply by contrast and volume 
     162    const double s = (sld-solvent_sld) * (length_a * length_b * length_c); 
    96163 
    97164    // Convert from [1e-12 A-1] to [cm-1] 
    98     return 1.0e-4 * square(volume * delrho * AP); 
     165    return 1.0e-4 * square(s * AP); 
    99166} 

  • sasmodels/models/rectangular_prism.py

    r0e55afe r71b751d  
    135135 
    136136source = ["lib/gauss76.c", "rectangular_prism.c"] 
     137have_Fq = True 
    137138 
    138139def ER(length_a, b2a_ratio, c2a_ratio): 

  • sasmodels/models/rpa.c

    reb63414 r71b751d  
    3939  double S31,S32,S33,S34,S41,S42,S43,S44; 
    4040  double Lad,Lbd,Lcd,Nav,Intg; 
    41    
     41 
    4242  // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc) 
    4343  //icase was shifted to N-1 from the original code 
    … …  
    5858    Kab = Kac = Kad = -0.0004; 
    5959  } 
    60   
     60 
    6161  // Set volume fraction of component D based on constraint that sum of vol frac =1 
    6262  Phi[3]=1.0-Phi[0]-Phi[1]-Phi[2]; 
    … …  
    8484  Phicd=sqrt(Phi[2]*Phi[3]); 
    8585 
    86   // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk  
     86  // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk 
    8787  Xa=q*q*b[0]*b[0]*N[0]/6.0; 
    8888  Xb=q*q*b[1]*b[1]*N[1]/6.0; 
    … …  
    311311  //Note that should multiply by Nav to get units of SLD which will become 
    312312  // Nav*2 in the next line (SLD^2) but then normalization in that line would 
    313   //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring.  
     313  //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring. 
    314314  Nav=6.022045e+23; 
    315315  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); 
    … …  
    320320 
    321321  //rescale for units of Lij^2 (fm^2 to cm^2) 
    322   Intg *= 1.0e-26;     
     322  Intg *= 1.0e-26; 
    323323 
    324324  return Intg; 

  • sasmodels/models/sc_paracrystal.c

    r108e70e r71b751d  
    8080} 
    8181 
    82  
    8382static double 
    8483Iqabc(double qa, double qb, double qc, 

  • sasmodels/models/sphere.py

    rc036ddb r71b751d  
    6767             ] 
    6868 
    69 source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"] 
     69source = ["lib/sas_3j1x_x.c"] 
     70have_Fq = True 
    7071 
    7172c_code = """ 
    7273static double form_volume(double radius) 
    7374{ 
    74     return sphere_volume(radius); 
     75    return M_4PI_3*cube(radius); 
    7576} 
    7677 
    77 static void Fq(double q, double *F1,double *F2, double sld, double solvent_sld, double radius) 
     78static void Fq(double q, double *f1, double *f2, double sld, double sld_solvent, double radius) 
    7879{ 
    79     const double fq = sas_3j1x_x(q*radius); 
    80     const double contrast = (sld - solvent_sld); 
    81     const double form = 1e-2 * contrast * sphere_volume(radius) * fq; 
    82     *F1 = form; 
    83     *F2 = form*form; 
     80    const double bes = sas_3j1x_x(q*radius); 
     81    const double contrast = (sld - sld_solvent); 
     82    const double form = contrast * form_volume(radius) * bes; 
     83    *f1 = 1.0e-2*form; 
     84    *f2 = 1.0e-4*form*form; 
    8485} 
    8586""" 
    86  
    87 # TODO: figure this out by inspection 
    88 have_Fq = True 
    8987 
    9088def ER(radius): 

  • sasmodels/models/spherical_sld.c

    r3330bb4 r71b751d  
    101101} 
    102102 
     103static void Fq( 
     104    double q, 
     105    double *F1, 
     106    double *F2, 
     107    double fp_n_shells, 
     108    double sld_solvent, 
     109    double sld[], 
     110    double thickness[], 
     111    double interface[], 
     112    double shape[], 
     113    double nu[], 
     114    double fp_n_steps) 
     115{ 
     116    // iteration for # of shells + core + solvent 
     117    int n_shells = (int)(fp_n_shells + 0.5); 
     118    int n_steps = (int)(fp_n_steps + 0.5); 
     119    double f=0.0; 
     120    double r=0.0; 
     121    for (int shell=0; shell<n_shells; shell++){ 
     122        const double sld_l = sld[shell]; 
     123 
     124        // uniform shell; r=0 => r^3=0 => f=0, so works for core as well. 
     125        f -= M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r); 
     126        r += thickness[shell]; 
     127        f += M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r); 
     128 
     129        // iterate over sub_shells in the interface 
     130        const double dr = interface[shell]/n_steps; 
     131        const double delta = (shell==n_shells-1 ? sld_solvent : sld[shell+1]) - sld_l; 
     132        const double nu_shell = fmax(fabs(nu[shell]), 1.e-14); 
     133        const int shape_shell = (int)(shape[shell]); 
     134 
     135        // if there is no interface the equations don't work 
     136        if (dr == 0.) continue; 
     137 
     138        double sld_in = sld_l; 
     139        for (int step=1; step <= n_steps; step++) { 
     140            // find sld_i at the outer boundary of sub-shell step 
     141            //const double z = (double)step/(double)n_steps; 
     142            const double z = (double)step/(double)n_steps; 
     143            const double fraction = blend(shape_shell, nu_shell, z); 
     144            const double sld_out = fraction*delta + sld_l; 
     145            // calculate slope 
     146            const double slope = (sld_out - sld_in)/dr; 
     147            const double contrast = sld_in - slope*r; 
     148 
     149            // iteration for the left and right boundary of the shells 
     150            f -= f_linear(q, r, contrast, slope); 
     151            r += dr; 
     152            f += f_linear(q, r, contrast, slope); 
     153            sld_in = sld_out; 
     154        } 
     155    } 
     156    // add in solvent effect 
     157    f -= M_4PI_3 * cube(r) * sld_solvent * sas_3j1x_x(q*r); 
     158 
     159    *F1 = 1e-2*f; 
     160    *F2 = 1e-4*f*f; 
     161} 

  • sasmodels/models/spherical_sld.py

    r2d81cfe r71b751d  
    209209source = ["lib/polevl.c", "lib/sas_erf.c", "lib/sas_3j1x_x.c", "spherical_sld.c"] 
    210210single = False  # TODO: fix low q behaviour 
     211have_Fq = True 
    211212 
    212213profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 

  • sasmodels/models/spinodal.py

    ref07e95 r71b751d  
    6060    :return:               Calculated intensity 
    6161    """ 
    62  
    6362    with errstate(divide='ignore'): 
    6463        x = q/q_0 

  • sasmodels/models/triaxial_ellipsoid.c

    r108e70e r71b751d  
    77} 
    88 
    9 static double 
    10 Iq(double q, 
     9 
     10static void 
     11Fq(double q, 
     12    double *F1, 
     13    double *F2, 
    1114    double sld, 
    1215    double sld_solvent, 
    … …  
    2023    const double zm = M_PI_4; 
    2124    const double zb = M_PI_4; 
    22     double outer = 0.0; 
     25    double outer_sum_F1 = 0.0; 
     26    double outer_sum_F2 = 0.0; 
    2327    for (int i=0;i<GAUSS_N;i++) { 
    2428        //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2; 
    … …  
    2630        const double pa_sinsq_phi = pa*square(sin(phi)); 
    2731 
    28         double inner = 0.0; 
     32        double inner_sum_F1 = 0.0; 
     33        double inner_sum_F2 = 0.0; 
    2934        const double um = 0.5; 
    3035        const double ub = 0.5; 
    … …  
    3439            const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 
    3540            const double fq = sas_3j1x_x(q*r); 
    36             inner += GAUSS_W[j] * fq * fq; 
     41            inner_sum_F1 += GAUSS_W[j] * fq; 
     42            inner_sum_F2 += GAUSS_W[j] * fq * fq; 
    3743        } 
    38         outer += GAUSS_W[i] * inner;  // correcting for dx later 
     44        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1;  // correcting for dx later 
     45        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2;  // correcting for dx later 
    3946    } 
    4047    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    41     const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    42     const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    43     const double drho = (sld - sld_solvent); 
    44     return 1.0e-4 * square(vol*drho) * fqsq; 
     48    outer_sum_F1 *= 0.25;  // = outer*um*zm*8.0/(4.0*M_PI); 
     49    outer_sum_F2 *= 0.25;  // = outer*um*zm*8.0/(4.0*M_PI); 
     50 
     51    const double volume = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     52    const double contrast = (sld - sld_solvent); 
     53    *F1 = 1.0e-2 * contrast * volume * outer_sum_F1; 
     54    *F2 = 1.0e-4 * square(contrast * volume) * outer_sum_F2; 
    4555} 
     56 
    4657 
    4758static double 

  • sasmodels/models/triaxial_ellipsoid.py

    r37f08d2 r71b751d  
    157157 
    158158source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 
     159have_Fq = True 
    159160 
    160161def ER(radius_equat_minor, radius_equat_major, radius_polar): 

  • sasmodels/models/two_lorentzian.py

    ref07e95 r71b751d  
    8787                                  power(q*lorentz_length_2, lorentz_exp_2)) 
    8888# pylint: enable=bad-whitespace 
    89  
    9089    return intensity 
    9190 

  • sasmodels/models/vesicle.c

    r925ad6e r71b751d  
    1 double form_volume(double radius, double thickness); 
    2  
    3 double Iq(double q,  
    4           double sld, double sld_solvent, double volfraction, 
    5           double radius, double thickness); 
    6  
    7 double form_volume(double radius, double thickness) 
     1static double 
     2form_volume(double radius, double thickness) 
    83{ 
    94    //note that for the vesicle model, the volume is ONLY the shell volume 
    … …  
    116} 
    127 
    13 double Iq(double q, 
     8 
     9static void 
     10Fq(double q, 
     11    double *F1, 
     12    double *F2, 
    1413    double sld, 
    1514    double sld_solvent, 
    … …  
    3130    vol = M_4PI_3*cube(radius); 
    3231    f = vol * sas_3j1x_x(q*radius) * contrast; 
    33   
     32 
    3433    //now the shell. No volume normalization as this is done by the caller 
    3534    contrast = sld-sld_solvent; 
    … …  
    3736    f += vol * sas_3j1x_x(q*(radius+thickness)) * contrast; 
    3837 
    39     //rescale to [cm-1].  
    40     f2 = volfraction * f*f*1.0e-4; 
    41      
    42     return f2; 
     38    //rescale to [cm-1]. 
     39    // With volume fraction as part of the model in the dilute limit need 
     40    // to return F2 = Vf <fq^2>.  In order for beta approx. to work correctly 
     41    // need F1^2/F2 equal to <fq>^2 / <fq^2>.  By returning F1 = sqrt(Vf) <fq> 
     42    // and F2 = Vf <fq^2> both conditions are satisfied. 
     43    // Since Vf is the volume fraction of vesicles of all radii, it is 
     44    // constant when averaging F1 and F2 over radii and so pops out of the 
     45    // polydispersity loop, so it is safe to apply it inside the model 
     46    // (albeit conceptually ugly). 
     47    *F1 = 1e-2 * sqrt(volfraction) * f; 
     48    *F2 = 1.0e-4 * volfraction * f * f; 
    4349} 

  • sasmodels/models/vesicle.py

    ref07e95 r71b751d  
    9494 
    9595source = ["lib/sas_3j1x_x.c", "vesicle.c"] 
     96have_Fq = True 
    9697 
    9798def ER(radius, thickness): 
Note: See TracChangeset for help on using the changeset viewer.