Changes in / [0159b5e:db472a5] in sasmodels


Ignore:
Files:
4 added
68 edited

Legend:

Unmodified
Added
Removed

  • explore/beta/sasfit_compare.py

    r7b0abf8 r2a12351b  
    208208    F1, F2 = np.zeros_like(q), np.zeros_like(q) 
    209209    for k, Rpk in enumerate(Rp_val): 
     210        print("ellipsoid cycle", k, "of", len(Rp_val)) 
    210211        for i, Rei in enumerate(Re_val): 
    211212            theory = ellipsoid_theta(q, Rpk, Rei, sld, sld_solvent) 
    … …  
    343344    #Sq = Iq/Pq 
    344345    #Iq = None#= Sq = None 
    345     r = I._kernel.results 
    346     return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r[1], Ibeta=Iq) 
     346    r = dict(I._kernel.results()) 
     347    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r["S_eff(Q)"], Ibeta=Iq) 
    347348 
    348349def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): 
    … …  
    401402        radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, 
    402403        background=0, 
    403         radius_polar_pd=.1, radius_polar_pd_type=pd_type, 
    404         radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type, 
     404        radius_polar_pd=0.1, radius_polar_pd_type=pd_type, 
     405        radius_equatorial_pd=0.1, radius_equatorial_pd_type=pd_type, 
    405406        volfraction=0.15, 
    406407        radius_effective=270.7543927018, 
    407408        #radius_effective=12.59921049894873, 
    408409        ) 
    409     target = sasmodels_theory(q, model, beta_mode=1, **pars) 
     410    target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
    410411    actual = ellipsoid_pe(q, norm='sasview', **pars) 
    411412    title = " ".join(("sasmodels", model, pd_type)) 

  • explore/precision.py

    r2a7e20e ree60aa7  
    357357) 
    358358add_function( 
    359     name="debye", 
     359    name="gauss_coil", 
    360360    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4, 
    361361    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4, 
    362362    ocl_function=make_ocl(""" 
    363363    const double qsq = q*q; 
    364     if (qsq < 1.0) { // Pade approximation 
     364    // For double: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide) 
     365    // For single: use O(7) Taylor with 0.8 cutoff (7 mad) 
     366    if (qsq < 0.0) { 
    365367        const double x = qsq; 
    366368        if (0) { // 0.36 single 
    … …  
    372374            const double B1=3./8., B2=3./56., B3=1./336.; 
    373375            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.); 
    374         } else if (1) { // 1.0 for single, 0.25 for double 
     376        } else if (0) { // 1.0 for single, 0.25 for double 
    375377            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
    376378            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
    … …  
    385387                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
    386388        } 
    387     } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double 
     389    } else if (qsq < 0.8) { 
    388390        const double x = qsq; 
    389391        const double C0 = +1.; 

  • sasmodels/core.py

    r71b751d ree60aa7  
    363363    model = load_model("cylinder@hardsphere*sphere") 
    364364    actual = [p.name for p in model.info.parameters.kernel_parameters] 
    365     target = ("sld sld_solvent radius length theta phi volfraction" 
    366               " beta_mode A_sld A_sld_solvent A_radius").split() 
     365    target = ("sld sld_solvent radius length theta phi" 
     366              " radius_effective volfraction " 
     367              " structure_factor_mode radius_effective_mode" 
     368              " A_sld A_sld_solvent A_radius").split() 
    367369    assert target == actual, "%s != %s"%(target, actual) 
    368370 

  • sasmodels/generate.py

    r6e45516 r5399809  
    802802                              for p in partable.kernel_parameters)) 
    803803    # Define the function calls 
     804    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, v) 0.0" 
    804805    if partable.form_volume_parameters: 
    805806        refs = _call_pars("_v.", partable.form_volume_parameters) 
    806807        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs)) 
     808        if model_info.effective_radius_type: 
     809            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, _v) effective_radius(mode, %s)"%(",".join(refs)) 
    807810    else: 
    808811        # Model doesn't have volume.  We could make the kernel run a little 
    … …  
    811814        call_volume = "#define CALL_VOLUME(v) 1.0" 
    812815    source.append(call_volume) 
     816    source.append(call_effective_radius) 
    813817    model_refs = _call_pars("_v.", partable.iq_parameters) 
    814818 

  • sasmodels/kernel.py

    r2d81cfe ree60aa7  
    4141    dtype = None  # type: np.dtype 
    4242 
    43     def __call__(self, call_details, values, cutoff, magnetic): 
     43    def Iq(self, call_details, values, cutoff, magnetic): 
    4444        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    45         raise NotImplementedError("need to implement __call__") 
     45        Pq, Reff = self.Pq_Reff(call_details, values, cutoff, magnetic, effective_radius_type=0) 
     46        return Pq 
     47    __call__ = Iq 
     48 
     49    def Pq_Reff(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     50        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     51        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     52        #print("returned",self.q_input.q, self.result) 
     53        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     54        total_weight = self.result[nout*self.q_input.nq + 0] 
     55        if total_weight == 0.: 
     56            total_weight = 1. 
     57        weighted_volume = self.result[nout*self.q_input.nq + 1] 
     58        weighted_radius = self.result[nout*self.q_input.nq + 2] 
     59        effective_radius = weighted_radius/total_weight 
     60        # compute I = scale*P + background 
     61        #           = scale*(sum(w*F^2)/sum w)/(sum (w*V)/sum w) + background 
     62        #           = scale/sum (w*V) * sum(w*F^2) + background 
     63        F2 = self.result[0:nout*self.q_input.nq:nout] 
     64        scale = values[0]/(weighted_volume if weighted_volume != 0.0 else 1.0) 
     65        background = values[1] 
     66        Pq = scale*F2 + background 
     67        #print("scale",scale,background) 
     68        return Pq, effective_radius 
     69 
     70    def beta(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     71        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     72        if self.dim == '2d': 
     73            raise NotImplementedError("beta not yet supported for 2D") 
     74        if not self.info.have_Fq: 
     75            raise NotImplementedError("beta not yet supported for "+self.info.id) 
     76        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     77        total_weight = self.result[2*self.q_input.nq + 0] 
     78        if total_weight == 0.: 
     79            total_weight = 1. 
     80        weighted_volume = self.result[2*self.q_input.nq + 1] 
     81        weighted_radius = self.result[2*self.q_input.nq + 2] 
     82        volume_average = weighted_volume/total_weight 
     83        effective_radius = weighted_radius/total_weight 
     84        F2 = self.result[0:2*self.q_input.nq:2]/total_weight 
     85        F1 = self.result[1:2*self.q_input.nq:2]/total_weight 
     86        return F1, F2, volume_average, effective_radius 
    4687 
    4788    def release(self): 
    4889        # type: () -> None 
    4990        pass 
     91 
     92    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     93        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     94        raise NotImplementedError() 

  • sasmodels/kernel_header.c

    r108e70e r296c52b  
    6868         #define NEED_EXPM1 
    6969         #define NEED_TGAMMA 
     70         #define NEED_CBRT 
    7071         // expf missing from windows? 
    7172         #define expf exp 
    … …  
    8586#  define pown(a,b) pow(a,b) 
    8687#endif // !USE_OPENCL 
     88 
     89#if defined(NEED_CBRT) 
     90   #define cbrt(_x) pow(_x, 0.33333333333333333333333) 
     91#endif 
    8792 
    8893#if defined(NEED_EXPM1) 

  • sasmodels/kernel_iq.c

    r70530778 r70530778  
    2727//      parameters in the parameter table. 
    2828//  CALL_VOLUME(table) : call the form volume function 
     29//  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function 
    2930//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
    3031//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
    … …  
    284285    global const double *q, // nq q values, with padding to boundary 
    285286    global double *result,  // nq+1 return values, again with padding 
    286     const double cutoff     // cutoff in the dispersity weight product 
     287    const double cutoff,     // cutoff in the dispersity weight product 
     288    int32_t effective_radius_type // which effective radius to compute 
    287289    ) 
    288290{ 
    … …  
    344346  #ifdef USE_OPENCL 
    345347    #if COMPUTE_F1_F2 
    346       double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    347       double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
    348       double this_F2 = (pd_start == 0 ? 0.0 : result[q_index]); 
    349       double this_F1 = (pd_start == 0 ? 0.0 : result[q_index+nq+1]); 
     348      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 
     349      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     350      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]); 
     351      double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]); 
     352      double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]); 
    350353    #else 
    351       double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     354      double weight_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     355      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]); 
     356      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]); 
    352357      double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
    353358    #endif 
    354359  #else // !USE_OPENCL 
    355360    #if COMPUTE_F1_F2 
    356       double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    357       double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     361      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 
     362      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     363      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]); 
    358364    #else 
    359       double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     365      double weight_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     366      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]); 
     367      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]); 
    360368    #endif 
    361369    if (pd_start == 0) { 
    … …  
    364372      #endif 
    365373      #if COMPUTE_F1_F2 
    366           for (int q_index=0; q_index < 2*nq+2; q_index++) result[q_index] = 0.0; 
     374          // 2*nq for F^2,F pairs 
     375          for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0; 
    367376      #else 
    368377          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
    369378      #endif 
    370379    } 
    371     //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
     380    //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_volume, result[0]); 
    372381#endif // !USE_OPENCL 
    373382 
    … …  
    684693    // Note: weight==0 must always be excluded 
    685694    if (weight > cutoff) { 
    686       pd_norm += weight * CALL_VOLUME(local_values.table); 
     695      weighted_volume += weight * CALL_VOLUME(local_values.table); 
    687696      #if COMPUTE_F1_F2 
    688697      weight_norm += weight; 
    689698      #endif 
     699      if (effective_radius_type != 0) { 
     700        weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 
     701      } 
    690702      BUILD_ROTATION(); 
    691703 
    … …  
    745757        #ifdef USE_OPENCL 
    746758          #if COMPUTE_F1_F2 
     759            this_F2 += weight * F2; 
    747760            this_F1 += weight * F1; 
    748             this_F2 += weight * F2; 
    749761          #else 
    750762            this_result += weight * scattering; 
    … …  
    752764        #else // !USE_OPENCL 
    753765          #if COMPUTE_F1_F2 
    754             result[q_index] += weight * F2; 
    755             result[q_index+nq+1] += weight * F1; 
     766            result[2*q_index+0] += weight * F2; 
     767            result[2*q_index+1] += weight * F1; 
    756768          #else 
    757769            result[q_index] += weight * scattering; 
    … …  
    782794#ifdef USE_OPENCL 
    783795  #if COMPUTE_F1_F2 
    784     result[q_index] = this_F2; 
    785     result[q_index+nq+1] = this_F1; 
     796    result[2*q_index+0] = this_F2; 
     797    result[2*q_index+1] = this_F1; 
    786798    if (q_index == 0) { 
    787       result[nq] = pd_norm; 
    788       result[2*nq+1] = weight_norm; 
     799      result[2*nq+0] = weight_norm; 
     800      result[2*nq+1] = weighted_volume; 
     801      result[2*nq+2] = weighted_radius; 
    789802    } 
    790803  #else 
    791804    result[q_index] = this_result; 
    792     if (q_index == 0) result[nq] = pd_norm; 
     805    if (q_index == 0) { 
     806      result[nq+0] = weight_norm; 
     807      result[nq+1] = weighted_volume; 
     808      result[nq+2] = weighted_radius; 
     809    } 
    793810  #endif 
    794811 
    795 //if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
     812//if (q_index == 0) printf("res: %g/%g\n", result[0], weigthed_volume); 
    796813#else // !USE_OPENCL 
    797814  #if COMPUTE_F1_F2 
    798     result[nq] = pd_norm; 
    799     result[2*nq+1] = weight_norm; 
     815    result[2*nq] = weight_norm; 
     816    result[2*nq+1] = weighted_volume; 
     817    result[2*nq+2] = weighted_radius; 
    800818  #else 
    801     result[nq] = pd_norm; 
     819    result[nq] = weight_norm; 
     820    result[nq+1] = weighted_volume; 
     821    result[nq+2] = weighted_radius; 
    802822  #endif 
    803 //printf("res: %g/%g\n", result[0], pd_norm); 
     823//printf("res: %g/%g\n", result[0], weighted_volume); 
    804824#endif // !USE_OPENCL 
    805825 

  • sasmodels/kernelcl.py

    rc036ddb r5399809  
    481481        # at this point, so instead using 32, which is good on the set of 
    482482        # architectures tested so far. 
     483        extra_q = 3  # total weight, weighted volume and weighted radius 
    483484        if self.is_2d: 
    484             # Note: 16 rather than 15 because result is 1 longer than input. 
    485             width = ((self.nq+16)//16)*16 
     485            width = ((self.nq+15+extra_q)//16)*16 
    486486            self.q = np.empty((width, 2), dtype=dtype) 
    487487            self.q[:self.nq, 0] = q_vectors[0] 
    488488            self.q[:self.nq, 1] = q_vectors[1] 
    489489        else: 
    490             # Note: 32 rather than 31 because result is 1 longer than input. 
    491             width = ((self.nq+32)//32)*32 
     490            width = ((self.nq+31+extra_q)//32)*32 
    492491            self.q = np.empty(width, dtype=dtype) 
    493492            self.q[:self.nq] = q_vectors[0] 
    … …  
    539538        self.dim = '2d' if q_input.is_2d else '1d' 
    540539        # leave room for f1/f2 results in case we need to compute beta for 1d models 
    541         num_returns = 1 if self.dim == '2d' else 2  # 
    542         # plus 1 for the normalization value 
    543         self.result = np.empty((q_input.nq+1)*num_returns, dtype) 
     540        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     541        # plus 3 weight, volume, radius 
     542        self.result = np.empty(q_input.nq*nout + 3, self.dtype) 
    544543 
    545544        # Inputs and outputs for each kernel call 
    … …  
    549548 
    550549        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    551                                   q_input.global_size[0] * num_returns * dtype.itemsize) 
     550                                  q_input.global_size[0] * nout * dtype.itemsize) 
    552551        self.q_input = q_input # allocated by GpuInput above 
    553552 
    … …  
    558557                     else np.float32)  # will never get here, so use np.float32 
    559558 
    560     def Iq(self, call_details, values, cutoff, magnetic): 
    561         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    562         self._call_kernel(call_details, values, cutoff, magnetic) 
    563         #print("returned",self.q_input.q, self.result) 
    564         pd_norm = self.result[self.q_input.nq] 
    565         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    566         background = values[1] 
    567         #print("scale",scale,background) 
    568         return scale*self.result[:self.q_input.nq] + background 
    569     __call__ = Iq 
    570  
    571     def beta(self, call_details, values, cutoff, magnetic): 
    572         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    573         if self.dim == '2d': 
    574             raise NotImplementedError("beta not yet supported for 2D") 
    575         self._call_kernel(call_details, values, cutoff, magnetic) 
    576         w_norm = self.result[2*self.q_input.nq + 1] 
    577         pd_norm = self.result[self.q_input.nq] 
    578         if w_norm == 0.: 
    579             w_norm = 1. 
    580         F2 = self.result[:self.q_input.nq]/w_norm 
    581         F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
    582         volume_avg = pd_norm/w_norm 
    583         return F1, F2, volume_avg 
    584  
    585     def _call_kernel(self, call_details, values, cutoff, magnetic): 
     559    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    586560        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    587561        context = self.queue.context 
    … …  
    597571            details_b, values_b, self.q_input.q_b, self.result_b, 
    598572            self.real(cutoff), 
     573            np.uint32(effective_radius_type), 
    599574        ] 
    600575        #print("Calling OpenCL") 

  • sasmodels/kerneldll.py

    r2f8cbb9 r6e7ba14  
    315315 
    316316        # int, int, int, int*, double*, double*, double*, double*, double 
    317         argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type] 
     317        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32] 
    318318        names = [generate.kernel_name(self.info, variant) 
    319319                 for variant in ("Iq", "Iqxy", "Imagnetic")] 
    … …  
    382382        self.dim = '2d' if q_input.is_2d else '1d' 
    383383        # leave room for f1/f2 results in case we need to compute beta for 1d models 
    384         num_returns = 1 if self.dim == '2d' else 2  # 
    385         # plus 1 for the normalization value 
    386         self.result = np.empty((q_input.nq+1)*num_returns, self.dtype) 
     384        nout = 2 if self.info.have_Fq else 1 
     385        # plus 3 weight, volume, radius 
     386        self.result = np.empty(q_input.nq*nout + 3, self.dtype) 
    387387        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    388388                     else np.float64 if self.q_input.dtype == generate.F64 
    389389                     else np.float128) 
    390390 
    391     def Iq(self, call_details, values, cutoff, magnetic): 
    392         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    393         self._call_kernel(call_details, values, cutoff, magnetic) 
    394         #print("returned",self.q_input.q, self.result) 
    395         pd_norm = self.result[self.q_input.nq] 
    396         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    397         background = values[1] 
    398         #print("scale",scale,background) 
    399         return scale*self.result[:self.q_input.nq] + background 
    400     __call__ = Iq 
    401  
    402     def beta(self, call_details, values, cutoff, magnetic): 
    403         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    404         if self.dim == '2d': 
    405             raise NotImplementedError("beta not yet supported for 2D") 
    406         self._call_kernel(call_details, values, cutoff, magnetic) 
    407         w_norm = self.result[2*self.q_input.nq + 1] 
    408         pd_norm = self.result[self.q_input.nq] 
    409         if w_norm == 0.: 
    410             w_norm = 1. 
    411         F2 = self.result[:self.q_input.nq]/w_norm 
    412         F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
    413         volume_avg = pd_norm/w_norm 
    414         return F1, F2, volume_avg 
    415  
    416     def _call_kernel(self, call_details, values, cutoff, magnetic): 
     391    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     392        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    417393        kernel = self.kernel[1 if magnetic else 0] 
    418394        args = [ 
    … …  
    425401            self.result.ctypes.data,   # results 
    426402            self.real(cutoff), # cutoff 
     403            effective_radius_type, # cutoff 
    427404        ] 
    428405        #print("Calling DLL") 

  • sasmodels/kernelpy.py

    r108e70e r5399809  
    1212 
    1313import numpy as np  # type: ignore 
     14from numpy import pi 
     15try: 
     16    from numpy import cbrt 
     17except ImportError: 
     18    def cbrt(x): return x ** (1.0/3.0) 
    1419 
    1520from .generate import F64 
    … …  
    158163 
    159164        # Generate a closure which calls the form_volume if it exists. 
    160         form_volume = model_info.form_volume 
    161         self._volume = ((lambda: form_volume(*volume_args)) if form_volume else 
    162                         (lambda: 1.0)) 
    163  
    164     def __call__(self, call_details, values, cutoff, magnetic): 
     165        self._volume_args = volume_args 
     166        volume = model_info.form_volume 
     167        radius = model_info.effective_radius 
     168        self._volume = ((lambda: volume(*volume_args)) if volume 
     169                        else (lambda: 1.0)) 
     170        self._radius = ((lambda mode: radius(mode, *volume_args)) if radius 
     171                        else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume 
     172                        else (lambda mode: 1.0)) 
     173 
     174 
     175 
     176    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    165177        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    166178        if magnetic: 
    … …  
    168180        #print("Calling python kernel") 
    169181        #call_details.show(values) 
    170         res = _loops(self._parameter_vector, self._form, self._volume, 
    171                      self.q_input.nq, call_details, values, cutoff) 
    172         return res 
     182        radius = ((lambda: 0.0) if effective_radius_type == 0 
     183                  else (lambda: self._radius(effective_radius_type))) 
     184        self.result = _loops( 
     185            self._parameter_vector, self._form, self._volume, radius, 
     186            self.q_input.nq, call_details, values, cutoff) 
    173187 
    174188    def release(self): 
    … …  
    183197           form,          # type: Callable[[], np.ndarray] 
    184198           form_volume,   # type: Callable[[], float] 
     199           form_radius,   # type: Callable[[], float] 
    185200           nq,            # type: int 
    186201           call_details,  # type: details.CallDetails 
    187202           values,        # type: np.ndarray 
    188            cutoff         # type: float 
     203           cutoff,        # type: float 
    189204          ): 
    190205    # type: (...) -> None 
    … …  
    198213    #                                                              # 
    199214    ################################################################ 
     215 
     216    # WARNING: Trickery ahead 
     217    # The parameters[] vector is embedded in the closures for form(), 
     218    # form_volume() and form_radius().  We set the initial vector from 
     219    # the values for the model parameters. As we loop through the polydispesity 
     220    # mesh, we update the components with the polydispersity values before 
     221    # calling the respective functions. 
    200222    n_pars = len(parameters) 
    201223    parameters[:] = values[2:n_pars+2] 
     224 
    202225    if call_details.num_active == 0: 
    203         pd_norm = float(form_volume()) 
    204         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    205         background = values[1] 
    206         return scale*form() + background 
    207  
    208     pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] 
    209     pd_weight = values[2+n_pars + call_details.num_weights:] 
    210  
    211     pd_norm = 0.0 
    212     partial_weight = np.NaN 
    213     weight = np.NaN 
    214  
    215     p0_par = call_details.pd_par[0] 
    216     p0_length = call_details.pd_length[0] 
    217     p0_index = p0_length 
    218     p0_offset = call_details.pd_offset[0] 
    219  
    220     pd_par = call_details.pd_par[:call_details.num_active] 
    221     pd_offset = call_details.pd_offset[:call_details.num_active] 
    222     pd_stride = call_details.pd_stride[:call_details.num_active] 
    223     pd_length = call_details.pd_length[:call_details.num_active] 
    224  
    225     total = np.zeros(nq, 'd') 
    226     for loop_index in range(call_details.num_eval): 
    227         # update polydispersity parameter values 
    228         if p0_index == p0_length: 
    229             pd_index = (loop_index//pd_stride)%pd_length 
    230             parameters[pd_par] = pd_value[pd_offset+pd_index] 
    231             partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
    232             p0_index = loop_index%p0_length 
    233  
    234         weight = partial_weight * pd_weight[p0_offset + p0_index] 
    235         parameters[p0_par] = pd_value[p0_offset + p0_index] 
    236         p0_index += 1 
    237         if weight > cutoff: 
    238             # Call the scattering function 
    239             # Assume that NaNs are only generated if the parameters are bad; 
    240             # exclude all q for that NaN.  Even better would be to have an 
    241             # INVALID expression like the C models, but that is too expensive. 
    242             Iq = np.asarray(form(), 'd') 
    243             if np.isnan(Iq).any(): 
    244                 continue 
    245  
    246             # update value and norm 
    247             total += weight * Iq 
    248             pd_norm += weight * form_volume() 
    249  
    250     scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    251     background = values[1] 
    252     return scale*total + background 
     226        total = form() 
     227        weight_norm = 1.0 
     228        weighted_volume = form_volume() 
     229        weighted_radius = form_radius() 
     230 
     231    else: 
     232        pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] 
     233        pd_weight = values[2+n_pars + call_details.num_weights:] 
     234 
     235        weight_norm = 0.0 
     236        weighted_volume = 0.0 
     237        weighted_radius = 0.0 
     238        partial_weight = np.NaN 
     239        weight = np.NaN 
     240 
     241        p0_par = call_details.pd_par[0] 
     242        p0_length = call_details.pd_length[0] 
     243        p0_index = p0_length 
     244        p0_offset = call_details.pd_offset[0] 
     245 
     246        pd_par = call_details.pd_par[:call_details.num_active] 
     247        pd_offset = call_details.pd_offset[:call_details.num_active] 
     248        pd_stride = call_details.pd_stride[:call_details.num_active] 
     249        pd_length = call_details.pd_length[:call_details.num_active] 
     250 
     251        total = np.zeros(nq, 'd') 
     252        for loop_index in range(call_details.num_eval): 
     253            # update polydispersity parameter values 
     254            if p0_index == p0_length: 
     255                pd_index = (loop_index//pd_stride)%pd_length 
     256                parameters[pd_par] = pd_value[pd_offset+pd_index] 
     257                partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
     258                p0_index = loop_index%p0_length 
     259 
     260            weight = partial_weight * pd_weight[p0_offset + p0_index] 
     261            parameters[p0_par] = pd_value[p0_offset + p0_index] 
     262            p0_index += 1 
     263            if weight > cutoff: 
     264                # Call the scattering function 
     265                # Assume that NaNs are only generated if the parameters are bad; 
     266                # exclude all q for that NaN.  Even better would be to have an 
     267                # INVALID expression like the C models, but that is expensive. 
     268                Iq = np.asarray(form(), 'd') 
     269                if np.isnan(Iq).any(): 
     270                    continue 
     271 
     272                # update value and norm 
     273                total += weight * Iq 
     274                weight_norm += weight 
     275                weighted_volume += weight * form_volume() 
     276                weighted_radius += weight * form_radius() 
     277 
     278    result = np.hstack((total, weight_norm, weighted_volume, weighted_radius)) 
     279    return result 
    253280 
    254281 

  • sasmodels/modelinfo.py

    r7b9e4dd r5399809  
    790790    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
    791791    # TODO: find Fq by inspection 
     792    info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None) 
    792793    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    793794    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    794795    info.source = getattr(kernel_module, 'source', []) 
    795796    info.c_code = getattr(kernel_module, 'c_code', None) 
     797    info.effective_radius = getattr(kernel_module, 'effective_radius', None) 
     798    info.ER = None  # CRUFT 
     799    info.VR = None  # CRUFT 
    796800    # TODO: check the structure of the tests 
    797801    info.tests = getattr(kernel_module, 'tests', []) 
    798     info.ER = getattr(kernel_module, 'ER', None) # type: ignore 
    799     info.VR = getattr(kernel_module, 'VR', None) # type: ignore 
    800802    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
    801803    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    … …  
    923925    #: use those functions.  Form factors are indicated by providing 
    924926    #: an :attr:`ER` function. 
     927    effective_radius_type = None   # type: List[str] 
     928    #: Returns the occupied volume and the total volume for each parameter set. 
     929    #: See :attr:`ER` for details on the parameters. 
    925930    source = None           # type: List[str] 
    926931    #: The set of tests that must pass.  The format of the tests is described 
    … …  
    943948    #: each parameter set.  Multiplicity parameters will be received as 
    944949    #: arrays, with one row per polydispersity condition. 
    945     ER = None               # type: Optional[Callable[[np.ndarray], np.ndarray]] 
    946     #: Returns the occupied volume and the total volume for each parameter set. 
    947     #: See :attr:`ER` for details on the parameters. 
    948     VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 
    949     #: Arbitrary C code containing supporting functions, etc., to be inserted 
    950     #: after everything in source.  This can include Iq and Iqxy functions with 
    951     #: the full function signature, including all parameters. 
    952950    c_code = None 
    953951    #: Returns the form volume for python-based models.  Form volume is needed 

  • sasmodels/models/barbell.c

    r71b751d ree60aa7  
    6262} 
    6363 
     64static double 
     65radius_from_volume(double radius_bell, double radius, double length) 
     66{ 
     67    const double vol_barbell = form_volume(radius_bell,radius,length); 
     68    return cbrt(0.75*vol_barbell/M_PI); 
     69} 
     70 
     71static double 
     72radius_from_totallength(double radius_bell, double radius, double length) 
     73{ 
     74    const double hdist = sqrt(square(radius_bell) - square(radius)); 
     75    return 0.5*length + hdist + radius_bell; 
     76} 
     77 
     78static double 
     79effective_radius(int mode, double radius_bell, double radius, double length) 
     80{ 
     81    switch (mode) { 
     82    case 1: // equivalent sphere 
     83        return radius_from_volume(radius_bell, radius , length); 
     84    case 2: // radius 
     85        return radius; 
     86    case 3: // half length 
     87        return 0.5*length; 
     88    case 4: // half total length 
     89        return radius_from_totallength(radius_bell,radius,length); 
     90    } 
     91} 
     92 
    6493static void 
    6594Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 

  • sasmodels/models/barbell.py

    r71b751d ree60aa7  
    117117source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 
    118118have_Fq = True 
     119effective_radius_type = [ 
     120    "equivalent sphere", "radius", "half length", "half total length", 
     121    ] 
    119122 
    120123def random(): 

  • sasmodels/models/capped_cylinder.c

    r71b751d ree60aa7  
    8484} 
    8585 
     86static double 
     87radius_from_volume(double radius, double radius_cap, double length) 
     88{ 
     89    const double vol_cappedcyl = form_volume(radius,radius_cap,length); 
     90    return cbrt(0.75*vol_cappedcyl/M_PI); 
     91} 
     92 
     93static double 
     94radius_from_totallength(double radius, double radius_cap, double length) 
     95{ 
     96    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); 
     97    return 0.5*length + hc; 
     98} 
     99 
     100static double 
     101effective_radius(int mode, double radius, double radius_cap, double length) 
     102{ 
     103    switch (mode) { 
     104    case 1: // equivalent sphere 
     105        return radius_from_volume(radius, radius_cap, length); 
     106    case 2: // radius 
     107        return radius; 
     108    case 3: // half length 
     109        return 0.5*length; 
     110    case 4: // half total length 
     111        return radius_from_totallength(radius, radius_cap,length); 
     112    } 
     113} 
     114 
    86115static void 
    87116Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 

  • sasmodels/models/capped_cylinder.py

    r71b751d ree60aa7  
    137137source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
    138138have_Fq = True 
     139effective_radius_type = [ 
     140    "equivalent sphere", "radius", "half length", "half total length", 
     141    ] 
    139142 
    140143def random(): 

  • sasmodels/models/core_multi_shell.c

    r71b751d ree60aa7  
    99 
    1010static double 
    11 form_volume(double core_radius, double fp_n, double thickness[]) 
     11outer_radius(double core_radius, double fp_n, double thickness[]) 
    1212{ 
    1313  double r = core_radius; 
    … …  
    1616    r += thickness[i]; 
    1717  } 
    18   return M_4PI_3 * cube(r); 
     18  return r; 
     19} 
     20 
     21static double 
     22form_volume(double core_radius, double fp_n, double thickness[]) 
     23{ 
     24  return M_4PI_3 * cube(outer_radius(core_radius, fp_n, thickness)); 
     25} 
     26 
     27static double 
     28effective_radius(int mode, double core_radius, double fp_n, double thickness[]) 
     29{ 
     30  switch (mode) { 
     31  case 1: // outer radius 
     32    return outer_radius(core_radius, fp_n, thickness); 
     33  case 2: // core radius 
     34    return core_radius; 
     35  } 
    1936} 
    2037 

  • sasmodels/models/core_multi_shell.py

    r71b751d ree60aa7  
    100100source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 
    101101have_Fq = True 
     102effective_radius_type = ["outer radius", "core radius"] 
    102103 
    103104def random(): 
    … …  
    144145    return np.asarray(z), np.asarray(rho) 
    145146 
    146 def ER(radius, n, thickness): 
    147     """Effective radius""" 
    148     n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse 
    149     return np.sum(thickness[:n], axis=0) + radius 
    150  
    151147demo = dict(sld_core=6.4, 
    152148            radius=60, 

  • sasmodels/models/core_shell_bicelle.c

    r71b751d ree60aa7  
    3434 
    3535    return t; 
     36} 
     37 
     38static double 
     39radius_from_volume(double radius, double thick_rim, double thick_face, double length) 
     40{ 
     41    const double volume_bicelle = form_volume(radius,thick_rim,thick_face,length); 
     42    return cbrt(0.75*volume_bicelle/M_PI); 
     43} 
     44 
     45static double 
     46radius_from_diagonal(double radius, double thick_rim, double thick_face, double length) 
     47{ 
     48    const double radius_tot = radius + thick_rim; 
     49    const double length_tot = length + 2.0*thick_face; 
     50    return sqrt(radius_tot*radius_tot + 0.25*length_tot*length_tot); 
     51} 
     52 
     53static double 
     54effective_radius(int mode, double radius, double thick_rim, double thick_face, double length) 
     55{ 
     56    switch (mode) { 
     57    case 1: // equivalent sphere 
     58        return radius_from_volume(radius, thick_rim, thick_face, length); 
     59    case 2: // outer rim radius 
     60        return radius + thick_rim; 
     61    case 3: // half outer thickness 
     62        return 0.5*length + thick_face; 
     63    case 4: // half diagonal 
     64        return radius_from_diagonal(radius,thick_rim,thick_face,length); 
     65    } 
    3666} 
    3767 

  • sasmodels/models/core_shell_bicelle.py

    r71b751d ree60aa7  
    155155          "core_shell_bicelle.c"] 
    156156have_Fq = True 
     157effective_radius_type = [ 
     158    "equivalent sphere", "outer rim radius", 
     159    "half outer thickness", "half diagonal", 
     160    ] 
    157161 
    158162def random(): 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r71b751d ree60aa7  
    88{ 
    99    return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face); 
     10} 
     11 
     12static double 
     13radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     14{ 
     15    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length); 
     16    return cbrt(0.75*volume_bicelle/M_PI); 
     17} 
     18 
     19static double 
     20radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     21{ 
     22    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor); 
     23    const double radius_max_tot = radius_max + thick_rim; 
     24    const double length_tot = length + 2.0*thick_face; 
     25    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot); 
     26} 
     27 
     28static double 
     29effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     30{ 
     31    switch (mode) { 
     32    case 1: // equivalent sphere 
     33        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     34    case 2: // outer rim average radius 
     35        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     36    case 3: // outer rim min radius 
     37        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     38    case 4: // outer max radius 
     39        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     40    case 5: // half outer thickness 
     41        return 0.5*length + thick_face; 
     42    case 6: // half diagonal 
     43        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     44    } 
    1045} 
    1146 

  • sasmodels/models/core_shell_bicelle_elliptical.py

    r71b751d ree60aa7  
    147147          "core_shell_bicelle_elliptical.c"] 
    148148have_Fq = True 
     149effective_radius_type = [ 
     150    "equivalent sphere", "outer rim average radius", "outer rim min radius", 
     151    "outer max radius", "half outer thickness", "half diagonal", 
     152    ] 
    149153 
    150154def random(): 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r71b751d rd299327  
    99    return M_PI*(  (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 
    1010                 square(r_minor)*x_core*2.0*thick_face  ); 
     11} 
     12 
     13static double 
     14radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     15{ 
     16    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length); 
     17    return cbrt(0.75*volume_bicelle/M_PI); 
     18} 
     19 
     20static double 
     21radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     22{ 
     23    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor); 
     24    const double radius_max_tot = radius_max + thick_rim; 
     25    const double length_tot = length + 2.0*thick_face; 
     26    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot); 
     27} 
     28 
     29static double 
     30effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     31{ 
     32    switch (mode) { 
     33    case 1: // equivalent sphere 
     34        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     35    case 2: // outer rim average radius 
     36        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     37    case 3: // outer rim min radius 
     38        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     39    case 4: // outer max radius 
     40        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     41    case 5: // half outer thickness 
     42        return 0.5*length + thick_face; 
     43    case 6: // half diagonal (this ignores the missing "corners", so may give unexpected answer if thick_face 
     44            //  or thick_rim is extremely large) 
     45        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     46    } 
    1147} 
    1248 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    r71b751d ree60aa7  
    160160          "core_shell_bicelle_elliptical_belt_rough.c"] 
    161161have_Fq = True 
     162effective_radius_type = [ 
     163    "equivalent sphere", "outer rim average radius", "outer rim min radius", 
     164    "outer max radius", "half outer thickness", "half diagonal", 
     165    ] 
    162166 
    163167demo = dict(scale=1, background=0, 

  • sasmodels/models/core_shell_cylinder.c

    r71b751d ree60aa7  
    1111{ 
    1212    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
     13} 
     14 
     15static double 
     16radius_from_volume(double radius, double thickness, double length) 
     17{ 
     18    const double volume_outer_cyl = form_volume(radius,thickness,length); 
     19    return cbrt(0.75*volume_outer_cyl/M_PI); 
     20} 
     21 
     22static double 
     23radius_from_diagonal(double radius, double thickness, double length) 
     24{ 
     25    const double radius_outer = radius + thickness; 
     26    const double length_outer = length + 2.0*thickness; 
     27    return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer); 
     28} 
     29 
     30static double 
     31effective_radius(int mode, double radius, double thickness, double length) 
     32{ 
     33    switch (mode) { 
     34    case 1: // equivalent sphere 
     35        return radius_from_volume(radius, thickness, length); 
     36    case 2: // outer radius 
     37        return radius + thickness; 
     38    case 3: // half outer length 
     39        return 0.5*length + thickness; 
     40    case 4: // half min outer length 
     41        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness); 
     42    case 5: // half max outer length 
     43        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness); 
     44    case 6: // half outer diagonal 
     45        return radius_from_diagonal(radius,thickness,length); 
     46    } 
    1347} 
    1448 

  • sasmodels/models/core_shell_cylinder.py

    re31b19a ree60aa7  
    132132source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 
    133133have_Fq = True 
    134  
    135 def ER(radius, thickness, length): 
    136     """ 
    137     Returns the effective radius used in the S*P calculation 
    138     """ 
    139     radius = radius + thickness 
    140     length = length + 2 * thickness 
    141     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    142     return 0.5 * (ddd) ** (1. / 3.) 
     134effective_radius_type = [ 
     135    "equivalent sphere", "outer radius", "half outer length", 
     136    "half min outer dimension", "half max outer dimension", 
     137    "half outer diagonal", 
     138    ] 
    143139 
    144140def random(): 

  • sasmodels/models/core_shell_ellipsoid.c

    r71b751d ree60aa7  
    3636    double vol = M_4PI_3*equat_shell*equat_shell*polar_shell; 
    3737    return vol; 
     38} 
     39 
     40static double 
     41radius_from_volume(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     42{ 
     43    const double volume_ellipsoid = form_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     44    return cbrt(0.75*volume_ellipsoid/M_PI); 
     45} 
     46 
     47static double 
     48radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     49{ 
     50    // Trivial cases 
     51    if (1.0 == x_core && 1.0 == x_polar_shell) return radius_equat_core + thick_shell; 
     52    if ((radius_equat_core + thick_shell)*(radius_equat_core*x_core + thick_shell*x_polar_shell) == 0.)  return 0.; 
     53 
     54    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     55    const double radius_equat_tot = radius_equat_core + thick_shell; 
     56    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 
     57    const double ratio = (radius_polar_tot < radius_equat_tot 
     58                          ? radius_polar_tot / radius_equat_tot 
     59                          : radius_equat_tot / radius_polar_tot); 
     60    const double e1 = sqrt(1.0 - ratio*ratio); 
     61    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     62    const double bL = (1.0 + e1) / (1.0 - e1); 
     63    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     64    const double delta = 0.75 * b1 * b2; 
     65    const double ddd = 2.0 * (delta + 1.0) * radius_polar_tot * radius_equat_tot * radius_equat_tot; 
     66    return 0.5 * cbrt(ddd); 
     67} 
     68 
     69static double 
     70effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     71{ 
     72    const double radius_equat_tot = radius_equat_core + thick_shell; 
     73    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 
     74 
     75    switch (mode) { 
     76    case 1: // equivalent sphere 
     77        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     78    case 2: // average outer curvature 
     79        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     80    case 3: // min outer radius 
     81        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     82    case 4: // max outer radius 
     83        return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     84    } 
    3885} 
    3986 

  • sasmodels/models/core_shell_ellipsoid.py

    r71b751d ree60aa7  
    146146source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    147147have_Fq = True 
    148  
    149 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 
    150     """ 
    151         Returns the effective radius used in the S*P calculation 
    152     """ 
    153     from .ellipsoid import ER as ellipsoid_ER 
    154     polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell 
    155     equat_outer = radius_equat_core + thick_shell 
    156     return ellipsoid_ER(polar_outer, equat_outer) 
     148effective_radius_type = [ 
     149    "equivalent sphere", "average outer curvature", 
     150    "min outer radius", "max outer radius", 
     151    ] 
    157152 
    158153def random(): 

  • sasmodels/models/core_shell_parallelepiped.c

    r71b751d ree60aa7  
    2525        2.0 * length_a * length_b * thick_rim_c; 
    2626#endif 
     27} 
     28 
     29static double 
     30radius_from_volume(double length_a, double length_b, double length_c, 
     31                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     32{ 
     33    const double volume = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     34    return cbrt(volume/M_4PI_3); 
     35} 
     36 
     37static double 
     38radius_from_crosssection(double length_a, double length_b, double thick_rim_a, double thick_rim_b) 
     39{ 
     40    const double area_xsec_paral = length_a*length_b + 2.0*thick_rim_a*length_b + 2.0*thick_rim_b*length_a; 
     41    return sqrt(area_xsec_paral/M_PI); 
     42} 
     43 
     44static double 
     45effective_radius(int mode, double length_a, double length_b, double length_c, 
     46                 double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     47{ 
     48    switch (mode) { 
     49    case 1: // equivalent sphere 
     50        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     51    case 2: // half outer length a 
     52        return 0.5 * length_a + thick_rim_a; 
     53    case 3: // half outer length b 
     54        return 0.5 * length_b + thick_rim_b; 
     55    case 4: // half outer length c 
     56        return 0.5 * length_c + thick_rim_c; 
     57    case 5: // equivalent circular cross-section 
     58        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b); 
     59    case 6: // half outer ab diagonal 
     60        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b)); 
     61    case 7: // half outer diagonal 
     62        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*thick_rim_c)); 
     63    } 
    2764} 
    2865 

  • sasmodels/models/core_shell_parallelepiped.py

    r71b751d ree60aa7  
    9898is the scattering length of the solvent. 
    9999 
    100 .. note::  
     100.. note:: 
    101101 
    102102   the code actually implements two substitutions: $d(cos\alpha)$ is 
    … …  
    227227source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] 
    228228have_Fq = True 
    229  
    230  
    231 def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c): 
    232     """ 
    233         Return equivalent radius (ER) 
    234     """ 
    235     from .parallelepiped import ER as ER_p 
    236  
    237     a = length_a + 2*thick_rim_a 
    238     b = length_b + 2*thick_rim_b 
    239     c = length_c + 2*thick_rim_c 
    240     return ER_p(a, b, c) 
    241  
    242 # VR defaults to 1.0 
     229effective_radius_type = [ 
     230    "equivalent sphere", 
     231    "half outer length_a", "half outer length_b", "half outer length_c", 
     232    "equivalent circular cross-section", 
     233    "half outer ab diagonal", "half outer diagonal", 
     234    ] 
    243235 
    244236def random(): 

  • sasmodels/models/core_shell_sphere.c

    r71b751d ree60aa7  
    33{ 
    44    return M_4PI_3 * cube(radius + thickness); 
     5} 
     6 
     7static double 
     8effective_radius(int mode, double radius, double thickness) 
     9{ 
     10    switch (mode) { 
     11    case 1: // outer radius 
     12        return radius + thickness; 
     13    case 2: // core radius 
     14        return radius; 
     15    } 
    516} 
    617 

  • sasmodels/models/core_shell_sphere.py

    rda1c8d1 ree60aa7  
    7878source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 
    7979have_Fq = True 
     80effective_radius_type = ["outer radius", "core radius"] 
    8081 
    8182demo = dict(scale=1, background=0, radius=60, thickness=10, 
    8283            sld_core=1.0, sld_shell=2.0, sld_solvent=0.0) 
    83  
    84 def ER(radius, thickness): 
    85     """ 
    86         Equivalent radius 
    87         @param radius: core radius 
    88         @param thickness: shell thickness 
    89     """ 
    90     return radius + thickness 
    9184 
    9285def random(): 
    … …  
    10396 
    10497tests = [ 
    105     [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
     98#    [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    10699 
    107100    # The SasView test result was 0.00169, with a background of 0.001 

  • sasmodels/models/cylinder.c

    r71b751d ree60aa7  
    1111{ 
    1212    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
     13} 
     14 
     15static double 
     16radius_from_volume(double radius, double length) 
     17{ 
     18    return cbrt(0.75*radius*radius*length); 
     19} 
     20 
     21static double 
     22radius_from_diagonal(double radius, double length) 
     23{ 
     24    return sqrt(radius*radius + 0.25*length*length); 
     25} 
     26 
     27static double 
     28effective_radius(int mode, double radius, double length) 
     29{ 
     30    switch (mode) { 
     31    case 1: 
     32        return radius_from_volume(radius, length); 
     33    case 2: 
     34        return radius; 
     35    case 3: 
     36        return 0.5*length; 
     37    case 4: 
     38        return (radius < 0.5*length ? radius : 0.5*length); 
     39    case 5: 
     40        return (radius > 0.5*length ? radius : 0.5*length); 
     41    case 6: 
     42        return radius_from_diagonal(radius,length); 
     43    } 
    1344} 
    1445 

  • sasmodels/models/cylinder.py

    r71b751d ree60aa7  
    139139source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
    140140have_Fq = True 
    141  
    142 def ER(radius, length): 
    143     """ 
    144         Return equivalent radius (ER) 
    145     """ 
    146     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    147     return 0.5 * (ddd) ** (1. / 3.) 
     141effective_radius_type = [ 
     142    "equivalent sphere", "radius", 
     143    "half length", "half min dimension", "half max dimension", "half diagonal", 
     144    ] 
    148145 
    149146def random(): 

  • sasmodels/models/ellipsoid.c

    r71b751d ree60aa7  
    44    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    55} 
     6 
     7static double 
     8radius_from_volume(double radius_polar, double radius_equatorial) 
     9{ 
     10    return cbrt(radius_polar*radius_equatorial*radius_equatorial); 
     11} 
     12 
     13static double 
     14radius_from_curvature(double radius_polar, double radius_equatorial) 
     15{ 
     16    // Trivial cases 
     17    if (radius_polar == radius_equatorial) return radius_polar; 
     18    if (radius_polar * radius_equatorial == 0.)  return 0.; 
     19 
     20    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     21    const double ratio = (radius_polar < radius_equatorial 
     22                          ? radius_polar / radius_equatorial 
     23                          : radius_equatorial / radius_polar); 
     24    const double e1 = sqrt(1.0 - ratio*ratio); 
     25    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     26    const double bL = (1.0 + e1) / (1.0 - e1); 
     27    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     28    const double delta = 0.75 * b1 * b2; 
     29    const double ddd = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial * radius_equatorial; 
     30    return 0.5 * cbrt(ddd); 
     31} 
     32 
     33static double 
     34effective_radius(int mode, double radius_polar, double radius_equatorial) 
     35{ 
     36    switch (mode) { 
     37    case 1: // equivalent sphere 
     38        return radius_from_volume(radius_polar, radius_equatorial); 
     39    case 2: // average curvature 
     40        return radius_from_curvature(radius_polar, radius_equatorial); 
     41    case 3: // min radius 
     42        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial); 
     43    case 4: // max radius 
     44        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial); 
     45    } 
     46} 
     47 
    648 
    749static void 

  • sasmodels/models/ellipsoid.py

    r0168844 ree60aa7  
    169169source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
    170170have_Fq = True 
    171  
    172 def ER(radius_polar, radius_equatorial): 
    173     # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
    174     ee = np.empty_like(radius_polar) 
    175     idx = radius_polar > radius_equatorial 
    176     ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2 
    177     idx = radius_polar < radius_equatorial 
    178     ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2 
    179     valid = (radius_polar * radius_equatorial != 0) & (radius_polar != radius_equatorial) 
    180     bd = 1.0 - ee[valid] 
    181     e1 = np.sqrt(ee[valid]) 
    182     b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd)) 
    183     bL = (1.0 + e1) / (1.0 - e1) 
    184     b2 = 1.0 + bd / 2 / e1 * np.log(bL) 
    185     delta = 0.75 * b1 * b2 
    186     ddd = 2.0 * (delta + 1.0) * (radius_polar * radius_equatorial**2)[valid] 
    187  
    188     r = np.zeros_like(radius_polar) 
    189     r[valid] = 0.5 * cbrt(ddd) 
    190     idx = radius_polar == radius_equatorial 
    191     r[idx] = radius_polar[idx] 
    192     return r 
     171effective_radius_type = [ 
     172    "equivalent sphere", "average curvature", "min radius", "max radius", 
     173    ] 
    193174 
    194175def random(): 

  • sasmodels/models/elliptical_cylinder.c

    r71b751d ree60aa7  
    33{ 
    44    return M_PI * radius_minor * radius_minor * r_ratio * length; 
     5} 
     6 
     7static double 
     8radius_from_volume(double radius_minor, double r_ratio, double length) 
     9{ 
     10    const double volume_ellcyl = form_volume(radius_minor,r_ratio,length); 
     11    return cbrt(0.75*volume_ellcyl/M_PI); 
     12} 
     13 
     14static double 
     15radius_from_min_dimension(double radius_minor, double r_ratio, double hlength) 
     16{ 
     17    const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     18    return (rad_min < hlength ? rad_min : hlength); 
     19} 
     20 
     21static double 
     22radius_from_max_dimension(double radius_minor, double r_ratio, double hlength) 
     23{ 
     24    const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     25    return (rad_max > hlength ? rad_max : hlength); 
     26} 
     27 
     28static double 
     29radius_from_diagonal(double radius_minor, double r_ratio, double length) 
     30{ 
     31    const double radius_max = (r_ratio > 1.0 ? radius_minor*r_ratio : radius_minor); 
     32    return sqrt(radius_max*radius_max + 0.25*length*length); 
     33} 
     34 
     35static double 
     36effective_radius(int mode, double radius_minor, double r_ratio, double length) 
     37{ 
     38    switch (mode) { 
     39    case 1: // equivalent sphere 
     40        return radius_from_volume(radius_minor, r_ratio, length); 
     41    case 2: // average radius 
     42        return 0.5*radius_minor*(1.0 + r_ratio); 
     43    case 3: // min radius 
     44        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     45    case 4: // max radius 
     46        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     47    case 5: // equivalent circular cross-section 
     48        return sqrt(radius_minor*radius_minor*r_ratio); 
     49    case 6: // half length 
     50        return 0.5*length; 
     51    case 7: // half min dimension 
     52        return radius_from_min_dimension(radius_minor,r_ratio,0.5*length); 
     53    case 8: // half max dimension 
     54        return radius_from_max_dimension(radius_minor,r_ratio,0.5*length); 
     55    case 9: // half diagonal 
     56        return radius_from_diagonal(radius_minor,r_ratio,length); 
     57    } 
    558} 
    659 

  • sasmodels/models/elliptical_cylinder.py

    r71b751d ree60aa7  
    123123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
    124124have_Fq = True 
     125effective_radius_type = [ 
     126    "equivalent sphere", "average radius", "min radius", "max radius", 
     127    "equivalent circular cross-section", 
     128    "half length", "half min dimension", "half max dimension", "half diagonal", 
     129    ] 
    125130 
    126131demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 
    127132            sld=4.0, sld_solvent=1.0, theta=10.0, phi=20, psi=30, 
    128133            theta_pd=10, phi_pd=2, psi_pd=3) 
    129  
    130 def ER(radius_minor, axis_ratio, length): 
    131     """ 
    132         Equivalent radius 
    133         @param radius_minor: Ellipse minor radius 
    134         @param axis_ratio: Ratio of major radius over minor radius 
    135         @param length: Length of the cylinder 
    136     """ 
    137     radius = sqrt(radius_minor * radius_minor * axis_ratio) 
    138     ddd = 0.75 * radius * (2 * radius * length 
    139                            + (length + radius) * (length + pi * radius)) 
    140     return 0.5 * (ddd) ** (1. / 3.) 
    141134 
    142135def random(): 
    … …  
    162155 
    163156tests = [ 
    164     [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], 
    165     [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], 
     157#    [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], 
     158#    [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], 
    166159 
    167160    # The SasView test result was 0.00169, with a background of 0.001 

  • sasmodels/models/fractal_core_shell.py

    reb3eb38 r2cc8aa2  
    134134    return radius + thickness 
    135135 
    136 tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    137  
     136#tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
     137tests = [ 
    138138#         # At some point the SasView 3.x test result was deemed incorrect. The 
    139139          #following tests were verified against NIST IGOR macros ver 7.850. 

  • sasmodels/models/fuzzy_sphere.py

    r71b751d ree60aa7  
    7878              ["sld_solvent", "1e-6/Ang^2",  3, [-inf, inf], "sld",    "Solvent scattering length density"], 
    7979              ["radius",      "Ang",        60, [0, inf],    "volume", "Sphere radius"], 
    80               ["fuzziness",   "Ang",        10, [0, inf],    "",       "std deviation of Gaussian convolution for interface (must be << radius)"], 
     80              ["fuzziness",   "Ang",        10, [0, inf],    "volume",       "std deviation of Gaussian convolution for interface (must be << radius)"], 
    8181             ] 
    8282# pylint: enable=bad-whitespace,line-too-long 
    8383 
    84 source = ["lib/sas_3j1x_x.c"] 
     84source = ["lib/sas_3j1x_x.c","fuzzy_sphere.c"] 
    8585have_Fq = True 
    86  
    87 c_code = """ 
    88 static double form_volume(double radius) 
    89 { 
    90     return M_4PI_3*cube(radius); 
    91 } 
    92  
    93 static void Fq(double q, double *F1, double *F2, double sld, double sld_solvent, 
    94                double radius, double fuzziness) 
    95 { 
    96     const double qr = q*radius; 
    97     const double bes = sas_3j1x_x(qr); 
    98     const double qf = exp(-0.5*square(q*fuzziness)); 
    99     const double contrast = (sld - sld_solvent); 
    100     const double form = contrast * form_volume(radius) * bes * qf; 
    101     *F1 = 1.0e-2*form; 
    102     *F2 = 1.0e-4*form*form; 
    103 } 
    104 """ 
    105  
    106 def ER(radius): 
    107     """ 
    108     Return radius 
    109     """ 
    110     return radius 
    111  
    112 # VR defaults to 1.0 
     86effective_radius_type = ["radius", "radius + fuzziness"] 
    11387 
    11488def random(): 

  • sasmodels/models/hollow_cylinder.c

    r71b751d ree60aa7  
    1515} 
    1616 
     17// TODO: interface to form_volume/shell_volume not yet settled 
     18static double 
     19shell_volume(double *total, double radius, double thickness, double length) 
     20{ 
     21    *total = M_PI*length*square(radius+thickness); 
     22    return *total - M_PI*length*radius*radius; 
     23} 
     24 
    1725static double 
    1826form_volume(double radius, double thickness, double length) 
    1927{ 
    20     double v_shell = M_PI*length*(square(radius+thickness) - radius*radius); 
    21     return v_shell; 
     28    double total; 
     29    return shell_volume(&total, radius, thickness, length); 
    2230} 
    2331 
     32static double 
     33radius_from_volume(double radius, double thickness, double length) 
     34{ 
     35    const double volume_outer_cyl = M_PI*square(radius + thickness)*length; 
     36    return cbrt(0.75*volume_outer_cyl/M_PI); 
     37} 
     38 
     39static double 
     40radius_from_diagonal(double radius, double thickness, double length) 
     41{ 
     42    return sqrt(square(radius + thickness) + 0.25*square(length)); 
     43} 
     44 
     45static double 
     46effective_radius(int mode, double radius, double thickness, double length) 
     47{ 
     48    switch (mode) { 
     49    case 1: // equivalent sphere 
     50        return radius_from_volume(radius, thickness, length); 
     51    case 2: // outer radius 
     52        return radius + thickness; 
     53    case 3: // half length 
     54        return 0.5*length; 
     55    case 4: // half outer min dimension 
     56        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length); 
     57    case 5: // half outer max dimension 
     58        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length); 
     59    case 6: // half outer diagonal 
     60        return radius_from_diagonal(radius,thickness,length); 
     61    } 
     62} 
    2463 
    2564static void 

  • sasmodels/models/hollow_cylinder.py

    r455aaa1 ree60aa7  
    100100source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    101101have_Fq = True 
    102  
    103 # pylint: disable=W0613 
    104 def ER(radius, thickness, length): 
    105     """ 
    106     :param radius:      Cylinder core radius 
    107     :param thickness:   Cylinder wall thickness 
    108     :param length:      Cylinder length 
    109     :return:            Effective radius 
    110     """ 
    111     router = radius + thickness 
    112     if router == 0 or length == 0: 
    113         return 0.0 
    114     len1 = router 
    115     len2 = length/2.0 
    116     term1 = len1*len1*2.0*len2/2.0 
    117     term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0) 
    118     ddd = 3.0*term1*term2 
    119     diam = pow(ddd, (1.0/3.0)) 
    120     return diam 
    121  
    122 def VR(radius, thickness, length): 
    123     """ 
    124     :param radius:      Cylinder radius 
    125     :param thickness:   Cylinder wall thickness 
    126     :param length:      Cylinder length 
    127     :return:            Volf ratio for P(q)*S(q) 
    128     """ 
    129     router = radius + thickness 
    130     vol_core = pi*radius*radius*length 
    131     vol_total = pi*router*router*length 
    132     vol_shell = vol_total - vol_core 
    133     return vol_total, vol_shell 
     102effective_radius_type = [ 
     103    "equivalent sphere", "outer radius", "half length", 
     104    "half outer min dimension", "half outer max dimension", 
     105    "half outer diagonal", 
     106    ] 
    134107 
    135108def random(): 
    … …  
    162135tests = [ 
    163136    [{}, 0.00005, 1764.926], 
    164     [{}, 'VR', 0.55555556], 
     137#    [{}, 'VR', 0.55555556], 
    165138    [{}, 0.001, 1756.76], 
    166139    [{}, (qx, qy), 2.36885476192], 

  • sasmodels/models/hollow_rectangular_prism.c

    r71b751d ree60aa7  
     1// TODO: interface to form_volume/shell_volume not yet settled 
    12static double 
    2 form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     3shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    34{ 
    45    double length_b = length_a * b2a_ratio; 
    … …  
    89    double c_core = length_c - 2.0*thickness; 
    910    double vol_core = a_core * b_core * c_core; 
    10     double vol_total = length_a * length_b * length_c; 
    11     double vol_shell = vol_total - vol_core; 
    12     return vol_shell; 
     11    *total = length_a * length_b * length_c; 
     12    return *total - vol_core; 
     13} 
     14 
     15static double 
     16form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     17{ 
     18    double total; 
     19    return shell_volume(&total, length_a, b2a_ratio, c2a_ratio, thickness); 
     20} 
     21 
     22static double 
     23effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     24//effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
     25//                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 
     26// NOTE length_a is external dimension! 
     27{ 
     28    switch (mode) { 
     29    case 1: // equivalent sphere 
     30        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
     31    case 2: // half length_a 
     32        return 0.5 * length_a; 
     33    case 3: // half length_b 
     34        return 0.5 * length_a*b2a_ratio; 
     35    case 4: // half length_c 
     36        return 0.5 * length_a*c2a_ratio; 
     37    case 5: // equivalent outer circular cross-section 
     38        return length_a*sqrt(b2a_ratio/M_PI); 
     39    case 6: // half ab diagonal 
     40        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     41    case 7: // half diagonal 
     42        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     43    } 
    1344} 
    1445 

  • sasmodels/models/hollow_rectangular_prism.py

    r455aaa1 ree60aa7  
    132132               "Solvent scattering length density"], 
    133133              ["length_a", "Ang", 35, [0, inf], "volume", 
    134                "Shorter side of the parallelepiped"], 
     134               "Shortest, external, size of the parallelepiped"], 
    135135              ["b2a_ratio", "Ang", 1, [0, inf], "volume", 
    136136               "Ratio sides b/a"], 
    … …  
    149149source = ["lib/gauss76.c", "hollow_rectangular_prism.c"] 
    150150have_Fq = True 
    151  
    152 def ER(length_a, b2a_ratio, c2a_ratio, thickness): 
    153     """ 
    154     Return equivalent radius (ER) 
    155     thickness parameter not used 
    156     """ 
    157     b_side = length_a * b2a_ratio 
    158     c_side = length_a * c2a_ratio 
    159  
    160     # surface average radius (rough approximation) 
    161     surf_rad = sqrt(length_a * b_side / pi) 
    162  
    163     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    164     return 0.5 * (ddd) ** (1. / 3.) 
    165  
    166 def VR(length_a, b2a_ratio, c2a_ratio, thickness): 
    167     """ 
    168     Return shell volume and total volume 
    169     """ 
    170     b_side = length_a * b2a_ratio 
    171     c_side = length_a * c2a_ratio 
    172     a_core = length_a - 2.0*thickness 
    173     b_core = b_side - 2.0*thickness 
    174     c_core = c_side - 2.0*thickness 
    175     vol_core = a_core * b_core * c_core 
    176     vol_total = length_a * b_side * c_side 
    177     vol_shell = vol_total - vol_core 
    178     return vol_total, vol_shell 
    179  
     151effective_radius_type = [ 
     152    "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     153    "equivalent outer circular cross-section", 
     154    "half ab diagonal", "half diagonal", 
     155    ] 
    180156 
    181157def random(): 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    r71b751d ree60aa7  
     1// TODO: interface to form_volume/shell_volume not yet settled 
    12static double 
    2 form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     3shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio) 
    34{ 
    45    double length_b = length_a * b2a_ratio; 
    56    double length_c = length_a * c2a_ratio; 
    67    double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
     8    *total = length_a * length_b * length_c; 
    79    return vol_shell; 
     10} 
     11 
     12static double 
     13form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     14{ 
     15    double total; 
     16    return shell_volume(&total, length_a, b2a_ratio, c2a_ratio); 
     17} 
     18 
     19 
     20static double 
     21effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
     22{ 
     23    switch (mode) { 
     24    case 1: // equivalent sphere 
     25        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
     26    case 2: // half length_a 
     27        return 0.5 * length_a; 
     28    case 3: // half length_b 
     29        return 0.5 * length_a*b2a_ratio; 
     30    case 4: // half length_c 
     31        return 0.5 * length_a*c2a_ratio; 
     32    case 5: // equivalent outer circular cross-section 
     33        return length_a*sqrt(b2a_ratio/M_PI); 
     34    case 6: // half ab diagonal 
     35        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     36    case 7: // half diagonal 
     37        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     38    } 
    839} 
    940 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    r6e7d7b6 ree60aa7  
    109109source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"] 
    110110have_Fq = True 
    111  
    112 def ER(length_a, b2a_ratio, c2a_ratio): 
    113     """ 
    114         Return equivalent radius (ER) 
    115     """ 
    116     b_side = length_a * b2a_ratio 
    117     c_side = length_a * c2a_ratio 
    118  
    119     # surface average radius (rough approximation) 
    120     surf_rad = sqrt(length_a * b_side / pi) 
    121  
    122     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    123     return 0.5 * (ddd) ** (1. / 3.) 
    124  
    125 def VR(length_a, b2a_ratio, c2a_ratio): 
    126     """ 
    127         Return shell volume and total volume 
    128     """ 
    129     b_side = length_a * b2a_ratio 
    130     c_side = length_a * c2a_ratio 
    131     vol_total = length_a * b_side * c_side 
    132     vol_shell = 2.0 * (length_a*b_side + length_a*c_side + b_side*c_side) 
    133     return vol_shell, vol_total 
     111effective_radius_type = [ 
     112    "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     113    "equivalent outer circular cross-section", 
     114    "half ab diagonal", "half diagonal", 
     115    ] 
    134116 
    135117 

  • sasmodels/models/mono_gauss_coil.py

    r2d81cfe ree60aa7  
    5454 
    5555import numpy as np 
    56 from numpy import inf, exp, errstate 
     56from numpy import inf 
    5757 
    5858name = "mono_gauss_coil" 
    … …  
    6969parameters = [ 
    7070    ["i_zero", "1/cm", 70.0, [0.0, inf], "", "Intensity at q=0"], 
    71     ["rg", "Ang", 75.0, [0.0, inf], "", "Radius of gyration"], 
     71    ["rg", "Ang", 75.0, [0.0, inf], "volume", "Radius of gyration"], 
    7272    ] 
    7373# pylint: enable=bad-whitespace, line-too-long 
    7474 
    75 # NB: Scale and Background are implicit parameters on every model 
    76 def Iq(q, i_zero, rg): 
    77     # pylint: disable = missing-docstring 
    78     z = (q * rg)**2 
     75source = ["mono_gauss_coil.c"] 
     76have_Fq = False 
     77effective_radius_type = ["R_g", "2R_g", "3R_g", "(5/3)^0.5*R_g"] 
    7978 
    80     with errstate(invalid='ignore'): 
    81         inten = (i_zero * 2.0) * (exp(-z) + z - 1.0)/z**2 
    82         inten[q == 0] = i_zero 
    83     return inten 
    84 Iq.vectorized = True # Iq accepts an array of q values 
    8579 
    8680def random(): 

  • sasmodels/models/multilayer_vesicle.c

    r71b751d ree60aa7  
    4747} 
    4848 
     49static double 
     50effective_radius(int mode, double radius, double thick_shell, double thick_solvent, double fp_n_shells) 
     51{ 
     52    // case 1: outer radius 
     53    return radius + fp_n_shells*thick_shell + (fp_n_shells - 1.0)*thick_solvent; 
     54} 
     55 
    4956static void 
    5057Fq(double q, 

  • sasmodels/models/multilayer_vesicle.py

    r71b751d ree60aa7  
    146146source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
    147147have_Fq = True 
    148  
    149 def ER(radius, thick_shell, thick_solvent, n_shells): 
    150     n_shells = int(n_shells+0.5) 
    151     return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent 
     148effective_radius_type = ["outer radius"] 
    152149 
    153150def random(): 

  • sasmodels/models/onion.c

    r71b751d ree60aa7  
    3030 
    3131static double 
    32 form_volume(double radius_core, double n_shells, double thickness[]) 
     32outer_radius(double radius_core, double n_shells, double thickness[]) 
    3333{ 
    3434  int n = (int)(n_shells+0.5); 
    … …  
    3737    r += thickness[i]; 
    3838  } 
    39   return M_4PI_3*cube(r); 
     39  return r; 
     40} 
     41 
     42static double 
     43form_volume(double radius_core, double n_shells, double thickness[]) 
     44{ 
     45  return M_4PI_3*cube(outer_radius(radius_core, n_shells, thickness)); 
     46} 
     47 
     48static double 
     49effective_radius(int mode, double radius_core, double n_shells, double thickness[]) 
     50{ 
     51  // case 1: outer radius 
     52  return outer_radius(radius_core, n_shells, thickness); 
    4053} 
    4154 

  • sasmodels/models/onion.py

    r71b751d ree60aa7  
    316316single = False 
    317317have_Fq = True 
     318effective_radius_type = ["outer radius"] 
    318319 
    319320profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    … …  
    365366 
    366367    return np.asarray(z), np.asarray(rho) 
    367  
    368 def ER(radius_core, n_shells, thickness): 
    369     """Effective radius""" 
    370     n = int(n_shells[0]+0.5) 
    371     return np.sum(thickness[:n], axis=0) + radius_core 
    372368 
    373369demo = { 

  • sasmodels/models/parallelepiped.c

    r71b751d ree60aa7  
    55} 
    66 
     7static double 
     8effective_radius(int mode, double length_a, double length_b, double length_c) 
     9{ 
     10    switch (mode) { 
     11    case 1: // equivalent sphere 
     12        return cbrt(0.75*length_a*length_b*length_c/M_PI); 
     13    case 2: // half length_a 
     14        return 0.5 * length_a; 
     15    case 3: // half length_b 
     16        return 0.5 * length_b; 
     17    case 4: // half length_c 
     18        return 0.5 * length_c; 
     19    case 5: // equivalent circular cross-section 
     20        return sqrt(length_a*length_b/M_PI); 
     21    case 6: // half ab diagonal 
     22        return 0.5*sqrt(length_a*length_a + length_b*length_b); 
     23    case 7: // half diagonal 
     24        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c); 
     25    } 
     26} 
    727 
    828static void 

  • sasmodels/models/parallelepiped.py

    r71b751d ree60aa7  
    140140   ensuring that the inequality $A < B < C$ is not violated,  The calculation 
    141141   will not report an error, but the results may be not correct. 
    142     
     142 
    143143.. _parallelepiped-orientation: 
    144144 
    … …  
    231231source = ["lib/gauss76.c", "parallelepiped.c"] 
    232232have_Fq = True 
    233  
    234 def ER(length_a, length_b, length_c): 
    235     """ 
    236     Return effective radius (ER) for P(q)*S(q) 
    237     """ 
    238     # now that axes can be in any size order, need to sort a,b,c 
    239     # where a~b and c is either much smaller or much larger 
    240     abc = np.vstack((length_a, length_b, length_c)) 
    241     abc = np.sort(abc, axis=0) 
    242     selector = (abc[1] - abc[0]) > (abc[2] - abc[1]) 
    243     length = np.where(selector, abc[0], abc[2]) 
    244     # surface average radius (rough approximation) 
    245     radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi) 
    246  
    247     ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius)) 
    248     return 0.5 * (ddd) ** (1. / 3.) 
    249  
    250 # VR defaults to 1.0 
    251  
     233effective_radius_type = [ 
     234    "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     235    "equivalent circular cross-section", "half ab diagonal", "half diagonal", 
     236    ] 
    252237 
    253238def random(): 

  • sasmodels/models/pearl_necklace.py

    r2d81cfe r2cc8aa2  
    140140            thick_string_pd=0.2, thick_string_pd_n=5, 
    141141           ) 
    142  
    143 tests = [[{}, 0.001, 17380.245], [{}, 'ER', 115.39502]] 
     142# ER function is not being used here, not that it is likely very sensible to  
     143# include an S(Q) with this model, the default in sasview 5.0 would be to the  
     144# "unconstrained" radius_effective. 
     145#tests = [[{}, 0.001, 17380.245], [{}, 'ER', 115.39502]] 
     146tests = [[{}, 0.001, 17380.245]] 

  • sasmodels/models/pringle.c

    r74768cb ree60aa7  
    104104} 
    105105 
     106static double 
     107effective_radius(int mode, double radius, double thickness, double alpha, double beta) 
     108{ 
     109    switch (mode) { 
     110    case 1: // equivalent sphere 
     111        return cbrt(0.75*radius*radius*thickness); 
     112    case 2: // radius 
     113        return radius; 
     114    } 
     115} 
     116 
    106117double Iq( 
    107118    double q, 

  • sasmodels/models/pringle.py

    ref07e95 ree60aa7  
    7272 
    7373 
    74 source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", \ 
     74source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", 
    7575          "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"] 
    76  
    77 def ER(radius, thickness, alpha, beta): 
    78     """ 
    79     Return equivalent radius (ER) 
    80     """ 
    81     ddd = 0.75 * radius * (2 * radius * thickness + (thickness + radius) \ 
    82                            * (thickness + pi * radius)) 
    83     return 0.5 * (ddd) ** (1. / 3.) 
     76effective_radius_type = ["equivalent sphere", "radius"] 
    8477 
    8578def random(): 

  • sasmodels/models/raspberry.c

    r71b751d ree60aa7  
    1 double form_volume(double radius_lg); 
     1double form_volume(double radius_lg, double radius_sm, double penetration); 
    22 
    33double Iq(double q, 
    … …  
    66          double radius_lg, double radius_sm, double penetration); 
    77 
    8 double form_volume(double radius_lg) 
     8double form_volume(double radius_lg, double radius_sm, double penetration) 
    99{ 
    1010    //Because of the complex structure, volume normalization must 
    … …  
    1212    double volume=1.0; 
    1313    return volume; 
     14} 
     15 
     16static double 
     17effective_radius(int mode, double radius_lg, double radius_sm, double penetration) 
     18{ 
     19    switch (mode) { 
     20    case 1: // radius_large 
     21        return radius_lg; 
     22    case 2: // radius_outer 
     23        return radius_lg + 2.0*radius_sm - penetration; 
     24    } 
    1425} 
    1526 

  • sasmodels/models/raspberry.py

    ref07e95 ree60aa7  
    145145              ["radius_lg", "Ang", 5000, [0, inf], "volume", 
    146146               "radius of large spheres"], 
    147               ["radius_sm", "Ang", 100, [0, inf], "", 
     147              ["radius_sm", "Ang", 100, [0, inf], "volume", 
    148148               "radius of small spheres"], 
    149               ["penetration", "Ang", 0, [-1, 1], "", 
     149              ["penetration", "Ang", 0, [-1, 1], "volume", 
    150150               "fractional penetration depth of small spheres into large sphere"], 
    151151             ] 
    152152 
    153153source = ["lib/sas_3j1x_x.c", "raspberry.c"] 
     154effective_radius_type = ["radius_large", "radius_outer"] 
    154155 
    155156def random(): 

  • sasmodels/models/rectangular_prism.c

    r71b751d ree60aa7  
    33{ 
    44    return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio); 
     5} 
     6 
     7static double 
     8effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
     9{ 
     10    switch (mode) { 
     11    case 1: // equivalent sphere 
     12        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
     13    case 2: // half length_a 
     14        return 0.5 * length_a; 
     15    case 3: // half length_b 
     16        return 0.5 * length_a*b2a_ratio; 
     17    case 4: // half length_c 
     18        return 0.5 * length_a*c2a_ratio; 
     19    case 5: // equivalent circular cross-section 
     20        return length_a*sqrt(b2a_ratio/M_PI); 
     21    case 6: // half ab diagonal 
     22        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     23    case 7: // half diagonal 
     24        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     25    } 
    526} 
    627 

  • sasmodels/models/rectangular_prism.py

    r71b751d ree60aa7  
    136136source = ["lib/gauss76.c", "rectangular_prism.c"] 
    137137have_Fq = True 
    138  
    139 def ER(length_a, b2a_ratio, c2a_ratio): 
    140     """ 
    141         Return equivalent radius (ER) 
    142     """ 
    143     b_side = length_a * b2a_ratio 
    144     c_side = length_a * c2a_ratio 
    145  
    146     # surface average radius (rough approximation) 
    147     surf_rad = sqrt(length_a * b_side / pi) 
    148  
    149     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    150     return 0.5 * (ddd) ** (1. / 3.) 
     138effective_radius_type = [ 
     139    "equivalent sphere", "half length_a", "half length_b", "half length_c", 
     140    "equivalent circular cross-section", "half ab diagonal", "half diagonal", 
     141    ] 
    151142 
    152143def random(): 

  • sasmodels/models/sphere.py

    r71b751d ree60aa7  
    6767             ] 
    6868 
    69 source = ["lib/sas_3j1x_x.c"] 
     69source = ["lib/sas_3j1x_x.c","sphere.c"] 
    7070have_Fq = True 
    71  
    72 c_code = """ 
    73 static double form_volume(double radius) 
    74 { 
    75     return M_4PI_3*cube(radius); 
    76 } 
    77  
    78 static void Fq(double q, double *f1, double *f2, double sld, double sld_solvent, double radius) 
    79 { 
    80     const double bes = sas_3j1x_x(q*radius); 
    81     const double contrast = (sld - sld_solvent); 
    82     const double form = contrast * form_volume(radius) * bes; 
    83     *f1 = 1.0e-2*form; 
    84     *f2 = 1.0e-4*form*form; 
    85 } 
    86 """ 
    87  
    88 def ER(radius): 
    89     """ 
    90     Return equivalent radius (ER) 
    91     """ 
    92     return radius 
    93  
    94 # VR defaults to 1.0 
     71effective_radius_type = ["radius"] 
    9572 
    9673def random(): 
    … …  
    10683      "radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, 
    10784     0.2, 0.228843], 
    108     [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.], 
    109     [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "VR", 1.], 
     85#    [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.], 
     86#    [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "VR", 1.], 
    11087] 

  • sasmodels/models/spherical_sld.c

    r71b751d ree60aa7  
    1 static double form_volume( 
    2     double fp_n_shells, 
    3     double thickness[], 
    4     double interface[]) 
     1static double 
     2outer_radius(double fp_n_shells, double thickness[], double interface[]) 
    53{ 
    64    int n_shells= (int)(fp_n_shells + 0.5); 
    … …  
    97        r += thickness[i] + interface[i]; 
    108    } 
    11     return M_4PI_3*cube(r); 
     9    return r; 
     10} 
     11 
     12static double form_volume( 
     13    double fp_n_shells, 
     14    double thickness[], 
     15    double interface[]) 
     16{ 
     17    return M_4PI_3*cube(outer_radius(fp_n_shells, thickness, interface)); 
     18} 
     19 
     20static double 
     21effective_radius(int mode, double fp_n_shells, double thickness[], double interface[]) 
     22{ 
     23    // case 1: outer radius 
     24    return outer_radius(fp_n_shells, thickness, interface); 
    1225} 
    1326 

  • sasmodels/models/spherical_sld.py

    r5601947 ree60aa7  
    2222 
    2323    0: erf($\nu z$) 
    24      
     24 
    2525    1: Rpow($z^\nu$) 
    26      
     26 
    2727    2: Lpow($z^\nu$) 
    28      
     28 
    2929    3: Rexp($-\nu z$) 
    30      
     30 
    3131    4: Lexp($-\nu z$) 
    3232 
    … …  
    217217              ["interface[n_shells]",  "Ang",        50.0,   [0, inf],       "volume", "thickness of the interface"], 
    218218              ["shape[n_shells]",      "",           0,      [SHAPES],       "", "interface shape"], 
    219               ["nu[n_shells]",         "",           2.5,    [0, inf],       "", "interface shape exponent"], 
     219              ["nu[n_shells]",         "",           2.5,    [1, inf],       "", "interface shape exponent"], 
    220220              ["n_steps",              "",           35,     [0, inf],       "", "number of steps in each interface (must be an odd integer)"], 
    221221             ] 
    … …  
    224224single = False  # TODO: fix low q behaviour 
    225225have_Fq = True 
     226effective_radius_type = ["outer radius"] 
    226227 
    227228profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    … …  
    269270 
    270271 
    271 def ER(n_shells, thickness, interface): 
    272     """Effective radius""" 
    273     n_shells = int(n_shells + 0.5) 
    274     total = (np.sum(thickness[:n_shells], axis=1) 
    275              + np.sum(interface[:n_shells], axis=1)) 
    276     return total 
    277  
    278  
    279272demo = { 
    280273    "n_shells": 5, 

  • sasmodels/models/triaxial_ellipsoid.c

    r71b751d ree60aa7  
    77} 
    88 
     9static double 
     10radius_from_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     11{ 
     12    return cbrt(radius_equat_minor*radius_equat_major*radius_polar); 
     13} 
     14 
     15static double 
     16radius_from_min_dimension(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     17{ 
     18    const double rad_equat_min = (radius_equat_minor < radius_equat_major ? radius_equat_minor : radius_equat_major); 
     19    return (rad_equat_min < radius_polar ? rad_equat_min : radius_polar); 
     20} 
     21 
     22static double 
     23radius_from_max_dimension(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     24{ 
     25    const double rad_equat_max = (radius_equat_minor < radius_equat_major ? radius_equat_major : radius_equat_minor); 
     26    return (rad_equat_max > radius_polar ? rad_equat_max : radius_polar); 
     27} 
     28 
     29static double 
     30effective_radius(int mode, double radius_equat_minor, double radius_equat_major, double radius_polar) 
     31{ 
     32    switch (mode) { 
     33    case 1: // equivalent sphere 
     34        return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar); 
     35    case 2: // min radius 
     36        return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
     37    case 3: // max radius 
     38        return radius_from_max_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
     39    } 
     40} 
    941 
    1042static void 

  • sasmodels/models/triaxial_ellipsoid.py

    r71b751d ree60aa7  
    158158source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 
    159159have_Fq = True 
    160  
    161 def ER(radius_equat_minor, radius_equat_major, radius_polar): 
    162     """ 
    163     Returns the effective radius used in the S*P calculation 
    164     """ 
    165     from .ellipsoid import ER as ellipsoid_ER 
    166  
    167     # now that radii can be in any size order, radii need sorting a,b,c 
    168     # where a~b and c is either much smaller or much larger 
    169     radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar)) 
    170     radii = np.sort(radii, axis=0) 
    171     selector = (radii[1] - radii[0]) > (radii[2] - radii[1]) 
    172     polar = np.where(selector, radii[0], radii[2]) 
    173     equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2])) 
    174     return ellipsoid_ER(polar, equatorial) 
     160effective_radius_type = ["equivalent sphere", "min radius", "max radius"] 
    175161 
    176162def random(): 

  • sasmodels/models/vesicle.c

    r71b751d ree60aa7  
     1// TODO: interface to form_volume/shell_volume not yet settled 
     2static double 
     3shell_volume(double *total, double radius, double thickness) 
     4{ 
     5    //note that for the vesicle model, the volume is ONLY the shell volume 
     6    *total = M_4PI_3 * cube(radius+thickness); 
     7    return *total - M_4PI_3*cube(radius); 
     8} 
     9 
    110static double 
    211form_volume(double radius, double thickness) 
    312{ 
    413    //note that for the vesicle model, the volume is ONLY the shell volume 
    5     return M_4PI_3*(cube(radius+thickness) - cube(radius)); 
     14    double total; 
     15    return shell_volume(&total, radius, thickness); 
    616} 
    717 
     18static double 
     19effective_radius(int mode, double radius, double thickness) 
     20{ 
     21    // case 1: outer radius 
     22    return radius + thickness; 
     23} 
    824 
    925static void 

  • sasmodels/models/vesicle.py

    rb477605 ree60aa7  
    100100source = ["lib/sas_3j1x_x.c", "vesicle.c"] 
    101101have_Fq = True 
    102  
    103 def ER(radius, thickness): 
    104     ''' 
    105     returns the effective radius used in the S*P calculation 
    106  
    107     :param radius: core radius 
    108     :param thickness: shell thickness 
    109     ''' 
    110     return radius + thickness 
    111  
    112 def VR(radius, thickness): 
    113     ''' 
    114     returns the volumes of the shell and of the whole sphere including the 
    115     core plus shell - is used to normalize when including polydispersity. 
    116  
    117     :param radius: core radius 
    118     :param thickness: shell thickness 
    119     :return whole: volume of core and shell 
    120     :return whole-core: volume of the shell 
    121     ''' 
    122  
    123     whole = 4./3. * pi * (radius + thickness)**3 
    124     core = 4./3. * pi * radius**3 
    125     return whole, whole - core 
     102effective_radius_type = ["outer radius"] 
    126103 
    127104def random(): 
    … …  
    151128         [{}, 0.100600200401, 1.77063682331], 
    152129         [{}, 0.5, 0.00355351388906], 
    153          [{}, 'ER', 130.], 
    154          [{}, 'VR', 0.54483386436], 
     130#         [{}, 'ER', 130.], 
     131#         [{}, 'VR', 0.54483386436], 
    155132        ] 

  • sasmodels/product.py

    rc0131d44 r33d7be3  
    2424# pylint: disable=unused-import 
    2525try: 
    26     from typing import Tuple, Callable 
     26    from typing import Tuple, Callable, Union 
    2727except ImportError: 
    2828    pass 
    … …  
    3737#] 
    3838 
    39 ER_ID = "radius_effective" 
    40 VF_ID = "volfraction" 
    41  
     39RADIUS_ID = "radius_effective" 
     40VOLFRAC_ID = "volfraction" 
    4241def make_extra_pars(p_info): 
    4342    pars = [] 
    … …  
    5150                "Structure factor calculation") 
    5251        pars.append(par) 
     52    if p_info.effective_radius_type is not None: 
     53        par = parse_parameter( 
     54                "radius_effective_mode", 
     55                "", 
     56                0, 
     57                [["unconstrained"] + p_info.effective_radius_type], 
     58                "", 
     59                "Effective radius calculation") 
     60        pars.append(par) 
    5361    return pars 
    5462 
    … …  
    6573    # have any magnetic parameters 
    6674    if not len(s_info.parameters.kernel_parameters) >= 2: 
    67         raise TypeError("S needs {} and {} as its first parameters".format(ER_ID, VF_ID)) 
    68     if not s_info.parameters.kernel_parameters[0].id == ER_ID: 
    69         raise TypeError("S needs {} as first parameter".format(ER_ID)) 
    70     if not s_info.parameters.kernel_parameters[1].id == VF_ID: 
    71         raise TypeError("S needs {} as second parameter".format(VF_ID)) 
     75        raise TypeError("S needs {} and {} as its first parameters".format(RADIUS_ID, VOLFRAC_ID)) 
     76    if not s_info.parameters.kernel_parameters[0].id == RADIUS_ID: 
     77        raise TypeError("S needs {} as first parameter".format(RADIUS_ID)) 
     78    if not s_info.parameters.kernel_parameters[1].id == VOLFRAC_ID: 
     79        raise TypeError("S needs {} as second parameter".format(VOLFRAC_ID)) 
    7280    if not s_info.parameters.magnetism_index == []: 
    7381        raise TypeError("S should not have SLD parameters") 
    … …  
    7583    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    7684 
    77     # Create list of parameters for the combined model.  Skip the first 
    78     # parameter of S, which we verified above is effective radius.  If there 
     85    # Create list of parameters for the combined model.  If there 
    7986    # are any names in P that overlap with those in S, modify the name in S 
    8087    # to distinguish it. 
    8188    p_set = set(p.id for p in p_pars.kernel_parameters) 
    8289    s_list = [(_tag_parameter(par) if par.id in p_set else par) 
    83               for par in s_pars.kernel_parameters[1:]] 
     90              for par in s_pars.kernel_parameters] 
    8491    # Check if still a collision after renaming.  This could happen if for 
    8592    # example S has volfrac and P has both volfrac and volfrac_S. 
    … …  
    8794        raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 
    8895 
     96    # make sure effective radius is not a polydisperse parameter in product 
     97    s_list[0] = copy(s_list[0]) 
     98    s_list[0].polydisperse = False 
     99 
    89100    translate_name = dict((old.id, new.id) for old, new 
    90                           in zip(s_pars.kernel_parameters[1:], s_list)) 
     101                          in zip(s_pars.kernel_parameters, s_list)) 
    91102    combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info) 
    92103    parameters = ParameterTable(combined_pars) 
    … …  
    94105    def random(): 
    95106        combined_pars = p_info.random() 
    96         s_names = set(par.id for par in s_pars.kernel_parameters[1:]) 
     107        s_names = set(par.id for par in s_pars.kernel_parameters) 
    97108        combined_pars.update((translate_name[k], v) 
    98109                             for k, v in s_info.random().items() 
    … …  
    157168    return par 
    158169 
    159 def _intermediates(P, S): 
    160     # type: (np.ndarray, np.ndarray) -> OrderedDict[str, np.ndarray] 
     170def _intermediates(P, S, effective_radius): 
     171    # type: (np.ndarray, np.ndarray, float) -> OrderedDict[str, np.ndarray] 
    161172    """ 
    162173    Returns intermediate results for standard product (P(Q)*S(Q)) 
    … …  
    165176        ("P(Q)", P), 
    166177        ("S(Q)", S), 
     178        ("effective_radius", effective_radius), 
    167179    )) 
    168180 
    169 def _intermediates_beta(F1, F2, S, scale, bg): 
    170     # type: (np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray) -> OrderedDict[str, np.ndarray] 
    171     """ 
    172     Returns intermediate results for beta approximation-enabled product 
     181def _intermediates_beta(F1,              # type: np.ndarray 
     182                        F2,              # type: np.ndarray 
     183                        S,               # type: np.ndarray 
     184                        scale,           # type: np.ndarray 
     185                        bg,              # type: np.ndarray 
     186                        effective_radius # type: float 
     187                        ): 
     188    # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] 
     189    """ 
     190    Returns intermediate results for beta approximation-enabled product. The result may be an array or a float. 
    173191    """ 
    174192    # TODO: 1. include calculated Q vector 
    … …  
    179197        ("beta(Q)", F1**2 / F2), 
    180198        ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), 
     199        ("effective_radius", effective_radius), 
    181200        # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), 
    182201    )) 
    … …  
    262281        p_values = np.hstack(p_values).astype(self.p_kernel.dtype) 
    263282 
    264         # Call ER and VR for P since these are needed for S. 
    265         p_er, p_vr = calc_er_vr(p_info, p_details, p_values) 
    266         s_vr = (volfrac/p_vr if p_vr != 0. else volfrac) 
    267         #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr)) 
    268  
    269283        # Construct the calling parameters for S. 
    270         # The  effective radius is not in the combined parameter list, so 
    271         # the number of 'S' parameters is one less than expected.  The 
    272         # computed effective radius needs to be added into the weights 
    273         # vector, especially since it is a polydisperse parameter in the 
    274         # stand-alone structure factor models.  We will added it at the 
    275         # end so the remaining offsets don't need to change. 
    276         s_npars = s_info.parameters.npars-1 
     284        s_npars = s_info.parameters.npars 
    277285        s_length = call_details.length[p_npars:p_npars+s_npars] 
    278286        s_offset = call_details.offset[p_npars:p_npars+s_npars] 
    … …  
    280288        s_offset = np.hstack((nweights, s_offset)) 
    281289        s_details = make_details(s_info, s_length, s_offset, nweights+1) 
    282         v, w = weights[:nweights], weights[nweights:] 
    283290        s_values = [ 
    284             # scale=1, background=0, radius_effective=p_er, volfraction=s_vr 
    285             [1., 0., p_er, s_vr], 
    286             # structure factor parameters start after scale, background and 
    287             # all the form factor parameters.  Skip the volfraction parameter 
    288             # as well, since it is computed elsewhere, and go to the end of the 
    289             # parameter list. 
    290             values[2+p_npars+1:2+p_npars+s_npars], 
    291             # no magnetism parameters to include for S 
    292             # add er into the (value, weights) pairs 
    293             v, [p_er], w, [1.0] 
     291            # scale=1, background=0, 
     292            [1., 0.], 
     293            values[2+p_npars:2+p_npars+s_npars], 
     294            weights, 
    294295        ] 
    295296        spacer = (32 - sum(len(v) for v in s_values)%32)%32 
    … …  
    298299 
    299300        # beta mode is the first parameter after the structure factor pars 
    300         beta_index = 2+p_npars+s_npars 
    301         beta_mode = values[beta_index] 
     301        extra_offset = 2+p_npars+s_npars 
     302        if p_info.have_Fq: 
     303            beta_mode = values[extra_offset] 
     304            extra_offset += 1 
     305        else: 
     306            beta_mode = 0 
     307        if p_info.effective_radius_type is not None: 
     308            effective_radius_type = int(values[extra_offset]) 
     309            extra_offset += 1 
     310        else: 
     311            effective_radius_type = 0 
    302312 
    303313        # Call the kernels 
    304         s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    305314        scale, background = values[0], values[1] 
    306315        if beta_mode: 
    307             F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic) 
     316            F1, F2, volume_avg, effective_radius = self.p_kernel.beta( 
     317                p_details, p_values, cutoff, magnetic, effective_radius_type) 
     318            if effective_radius_type > 0: 
     319                # Plug R_eff from p into S model (initial value and pd value) 
     320                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 
     321            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    308322            combined_scale = scale*volfrac/volume_avg 
    309323 
    310             self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background) 
     324            self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background, effective_radius) 
    311325            final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background 
    312326 
    313327        else: 
    314             p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic) 
    315  
    316             self.results = lambda: _intermediates(p_result, s_result) 
     328            p_result, effective_radius = self.p_kernel.Pq_Reff( 
     329                p_details, p_values, cutoff, magnetic, effective_radius_type) 
     330            if effective_radius_type > 0: 
     331                # Plug R_eff from p into S model (initial value and pd value) 
     332                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 
     333            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
     334            # remember the parts for plotting later 
     335            self.results = lambda: _intermediates(p_result, s_result, effective_radius) 
    317336            final_result = scale*(p_result*s_result) + background 
    318337 

  • sasmodels/sasview_model.py

    r84f2962 rc952e59  
    272272    attrs['category'] = model_info.category 
    273273    attrs['is_structure_factor'] = model_info.structure_factor 
    274     attrs['is_form_factor'] = model_info.ER is not None 
     274    attrs['is_form_factor'] = model_info.effective_radius_type is not None 
    275275    attrs['is_multiplicity_model'] = variants[0] > 1 
    276276    attrs['multiplicity_info'] = variants 
Note: See TracChangeset for help on using the changeset viewer.