Changeset c4e7a5f in sasmodels


Ignore:
Timestamp:
Apr 6, 2016 1:39:39 PM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
3a45c2c, 0278e3f
Parents:
4d76711
Message:

fix spherical correction in old style kernels

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/kernel_template.c

    rcf85329 rc4e7a5f  
    158158#ifdef IQ_OPEN_LOOPS 
    159159    double ret=0.0, norm=0.0; 
    160   #ifdef VOLUME_PARAMETERS 
    161     double vol=0.0, norm_vol=0.0; 
    162   #endif 
    163160    IQ_OPEN_LOOPS 
    164161    //for (int radius_i=0; radius_i < Nradius; radius_i++) { 
     
    172169      if (!isnan(scattering)) { 
    173170        ret += weight*scattering; 
     171      #ifdef VOLUME_PARAMETERS 
     172        norm += weight * form_volume(VOLUME_PARAMETERS); 
     173      #else 
    174174        norm += weight; 
    175       #ifdef VOLUME_PARAMETERS 
    176         const double vol_weight = VOLUME_WEIGHT_PRODUCT; 
    177         vol += vol_weight*form_volume(VOLUME_PARAMETERS); 
    178         norm_vol += vol_weight; 
    179175      #endif 
    180176      } 
     
    182178    } 
    183179    IQ_CLOSE_LOOPS 
    184   #ifdef VOLUME_PARAMETERS 
    185     if (vol*norm_vol != 0.0) { 
    186       ret *= norm_vol/vol; 
    187     } 
    188   #endif 
    189     result[i] = scale*ret/norm+background; 
     180    // norm can only be zero if volume is zero, so no scattering 
     181    result[i] = (norm > 0. ? scale*ret/norm + background : background); 
    190182#else 
    191183    result[i] = scale*Iq(qi, IQ_PARAMETERS) + background; 
     
    241233    //  const double radius = loops[2*(radius_i)]; 
    242234    //  const double radius_w = loops[2*(radius_i)+1]; 
    243  
    244     const double weight = IQXY_WEIGHT_PRODUCT; 
     235    double weight = IQXY_WEIGHT_PRODUCT; 
    245236    if (weight > cutoff) { 
    246237 
    247238      const double scattering = Iqxy(qxi, qyi, IQXY_PARAMETERS); 
    248239      if (!isnan(scattering)) { // if scattering is bad, exclude it from sum 
    249       //if (scattering >= 0.0) { // scattering cannot be negative 
    250         // TODO: use correct angle for spherical correction 
    251         // Definition of theta and phi are probably reversed relative to the 
    252         // equation which gave rise to this correction, leading to an 
    253         // attenuation of the pattern as theta moves through pi/2.  Either 
    254         // reverse the meanings of phi and theta in the forms, or use phi 
    255         // rather than theta in this correction.  Current code uses cos(theta) 
    256         // so that values match those of sasview. 
    257       #if defined(IQXY_HAS_THETA) // && 0 
    258         const double spherical_correction 
    259           = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2:1.0); 
    260         const double next = spherical_correction * weight * scattering; 
    261       #else 
    262         const double next = weight * scattering; 
    263       #endif 
     240      #if defined(IQXY_HAS_THETA) 
     241        // Force a nominal value for the spherical correction even when 
     242        // theta is +90/-90 so that there are no divide by zero problems. 
     243        // For cos(theta) fixed at 90, we effectively multiply top and bottom 
     244        // by 1e-6, so the effect cancels. 
     245        const double spherical_correction = fmax(fabs(cos(M_PI_180*theta)), 1e-6); 
     246        weight *= spherical_correction; 
     247      #endif 
     248      const double next = weight * scattering; 
    264249      #if USE_KAHAN_SUMMATION 
    265250        const double y = next - accumulated_error; 
     
    270255        ret += next; 
    271256      #endif 
     257      #ifdef VOLUME_PARAMETERS 
     258        norm += weight*form_volume(VOLUME_PARAMETERS); 
     259      #else 
    272260        norm += weight; 
    273       #ifdef VOLUME_PARAMETERS 
    274         const double vol_weight = VOLUME_WEIGHT_PRODUCT; 
    275         vol += vol_weight*form_volume(VOLUME_PARAMETERS); 
    276         norm_vol += vol_weight; 
    277261      #endif 
    278262      } 
     
    280264    } 
    281265    IQXY_CLOSE_LOOPS 
    282   #ifdef VOLUME_PARAMETERS 
    283     if (vol*norm_vol != 0.0) { 
    284       ret *= norm_vol/vol; 
    285     } 
    286   #endif 
    287     result[i] = scale*ret/norm+background; 
     266    // norm can only be zero if volume is zero, so no scattering 
     267    result[i] = (norm>0. ? scale*ret/norm + background : background); 
    288268#else 
    289269    result[i] = scale*Iqxy(qxi, qyi, IQXY_PARAMETERS) + background; 
Note: See TracChangeset for help on using the changeset viewer.