Changes in / [c88f983:2773c66] in sasmodels


Ignore:
Files:
3 added
65 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)) 

  • sasmodels/core.py

    r71b751d r5399809  
    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 structure_factor_mode" 
     367              " A_sld A_sld_solvent A_radius").split() 
    367368    assert target == actual, "%s != %s"%(target, actual) 
    368369 

  • 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 r5399809  
    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 NotImpmentedError() 

  • 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 rd277229  
    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    if (mode == 1) { 
     82        return radius_from_volume(radius_bell, radius , length); 
     83    } else if (mode == 2) { 
     84        return radius; 
     85    } else if (mode == 3) { 
     86        return 0.5*length; 
     87    } else { 
     88        return radius_from_totallength(radius_bell,radius,length); 
     89    } 
     90} 
     91 
    6492static void 
    6593Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 

  • sasmodels/models/barbell.py

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

  • sasmodels/models/capped_cylinder.c

    r71b751d rd277229  
    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    if (mode == 1) { 
     104        return radius_from_volume(radius, radius_cap, length); 
     105    } else if (mode == 2) { 
     106        return radius; 
     107    } else if (mode == 3) { 
     108        return 0.5*length; 
     109    } else { 
     110        return radius_from_totallength(radius, radius_cap,length); 
     111    } 
     112} 
     113 
    86114static void 
    87115Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 

  • sasmodels/models/capped_cylinder.py

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

  • sasmodels/models/core_multi_shell.c

    r71b751d rd277229  
    1919} 
    2020 
     21static double 
     22outer_radius(double core_radius, double fp_n, double thickness[]) 
     23{ 
     24  double r = core_radius; 
     25  int n = (int)(fp_n+0.5); 
     26  for (int i=0; i < n; i++) { 
     27    r += thickness[i]; 
     28  } 
     29  return r; 
     30} 
     31 
     32static double 
     33effective_radius(int mode, double core_radius, double fp_n, double thickness[]) 
     34{ 
     35    if (mode == 1) { 
     36        double r = core_radius; 
     37        int n = (int)(fp_n+0.5); 
     38        for (int i=0; i < n; i++) { 
     39            r += thickness[i]; 
     40        } 
     41        return r; 
     42        //return outer_radius(core_radius,fp_n,thickness); 
     43    } else { 
     44        return core_radius; 
     45    } 
     46} 
     47 
    2148static void 
    2249Fq(double q, double *F1, double *F2, double core_sld, double core_radius, 

  • sasmodels/models/core_multi_shell.py

    r71b751d rd277229  
    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 
     147#def ER(radius, n, thickness): 
     148#    """Effective radius""" 
     149#    n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse 
     150#    return np.sum(thickness[:n], axis=0) + radius 
    150151 
    151152demo = dict(sld_core=6.4, 

  • sasmodels/models/core_shell_bicelle.c

    r71b751d rd277229  
    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    if (mode == 1) { 
     57        return radius_from_volume(radius, thick_rim, thick_face, length); 
     58    } else if (mode == 2) { 
     59        return radius + thick_rim; 
     60    } else if (mode == 3) { 
     61        return 0.5*length + thick_face; 
     62    } else { 
     63        return radius_from_diagonal(radius,thick_rim,thick_face,length); 
     64    } 
    3665} 
    3766 

  • sasmodels/models/core_shell_bicelle.py

    r71b751d rd277229  
    155155          "core_shell_bicelle.c"] 
    156156have_Fq = True 
     157effective_radius_type = ["equivalent sphere","outer rim radius", "half outer thickness","half diagonal"] 
    157158 
    158159def random(): 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r71b751d rd277229  
    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    if (mode == 1) { 
     32        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     33    } else if (mode == 2) { 
     34        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     35    } else if (mode == 3) { 
     36        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     37    } else if (mode == 4) { 
     38        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     39    } else if (mode ==5) { 
     40        return 0.5*length + thick_face; 
     41    } else { 
     42        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     43    } 
    1044} 
    1145 

  • sasmodels/models/core_shell_bicelle_elliptical.py

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

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r71b751d rd277229  
    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    if (mode == 1) { 
     33        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     34    } else if (mode == 2) { 
     35        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     36    } else if (mode == 3) { 
     37        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     38    } else if (mode == 4) { 
     39        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     40    } else if (mode ==5) { 
     41        return 0.5*length + thick_face; 
     42    } else { 
     43        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     44    } 
    1145} 
    1246 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

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

  • sasmodels/models/core_shell_cylinder.c

    r71b751d rd277229  
    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 + 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    if (mode == 1) { 
     34        return radius_from_volume(radius, thickness, length); 
     35    } else if (mode == 2) { 
     36        return radius + thickness; 
     37    } else if (mode == 3) { 
     38        return 0.5*length + thickness; 
     39    } else if (mode == 4) { 
     40        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness); 
     41    } else if (mode == 5) { 
     42        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness); 
     43    } else { 
     44        return radius_from_diagonal(radius,thickness,length); 
     45    } 
    1346} 
    1447 

  • sasmodels/models/core_shell_cylinder.py

    r71b751d rd277229  
    126126source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 
    127127have_Fq = True 
     128effective_radius_type = ["equivalent sphere","outer radius","half outer length","half min outer dimension", 
     129                         "half max outer dimension","half outer diagonal"] 
    128130 
    129 def ER(radius, thickness, length): 
    130     """ 
    131     Returns the effective radius used in the S*P calculation 
    132     """ 
    133     radius = radius + thickness 
    134     length = length + 2 * thickness 
    135     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    136     return 0.5 * (ddd) ** (1. / 3.) 
     131#def ER(radius, thickness, length): 
     132#    """ 
     133#    Returns the effective radius used in the S*P calculation 
     134#    """ 
     135#    radius = radius + thickness 
     136#    length = length + 2 * thickness 
     137#    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
     138#    return 0.5 * (ddd) ** (1. / 3.) 
    137139 
    138140def VR(radius, thickness, length): 

  • sasmodels/models/core_shell_ellipsoid.c

    r71b751d r3c60146  
    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    if (mode == 1) { 
     76        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     77    } else if (mode == 2) { 
     78        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     79    } else if (mode == 3) { 
     80        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     81    } else { 
     82        return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     83    } 
    3884} 
    3985 

  • sasmodels/models/core_shell_ellipsoid.py

    r71b751d rd277229  
    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 = ["equivalent sphere","average outer curvature", "min outer radius", "max outer radius"] 
     149 
     150#def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 
     151#    """ 
     152#        Returns the effective radius used in the S*P calculation 
     153#    """ 
     154#    from .ellipsoid import ER as ellipsoid_ER 
     155#    polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell 
     156#    equat_outer = radius_equat_core + thick_shell 
     157#    return ellipsoid_ER(polar_outer, equat_outer) 
    157158 
    158159def random(): 

  • sasmodels/models/core_shell_parallelepiped.c

    r71b751d rd277229  
    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_paral = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     34    return cbrt(0.75*volume_paral/M_PI); 
     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    if (mode == 1) { 
     49        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     50    } else if (mode == 2) { 
     51        return 0.5 * (length_a + thick_rim_a); 
     52    } else if (mode == 3) { 
     53        return 0.5 * (length_b + thick_rim_b); 
     54    } else if (mode == 4) { 
     55        return 0.5 * (length_c + thick_rim_c); 
     56    } else if (mode == 5) { 
     57        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b); 
     58    } else if (mode == 6) { 
     59        return 0.5*sqrt(square(length_a+thick_rim_a) + square(length_b+thick_rim_b)); 
     60    } else { 
     61        return 0.5*sqrt(square(length_a+thick_rim_a) + square(length_b+thick_rim_b) + square(length_c+thick_rim_c)); 
     62    } 
    2763} 
    2864 

  • sasmodels/models/core_shell_parallelepiped.py

    r71b751d rd277229  
    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) 
     229effective_radius_type = ["equivalent sphere","half outer length_a", "half outer length_b", "half outer length_c", 
     230                         "equivalent circular cross-section","half outer ab diagonal","half outer diagonal"] 
     231 
     232#def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c): 
     233#    """ 
     234#        Return equivalent radius (ER) 
     235#    """ 
     236#    from .parallelepiped import ER as ER_p 
     237# 
     238#    a = length_a + 2*thick_rim_a 
     239#    b = length_b + 2*thick_rim_b 
     240#    c = length_c + 2*thick_rim_c 
     241#    return ER_p(a, b, c) 
    241242 
    242243# VR defaults to 1.0 

  • sasmodels/models/core_shell_sphere.c

    r71b751d rd277229  
    33{ 
    44    return M_4PI_3 * cube(radius + thickness); 
     5} 
     6 
     7static double 
     8effective_radius(int mode, double radius, double thickness) 
     9{ 
     10    if (mode == 1) { 
     11        return radius + thickness; 
     12    } else { 
     13        return radius; 
     14    } 
    515} 
    616 

  • sasmodels/models/core_shell_sphere.py

    r71b751d rd277229  
    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) 
    8384 
    84 def ER(radius, thickness): 
    85     """ 
    86         Equivalent radius 
    87         @param radius: core radius 
    88         @param thickness: shell thickness 
    89     """ 
    90     return radius + thickness 
     85#def ER(radius, thickness): 
     86#    """ 
     87#        Equivalent radius 
     88#        @param radius: core radius 
     89#        @param thickness: shell thickness 
     90#    """ 
     91#    return radius + thickness 
    9192 
    9293def VR(radius, thickness): 

  • sasmodels/models/cylinder.c

    r71b751d rd277229  
    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    if (mode == 1) { 
     31        return radius_from_volume(radius, length); 
     32    } else if (mode == 2) { 
     33        return radius; 
     34    } else if (mode == 3) { 
     35        return 0.5*length; 
     36    } else if (mode == 4) { 
     37        return (radius < 0.5*length ? radius : 0.5*length); 
     38    } else if (mode == 5) { 
     39        return (radius > 0.5*length ? radius : 0.5*length); 
     40    } else { 
     41        return radius_from_diagonal(radius,length); 
     42    } 
    1343} 
    1444 

  • sasmodels/models/cylinder.py

    r71b751d rd277229  
    139139source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
    140140have_Fq = True 
     141effective_radius_type = ["equivalent sphere","radius","half length","half min dimension","half max dimension","half diagonal"] 
    141142 
    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.) 
     143#def ER(radius, length): 
     144#    """ 
     145#        Return equivalent radius (ER) 
     146#    """ 
     147#    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
     148#    return 0.5 * (ddd) ** (1. / 3.) 
    148149 
    149150def random(): 

  • sasmodels/models/ellipsoid.c

    r71b751d rd277229  
    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    if (mode == 1) { 
     37        return radius_from_volume(radius_polar, radius_equatorial); 
     38    } else if (mode == 2) { 
     39        return radius_from_curvature(radius_polar, radius_equatorial); 
     40    } else if (mode == 3) { 
     41        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial); 
     42    } else { 
     43        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial); 
     44    } 
     45} 
     46 
    647 
    748static void 

  • sasmodels/models/ellipsoid.py

    r0168844 rd277229  
    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 = ["equivalent sphere","average curvature", "min radius", "max radius"] 
    193172 
    194173def random(): 

  • sasmodels/models/elliptical_cylinder.c

    r71b751d rd277229  
    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 length) 
     16{ 
     17    const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     18    return (rad_min < length ? rad_min : length); 
     19} 
     20 
     21static double 
     22radius_from_max_dimension(double radius_minor, double r_ratio, double length) 
     23{ 
     24    const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     25    return (rad_max > length ? rad_max : length); 
     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    if (mode == 1) { 
     39        return radius_from_volume(radius_minor, r_ratio, length); 
     40    } else if (mode == 2) { 
     41        return 0.5*radius_minor*(1.0 + r_ratio); 
     42    } else if (mode == 3) { 
     43        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     44    } else if (mode == 4) { 
     45        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     46    } else if (mode == 5) { 
     47        return sqrt(radius_minor*radius_minor*r_ratio); 
     48    } else if (mode == 6) { 
     49        return 0.5*length; 
     50    } else if (mode == 7) { 
     51        return radius_from_min_dimension(radius_minor,r_ratio,length); 
     52    } else if (mode == 8) { 
     53        return radius_from_max_dimension(radius_minor,r_ratio,length); 
     54    } else { 
     55        return radius_from_diagonal(radius_minor,r_ratio,length); 
     56    } 
    557} 
    658 

  • sasmodels/models/elliptical_cylinder.py

    r71b751d rd277229  
    123123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
    124124have_Fq = True 
     125effective_radius_type = ["equivalent sphere","average radius","min radius","max radius", 
     126                         "equivalent circular cross-section","half length","half min dimension","half max dimension","half diagonal"] 
    125127 
    126128demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 
    … …  
    128130            theta_pd=10, phi_pd=2, psi_pd=3) 
    129131 
    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.) 
     132#def ER(radius_minor, axis_ratio, length): 
     133#    """ 
     134#        Equivalent radius 
     135#        @param radius_minor: Ellipse minor radius 
     136#        @param axis_ratio: Ratio of major radius over minor radius 
     137#        @param length: Length of the cylinder 
     138#    """ 
     139#    radius = sqrt(radius_minor * radius_minor * axis_ratio) 
     140#    ddd = 0.75 * radius * (2 * radius * length 
     141#                           + (length + radius) * (length + pi * radius)) 
     142#    return 0.5 * (ddd) ** (1. / 3.) 
    141143 
    142144def random(): 

  • sasmodels/models/fuzzy_sphere.py

    r71b751d rd277229  
    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 
     86effective_radius_type = ["radius","radius + fuzziness"] 
    8687 
    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 
     88#def ER(radius): 
     89#    """ 
     90#    Return radius 
     91#    """ 
     92#    return radius 
    11193 
    11294# VR defaults to 1.0 

  • sasmodels/models/hollow_cylinder.c

    r71b751d rd277229  
    2222} 
    2323 
     24static double 
     25radius_from_volume(double radius, double thickness, double length) 
     26{ 
     27    const double volume_outer_cyl = M_PI*square(radius + thickness)*length; 
     28    return cbrt(0.75*volume_outer_cyl/M_PI); 
     29} 
     30 
     31static double 
     32radius_from_diagonal(double radius, double thickness, double length) 
     33{ 
     34    return sqrt(square(radius + thickness) + 0.25*square(length)); 
     35} 
     36 
     37static double 
     38effective_radius(int mode, double radius, double thickness, double length) 
     39{ 
     40    if (mode == 1) { 
     41        return radius_from_volume(radius, thickness, length); 
     42    } else if (mode == 2) { 
     43        return radius + thickness; 
     44    } else if (mode == 3) { 
     45        return 0.5*length; 
     46    } else if (mode == 4) { 
     47        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length); 
     48    } else if (mode == 5) { 
     49        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length); 
     50    } else { 
     51        return radius_from_diagonal(radius,thickness,length); 
     52    } 
     53} 
    2454 
    2555static void 

  • sasmodels/models/hollow_cylinder.py

    r71b751d rd277229  
    9090source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    9191have_Fq = True 
     92effective_radius_type = ["equivalent sphere","outer radius","half length", 
     93                         "half outer min dimension","half outer max dimension","half outer diagonal"] 
    9294 
    9395# pylint: disable=W0613 
    94 def ER(radius, thickness, length): 
    95     """ 
    96     :param radius:      Cylinder core radius 
    97     :param thickness:   Cylinder wall thickness 
    98     :param length:      Cylinder length 
    99     :return:            Effective radius 
    100     """ 
    101     router = radius + thickness 
    102     if router == 0 or length == 0: 
    103         return 0.0 
    104     len1 = router 
    105     len2 = length/2.0 
    106     term1 = len1*len1*2.0*len2/2.0 
    107     term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0) 
    108     ddd = 3.0*term1*term2 
    109     diam = pow(ddd, (1.0/3.0)) 
    110     return diam 
     96#def ER(radius, thickness, length): 
     97#    """ 
     98#    :param radius:      Cylinder core radius 
     99#    :param thickness:   Cylinder wall thickness 
     100#    :param length:      Cylinder length 
     101#    :return:            Effective radius 
     102#    """ 
     103#    router = radius + thickness 
     104#    if router == 0 or length == 0: 
     105#        return 0.0 
     106#    len1 = router 
     107#    len2 = length/2.0 
     108#    term1 = len1*len1*2.0*len2/2.0 
     109#    term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0) 
     110#    ddd = 3.0*term1*term2 
     111#    diam = pow(ddd, (1.0/3.0)) 
     112#    return diam 
    111113 
    112114def VR(radius, thickness, length): 

  • sasmodels/models/hollow_rectangular_prism.c

    r71b751d rd277229  
    1111    double vol_shell = vol_total - vol_core; 
    1212    return vol_shell; 
     13} 
     14 
     15static double 
     16effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     17{ 
     18    if (mode == 1) { 
     19        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
     20    } else if (mode == 2) { 
     21        return 0.5 * length_a; 
     22    } else if (mode == 3) { 
     23        return 0.5 * length_a*b2a_ratio; 
     24    } else if (mode == 4) { 
     25        return 0.5 * length_a*c2a_ratio; 
     26    } else if (mode == 5) { 
     27        return length_a*sqrt(b2a_ratio/M_PI); 
     28    } else if (mode == 6) { 
     29        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     30    } else { 
     31        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     32    } 
    1333} 
    1434 

  • sasmodels/models/hollow_rectangular_prism.py

    r71b751d rd277229  
    143143source = ["lib/gauss76.c", "hollow_rectangular_prism.c"] 
    144144have_Fq = True 
    145  
    146 def ER(length_a, b2a_ratio, c2a_ratio, thickness): 
    147     """ 
    148     Return equivalent radius (ER) 
    149     thickness parameter not used 
    150     """ 
    151     b_side = length_a * b2a_ratio 
    152     c_side = length_a * c2a_ratio 
    153  
    154     # surface average radius (rough approximation) 
    155     surf_rad = sqrt(length_a * b_side / pi) 
    156  
    157     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    158     return 0.5 * (ddd) ** (1. / 3.) 
     145effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
     146                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 
     147 
     148#def ER(length_a, b2a_ratio, c2a_ratio, thickness): 
     149#    """ 
     150#    Return equivalent radius (ER) 
     151#    thickness parameter not used 
     152#    """ 
     153#    b_side = length_a * b2a_ratio 
     154#    c_side = length_a * c2a_ratio 
     155# 
     156#    # surface average radius (rough approximation) 
     157#    surf_rad = sqrt(length_a * b_side / pi) 
     158# 
     159#    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
     160#    return 0.5 * (ddd) ** (1. / 3.) 
    159161 
    160162def VR(length_a, b2a_ratio, c2a_ratio, thickness): 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    r71b751d rd277229  
    66    double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
    77    return vol_shell; 
     8} 
     9 
     10static double 
     11effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio) 
     12{ 
     13    if (mode == 1) { 
     14        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
     15    } else if (mode == 2) { 
     16        return 0.5 * length_a; 
     17    } else if (mode == 3) { 
     18        return 0.5 * length_a*b2a_ratio; 
     19    } else if (mode == 4) { 
     20        return 0.5 * length_a*c2a_ratio; 
     21    } else if (mode == 5) { 
     22        return length_a*sqrt(b2a_ratio/M_PI); 
     23    } else if (mode == 6) { 
     24        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     25    } else { 
     26        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     27    } 
    828} 
    929 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    r71b751d rd277229  
    103103source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"] 
    104104have_Fq = True 
     105effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
     106                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 
    105107 
    106 def ER(length_a, b2a_ratio, c2a_ratio): 
    107     """ 
    108         Return equivalent radius (ER) 
    109     """ 
    110     b_side = length_a * b2a_ratio 
    111     c_side = length_a * c2a_ratio 
    112  
    113     # surface average radius (rough approximation) 
    114     surf_rad = sqrt(length_a * b_side / pi) 
    115  
    116     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    117     return 0.5 * (ddd) ** (1. / 3.) 
     108#def ER(length_a, b2a_ratio, c2a_ratio): 
     109#    """ 
     110#        Return equivalent radius (ER) 
     111#    """ 
     112#    b_side = length_a * b2a_ratio 
     113#    c_side = length_a * c2a_ratio 
     114# 
     115#    # surface average radius (rough approximation) 
     116#    surf_rad = sqrt(length_a * b_side / pi) 
     117# 
     118#    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
     119#    return 0.5 * (ddd) ** (1. / 3.) 
    118120 
    119121def VR(length_a, b2a_ratio, c2a_ratio): 

  • sasmodels/models/mono_gauss_coil.py

    r2d81cfe rd277229  
    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    ] 
     73 
     74source = ["mono_gauss_coil.c"] 
     75have_Fq = False 
     76effective_radius_type = ["R_g","2R_g","3R_g","(5/3)^0.5*R_g"] 
     77 
     78 
    7379# pylint: enable=bad-whitespace, line-too-long 
    7480 
    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 
    79  
    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 
     81## NB: Scale and Background are implicit parameters on every model 
     82#def Iq(q, i_zero, rg): 
     83#    # pylint: disable = missing-docstring 
     84#    z = (q * rg)**2 
     85# 
     86#    with errstate(invalid='ignore'): 
     87#        inten = (i_zero * 2.0) * (exp(-z) + z - 1.0)/z**2 
     88#        inten[q == 0] = i_zero 
     89#    return inten 
     90#Iq.vectorized = True # Iq accepts an array of q values 
    8591 
    8692def random(): 

  • sasmodels/models/multilayer_vesicle.c

    r71b751d rd277229  
    4747} 
    4848 
     49static double 
     50effective_radius(int mode, double radius, double thick_shell, double thick_solvent, double fp_n_shells) 
     51{ 
     52    return radius + fp_n_shells*thick_shell + (fp_n_shells - 1.0)*thick_solvent; 
     53} 
     54 
    4955static void 
    5056Fq(double q, 

  • sasmodels/models/multilayer_vesicle.py

    r71b751d rd277229  
    146146source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
    147147have_Fq = True 
     148effective_radius_type = ["outer radius"] 
    148149 
    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 
     150#def ER(radius, thick_shell, thick_solvent, n_shells): 
     151#    n_shells = int(n_shells+0.5) 
     152#    return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent 
    152153 
    153154def random(): 

  • sasmodels/models/onion.c

    r71b751d rd277229  
    4040} 
    4141 
     42static double 
     43effective_radius(int mode, double radius_core, double n_shells, double thickness[]) 
     44{ 
     45  int n = (int)(n_shells+0.5); 
     46  double r = radius_core; 
     47  for (int i=0; i < n; i++) { 
     48    r += thickness[i]; 
     49  } 
     50  return r; 
     51} 
     52 
    4253static void 
    4354Fq(double q, double *F1, double *F2, double sld_core, double radius_core, double sld_solvent, 

  • sasmodels/models/onion.py

    r71b751d rd277229  
    316316single = False 
    317317have_Fq = True 
     318effective_radius_type = ["outer radius"] 
    318319 
    319320profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    … …  
    366367    return np.asarray(z), np.asarray(rho) 
    367368 
    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 
     369#def ER(radius_core, n_shells, thickness): 
     370#    """Effective radius""" 
     371#    n = int(n_shells[0]+0.5) 
     372#    return np.sum(thickness[:n], axis=0) + radius_core 
    372373 
    373374demo = { 

  • sasmodels/models/parallelepiped.c

    r71b751d rd277229  
    33{ 
    44    return length_a * length_b * length_c; 
     5} 
     6 
     7static double 
     8effective_radius(int mode, double length_a, double length_b, double length_c) 
     9{ 
     10    if (mode == 1) { 
     11        return cbrt(0.75*length_a*length_b*length_c/M_PI); 
     12    } else if (mode == 2) { 
     13        return 0.5 * length_a; 
     14    } else if (mode == 3) { 
     15        return 0.5 * length_b; 
     16    } else if (mode == 4) { 
     17        return 0.5 * length_c; 
     18    } else if (mode == 5) { 
     19        return sqrt(length_a*length_b/M_PI); 
     20    } else if (mode == 6) { 
     21        return 0.5*sqrt(length_a*length_a + length_b*length_b); 
     22    } else { 
     23        return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c); 
     24    } 
    525} 
    626 

  • sasmodels/models/parallelepiped.py

    r71b751d rd277229  
    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.) 
     233effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
     234                         "equivalent circular cross-section","half ab diagonal","half diagonal"] 
     235 
     236#def ER(length_a, length_b, length_c): 
     237#    """ 
     238#    Return effective radius (ER) for P(q)*S(q) 
     239#    """ 
     240#    # now that axes can be in any size order, need to sort a,b,c 
     241#    # where a~b and c is either much smaller or much larger 
     242#    abc = np.vstack((length_a, length_b, length_c)) 
     243#    abc = np.sort(abc, axis=0) 
     244#    selector = (abc[1] - abc[0]) > (abc[2] - abc[1]) 
     245#    length = np.where(selector, abc[0], abc[2]) 
     246#    # surface average radius (rough approximation) 
     247#    radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi) 
     248# 
     249#    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius)) 
     250#    return 0.5 * (ddd) ** (1. / 3.) 
    249251 
    250252# VR defaults to 1.0 

  • sasmodels/models/pringle.c

    r74768cb rd277229  
    104104} 
    105105 
     106static double 
     107effective_radius(int mode, double radius, double thickness, double alpha, double beta) 
     108{ 
     109    if (mode == 1) { 
     110        return cbrt(0.75*radius*radius*thickness); 
     111    } else { 
     112        return radius; 
     113    } 
     114} 
     115 
    106116double Iq( 
    107117    double q, 

  • sasmodels/models/pringle.py

    ref07e95 rd277229  
    7474source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", \ 
    7575          "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"] 
     76effective_radius_type = ["equivalent sphere","radius"] 
    7677 
    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.) 
     78#def ER(radius, thickness, alpha, beta): 
     79#    """ 
     80#    Return equivalent radius (ER) 
     81#    """ 
     82#    ddd = 0.75 * radius * (2 * radius * thickness + (thickness + radius) \ 
     83#                           * (thickness + pi * radius)) 
     84#    return 0.5 * (ddd) ** (1. / 3.) 
    8485 
    8586def random(): 

  • sasmodels/models/raspberry.c

    r71b751d rd277229  
    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    if (mode == 1) { 
     20        return radius_lg; 
     21    } else { 
     22        return radius_lg + 2.0*radius_sm - penetration; 
     23    } 
    1424} 
    1525 

  • sasmodels/models/raspberry.py

    ref07e95 rd277229  
    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 rd277229  
    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    if (mode == 1) { 
     11        return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 
     12    } else if (mode == 2) { 
     13        return 0.5 * length_a; 
     14    } else if (mode == 3) { 
     15        return 0.5 * length_a*b2a_ratio; 
     16    } else if (mode == 4) { 
     17        return 0.5 * length_a*c2a_ratio; 
     18    } else if (mode == 5) { 
     19        return length_a*sqrt(b2a_ratio/M_PI); 
     20    } else if (mode == 6) { 
     21        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     22    } else { 
     23        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     24    } 
    525} 
    626 

  • sasmodels/models/rectangular_prism.py

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

  • sasmodels/models/sphere.py

    r71b751d rd277229  
    6767             ] 
    6868 
    69 source = ["lib/sas_3j1x_x.c"] 
     69source = ["lib/sas_3j1x_x.c","sphere.c"] 
    7070have_Fq = True 
     71effective_radius_type = ["radius"] 
    7172 
    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 
     73#def ER(radius): 
     74#    """ 
     75#    Return equivalent radius (ER) 
     76#    """ 
     77#    return radius 
    9378 
    9479# VR defaults to 1.0 

  • sasmodels/models/spherical_sld.c

    r71b751d rd277229  
    1010    } 
    1111    return M_4PI_3*cube(r); 
     12} 
     13 
     14static double 
     15effective_radius(int mode, double fp_n_shells, double thickness[], double interface[]) 
     16{ 
     17    int n_shells= (int)(fp_n_shells + 0.5); 
     18    double r = 0.0; 
     19    for (int i=0; i < n_shells; i++) { 
     20        r += thickness[i] + interface[i]; 
     21    } 
     22    return r; 
    1223} 
    1324 

  • sasmodels/models/spherical_sld.py

    r71b751d rd277229  
    210210single = False  # TODO: fix low q behaviour 
    211211have_Fq = True 
     212effective_radius_type = ["outer radius"] 
    212213 
    213214profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    … …  
    255256 
    256257 
    257 def ER(n_shells, thickness, interface): 
    258     """Effective radius""" 
    259     n_shells = int(n_shells + 0.5) 
    260     total = (np.sum(thickness[:n_shells], axis=1) 
    261              + np.sum(interface[:n_shells], axis=1)) 
    262     return total 
     258#def ER(n_shells, thickness, interface): 
     259#    """Effective radius""" 
     260#    n_shells = int(n_shells + 0.5) 
     261#    total = (np.sum(thickness[:n_shells], axis=1) 
     262#             + np.sum(interface[:n_shells], axis=1)) 
     263#    return total 
    263264 
    264265 

  • sasmodels/models/triaxial_ellipsoid.c

    r71b751d rd277229  
    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    if (mode == 1) { 
     33        return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar); 
     34    } else if (mode == 2) { 
     35        return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
     36    } else { 
     37        return radius_from_max_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
     38    } 
     39} 
    940 
    1041static void 

  • sasmodels/models/triaxial_ellipsoid.py

    r71b751d rd277229  
    158158source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 
    159159have_Fq = True 
     160effective_radius_type = ["equivalent sphere","min radius", "max radius"] 
    160161 
    161162def ER(radius_equat_minor, radius_equat_major, radius_polar): 

  • sasmodels/models/vesicle.c

    r71b751d rd277229  
    66} 
    77 
     8static double 
     9effective_radius(int mode, double radius, double thickness) 
     10{ 
     11    return radius + thickness; 
     12} 
    813 
    914static void 

  • sasmodels/models/vesicle.py

    r71b751d rd277229  
    9595source = ["lib/sas_3j1x_x.c", "vesicle.c"] 
    9696have_Fq = True 
     97effective_radius_type = ["outer radius"] 
    9798 
    98 def ER(radius, thickness): 
    99     ''' 
    100     returns the effective radius used in the S*P calculation 
    101  
    102     :param radius: core radius 
    103     :param thickness: shell thickness 
    104     ''' 
    105     return radius + thickness 
     99#def ER(radius, thickness): 
     100#    ''' 
     101#    returns the effective radius used in the S*P calculation 
     102# 
     103#    :param radius: core radius 
     104#    :param thickness: shell thickness 
     105#    ''' 
     106#    return radius + thickness 
    106107 
    107108def VR(radius, thickness): 

  • 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 r84f2962  
    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 
    … …  
    381381                continue 
    382382            self.params[p.name] = p.default 
    383             self.details[p.id] = [p.units, p.limits[0], p.limits[1]] 
     383            if p.limits and type(p.limits) is list and len(p.limits) > 1: 
     384                self.details[p.id] = [p.units if not p.choices else p.choices, p.limits[0], p.limits[1]] 
     385            else: 
     386                self.details[p.id] = [p.units if not p.choices else p.choices, None, None] 
    384387            if p.polydisperse: 
    385388                self.details[p.id+".width"] = [ 
Note: See TracChangeset for help on using the changeset viewer.