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

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 0b12abb5 was 13eb1c4, checked in by Jae Cho <jhjcho@…>, 15 years ago

updated sld parameters according to the changes in NIST lib: missed it in the last commit

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