Changeset 3fe701a in sasview for sansmodels/src/sans
- Timestamp:
- Jun 11, 2009 3:05:02 PM (15 years ago)
- 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
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/sans/models/c_extensions/elliptical_cylinder.c
rae3ce4e r3fe701a 18 18 double elliptical_cylinder_analytical_1D(EllipticalCylinderParameters *pars, double q) { 19 19 double dp[6]; 20 20 21 21 // Fill paramater array 22 22 dp[0] = pars->scale; … … 26 26 dp[4] = pars->contrast; 27 27 dp[5] = pars->background; 28 28 29 29 // Call library function to evaluate model 30 return EllipCyl20(dp, q); 30 return EllipCyl20(dp, q); 31 31 } 32 32 33 double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double psi ) {33 double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double psi, double nu) { 34 34 double qr; 35 35 double qL; 36 36 double r_major; 37 37 double kernel; 38 38 39 39 r_major = pars->r_ratio * pars->r_minor; 40 40 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) ); 42 42 qL = q*pars->length*cos(alpha)/2.0; 43 43 44 44 kernel = 2.0*NR_BessJ1(qr)/qr * sin(qL)/qL; 45 45 return kernel*kernel; … … 56 56 q = sqrt(qx*qx+qy*qy); 57 57 return elliptical_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q); 58 } 58 } 59 59 60 60 /** … … 62 62 * @param pars: parameters of the cylinder 63 63 * @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) 65 66 * @return: function value 66 67 */ 67 68 double elliptical_cylinder_analytical_2D(EllipticalCylinderParameters *pars, double q, double phi) { 68 69 return elliptical_cylinder_analytical_2D_scaled(pars, q, cos(phi), sin(phi)); 69 } 70 } 70 71 71 72 /** … … 79 80 double elliptical_cylinder_analytical_2D_scaled(EllipticalCylinderParameters *pars, double q, double q_x, double q_y) { 80 81 double cyl_x, cyl_y, cyl_z; 82 double ell_x, ell_y; 81 83 double q_z; 82 84 double alpha, vol, cos_val; 85 double nu, cos_nu; 83 86 double answer; 84 85 // 87 88 //Cylinder orientation 86 89 cyl_x = sin(pars->cyl_theta) * cos(pars->cyl_phi); 87 90 cyl_y = sin(pars->cyl_theta) * sin(pars->cyl_phi); 88 91 cyl_z = cos(pars->cyl_theta); 89 92 90 93 // q vector 91 94 q_z = 0; 92 95 93 96 // Compute the angle btw vector q and the 94 97 // axis of the cylinder 95 98 cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; 96 99 97 100 // The following test should always pass 98 101 if (fabs(cos_val)>1.0) { … … 100 103 return 0; 101 104 } 102 105 103 106 // Note: cos(alpha) = 0 and 1 will get an 104 107 // undefined value from CylKernel 105 108 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 109 133 // Multiply by contrast^2 110 134 answer *= pars->contrast*pars->contrast; 111 135 112 136 //normalize by cylinder volume 113 137 //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 114 138 vol = acos(-1.0) * pars->r_minor * pars->r_minor * pars->r_ratio * pars->length; 115 139 answer *= vol; 116 140 117 141 //convert to [cm-1] 118 142 answer *= 1.0e8; 119 143 120 144 //Scale 121 145 answer *= pars->scale; 122 146 123 147 // add in the background 124 148 answer += pars->background; 125 149 126 150 return answer; 127 151 } 128 152
Note: See TracChangeset
for help on using the changeset viewer.