source: sasview/sansmodels/src/sans/models/c_extensions/triaxial_ellipsoid.c @ a24f530

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 a24f530 was a24f530, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

Re #4 A few more warnings

  • Property mode set to 100644
File size: 4.3 KB
Line 
1/**
2 * Scattering model for a cylinder
3 * @author: Mathieu Doucet / UTK
4 */
5
6#include "triaxial_ellipsoid.h"
7#include <math.h>
8#include "libCylinder.h"
9#include <stdio.h>
10#include <stdlib.h>
11
12
13/**
14 * Function to evaluate 1D scattering function
15 * @param pars: parameters of the triaxial ellipsoid
16 * @param q: q-value
17 * @return: function value
18 */
19double triaxial_ellipsoid_analytical_1D(TriaxialEllipsoidParameters *pars, double q) {
20        double dp[7];
21
22        // Fill paramater array
23        dp[0] = pars->scale;
24        dp[1] = pars->semi_axisA;
25        dp[2] = pars->semi_axisB;
26        dp[3] = pars->semi_axisC;
27        dp[4] = pars->sldEll;
28        dp[5] = pars->sldSolv;
29        dp[6] = pars->background;
30
31        // Call library function to evaluate model
32        return TriaxialEllipsoid(dp, q);
33}
34
35double triaxial_ellipsoid_kernel(TriaxialEllipsoidParameters *pars, double q, double alpha, double nu) {
36        double t,a,b,c;
37        double kernel;
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.0){
45                kernel  = 1.0;
46        }else{
47                kernel  = 3.0*(sin(t)-t*cos(t))/(t*t*t);
48        }
49        return kernel*kernel;
50}
51
52
53/**
54 * Function to evaluate 2D scattering function
55 * @param pars: parameters of the triaxial ellipsoid
56 * @param q: q-value
57 * @return: function value
58 */
59double triaxial_ellipsoid_analytical_2DXY(TriaxialEllipsoidParameters *pars, double qx, double qy) {
60        double q;
61        q = sqrt(qx*qx+qy*qy);
62    return triaxial_ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q);
63}
64
65
66/**
67 * Function to evaluate 2D scattering function
68 * @param pars: parameters of the triaxial ellipsoid
69 * @param q: q-value
70 * @param phi: angle phi
71 * @return: function value
72 */
73double triaxial_ellipsoid_analytical_2D(TriaxialEllipsoidParameters *pars, double q, double phi) {
74    return triaxial_ellipsoid_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
75}
76
77/**
78 * Function to evaluate 2D scattering function
79 * @param pars: parameters of the triaxial ellipsoid
80 * @param q: q-value
81 * @param q_x: q_x / q
82 * @param q_y: q_y / q
83 * @return: function value
84 */
85double triaxial_ellipsoid_analytical_2D_scaled(TriaxialEllipsoidParameters *pars, double q, double q_x, double q_y) {
86        double cyl_x, cyl_y, cyl_z, ell_x, ell_y;
87        double q_z;
88        double cos_nu,nu;
89        double alpha, vol, cos_val;
90        double answer;
91    double pi = 4.0*atan(1.0);
92
93        //convert angle degree to radian
94        double theta = pars->axis_theta * pi/180.0;
95        double phi = pars->axis_phi * pi/180.0;
96        double psi = pars->axis_psi * pi/180.0;
97
98    // Cylinder orientation
99    cyl_x = sin(theta) * cos(phi);
100    cyl_y = sin(theta) * sin(phi);
101    cyl_z = cos(theta);
102
103    // q vector
104    q_z = 0.0;
105
106        //dx = 1.0;
107        //dy = 1.0;
108    // Compute the angle btw vector q and the
109    // axis of the cylinder
110    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
111
112    // The following test should always pass
113    if (fabs(cos_val)>1.0) {
114        printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
115        return 0;
116    }
117
118    // Note: cos(alpha) = 0 and 1 will get an
119    // undefined value from CylKernel
120        alpha = acos( cos_val );
121
122    //ellipse orientation:
123        // the elliptical corss section was transformed and projected
124        // into the detector plane already through sin(alpha)and furthermore psi remains as same
125        // on the detector plane.
126        // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt
127        // the wave vector q.
128
129        //x- y- component on the detector plane.
130    ell_x =  cos(psi);
131    ell_y =  sin(psi);
132
133    // calculate the axis of the ellipse wrt q-coord.
134    cos_nu = ell_x*q_x + ell_y*q_y;
135    nu = acos(cos_nu);
136
137        // Call the IGOR library function to get the kernel
138        answer = triaxial_ellipsoid_kernel(pars, q, alpha, nu);
139
140        // Multiply by contrast^2
141        answer *= (pars->sldEll- pars->sldSolv)*(pars->sldEll- pars->sldSolv);
142
143        //normalize by cylinder volume
144        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
145    vol = 4.0* pi/3.0  * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC;
146        answer *= vol;
147        //convert to [cm-1]
148        answer *= 1.0e8;
149        //Scale
150        answer *= pars->scale;
151
152        // add in the background
153        answer += pars->background;
154
155        return answer;
156}
Note: See TracBrowser for help on using the repository browser.