Changeset 885753a in sasmodels


Ignore:
Timestamp:
Apr 6, 2018 10:35:22 AM (4 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:
dc6f601
Parents:
5770493
Message:

correct magnetic calculations; fractional polarization no longer matches 3.1.2

Location:
sasmodels
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/details.py

    r108e70e r885753a  
    243243    offset = np.cumsum(np.hstack((0, length))) 
    244244    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
    245     # Pad value array to a 32 value boundaryd 
     245    # Pad value array to a 32 value boundary 
    246246    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    247247    extra = (32 - data_len%32)%32 
     
    250250    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
    251251    #call_details.show() 
     252    #print("data", data) 
    252253    return call_details, data, is_magnetic 
    253254 
     
    296297    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 
    297298    mag = mag.reshape(-1, 3) 
    298     scale = mag[:, 0] 
    299     if np.any(scale): 
     299    if np.any(mag[:, 0] != 0.0): 
     300        M0 = mag[:, 0].copy() 
    300301        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    301         cos_theta = cos(theta) 
    302         mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
    303         mag[:, 1] = scale*sin(theta)  # my 
    304         mag[:, 2] = -scale*cos_theta*sin(phi)  # mz 
     302        mag[:, 0] = +M0*cos(theta)*cos(phi)  # mx 
     303        mag[:, 1] = +M0*sin(theta) # my 
     304        mag[:, 2] = -M0*cos(theta)*sin(phi)  # mz 
    305305        return True 
    306306    else: 
  • sasmodels/kernel_iq.c

    rd86f0fc r885753a  
    7979//     du * (m_sigma_y + 1j*m_sigma_z); 
    8080// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag 
    81 static void set_spin_weights(double in_spin, double out_spin, double spins[4]) 
     81static void set_spin_weights(double in_spin, double out_spin, double weight[6]) 
    8282{ 
    8383  in_spin = clip(in_spin, 0.0, 1.0); 
    8484  out_spin = clip(out_spin, 0.0, 1.0); 
    85   spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd 
    86   spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real 
    87   spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real 
    88   spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu 
    89   spins[4] = spins[1]; // du imag 
    90   spins[5] = spins[2]; // ud imag 
     85  // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that 
     86  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
     87  // which is likely to be the case for simple models. 
     88  weight[0] = ((1.0-in_spin) * (1.0-out_spin)); // dd 
     89  weight[1] = ((1.0-in_spin) * out_spin);       // du.real 
     90  weight[2] = (in_spin * (1.0-out_spin));       // ud.real 
     91  weight[3] = (in_spin * out_spin);             // uu 
     92  weight[4] = weight[1]; // du.imag 
     93  weight[5] = weight[2]; // ud.imag 
    9194} 
    9295 
    9396// Compute the magnetic sld 
    9497static double mag_sld( 
    95   const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
     98  const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag 
    9699  const double qx, const double qy, 
    97100  const double px, const double py, 
     
    106109      case 0: // uu => sld - D M_perpx 
    107110          return sld - px*perp; 
    108       case 1: // ud real => -D M_perpy 
     111      case 1: // ud.real => -D M_perpy 
    109112          return py*perp; 
    110       case 2: // du real => -D M_perpy 
     113      case 2: // du.real => -D M_perpy 
    111114          return py*perp; 
    112       case 3: // dd real => sld + D M_perpx 
     115      case 3: // dd => sld + D M_perpx 
    113116          return sld + px*perp; 
    114117    } 
    115118  } else { 
    116119    if (xs== 4) { 
    117       return -mz;  // ud imag => -D M_perpz 
     120      return -mz;  // du.imag => +D M_perpz 
    118121    } else { // index == 5 
    119       return mz;   // du imag => D M_perpz 
     122      return +mz;  // ud.imag => -D M_perpz 
    120123    } 
    121124  } 
     
    312315  //     up_angle = values[NUM_PARS+4]; 
    313316  // TODO: could precompute more magnetism parameters before calling the kernel. 
    314   double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
     317  double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
    315318  double cos_mspin, sin_mspin; 
    316   set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
     319  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights); 
    317320  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    318321#endif // MAGNETIC 
     
    661664            // loop over uu, ud real, du real, dd, ud imag, du imag 
    662665            for (unsigned int xs=0; xs<6; xs++) { 
    663               const double xs_weight = spins[xs]; 
     666              const double xs_weight = xs_weights[xs]; 
    664667              if (xs_weight > 1.e-8) { 
    665668                // Since the cross section weight is significant, set the slds 
     
    674677                  local_values.vector[sld_index] = 
    675678                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz); 
     679//if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n", 
     680//q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz); 
    676681                } 
    677682                scattering += xs_weight * CALL_KERNEL(); 
Note: See TracChangeset for help on using the changeset viewer.