Changeset 767dca8 in sasmodels
 Timestamp:
 Oct 27, 2017 12:00:07 PM (22 months ago)
 Branches:
 master, core_shell_microgels, magnetic_model, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 da5536f
 Parents:
 a9bc435
 Location:
 sasmodels
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

sasmodels/details.py
r2bccb5a r767dca8 234 234 # skipping scale and background when building values and weights 235 235 values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 236 weights = correct_theta_weights(kernel.info.parameters, dispersity, weights)236 #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights) 237 237 length = np.array([len(w) for w in weights]) 238 238 offset = np.cumsum(np.hstack((0, length))) 
sasmodels/kernel_iq.c
r2c108a3 r767dca8 32 32 // INVALID(table) : test if the current point is feesible to calculate. This 33 33 // will be defined in the kernel definition file. 34 35 // Set SCALED_PHI to 1 if dphi represents distance rather than longitude. 36 // That is, in order to make an equally spaced mesh on a curved surface, 37 // you will need to scale longitude displacement by cos(latitude) to account 38 // for the reduction in distance per degree latitude as you move away from 39 // the equator. Points on the mesh with scaled dphi values more than 180 40 // degrees are not included in the dispersity calculation. This should make 41 // the integral over the entire surface work by sampling fewer points near 42 // the poles. 43 # define USE_SCALED_DPHI 1 34 44 35 45 #ifndef _PAR_BLOCK_ // protected block so we can include this code twice. … … 450 460 // unoriented 1D 451 461 double qk; 462 #define SCALE_DPHI_AND_SET_WEIGHT() const double weight=weight0 452 463 #define FETCH_Q() do { qk = q[q_index]; } while (0) 453 464 #define BUILD_ROTATION() do {} while(0) … … 458 469 // unoriented 2D 459 470 double qx, qy; 471 #define SCALE_DPHI_AND_SET_WEIGHT() const double weight=weight0 460 472 #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 461 473 #define BUILD_ROTATION() do {} while(0) … … 475 487 const double theta = values[details>theta_par+2]; 476 488 const double phi = values[details>theta_par+3]; 477 #define BUILD_ROTATION() qac_rotation(&rotation, \ 478 theta, phi, local_values.table.theta, local_values.table.phi) 489 double dtheta, dphi, weight; 490 #if USE_SCALED_DPHI 491 #define SCALE_DPHI_AND_SET_WEIGHT() do { \ 492 dtheta = local_values.table.theta; \ 493 dphi = local_values.table.phi; \ 494 weight = weight0; \ 495 if (dtheta != 90.0) dphi /= fabs(cos(dtheta*M_PI_180)); \ 496 else if (dphi != 0.0) weight = 0.; \ 497 if (fabs(dphi) >= 180.) weight = 0.; \ 498 } while (0) 499 #else 500 #define SCALE_DPHI_AND_SET_WEIGHT() do { \ 501 dtheta = local_values.table.theta; \ 502 dphi = local_values.table.phi; \ 503 weight = fabs(cos(dtheta*M_PI_180)) * weight0; \ 504 } while (0) 505 #endif 506 #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi); 479 507 #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc) 480 508 #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table) … … 493 521 const double phi = values[details>theta_par+3]; 494 522 const double psi = values[details>theta_par+4]; 495 #define BUILD_ROTATION() qabc_rotation(&rotation, \ 496 theta, phi, psi, local_values.table.theta, \ 497 local_values.table.phi, local_values.table.psi) 523 double dtheta, dphi, dpsi, weight; 524 #if USE_SCALED_DPHI 525 #define SCALE_DPHI_AND_SET_WEIGHT() do { \ 526 dtheta = local_values.table.theta; \ 527 dphi = local_values.table.phi; \ 528 dpsi = local_values.table.psi; \ 529 weight = weight0; \ 530 if (dtheta != 90.0) dphi /= fabs(cos(dtheta*M_PI_180)); \ 531 else if (dphi != 0.0) weight = 0.; \ 532 if (fabs(dphi) >= 180.) weight = 0.; \ 533 } while (0) 534 #else 535 #define SCALE_DPHI_AND_SET_WEIGHT() do { \ 536 dtheta = local_values.table.theta; \ 537 dphi = local_values.table.phi; \ 538 dpsi = local_values.table.psi; \ 539 weight = fabs(cos(dtheta*M_PI_180)) * weight0; \ 540 } while (0) 541 #endif 542 #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, dpsi) 498 543 #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 499 544 #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) … … 551 596 552 597 //if (q_index==0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n");} 598 SCALE_DPHI_AND_SET_WEIGHT(); 553 599 554 600 // ====== loop body ======= … … 559 605 // Accumulate I(q) 560 606 // Note: weight==0 must always be excluded 561 if (weight 0> cutoff) {562 pd_norm += weight 0* CALL_VOLUME(local_values.table);607 if (weight > cutoff) { 608 pd_norm += weight * CALL_VOLUME(local_values.table); 563 609 BUILD_ROTATION(); 564 610 … … 608 654 const double scattering = CALL_KERNEL(); 609 655 #endif // !MAGNETIC 610 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight , weight0);656 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 611 657 612 658 #ifdef USE_OPENCL 613 this_result += weight 0* scattering;659 this_result += weight * scattering; 614 660 #else // !USE_OPENCL 615 result[q_index] += weight 0* scattering;661 result[q_index] += weight * scattering; 616 662 #endif // !USE_OPENCL 617 663 } … … 651 697 #undef PD_OPEN 652 698 #undef PD_CLOSE 699 #undef SCALE_DPHI_AND_SET_WEIGHT 653 700 #undef FETCH_Q 654 701 #undef BUILD_ROTATION
Note: See TracChangeset
for help on using the changeset viewer.