Changeset 767dca8 in sasmodels


Ignore:
Timestamp:
Oct 27, 2017 12:00:07 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
da5536f
Parents:
a9bc435
Message:

use alternative jitter, which maintains constant arc length for the phi distribution

Location:
sasmodels
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/details.py

    r2bccb5a r767dca8  
    234234    # skipping scale and background when building values and weights 
    235235    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) 
    237237    length = np.array([len(w) for w in weights]) 
    238238    offset = np.cumsum(np.hstack((0, length))) 
  • sasmodels/kernel_iq.c

    r2c108a3 r767dca8  
    3232//  INVALID(table) : test if the current point is feesible to calculate.  This 
    3333//      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 
    3444 
    3545#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
     
    450460  // unoriented 1D 
    451461  double qk; 
     462  #define SCALE_DPHI_AND_SET_WEIGHT() const double weight=weight0 
    452463  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    453464  #define BUILD_ROTATION() do {} while(0) 
     
    458469  // unoriented 2D 
    459470  double qx, qy; 
     471  #define SCALE_DPHI_AND_SET_WEIGHT() const double weight=weight0 
    460472  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
    461473  #define BUILD_ROTATION() do {} while(0) 
     
    475487  const double theta = values[details->theta_par+2]; 
    476488  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); 
    479507  #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc) 
    480508  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table) 
     
    493521  const double phi = values[details->theta_par+3]; 
    494522  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) 
    498543  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
    499544  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     
    551596 
    552597//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(); 
    553599 
    554600  // ====== loop body ======= 
     
    559605    // Accumulate I(q) 
    560606    // Note: weight==0 must always be excluded 
    561     if (weight0 > cutoff) { 
    562       pd_norm += weight0 * CALL_VOLUME(local_values.table); 
     607    if (weight > cutoff) { 
     608      pd_norm += weight * CALL_VOLUME(local_values.table); 
    563609      BUILD_ROTATION(); 
    564610 
     
    608654          const double scattering = CALL_KERNEL(); 
    609655        #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); 
    611657 
    612658        #ifdef USE_OPENCL 
    613           this_result += weight0 * scattering; 
     659          this_result += weight * scattering; 
    614660        #else // !USE_OPENCL 
    615           result[q_index] += weight0 * scattering; 
     661          result[q_index] += weight * scattering; 
    616662        #endif // !USE_OPENCL 
    617663      } 
     
    651697#undef PD_OPEN 
    652698#undef PD_CLOSE 
     699#undef SCALE_DPHI_AND_SET_WEIGHT 
    653700#undef FETCH_Q 
    654701#undef BUILD_ROTATION 
Note: See TracChangeset for help on using the changeset viewer.