Ignore:
Timestamp:
Aug 14, 2009 3:58:58 PM (16 years ago)
Author:
Jae Cho <jhjcho@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
d5bd424
Parents:
cad821b
Message:

working on 2D models. Still need smore corrections and unit tests.

Location:
sansmodels/src/sans/models/c_extensions
Files:
2 added
11 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/sans/models/c_extensions/binaryHs.h

    r3d25331f r975ec8e  
    77        [PYTHONCLASS] = BinaryHSModel 
    88        [DISP_PARAMS] = l_radius,s_radius 
    9         [DESCRIPTION] =<text> 
    10                                                 Model parameters: 
    11                                                  
    12                                                 l_radius : large radius of the binary hard sphere 
    13                                                 s_radius : small radius of the binary hard sphere 
    14                                                 vol_frac_ls : volume fraction of large spheres  
    15                                                 vol_frac_ss : volume fraction of small spheres 
    16                                                 ls_sld: large sphere  scattering length density 
    17                                                 ss_sld: small sphere scattering length density 
    18                                                 solvent_sld: solvent scattering length density 
    19                                                 background: incoherent background 
     9        [DESCRIPTION] =<text> Model parameters: l_radius : large radius of binary hard sphere 
     10                        s_radius : small radius of binary hard sphere 
     11                        vol_frac_ls : volume fraction of large spheres 
     12                        vol_frac_ss : volume fraction of small spheres 
     13                        ls_sld: large sphere  scattering length density 
     14                        ss_sld: small sphere scattering length density 
     15                        solvent_sld: solvent scattering length density 
     16                        background: incoherent background 
    2017               </text> 
    2118        [FIXED]= l_radius.width;s_radius.width 
    … …  
    2320 */ 
    2421typedef struct { 
    25      
     22 
    2623        ///     large radius of the binary hard sphere [A] 
    27     //  [DEFAULT]=l_radius= 160.0 [A] 
     24    //  [DEFAULT]=l_radius= 100.0 [A] 
    2825    double l_radius; 
    2926 
    … …  
    3229    double s_radius; 
    3330 
    34         ///     volume fraction of large spheres  
    35     //  [DEFAULT]=vol_frac_ls= 0.2  
     31        ///     volume fraction of large spheres 
     32    //  [DEFAULT]=vol_frac_ls= 0.1 
    3633    double vol_frac_ls; 
    3734 
    38         ///     volume fraction of small spheres  
    39     //  [DEFAULT]=vol_frac_ss= 0.2  
     35        ///     volume fraction of small spheres 
     36    //  [DEFAULT]=vol_frac_ss= 0.2 
    4037    double vol_frac_ss; 
    4138 

  • sansmodels/src/sans/models/c_extensions/elliptical_cylinder.c

    r3fe701a r975ec8e  
    3131} 
    3232 
    33 double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double psi, double nu) { 
     33double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double nu) { 
    3434        double qr; 
    3535        double qL; 
     36        double Be,Si; 
    3637        double r_major; 
    3738        double kernel; 
    … …  
    4243        qL = q*pars->length*cos(alpha)/2.0; 
    4344 
    44         kernel = 2.0*NR_BessJ1(qr)/qr * sin(qL)/qL; 
     45        if (qr==0){ 
     46                Be = 0.5; 
     47        }else{ 
     48                Be = NR_BessJ1(qr)/qr; 
     49        } 
     50        if (qL==0){ 
     51                Si = 1.0; 
     52        }else{ 
     53                Si = sin(qL)/qL; 
     54        } 
     55 
     56 
     57        kernel = 2.0*Be * Si; 
    4558        return kernel*kernel; 
    4659} 
    … …  
    129142    } 
    130143 
    131         answer = elliptical_cylinder_kernel(pars, q, alpha, pars->cyl_psi,nu); 
     144        answer = elliptical_cylinder_kernel(pars, q, alpha,nu); 
    132145 
    133146        // Multiply by contrast^2 

  • sansmodels/src/sans/models/c_extensions/lamellar.c

    r42f193a r975ec8e  
    2222        // Fill paramater array 
    2323        dp[0] = pars->scale; 
    24         dp[1] = pars->delta; 
    25         dp[2] = pars->sigma; 
    26         dp[3] = pars->contrast; 
     24        dp[1] = pars->bi_thick; 
     25        dp[2] = pars->sld_bi; 
     26        dp[3] = pars->sld_sol; 
    2727        dp[4] = pars->background; 
    2828 
    2929 
    3030        // Call library function to evaluate model 
    31         return LamellarFF(dp, q); 
     31        return lamellar_kernel(dp, q); 
    3232} 
     33 
     34 
    3335/** 
    3436 * Function to evaluate 2D scattering function 

  • sansmodels/src/sans/models/c_extensions/lamellar.h

    r42f193a r975ec8e  
    33/** Structure definition for lamellar parameters 
    44 * [PYTHONCLASS] = LamellarModel 
     5 * [DISP_PARAMS] = bi_thick 
    56   [DESCRIPTION] = <text>[Dilute Lamellar Form Factor](from a lyotropic lamellar phase) 
    6            I(q)= 2*pi*P(q)/(delta *q^(2)), where 
    7                 P(q)=2*(contrast/q)^(2)*(1-cos(q*delta) 
    8                 *e^(1/2*(q*sigma)^(2)). 
    9                 delta = bilayer thickness 
    10                 sigma = variation in bilayer thickness 
    11                         = delta*polydispersity 
    12                 contrast = SLD_solvent - SLD_bilayer 
    13         Note: the polydispersity in delta is included. 
     7                I(q)= 2*pi*P(q)/(delta *q^(2)), where 
     8                P(q)=2*(contrast/q)^(2)*(1-cos(q*delta))^(2)) 
     9                bi_thick = bilayer thickness 
     10                sld_bi = SLD of bilayer 
     11                sld_sol = SLD of solvent 
     12                background = Incoherent background 
     13                scale = scale factor 
     14 
    1415 </text> 
     16[FIXED]= <text>bi_thick.width</text> 
    1517 **/ 
    1618typedef struct { 
    … …  
    1921    double scale; 
    2022    /// delta bilayer thickness [A] 
    21     //  [DEFAULT]=delta=50.0 [A] 
    22     double delta; 
    23     /// variation in bilayer thickness 
    24     //  [DEFAULT]=sigma=0.15 
    25     double sigma; 
    26     /// Contrast [1/A²] 
    27     //  [DEFAULT]=contrast=5.3e-6 [1/A²] 
    28     double contrast; 
     23    //  [DEFAULT]=bi_thick=50.0 [A] 
     24    double bi_thick; 
     25    /// SLD of bilayer [1/A²] 
     26    //  [DEFAULT]=sld_bi=1.0e-6 [1/A²] 
     27    double sld_bi; 
     28    /// SLD of solvent [1/A²] 
     29    //  [DEFAULT]=sld_sol=6.3e-6 [1/A²] 
     30    double sld_sol; 
    2931        /// Incoherent Background [1/cm] 0.00 
    3032        //  [DEFAULT]=background=0.0 [1/cm] 

  • sansmodels/src/sans/models/c_extensions/oblate.c

    r96b59384 r975ec8e  
    6666 * @return: function value 
    6767 */ 
    68 /* 
     68 
    6969double oblate_analytical_2D_scaled(OblateParameters *pars, double q, double q_x, double q_y) { 
    7070 
    7171        return 1.0; 
    7272} 
    73 */ 
     73 

  • sansmodels/src/sans/models/c_extensions/oblate.h

    r96b59384 r975ec8e  
    2424 
    2525   [FIXED] = <text>major_core.width;minor_core.width; major_shell.width; minor_shell.width</text> 
    26    [ORIENTATION_PARAMS] = 
     26   [ORIENTATION_PARAMS]= <text>axis_phi; axis_theta; axis_phi.width; axis_theta.width</text> 
    2727 
    2828 **/ 
    … …  
    5252        //  [DEFAULT]=background=0.001 [1/cm] 
    5353        double background; 
    54         /*//Disable for now 
     54        //Disable for now 
    5555    /// Orientation of the oblate axis w/respect incoming beam [rad] 
    5656    //  [DEFAULT]=axis_theta=1.0 [rad] 
    … …  
    5959    //  [DEFAULT]=axis_phi=1.0 [rad] 
    6060    double axis_phi; 
    61         */ 
     61 
    6262} OblateParameters; 
    6363 

  • sansmodels/src/sans/models/c_extensions/parallelepiped.c

    r5068697 r975ec8e  
    1919double parallelepiped_analytical_1D(ParallelepipedParameters *pars, double q) { 
    2020        double dp[6]; 
    21          
     21 
    2222        // Fill paramater array 
    2323        dp[0] = pars->scale; 
    … …  
    2929 
    3030        // Call library function to evaluate model 
    31         return Parallelepiped(dp, q);    
     31        return Parallelepiped(dp, q); 
    3232} 
     33 
     34 
     35double pkernel(double a, double b,double c, double ala, double alb, double alc){ 
     36    // mu passed in is really mu*sqrt(1-sig^2) 
     37    double argA,argB,argC,tmp1,tmp2,tmp3;                       //local variables 
     38 
     39    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
     40    argA = a*ala; 
     41    argB = b*alb; 
     42    argC = c*alc; 
     43    if(argA==0.0) { 
     44                tmp1 = 1.0; 
     45        } else { 
     46                tmp1 = sin(argA)*sin(argA)/argA/argA; 
     47    } 
     48    if (argB==0.0) { 
     49                tmp2 = 1.0; 
     50        } else { 
     51                tmp2 = sin(argB)*sin(argB)/argB/argB; 
     52    } 
     53 
     54    if (argC==0.0) { 
     55                tmp3 = 1.0; 
     56        } else { 
     57                tmp3 = sin(argC)*sin(argC)/argC/argC; 
     58    } 
     59 
     60    return (tmp1*tmp2*tmp3); 
     61 
     62}//Function pkernel() 
     63 
     64 
     65 
     66 
    3367/** 
    3468 * Function to evaluate 2D scattering function 
    … …  
    4175        q = sqrt(qx*qx+qy*qy); 
    4276    return parallelepiped_analytical_2D_scaled(pars, q, qx/q, qy/q); 
    43 }  
     77} 
    4478 
    4579 
    … …  
    5387double parallelepiped_analytical_2D(ParallelepipedParameters *pars, double q, double phi) { 
    5488    return parallelepiped_analytical_2D_scaled(pars, q, cos(phi), sin(phi)); 
    55 }  
    56          
     89} 
     90 
    5791/** 
    5892 * Function to evaluate 2D scattering function 
    … …  
    6498 */ 
    6599double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y) { 
    66         double parallel_x, parallel_y, parallel_z; 
     100        double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y, parallel_z; 
    67101        double q_z; 
    68         double alpha, vol, cos_val; 
    69         double aa, mu, uu; 
     102        double alpha, vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC; 
     103 
    70104        double answer; 
    71          
    72          
     105        double pi = 4.0*atan(1.0); 
    73106 
    74     // parallelepiped orientation 
    75     parallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi); 
    76     parallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi); 
    77     parallel_z = cos(pars->parallel_theta); 
    78       
     107        edgeA = pars->short_edgeA; 
     108        edgeB = pars->longer_edgeB; 
     109        edgeC = pars->longuest_edgeC; 
     110 
     111 
     112    // parallelepiped c axis orientation 
     113    cparallel_x = sin(pars->parallel_theta) * cos(pars->parallel_phi); 
     114    cparallel_y = sin(pars->parallel_theta) * sin(pars->parallel_phi); 
     115    cparallel_z = cos(pars->parallel_theta); 
     116 
    79117    // q vector 
    80118    q_z = 0; 
    81          
     119 
    82120    // Compute the angle btw vector q and the 
    83121    // axis of the parallelepiped 
    84     cos_val = parallel_x*q_x + parallel_y*q_y + parallel_z*q_z; 
    85      
     122    cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z; 
     123    alpha = acos(cos_val_c); 
     124 
     125    // parallelepiped a axis orientation 
     126    parallel_x = (1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi); 
     127    parallel_y = (1-sin(pars->parallel_theta)*sin(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi); 
     128    cos_val_a = (parallel_x*q_x + parallel_y*q_y); 
     129 
     130 
     131 
     132    // parallelepiped b axis orientation 
     133    bparallel_x = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*cos(pars->parallel_psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi); 
     134    bparallel_y = (1-sin(pars->parallel_theta)*cos(pars->parallel_phi))*sin(pars->parallel_psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi); 
     135    // axis of the parallelepiped 
     136    cos_val_b = (bparallel_x*q_x + bparallel_y*q_y) ; 
     137 
     138 
     139 
    86140    // The following test should always pass 
    87     if (fabs(cos_val)>1.0) { 
     141    if (fabs(cos_val_c)>1.0) { 
    88142        printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    89143        return 0; 
    90144    } 
    91      
    92     // Note: cos(alpha) = 0 and 1 will get an 
    93     // undefined value from PPKernel 
    94         alpha = acos( cos_val ); 
    95145 
    96         aa = pars->short_edgeA/pars->longer_edgeB; 
    97         mu = 1.0; 
    98         uu = 1.0; 
    99          
    100146        // Call the IGOR library function to get the kernel 
    101         answer = PPKernel( aa, mu, uu); 
    102          
     147        answer = pkernel( q*edgeA, q*edgeB, q*edgeC, cos_val_a,cos_val_b,cos_val_c); 
     148 
    103149        // Multiply by contrast^2 
    104150        answer *= pars->contrast*pars->contrast; 
    105          
     151 
    106152        //normalize by cylinder volume 
    107153        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel 
    108     vol = pars->short_edgeA* pars->longer_edgeB * pars->longuest_edgeC; 
     154    vol = 8*edgeA* edgeB * edgeC; 
    109155        answer *= vol; 
    110          
     156 
    111157        //convert to [cm-1] 
    112158        answer *= 1.0e8; 
    113          
     159 
    114160        //Scale 
    115161        answer *= pars->scale; 
    116          
     162 
    117163        // add in the background 
    118164        answer += pars->background; 
    119          
     165 
    120166        return answer; 
    121167} 
    122      
    123168 

  • sansmodels/src/sans/models/c_extensions/parallelepiped.h

    r5068697 r975ec8e  
    77/** Structure definition for Parallelepiped parameters 
    88 * [PYTHONCLASS] = ParallelepipedModel 
    9  * [DISP_PARAMS] = short_edgeA, longer_edgeB, longuest_edgeC,parallel_phi, parallel_theta 
     9 * [DISP_PARAMS] = short_edgeA, longer_edgeB, longuest_edgeC,parallel_phi,parallel_psi, parallel_theta 
    1010   [DESCRIPTION] = <text> Calculates the form factor for a rectangular solid with uniform scattering length density. 
    11                  
     11 
    1212                scale:Scale factor 
    1313                short_edgeA: Shortest edge of the parallelepiped [A] 
    … …  
    1515                longuest_edgeC: Longuest edge of the parallelepiped [A] 
    1616                constrast: particle_sld - solvent_sld 
    17                 background:Incoherent Background [1/cm]  
     17                background:Incoherent Background [1/cm] 
    1818                </text> 
    19         [FIXED]= <text>short_edgeA.width; longer_edgeB.width; longuest_edgeC.width;parallel_phi.width; parallel_theta.width</text> 
    20         [ORIENTATION_PARAMS]= <text>parallel_phi; parallel_theta; parallel_phi.width; parallel_theta.width</text> 
     19        [FIXED]= <text>short_edgeA.width; longer_edgeB.width; longuest_edgeC.width;parallel_phi.width;parallel_psi.width; parallel_theta.width</text> 
     20        [ORIENTATION_PARAMS]= <text>parallel_phi;parallel_psi; parallel_theta; parallel_phi.width;parallel_psi.width; parallel_theta.width</text> 
    2121 
    2222 
    … …  
    3838    //  [DEFAULT]=contrast=5.3e-6 [1/A²] 
    3939    double contrast; 
    40         /// Incoherent Background [1/cm]  
     40        /// Incoherent Background [1/cm] 
    4141        //  [DEFAULT]=background=0.0 [1/cm] 
    4242        double background; 
    4343    /// Orientation of the parallelepiped axis w/respect incoming beam [rad] 
    44     //  [DEFAULT]=parallel_theta=1.0 [rad] 
     44    //  [DEFAULT]=parallel_theta=0.0 [rad] 
    4545    double parallel_theta; 
    46     /// Orientation of the parallelepiped in the plane of the detector [rad] 
    47     //  [DEFAULT]=parallel_phi=1.0 [rad] 
     46    /// Orientation of the longitudinal axis of the parallelepiped in the plane of the detector [rad] 
     47    //  [DEFAULT]=parallel_phi=0.0 [rad] 
    4848    double parallel_phi; 
     49    /// Orientation of the cross-sectional minor axis of the parallelepiped in the plane of the detector [rad] 
     50    //  [DEFAULT]=parallel_psi=0.0 [rad] 
     51    double parallel_psi; 
    4952 
    5053 
    … …  
    6063double parallelepiped_analytical_2DXY(ParallelepipedParameters *pars, double qx, double qy); 
    6164double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y); 
    62  
    6365#endif 

  • sansmodels/src/sans/models/c_extensions/stacked_disks.h

    r5068697 r975ec8e  
    77/** Structure definition for stacked disks parameters 
    88 * [PYTHONCLASS] = StackedDisksModel 
    9  * [DISP_PARAMS] = length, radius, axis_theta, axis_phi 
     9 * [DISP_PARAMS] = core_thick, layer_thick, radius, axis_theta, axis_phi 
    1010   [DESCRIPTION] = <text> 
    1111                </text> 
    12         [FIXED]= <text>length.width; radius.width; axis_theta.width; axis_phi.width</text> 
     12        [FIXED]= <text>core_thick.width;layer_thick.width; radius.width; axis_theta.width; axis_phi.width</text> 
    1313        [ORIENTATION_PARAMS]= <text>axis_phi; axis_theta; axis_phi.width; axis_theta.width</text> 
    1414 
    … …  
    2222    //  [DEFAULT]=radius=3000 [A] 
    2323    double radius; 
    24     /// Length of the staked disk [A] 
    25     //  [DEFAULT]=length=10.0 [A] 
    26     double length; 
     24    /// Thickness of the core disk [A] 
     25    //  [DEFAULT]=core_thick=10.0 [A] 
     26    double core_thick; 
    2727        /// Thickness of the staked disk [A] 
    28     //  [DEFAULT]=thickness=15.0 [A] 
    29     double thickness; 
     28    //  [DEFAULT]=layer_thick=15.0 [A] 
     29    double layer_thick; 
    3030        /// Core scattering length density[1/A²] 
    3131    //  [DEFAULT]=core_sld=4e-6 [1/A²] 
    … …  
    3737    //  [DEFAULT]=solvent_sld=5.0e-6 [1/A²] 
    3838    double solvent_sld; 
    39     /// number of layers 
    40     //  [DEFAULT]=nlayers=1 
    41         double nlayers; 
    42     /// GSD of disks spacing 
    43     //  [DEFAULT]=spacing=0 
    44     double spacing; 
    45         /// Incoherent Background [1/cm]  
     39    /// number of stacking 
     40    //  [DEFAULT]=n_stacking=1 
     41        double n_stacking; 
     42    /// GSD of disks sigma_d 
     43    //  [DEFAULT]=sigma_d=0 
     44    double sigma_d; 
     45        /// Incoherent Background [1/cm] 
    4646        //  [DEFAULT]=background=0.001 [1/cm] 
    4747        double background; 
    4848    /// Orientation of the staked disk axis w/respect incoming beam [rad] 
    49     //  [DEFAULT]=axis_theta=1.0 [rad] 
     49    //  [DEFAULT]=axis_theta=0.0 [rad] 
    5050    double axis_theta; 
    5151    /// Orientation of the  staked disk in the plane of the detector [rad] 
    52     //  [DEFAULT]=axis_phi=1.0 [rad] 
     52    //  [DEFAULT]=axis_phi=0.0 [rad] 
    5353    double axis_phi; 
    5454 

  • sansmodels/src/sans/models/c_extensions/triaxial_ellipsoid.c

    r34c3020 r975ec8e  
    1919double triaxial_ellipsoid_analytical_1D(TriaxialEllipsoidParameters *pars, double q) { 
    2020        double dp[6]; 
    21          
     21 
    2222        // Fill paramater array 
    2323        dp[0] = pars->scale; 
    … …  
    2727        dp[4] = pars->contrast; 
    2828        dp[5] = pars->background; 
    29          
     29 
    3030        // Call library function to evaluate model 
    31         return TriaxialEllipsoid(dp, q);         
     31        return TriaxialEllipsoid(dp, q); 
    3232} 
     33 
     34double triaxial_ellipsoid_kernel(TriaxialEllipsoidParameters *pars, double q, double alpha, double nu) { 
     35        double t,a,b,c; 
     36        double kernel; 
     37        double pi = acos(-1.0); 
     38 
     39        a = pars->semi_axisA ; 
     40        b = pars->semi_axisB ; 
     41        c = pars->semi_axisC ; 
     42 
     43        t = q * sqrt(a*a*cos(nu)*cos(nu)+b*b*sin(nu)*sin(nu)*sin(alpha)*sin(alpha)+c*c*cos(alpha)*cos(alpha)); 
     44        if (t==0){ 
     45                kernel  = 1.0; 
     46        }else{ 
     47                kernel  = 3*(sin(t)-t*cos(t))/(t*t*t); 
     48        } 
     49        return kernel*kernel; 
     50} 
     51 
    3352 
    3453/** 
    … …  
    4261        q = sqrt(qx*qx+qy*qy); 
    4362    return triaxial_ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q); 
    44 }  
     63} 
    4564 
    4665 
    … …  
    5473double triaxial_ellipsoid_analytical_2D(TriaxialEllipsoidParameters *pars, double q, double phi) { 
    5574    return triaxial_ellipsoid_analytical_2D_scaled(pars, q, cos(phi), sin(phi)); 
    56 }  
    57          
     75} 
     76 
    5877/** 
    5978 * Function to evaluate 2D scattering function 
    … …  
    6584 */ 
    6685double triaxial_ellipsoid_analytical_2D_scaled(TriaxialEllipsoidParameters *pars, double q, double q_x, double q_y) { 
    67         double cyl_x, cyl_y, cyl_z; 
     86        double cyl_x, cyl_y, cyl_z, ell_x, ell_y; 
    6887        double q_z; 
    69         double dx,  dy; 
     88        double cos_nu,nu; 
    7089        double alpha, vol, cos_val; 
    7190        double answer; 
    … …  
    7594    cyl_y = sin(pars->axis_theta) * sin(pars->axis_phi); 
    7695    cyl_z = cos(pars->axis_theta); 
    77       
     96 
    7897    // q vector 
    7998    q_z = 0; 
    80        
    81         dx = 1.0; 
    82         dy = 1.0; 
     99 
     100        //dx = 1.0; 
     101        //dy = 1.0; 
    83102    // Compute the angle btw vector q and the 
    84103    // axis of the cylinder 
    85104    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; 
    86      
     105 
    87106    // The following test should always pass 
    88107    if (fabs(cos_val)>1.0) { 
    … …  
    90109        return 0; 
    91110    } 
    92      
     111 
    93112    // Note: cos(alpha) = 0 and 1 will get an 
    94113    // undefined value from CylKernel 
    95114        alpha = acos( cos_val ); 
    96          
     115 
     116    //ellipse orientation: 
     117        // the elliptical corss section was transformed and projected 
     118        // into the detector plane already through sin(alpha)and furthermore psi remains as same 
     119        // on the detector plane. 
     120        // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt 
     121        // the wave vector q. 
     122 
     123        //x- y- component on the detector plane. 
     124    ell_x =  cos(pars->axis_psi); 
     125    ell_y =  sin(pars->axis_psi); 
     126 
     127    // calculate the axis of the ellipse wrt q-coord. 
     128    cos_nu = ell_x*q_x + ell_y*q_y; 
     129    nu = acos(cos_nu); 
     130 
    97131        // Call the IGOR library function to get the kernel 
    98         answer = TriaxialKernel(q,pars->semi_axisA, pars->semi_axisB, pars->semi_axisC, dx, dy); 
    99          
     132        answer = triaxial_ellipsoid_kernel(pars, q, alpha, nu); 
     133 
    100134        // Multiply by contrast^2 
    101135        answer *= pars->contrast*pars->contrast; 
    102          
     136 
    103137        //normalize by cylinder volume 
    104138        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    105139    vol = 4/3 * pi * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC; 
    106140        answer *= vol; 
    107          
     141 
    108142        //convert to [cm-1] 
    109143        answer *= 1.0e8; 
    110          
     144 
    111145        //Scale 
    112146        answer *= pars->scale; 
    113          
     147 
    114148        // add in the background 
    115149        answer += pars->background; 
    116          
     150 
    117151        return answer; 
    118152} 
    119      

  • sansmodels/src/sans/models/c_extensions/triaxial_ellipsoid.h

    r5068697 r975ec8e  
    33/** Structure definition for cylinder parameters 
    44 * [PYTHONCLASS] = TriaxialEllipsoidModel 
    5  * [DISP_PARAMS] = axis_theta, axis_phi 
     5 * [DISP_PARAMS] = axis_theta, axis_phi, axis_psi 
    66   [DESCRIPTION] = <text> Note: 
    77                                        Constraints must be applied during fitting to ensure that the inequality a<b<c is not 
    88                                                  violated. The calculation will not report an error, but the results will not be correct. 
    99                                   </text> 
    10         [FIXED]= <text>axis_phi.width; axis_theta.width; semi_axisA.width; semi_axisB.width; semi_axisC.width </text> 
    11         [ORIENTATION_PARAMS]= <text>axis_phi; axis_theta; axis_phi.width; axis_theta.width</text> 
     10        [FIXED]= <text>axis_psi.width; axis_phi.width; axis_theta.width; semi_axisA.width; semi_axisB.width; semi_axisC.width </text> 
     11        [ORIENTATION_PARAMS]= <text>axis_psi; axis_phi; axis_theta; axis_psi.width; axis_phi.width; axis_theta.width</text> 
    1212 
    1313 
    … …  
    3333        double background; 
    3434    /// Orientation of the triaxial_ellipsoid axis w/respect incoming beam [rad] 
    35     //  [DEFAULT]=axis_theta=1.0 [rad] 
     35    //  [DEFAULT]=axis_theta=0.0 [rad] 
    3636    double axis_theta; 
    3737    /// Orientation of the triaxial_ellipsoid in the plane of the detector [rad] 
    38     //  [DEFAULT]=axis_phi=1.0 [rad] 
     38    //  [DEFAULT]=axis_phi=0.0 [rad] 
    3939    double axis_phi; 
     40    /// Orientation of the cross section of the triaxial_ellipsoid in the plane of the detector [rad] 
     41    //  [DEFAULT]=axis_psi=0.0 [rad] 
     42    double axis_psi; 
    4043 
    4144} TriaxialEllipsoidParameters; 
Note: See TracChangeset for help on using the changeset viewer.