Changeset 3fe701a in sasview for sansmodels


Ignore:
Timestamp:
Jun 11, 2009 3:05:02 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:
0d88a09
Parents:
240b9966
Message:

think 2d bug fixed

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    rae3ce4e r3fe701a  
    1818double elliptical_cylinder_analytical_1D(EllipticalCylinderParameters *pars, double q) { 
    1919        double dp[6]; 
    20          
     20 
    2121        // Fill paramater array 
    2222        dp[0] = pars->scale; 
    … …  
    2626        dp[4] = pars->contrast; 
    2727        dp[5] = pars->background; 
    28          
     28 
    2929        // Call library function to evaluate model 
    30         return EllipCyl20(dp, q);        
     30        return EllipCyl20(dp, q); 
    3131} 
    3232 
    33 double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double psi) { 
     33double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double psi, double nu) { 
    3434        double qr; 
    3535        double qL; 
    3636        double r_major; 
    3737        double kernel; 
    38          
     38 
    3939        r_major = pars->r_ratio * pars->r_minor; 
    4040 
    41         qr = q*sin(alpha)*sqrt( r_major*r_major*sin(psi)*sin(psi) + pars->r_minor*pars->r_minor*cos(psi)*cos(psi) ); 
     41        qr = q*sin(alpha)*sqrt( r_major*r_major*sin(nu)*sin(nu) + pars->r_minor*pars->r_minor*cos(nu)*cos(nu) ); 
    4242        qL = q*pars->length*cos(alpha)/2.0; 
    43          
     43 
    4444        kernel = 2.0*NR_BessJ1(qr)/qr * sin(qL)/qL; 
    4545        return kernel*kernel; 
    … …  
    5656        q = sqrt(qx*qx+qy*qy); 
    5757    return elliptical_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q); 
    58 }  
     58} 
    5959 
    6060/** 
    … …  
    6262 * @param pars: parameters of the cylinder 
    6363 * @param q: q-value 
    64  * @param phi: angle phi 
     64 * @param theta: angle theta = angle wrt z axis 
     65 * @param phi: angle phi = angle around y axis (starting from the x+-direction as phi = 0) 
    6566 * @return: function value 
    6667 */ 
    6768double elliptical_cylinder_analytical_2D(EllipticalCylinderParameters *pars, double q, double phi) { 
    6869    return elliptical_cylinder_analytical_2D_scaled(pars, q, cos(phi), sin(phi)); 
    69 }  
     70} 
    7071 
    7172/** 
    … …  
    7980double elliptical_cylinder_analytical_2D_scaled(EllipticalCylinderParameters *pars, double q, double q_x, double q_y) { 
    8081        double cyl_x, cyl_y, cyl_z; 
     82        double ell_x, ell_y; 
    8183        double q_z; 
    8284        double alpha, vol, cos_val; 
     85        double nu, cos_nu; 
    8386        double answer; 
    84          
    85     // Cylinder orientation 
     87 
     88    //Cylinder orientation 
    8689    cyl_x = sin(pars->cyl_theta) * cos(pars->cyl_phi); 
    8790    cyl_y = sin(pars->cyl_theta) * sin(pars->cyl_phi); 
    8891    cyl_z = cos(pars->cyl_theta); 
    89       
     92 
    9093    // q vector 
    9194    q_z = 0; 
    92          
     95 
    9396    // Compute the angle btw vector q and the 
    9497    // axis of the cylinder 
    9598    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; 
    96      
     99 
    97100    // The following test should always pass 
    98101    if (fabs(cos_val)>1.0) { 
    … …  
    100103        return 0; 
    101104    } 
    102      
     105 
    103106    // Note: cos(alpha) = 0 and 1 will get an 
    104107    // undefined value from CylKernel 
    105108        alpha = acos( cos_val ); 
    106          
    107         answer = elliptical_cylinder_kernel(pars, q, alpha, pars->cyl_psi); 
    108          
     109 
     110    //ellipse orientation: 
     111        // the elliptical corss section was transformed and projected 
     112        // into the detector plane already through sin(alpha)and furthermore psi remains as same 
     113        // on the detector plane. 
     114        // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt 
     115        // the wave vector q. 
     116 
     117        //x- y- component on the detector plane. 
     118    ell_x =  cos(pars->cyl_psi); 
     119    ell_y =  sin(pars->cyl_psi); 
     120 
     121    // calculate the axis of the ellipse wrt q-coord. 
     122    cos_nu = ell_x*q_x + ell_y*q_y; 
     123    nu = acos(cos_nu); 
     124 
     125    // The following test should always pass 
     126    if (fabs(cos_nu)>1.0) { 
     127        printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n"); 
     128        return 0; 
     129    } 
     130 
     131        answer = elliptical_cylinder_kernel(pars, q, alpha, pars->cyl_psi,nu); 
     132 
    109133        // Multiply by contrast^2 
    110134        answer *= pars->contrast*pars->contrast; 
    111          
     135 
    112136        //normalize by cylinder volume 
    113137        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    114138    vol = acos(-1.0) * pars->r_minor * pars->r_minor * pars->r_ratio * pars->length; 
    115139        answer *= vol; 
    116          
     140 
    117141        //convert to [cm-1] 
    118142        answer *= 1.0e8; 
    119          
     143 
    120144        //Scale 
    121145        answer *= pars->scale; 
    122          
     146 
    123147        // add in the background 
    124148        answer += pars->background; 
    125          
     149 
    126150        return answer; 
    127151} 
    128      
     152 
Note: See TracChangeset for help on using the changeset viewer.