[8a48713] | 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 |
---|
[5068697] | 15 | * @param pars: parameters of the parallelepiped |
---|
[8a48713] | 16 | * @param q: q-value |
---|
| 17 | * @return: function value |
---|
| 18 | */ |
---|
| 19 | double parallelepiped_analytical_1D(ParallelepipedParameters *pars, double q) { |
---|
[f10063e] | 20 | double dp[7]; |
---|
[975ec8e] | 21 | |
---|
[8a48713] | 22 | // Fill paramater array |
---|
| 23 | dp[0] = pars->scale; |
---|
[2cb89e7] | 24 | dp[1] = pars->short_a; |
---|
[8e36cdd] | 25 | dp[2] = pars->short_b; |
---|
| 26 | dp[3] = pars->long_c; |
---|
[f10063e] | 27 | dp[4] = pars->sldPipe; |
---|
| 28 | dp[5] = pars->sldSolv; |
---|
| 29 | dp[6] = pars->background; |
---|
[8a48713] | 30 | |
---|
| 31 | // Call library function to evaluate model |
---|
[975ec8e] | 32 | return Parallelepiped(dp, q); |
---|
[8a48713] | 33 | } |
---|
[975ec8e] | 34 | |
---|
| 35 | |
---|
| 36 | double 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 |
---|
[a1ecfc3] | 41 | argA = a*ala/2.0; |
---|
| 42 | argB = b*alb/2.0; |
---|
| 43 | argC = c*alc/2.0; |
---|
[975ec8e] | 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 | |
---|
[8a48713] | 68 | /** |
---|
| 69 | * Function to evaluate 2D scattering function |
---|
[5068697] | 70 | * @param pars: parameters of the parallelepiped |
---|
[8a48713] | 71 | * @param q: q-value |
---|
| 72 | * @return: function value |
---|
| 73 | */ |
---|
| 74 | double 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); |
---|
[975ec8e] | 78 | } |
---|
[8a48713] | 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 | */ |
---|
| 88 | double parallelepiped_analytical_2D(ParallelepipedParameters *pars, double q, double phi) { |
---|
| 89 | return parallelepiped_analytical_2D_scaled(pars, q, cos(phi), sin(phi)); |
---|
[975ec8e] | 90 | } |
---|
| 91 | |
---|
[8a48713] | 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 | */ |
---|
| 100 | double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y) { |
---|
[74c9a30] | 101 | double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y; |
---|
[8a48713] | 102 | double q_z; |
---|
[975ec8e] | 103 | double alpha, vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC; |
---|
| 104 | |
---|
[8a48713] | 105 | double answer; |
---|
[975ec8e] | 106 | double pi = 4.0*atan(1.0); |
---|
| 107 | |
---|
[890ac7f1] | 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 | |
---|
[2cb89e7] | 113 | edgeA = pars->short_a; |
---|
[8e36cdd] | 114 | edgeB = pars->short_b; |
---|
| 115 | edgeC = pars->long_c; |
---|
[975ec8e] | 116 | |
---|
| 117 | |
---|
| 118 | // parallelepiped c axis orientation |
---|
[4628e31] | 119 | cparallel_x = sin(theta) * cos(phi); |
---|
| 120 | cparallel_y = sin(theta) * sin(phi); |
---|
| 121 | cparallel_z = cos(theta); |
---|
[975ec8e] | 122 | |
---|
[8a48713] | 123 | // q vector |
---|
[a1ecfc3] | 124 | q_z = 0.0; |
---|
[975ec8e] | 125 | |
---|
[8a48713] | 126 | // Compute the angle btw vector q and the |
---|
| 127 | // axis of the parallelepiped |
---|
[8e36cdd] | 128 | cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z; |
---|
[975ec8e] | 129 | alpha = acos(cos_val_c); |
---|
| 130 | |
---|
| 131 | // parallelepiped a axis orientation |
---|
[4628e31] | 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); |
---|
[3c102d4] | 134 | |
---|
[8e36cdd] | 135 | cos_val_a = parallel_x*q_x + parallel_y*q_y; |
---|
[975ec8e] | 136 | |
---|
| 137 | |
---|
| 138 | |
---|
| 139 | // parallelepiped b axis orientation |
---|
[4628e31] | 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); |
---|
[975ec8e] | 142 | // axis of the parallelepiped |
---|
[8e36cdd] | 143 | cos_val_b = sin(acos(cos_val_a)) ; |
---|
[975ec8e] | 144 | |
---|
| 145 | |
---|
| 146 | |
---|
[8a48713] | 147 | // The following test should always pass |
---|
[975ec8e] | 148 | if (fabs(cos_val_c)>1.0) { |
---|
[8a48713] | 149 | printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); |
---|
| 150 | return 0; |
---|
| 151 | } |
---|
[975ec8e] | 152 | |
---|
[8a48713] | 153 | // Call the IGOR library function to get the kernel |
---|
[8e36cdd] | 154 | answer = pkernel( q*edgeA, q*edgeB, q*edgeC, sin(alpha)*cos_val_a,sin(alpha)*cos_val_b,cos_val_c); |
---|
[975ec8e] | 155 | |
---|
[8a48713] | 156 | // Multiply by contrast^2 |
---|
[f10063e] | 157 | answer *= (pars->sldPipe - pars->sldSolv) * (pars->sldPipe - pars->sldSolv); |
---|
[975ec8e] | 158 | |
---|
[8a48713] | 159 | //normalize by cylinder volume |
---|
| 160 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel |
---|
[2cb89e7] | 161 | vol = edgeA* edgeB * edgeC; |
---|
[8a48713] | 162 | answer *= vol; |
---|
[975ec8e] | 163 | |
---|
[8a48713] | 164 | //convert to [cm-1] |
---|
| 165 | answer *= 1.0e8; |
---|
[975ec8e] | 166 | |
---|
[8a48713] | 167 | //Scale |
---|
| 168 | answer *= pars->scale; |
---|
[975ec8e] | 169 | |
---|
[8a48713] | 170 | // add in the background |
---|
| 171 | answer += pars->background; |
---|
[975ec8e] | 172 | |
---|
[8a48713] | 173 | return answer; |
---|
| 174 | } |
---|
| 175 | |
---|