source: sasview/sansmodels/src/sans/models/c_extensions/elliptical_cylinder.c @ d6da3b1

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since d6da3b1 was 4628e31, checked in by Jae Cho <jhjcho@…>, 14 years ago

changed the unit of angles into degrees

  • Property mode set to 100644
File size: 4.5 KB
Line 
1/**
2 * Scattering model for a cylinder with elliptical cross-section
3 */
4
5#include "elliptical_cylinder.h"
6#include <math.h>
7#include "libCylinder.h"
8#include <stdio.h>
9#include <stdlib.h>
10
11
12/**
13 * Function to evaluate 1D scattering function
14 * @param pars: parameters of the cylinder
15 * @param q: q-value
16 * @return: function value
17 */
18double elliptical_cylinder_analytical_1D(EllipticalCylinderParameters *pars, double q) {
19        double dp[7];
20
21        // Fill paramater array
22        dp[0] = pars->scale;
23        dp[1] = pars->r_minor;
24        dp[2] = pars->r_ratio;
25        dp[3] = pars->length;
26        dp[4] = pars->sldCyl;
27        dp[5] = pars->sldSolv;
28        dp[6] = pars->background;
29
30        // Call library function to evaluate model
31        return EllipCyl20(dp, q);
32}
33
34double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double nu) {
35        double qr;
36        double qL;
37        double Be,Si;
38        double r_major;
39        double kernel;
40
41        r_major = pars->r_ratio * pars->r_minor;
42
43        qr = q*sin(alpha)*sqrt( r_major*r_major*sin(nu)*sin(nu) + pars->r_minor*pars->r_minor*cos(nu)*cos(nu) );
44        qL = q*pars->length*cos(alpha)/2.0;
45
46        if (qr==0){
47                Be = 0.5;
48        }else{
49                Be = NR_BessJ1(qr)/qr;
50        }
51        if (qL==0){
52                Si = 1.0;
53        }else{
54                Si = sin(qL)/qL;
55        }
56
57
58        kernel = 2.0*Be * Si;
59        return kernel*kernel;
60}
61
62/**
63 * Function to evaluate 2D scattering function
64 * @param pars: parameters of the cylinder
65 * @param q: q-value
66 * @return: function value
67 */
68double elliptical_cylinder_analytical_2DXY(EllipticalCylinderParameters *pars, double qx, double qy) {
69        double q;
70        q = sqrt(qx*qx+qy*qy);
71    return elliptical_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
72}
73
74/**
75 * Function to evaluate 2D scattering function
76 * @param pars: parameters of the cylinder
77 * @param q: q-value
78 * @param theta: angle theta = angle wrt z axis
79 * @param phi: angle phi = angle around y axis (starting from the x+-direction as phi = 0)
80 * @return: function value
81 */
82double elliptical_cylinder_analytical_2D(EllipticalCylinderParameters *pars, double q, double phi) {
83    return elliptical_cylinder_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
84}
85
86/**
87 * Function to evaluate 2D scattering function
88 * @param pars: parameters of the cylinder
89 * @param q: q-value
90 * @param q_x: q_x / q
91 * @param q_y: q_y / q
92 * @return: function value
93 */
94double elliptical_cylinder_analytical_2D_scaled(EllipticalCylinderParameters *pars, double q, double q_x, double q_y) {
95        double cyl_x, cyl_y, cyl_z;
96        double ell_x, ell_y;
97        double q_z;
98        double alpha, vol, cos_val;
99        double nu, cos_nu;
100        double answer;
101        //convert angle degree to radian
102        double pi = 4.0*atan(1.0);
103        double theta = pars->cyl_theta * pi/180.0;
104        double phi = pars->cyl_phi * pi/180.0;
105        double psi = pars->cyl_psi * pi/180.0;
106
107    //Cylinder orientation
108    cyl_x = sin(theta) * cos(phi);
109    cyl_y = sin(theta) * sin(phi);
110    cyl_z = cos(theta);
111
112    // q vector
113    q_z = 0;
114
115    // Compute the angle btw vector q and the
116    // axis of the cylinder
117    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
118
119    // The following test should always pass
120    if (fabs(cos_val)>1.0) {
121        printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
122        return 0;
123    }
124
125    // Note: cos(alpha) = 0 and 1 will get an
126    // undefined value from CylKernel
127        alpha = acos( cos_val );
128
129    //ellipse orientation:
130        // the elliptical corss section was transformed and projected
131        // into the detector plane already through sin(alpha)and furthermore psi remains as same
132        // on the detector plane.
133        // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt
134        // the wave vector q.
135
136        //x- y- component on the detector plane.
137    ell_x =  cos(psi);
138    ell_y =  sin(psi);
139
140    // calculate the axis of the ellipse wrt q-coord.
141    cos_nu = ell_x*q_x + ell_y*q_y;
142    nu = acos(cos_nu);
143
144    // The following test should always pass
145    if (fabs(cos_nu)>1.0) {
146        printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
147        return 0;
148    }
149
150        answer = elliptical_cylinder_kernel(pars, q, alpha,nu);
151
152        // Multiply by contrast^2
153        answer *= (pars->sldCyl - pars->sldSolv) * (pars->sldCyl - pars->sldSolv);
154
155        //normalize by cylinder volume
156        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
157    vol = acos(-1.0) * pars->r_minor * pars->r_minor * pars->r_ratio * pars->length;
158        answer *= vol;
159
160        //convert to [cm-1]
161        answer *= 1.0e8;
162
163        //Scale
164        answer *= pars->scale;
165
166        // add in the background
167        answer += pars->background;
168
169        return answer;
170}
171
Note: See TracBrowser for help on using the repository browser.