source: sasview/sansmodels/src/c_extensions/sc.c @ 31af069

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

Reorganizing models in preparation of cpp cleanup

  • Property mode set to 100644
File size: 3.6 KB
Line 
1/*
2 * Scattering model for a SC_ParaCrystal
3 */
4#include "sc.h"
5#include "libSphere.h"
6#include <math.h>
7#include <stdio.h>
8
9
10/**
11 * Function to evaluate 1D scattering function
12 * @param pars: parameters of the SC_ParaCrystal
13 * @param q: q-value
14 * @return: function value
15 */
16double sc_analytical_1D(SCParameters *pars, double q) {
17        double dp[7];
18        double result;
19
20        dp[0] = pars->scale;
21        dp[1] = pars->dnn;
22        dp[2] = pars->d_factor;
23        dp[3] = pars->radius;
24        dp[4] = pars->sldSph;
25        dp[5] = pars->sldSolv;
26        dp[6] = pars->background;
27
28        result = SC_ParaCrystal(dp, q);
29        // This FIXES a singualrity the kernel in libigor.
30        if ( result == INFINITY || result == NAN){
31                result = pars->background;
32        }
33        return result;
34}
35
36/**
37 * Function to evaluate 2D scattering function
38 * @param pars: parameters of the SC_ParaCrystal
39 * @param q: q-value
40 * @return: function value
41 */
42double sc_analytical_2DXY(SCParameters *pars, double qx, double qy){
43        double q;
44        q = sqrt(qx*qx+qy*qy);
45        return sc_analytical_2D_scaled(pars, q, qx/q, qy/q);
46}
47
48double sc_analytical_2D(SCParameters *pars, double q, double phi) {
49        return sc_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
50}
51
52/**
53 * Function to evaluate 2D scattering function
54 * @param pars: parameters of the SCCrystalModel
55 * @param q: q-value
56 * @param q_x: q_x / q
57 * @param q_y: q_y / q
58 * @return: function value
59 */
60double sc_analytical_2D_scaled(SCParameters *pars, double q, double q_x, double q_y) {
61        double a3_x, a3_y, a3_z, a2_x, a2_y, a1_x, a1_y;
62        double q_z;
63        double alpha, cos_val_a3, cos_val_a2, cos_val_a1;
64        double a1_dot_q, a2_dot_q,a3_dot_q;
65        double answer;
66        double Pi = 4.0*atan(1.0);
67        double aa, Da, qDa_2, latticeScale, Zq;
68
69        double dp[5];
70  //convert angle degree to radian
71  double theta = pars->theta * Pi/180.0;
72  double phi = pars->phi * Pi/180.0;
73  double psi = pars->psi * Pi/180.0;
74        dp[0] = 1.0;
75        dp[1] = pars->radius;
76        dp[2] = pars->sldSph;
77        dp[3] = pars->sldSolv;
78        dp[4] = 0.0;
79
80
81        aa = pars->dnn;
82        Da = pars->d_factor*aa;
83        qDa_2 = pow(q*Da,2.0);
84
85        latticeScale = (4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa,3.0);
86        /// Angles here are respect to detector coordinate instead of against q coordinate(PRB 36, 3, 1754)
87    // a3 axis orientation
88    a3_x = sin(theta) * cos(phi);//negative sign here???
89    a3_y = sin(theta) * sin(phi);
90    a3_z = cos(theta);
91
92    // q vector
93    q_z = 0.0;
94
95    // Compute the angle btw vector q and the a3 axis
96    cos_val_a3 = a3_x*q_x + a3_y*q_y + a3_z*q_z;
97    alpha = acos(cos_val_a3);
98    //alpha = acos(cos_val_a3);
99    a3_dot_q = aa*q*cos_val_a3;
100    // a1 axis orientation
101    a1_x = sin(psi);
102    a1_y = cos(psi);
103
104    cos_val_a1 = a1_x*q_x + a1_y*q_y;
105    a1_dot_q = aa*q*cos_val_a1*sin(alpha);
106
107    // a2 axis orientation
108    a2_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi);
109    a2_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi);
110    // a2 axis
111    cos_val_a2 =  sin(acos(cos_val_a1));//a2_x*q_x + a2_y*q_y;
112    a2_dot_q = aa*q*cos_val_a2*sin(alpha);
113
114    // The following test should always pass
115    if (fabs(cos_val_a3)>1.0) {
116        printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
117        return 0;
118    }
119    // Call Zq=Z1*Z2*Z3
120    Zq = (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a1_dot_q)+exp(-qDa_2));
121    Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a2_dot_q)+exp(-qDa_2));
122    Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a3_dot_q)+exp(-qDa_2));
123
124        // Use SphereForm directly from libigor
125        answer = SphereForm(dp,q)*Zq;
126
127        //consider scales
128        answer *= latticeScale * pars->scale;
129
130        // This FIXES a singualrity the kernel in libigor.
131        if ( answer == INFINITY || answer == NAN){
132                answer = 0.0;
133        }
134
135        // add background
136        answer += pars->background;
137
138        return answer;
139}
Note: See TracBrowser for help on using the repository browser.