[ae3ce4e] | 1 | #include "simcylinder.h" |
---|
| 2 | #include <stdio.h> |
---|
| 3 | #include <stdlib.h> |
---|
| 4 | #include <math.h> |
---|
| 5 | #include <time.h> |
---|
| 6 | #include <memory.h> |
---|
| 7 | #include "modelCalculations.h" |
---|
| 8 | #include "libCylinder.h" |
---|
| 9 | |
---|
| 10 | double test_analytical_2D(SimCylinderParameters *pars, double q, double phi) { |
---|
| 11 | double cyl_x, cyl_y, cyl_z; |
---|
| 12 | double q_x, q_y, q_z, lenq; |
---|
| 13 | double theta, alpha, f, vol, sin_val, cos_val; |
---|
| 14 | double answer; |
---|
| 15 | |
---|
| 16 | // Cylinder orientation |
---|
| 17 | cyl_x = sin(pars->theta) * cos(pars->phi); |
---|
| 18 | cyl_y = sin(pars->theta) * sin(pars->phi); |
---|
| 19 | cyl_z = cos(pars->theta); |
---|
| 20 | |
---|
| 21 | // q orientation |
---|
| 22 | |
---|
| 23 | q_x = cos(phi); |
---|
| 24 | q_y = sin(phi); |
---|
| 25 | q_z = 0; |
---|
| 26 | |
---|
| 27 | // Length of q vector |
---|
| 28 | //lenq = sqrt(q_x*q_x + q_y*q_y + q_z*q_z); |
---|
| 29 | |
---|
| 30 | // Normalize unit vector q |
---|
| 31 | //q_x = q_x/lenq; |
---|
| 32 | //q_y = q_y/lenq; |
---|
| 33 | //q_z = q_z/lenq; |
---|
| 34 | |
---|
| 35 | // Compute the angle btw vector q and the |
---|
| 36 | // axis of the cylinder |
---|
| 37 | cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z; |
---|
| 38 | |
---|
| 39 | // The following test should always pass |
---|
| 40 | if (fabs(cos_val)>1.0) { |
---|
| 41 | printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n"); |
---|
| 42 | return 0; |
---|
| 43 | } |
---|
| 44 | |
---|
| 45 | alpha = acos( cos_val ); |
---|
| 46 | |
---|
| 47 | // Call the IGOR library function to get the kernel |
---|
| 48 | answer = acos(-1.0)/2.0 * CylKernel(q, pars->radius, pars->length/2.0, alpha)/sin(alpha); |
---|
| 49 | |
---|
| 50 | //normalize by cylinder volume |
---|
| 51 | //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 52 | vol = acos(-1.0) * pars->radius * pars->radius * pars->length; |
---|
| 53 | answer *= vol; |
---|
| 54 | |
---|
| 55 | //Scale |
---|
| 56 | answer *= pars->scale; |
---|
| 57 | |
---|
| 58 | //convert to [cm-1] |
---|
| 59 | //answer *= 1.0e8; |
---|
| 60 | |
---|
| 61 | return answer; |
---|
| 62 | } |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | /// 1D scattering function |
---|
| 66 | double simcylinder_analytical_1D(SimCylinderParameters *pars, double q) { |
---|
| 67 | /*** |
---|
| 68 | * Things to keep here: |
---|
| 69 | * - volume calc |
---|
| 70 | * - point generation |
---|
| 71 | * |
---|
| 72 | */ |
---|
| 73 | |
---|
| 74 | // Check if Rho array is available |
---|
| 75 | int volume_points; |
---|
| 76 | int r_points; |
---|
| 77 | int ptsGenerated; |
---|
| 78 | |
---|
| 79 | double r_step; |
---|
| 80 | double vol; |
---|
| 81 | double bin_width; |
---|
| 82 | double twopiq; |
---|
| 83 | double tmp; |
---|
| 84 | |
---|
| 85 | //printf("returning %g %g\n",test_analytical_2D(pars, q, 0.0), pars->length); |
---|
| 86 | return test_analytical_2D(pars, q, 0.0); |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | // These should be parameters |
---|
| 90 | vol = 1.0*acos(-1.0)*pars->radius*pars->radius*pars->length; |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | r_points = pars->calcPars.r_points; |
---|
| 94 | volume_points = pars->calcPars.volume_points; |
---|
| 95 | |
---|
| 96 | if(pars->calcPars.isPointMemAllocated==0) { |
---|
| 97 | |
---|
| 98 | // Call modelCalc function here to init_volume |
---|
| 99 | |
---|
| 100 | //twopiq = (2*acos(-1.0)*pars->qmax); |
---|
| 101 | //tmp = twopiq*twopiq*twopiq*vol; |
---|
| 102 | //volume_points = (int) floor(tmp); |
---|
| 103 | //r_points = (int) floor(10.0*pow(tmp,0.3333)); |
---|
| 104 | //r_points = (int) floor(pow(tmp,0.3333)); |
---|
| 105 | //printf("v, r = %d, %d\n",volume_points, r_points); |
---|
| 106 | |
---|
| 107 | // TEST |
---|
| 108 | volume_points = 1000; |
---|
| 109 | r_points = 1000; |
---|
| 110 | |
---|
| 111 | pars->calcPars.volume_points = volume_points; |
---|
| 112 | pars->calcPars.r_points = r_points; |
---|
| 113 | |
---|
| 114 | // Memory allocation |
---|
| 115 | pars->calcPars.points = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint)); |
---|
| 116 | if(pars->calcPars.points==NULL) { |
---|
| 117 | printf("Problem allocating memory for 1D volume points\n"); |
---|
| 118 | return -1.0; |
---|
| 119 | } |
---|
| 120 | pars->calcPars.isPointMemAllocated=1; |
---|
| 121 | } |
---|
| 122 | |
---|
| 123 | r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points; |
---|
| 124 | |
---|
| 125 | //printf("step=%g r=%g l=%g\n", r_step, pars->radius, pars->length); |
---|
| 126 | if(pars->calcPars.isRhoAvailable==0) { |
---|
| 127 | // Generate random points accross the volume |
---|
| 128 | ptsGenerated = simcylinder_generatePoints(pars->calcPars.points, volume_points, pars->radius, pars->length); |
---|
| 129 | |
---|
| 130 | // Consistency check |
---|
| 131 | if(ptsGenerated <= 0) { |
---|
| 132 | // Set error code here |
---|
| 133 | return 0; |
---|
| 134 | } |
---|
| 135 | |
---|
| 136 | // Allocate memory |
---|
| 137 | |
---|
| 138 | pars->calcPars.rho = (double*) malloc(r_points*sizeof(double)); |
---|
| 139 | if(pars->calcPars.rho==NULL){ |
---|
| 140 | printf("Problem allocating memory for 1D correlation points\n"); |
---|
| 141 | return 0; |
---|
| 142 | } |
---|
| 143 | |
---|
| 144 | |
---|
| 145 | if(modelcalculations_calculatePairCorrelation_1D(pars->calcPars.points, volume_points, pars->calcPars.rho, r_points, r_step)<=0){ |
---|
| 146 | printf("Error occured!\n"); |
---|
| 147 | return 0; |
---|
| 148 | }; |
---|
| 149 | pars->calcPars.isRhoAvailable=1; |
---|
| 150 | } |
---|
| 151 | |
---|
| 152 | // Calculate I(q,phi) and return that value |
---|
| 153 | |
---|
| 154 | return acos(-1.0)/2.0*modelcalculations_calculateIq_1D(pars->calcPars.rho, r_points, r_step, q) * vol; |
---|
| 155 | |
---|
| 156 | } |
---|
| 157 | |
---|
| 158 | /// 1D scattering function |
---|
| 159 | double simcylinder_analytical_2D(SimCylinderParameters *pars, double q, double phi) { |
---|
| 160 | // Check if Rho array is available |
---|
| 161 | int volume_points; |
---|
| 162 | int r_points; |
---|
| 163 | int ptsGenerated; |
---|
| 164 | double bin_width; |
---|
| 165 | double r_step; |
---|
| 166 | double vol; |
---|
| 167 | double retval; |
---|
| 168 | |
---|
| 169 | // These should be parameters |
---|
| 170 | //r_points = pars->calcPars.r_points; |
---|
| 171 | //volume_points = pars->calcPars.volume_points; |
---|
| 172 | volume_points = 50000; |
---|
| 173 | r_points = 100; |
---|
| 174 | if(pars->calcPars.isPointMemAllocated_2D==0) { |
---|
| 175 | //volume_points = 1000; |
---|
| 176 | //r_points = 100; |
---|
| 177 | pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint)); |
---|
| 178 | if(pars->calcPars.points_2D==NULL) { |
---|
| 179 | printf("Problem allocating memory for 2D volume points\n"); |
---|
| 180 | return -1.0; |
---|
| 181 | } |
---|
| 182 | pars->calcPars.isPointMemAllocated_2D=1; |
---|
| 183 | } |
---|
| 184 | |
---|
| 185 | r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points; |
---|
| 186 | if(pars->calcPars.isRhoAvailable_2D==0) { |
---|
| 187 | // Initialize random number generator |
---|
| 188 | int seed; |
---|
| 189 | seed = 10000; |
---|
| 190 | srand(seed); |
---|
| 191 | |
---|
| 192 | // Generate random points accross the volume |
---|
| 193 | ptsGenerated = simcylinder_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius, pars->length); |
---|
| 194 | |
---|
| 195 | // Calculate correlation function |
---|
| 196 | pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*sizeof(float)); |
---|
| 197 | if(pars->calcPars.rho_2D==NULL){ |
---|
| 198 | printf("Problem allocating memory for 2D correlations points\n"); |
---|
| 199 | return -1.0; |
---|
| 200 | } |
---|
| 201 | //if(modelcalculations_calculatePairCorrelation_2D(pars->calcPars.points_2D, volume_points, pars->calcPars.rho_2D, r_points, r_step)==NULL){ |
---|
| 202 | if(modelcalculations_calculatePairCorrelation_2D_vector(pars->calcPars.points_2D, |
---|
| 203 | volume_points, pars->calcPars.rho_2D, r_points, r_step, |
---|
| 204 | pars->theta, pars->phi,0.0)==0){ |
---|
| 205 | return 0; |
---|
| 206 | }; |
---|
| 207 | pars->calcPars.isRhoAvailable_2D=1; |
---|
| 208 | } |
---|
| 209 | // Calculate I(q,phi) and return that value |
---|
| 210 | vol = acos(-1.0)*pars->radius*pars->radius*pars->length; |
---|
| 211 | |
---|
| 212 | //printf("in ana_2D %f %f\n",q, phi); |
---|
| 213 | retval = modelcalculations_calculateIq_2D(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol; |
---|
| 214 | //printf("I=%g %f %f\n",retval, q, pars->theta); |
---|
| 215 | return acos(-1.0)/2.0*retval; |
---|
| 216 | } |
---|
| 217 | |
---|
| 218 | /// 1D scattering function |
---|
| 219 | double simcylinder_analytical_2D_3Darray(SimCylinderParameters *pars, double q, double phi) { |
---|
| 220 | // Check if Rho array is available |
---|
| 221 | int volume_points; |
---|
| 222 | int r_points; |
---|
| 223 | int ptsGenerated; |
---|
| 224 | double bin_width; |
---|
| 225 | double r_step; |
---|
| 226 | double vol; |
---|
| 227 | |
---|
| 228 | // These should be parameters |
---|
| 229 | volume_points = 5000; |
---|
| 230 | r_points = 100; |
---|
| 231 | |
---|
| 232 | if(pars->calcPars.isPointMemAllocated_2D==0) { |
---|
| 233 | pars->calcPars.points_2D = (SpacePoint*)malloc(volume_points*sizeof(SpacePoint)); |
---|
| 234 | if(pars->calcPars.points_2D==NULL) { |
---|
| 235 | printf("Problem allocating memory for 2D volume points\n"); |
---|
| 236 | return -1.0; |
---|
| 237 | } |
---|
| 238 | pars->calcPars.isPointMemAllocated_2D=1; |
---|
| 239 | } |
---|
| 240 | |
---|
| 241 | r_step = sqrt(pars->radius*pars->radius*4.0+pars->length*pars->length)/(double)r_points; |
---|
| 242 | |
---|
| 243 | if(pars->calcPars.isRhoAvailable_2D==0) { |
---|
| 244 | // Initialize random number generator |
---|
| 245 | int seed; |
---|
| 246 | seed = 10000; |
---|
| 247 | srand(seed); |
---|
| 248 | |
---|
| 249 | // Generate random points accross the volume |
---|
| 250 | ptsGenerated = simcylinder_generatePoints(pars->calcPars.points_2D, volume_points, pars->radius, pars->length); |
---|
| 251 | |
---|
| 252 | // Calculate correlation function |
---|
| 253 | // TODO: Check memory leaks |
---|
| 254 | printf("A\n"); |
---|
| 255 | pars->calcPars.rho_2D = (float*) malloc(r_points*r_points*r_points*sizeof(float)); |
---|
| 256 | if(pars->calcPars.rho_2D==NULL){ |
---|
| 257 | printf("Problem allocating memory for 2D correlations points\n"); |
---|
| 258 | return -1.0; |
---|
| 259 | } |
---|
| 260 | printf("B\n"); |
---|
| 261 | if(modelcalculations_calculatePairCorrelation_2D_3Darray(pars->calcPars.points_2D, |
---|
| 262 | volume_points, pars->calcPars.rho_2D, r_points, r_step)<0){ |
---|
| 263 | return 0; |
---|
| 264 | }; |
---|
| 265 | printf("C\n"); |
---|
| 266 | pars->calcPars.isRhoAvailable_2D=1; |
---|
| 267 | } |
---|
| 268 | // Calculate I(q,phi) and return that value |
---|
| 269 | vol = acos(-1.0)*pars->radius*pars->radius*pars->length; |
---|
| 270 | |
---|
| 271 | printf("in ana_2D %f %f\n",q, phi); |
---|
| 272 | return modelcalculations_calculateIq_2D_3Darray(pars->calcPars.rho_2D, r_points, r_step, q, phi)*vol; |
---|
| 273 | |
---|
| 274 | |
---|
| 275 | } |
---|
| 276 | |
---|
| 277 | |
---|
| 278 | /** |
---|
| 279 | * Generate points randomly accross the volume |
---|
| 280 | * @param points [SpacePoint*] Array of 3D points to be filled |
---|
| 281 | * @param n [int] Number of points to generat |
---|
| 282 | * @param radius [double] Radius of the sphere |
---|
| 283 | * @return Number of points generated |
---|
| 284 | */ |
---|
| 285 | int simcylinder_generatePoints(SpacePoint * points, int n, double radius, double length) { |
---|
| 286 | int i; |
---|
| 287 | int testcounter; |
---|
| 288 | double x, y, z; |
---|
| 289 | |
---|
| 290 | // Create points |
---|
| 291 | // To have a uniform density, you want to generate |
---|
| 292 | // random points in a box and keep only those that are |
---|
| 293 | // within the volume. |
---|
| 294 | |
---|
| 295 | // Initialize random number generator |
---|
| 296 | int seed; |
---|
| 297 | time_t now; |
---|
| 298 | |
---|
| 299 | time(&now); |
---|
| 300 | //seed = 10000; |
---|
| 301 | |
---|
| 302 | seed = (int)floor(fmod(now,10000)); |
---|
| 303 | //seed = 10009; |
---|
| 304 | srand(seed); |
---|
| 305 | printf("Seed = %i\n", seed); |
---|
| 306 | |
---|
| 307 | testcounter = 0; |
---|
| 308 | |
---|
| 309 | memset(points,0,n*sizeof(SpacePoint)); |
---|
| 310 | for(i=0;i<n;i++) { |
---|
| 311 | // Generate in a box centered around zero |
---|
| 312 | x = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius; |
---|
| 313 | y = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * radius; |
---|
| 314 | z = (2.0*((double)rand())/((double)(RAND_MAX)+(double)(1))-1.0) * length/2.0; |
---|
| 315 | |
---|
| 316 | // reject those that are not within the volume |
---|
| 317 | if( sqrt(x*x+y*y) < radius && fabs(z)<length/2.0) { |
---|
| 318 | points[i].x = x; |
---|
| 319 | points[i].y = y; |
---|
| 320 | points[i].z = z; |
---|
| 321 | testcounter++; |
---|
| 322 | } else { |
---|
| 323 | i--; |
---|
| 324 | } |
---|
| 325 | } |
---|
| 326 | //printf("test counter = %d\n", testcounter); |
---|
| 327 | |
---|
| 328 | // Consistency check |
---|
| 329 | if(testcounter != n) { |
---|
| 330 | return -1; |
---|
| 331 | } |
---|
| 332 | |
---|
| 333 | return testcounter; |
---|
| 334 | } |
---|