source: sasview/sansmodels/src/sans/models/c_extensions/parallelepiped.c @ 01de557

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

Re #5 fixing samsmodels compilation on MSVC

  • Property mode set to 100644
File size: 4.7 KB
Line 
1/**
2 * Scattering model for a parallelepiped
3 * TODO: Add 2D analysis
4 */
5
6#include "parallelepiped.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 parallelepiped
16 * @param q: q-value
17 * @return: function value
18 */
19double parallelepiped_analytical_1D(ParallelepipedParameters *pars, double q) {
20        double dp[7];
21
22        // Fill paramater array
23        dp[0] = pars->scale;
24        dp[1] = pars->short_a;
25        dp[2] = pars->short_b;
26        dp[3] = pars->long_c;
27        dp[4] = pars->sldPipe;
28        dp[5] = pars->sldSolv;
29        dp[6] = pars->background;
30
31        // Call library function to evaluate model
32        return Parallelepiped(dp, q);
33}
34
35
36double pkernel(double a, double b,double c, double ala, double alb, double alc){
37    // mu passed in is really mu*sqrt(1-sig^2)
38    double argA,argB,argC,tmp1,tmp2,tmp3;                       //local variables
39
40    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
41    argA = a*ala/2.0;
42    argB = b*alb/2.0;
43    argC = c*alc/2.0;
44    if(argA==0.0) {
45                tmp1 = 1.0;
46        } else {
47                tmp1 = sin(argA)*sin(argA)/argA/argA;
48    }
49    if (argB==0.0) {
50                tmp2 = 1.0;
51        } else {
52                tmp2 = sin(argB)*sin(argB)/argB/argB;
53    }
54
55    if (argC==0.0) {
56                tmp3 = 1.0;
57        } else {
58                tmp3 = sin(argC)*sin(argC)/argC/argC;
59    }
60
61    return (tmp1*tmp2*tmp3);
62
63}//Function pkernel()
64
65
66
67
68/**
69 * Function to evaluate 2D scattering function
70 * @param pars: parameters of the parallelepiped
71 * @param q: q-value
72 * @return: function value
73 */
74double parallelepiped_analytical_2DXY(ParallelepipedParameters *pars, double qx, double qy) {
75        double q;
76        q = sqrt(qx*qx+qy*qy);
77    return parallelepiped_analytical_2D_scaled(pars, q, qx/q, qy/q);
78}
79
80
81/**
82 * Function to evaluate 2D scattering function
83 * @param pars: parameters of the Parallelepiped
84 * @param q: q-value
85 * @param phi: angle phi
86 * @return: function value
87 */
88double parallelepiped_analytical_2D(ParallelepipedParameters *pars, double q, double phi) {
89    return parallelepiped_analytical_2D_scaled(pars, q, cos(phi), sin(phi));
90}
91
92/**
93 * Function to evaluate 2D scattering function
94 * @param pars: parameters of the parallelepiped
95 * @param q: q-value
96 * @param q_x: q_x / q
97 * @param q_y: q_y / q
98 * @return: function value
99 */
100double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y) {
101        double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y, parallel_z;
102        double q_z;
103        double alpha, vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC;
104
105        double answer;
106        double pi = 4.0*atan(1.0);
107
108  //convert angle degree to radian
109  double theta = pars->parallel_theta * pi/180.0;
110  double phi = pars->parallel_phi * pi/180.0;
111  double psi = pars->parallel_psi * pi/180.0;
112
113        edgeA = pars->short_a;
114        edgeB = pars->short_b;
115        edgeC = pars->long_c;
116
117
118    // parallelepiped c axis orientation
119    cparallel_x = sin(theta) * cos(phi);
120    cparallel_y = sin(theta) * sin(phi);
121    cparallel_z = cos(theta);
122
123    // q vector
124    q_z = 0.0;
125
126    // Compute the angle btw vector q and the
127    // axis of the parallelepiped
128    cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z;
129    alpha = acos(cos_val_c);
130
131    // parallelepiped a axis orientation
132    parallel_x = sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi);
133    parallel_y = cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi);
134
135    cos_val_a = parallel_x*q_x + parallel_y*q_y;
136
137
138
139    // parallelepiped b axis orientation
140    bparallel_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi);
141    bparallel_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi);
142    // axis of the parallelepiped
143    cos_val_b = sin(acos(cos_val_a)) ;
144
145
146
147    // The following test should always pass
148    if (fabs(cos_val_c)>1.0) {
149        printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
150        return 0;
151    }
152
153        // Call the IGOR library function to get the kernel
154        answer = pkernel( q*edgeA, q*edgeB, q*edgeC, sin(alpha)*cos_val_a,sin(alpha)*cos_val_b,cos_val_c);
155
156        // Multiply by contrast^2
157        answer *= (pars->sldPipe - pars->sldSolv) * (pars->sldPipe - pars->sldSolv);
158
159        //normalize by cylinder volume
160        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel
161    vol = edgeA* edgeB * edgeC;
162        answer *= vol;
163
164        //convert to [cm-1]
165        answer *= 1.0e8;
166
167        //Scale
168        answer *= pars->scale;
169
170        // add in the background
171        answer += pars->background;
172
173        return answer;
174}
175
Note: See TracBrowser for help on using the repository browser.