Changes in / [2773c66:c88f983] in sasmodels


Ignore:
Files:
3 deleted
65 edited

Legend:

Unmodified
Added
Removed

  • explore/beta/sasfit_compare.py

    r2a12351b r7b0abf8  
    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)) 
    211210        for i, Rei in enumerate(Re_val): 
    212211            theory = ellipsoid_theta(q, Rpk, Rei, sld, sld_solvent) 
    … …  
    344343    #Sq = Iq/Pq 
    345344    #Iq = None#= Sq = None 
    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) 
     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) 
    348347 
    349348def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): 
    … …  
    402401        radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, 
    403402        background=0, 
    404         radius_polar_pd=0.1, radius_polar_pd_type=pd_type, 
    405         radius_equatorial_pd=0.1, radius_equatorial_pd_type=pd_type, 
     403        radius_polar_pd=.1, radius_polar_pd_type=pd_type, 
     404        radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type, 
    406405        volfraction=0.15, 
    407406        radius_effective=270.7543927018, 
    408407        #radius_effective=12.59921049894873, 
    409408        ) 
    410     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     409    target = sasmodels_theory(q, model, beta_mode=1, **pars) 
    411410    actual = ellipsoid_pe(q, norm='sasview', **pars) 
    412411    title = " ".join(("sasmodels", model, pd_type)) 

  • sasmodels/core.py

    r5399809 r71b751d  
    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" 
    366               " radius_effective volfraction structure_factor_mode" 
    367               " A_sld A_sld_solvent A_radius").split() 
     365    target = ("sld sld_solvent radius length theta phi volfraction" 
     366              " beta_mode A_sld A_sld_solvent A_radius").split() 
    368367    assert target == actual, "%s != %s"%(target, actual) 
    369368 

  • sasmodels/generate.py

    r5399809 r6e45516  
    802802                              for p in partable.kernel_parameters)) 
    803803    # Define the function calls 
    804     call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, v) 0.0" 
    805804    if partable.form_volume_parameters: 
    806805        refs = _call_pars("_v.", partable.form_volume_parameters) 
    807806        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)) 
    810807    else: 
    811808        # Model doesn't have volume.  We could make the kernel run a little 
    … …  
    814811        call_volume = "#define CALL_VOLUME(v) 1.0" 
    815812    source.append(call_volume) 
    816     source.append(call_effective_radius) 
    817813    model_refs = _call_pars("_v.", partable.iq_parameters) 
    818814 

  • sasmodels/kernel.py

    r5399809 r2d81cfe  
    4141    dtype = None  # type: np.dtype 
    4242 
    43     def Iq(self, call_details, values, cutoff, magnetic): 
     43    def __call__(self, call_details, values, cutoff, magnetic): 
    4444        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    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 
     45        raise NotImplementedError("need to implement __call__") 
    8746 
    8847    def release(self): 
    8948        # type: () -> None 
    9049        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

    r296c52b r108e70e  
    6868         #define NEED_EXPM1 
    6969         #define NEED_TGAMMA 
    70          #define NEED_CBRT 
    7170         // expf missing from windows? 
    7271         #define expf exp 
    … …  
    8685#  define pown(a,b) pow(a,b) 
    8786#endif // !USE_OPENCL 
    88  
    89 #if defined(NEED_CBRT) 
    90    #define cbrt(_x) pow(_x, 0.33333333333333333333333) 
    91 #endif 
    9287 
    9388#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 
    3029//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
    3130//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
    … …  
    285284    global const double *q, // nq q values, with padding to boundary 
    286285    global double *result,  // nq+1 return values, again with padding 
    287     const double cutoff,     // cutoff in the dispersity weight product 
    288     int32_t effective_radius_type // which effective radius to compute 
     286    const double cutoff     // cutoff in the dispersity weight product 
    289287    ) 
    290288{ 
    … …  
    346344  #ifdef USE_OPENCL 
    347345    #if COMPUTE_F1_F2 
    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]); 
     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]); 
    353350    #else 
    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]); 
     351      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    357352      double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
    358353    #endif 
    359354  #else // !USE_OPENCL 
    360355    #if COMPUTE_F1_F2 
    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]); 
     356      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     357      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
    364358    #else 
    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]); 
     359      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    368360    #endif 
    369361    if (pd_start == 0) { 
    … …  
    372364      #endif 
    373365      #if COMPUTE_F1_F2 
    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; 
     366          for (int q_index=0; q_index < 2*nq+2; q_index++) result[q_index] = 0.0; 
    376367      #else 
    377368          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
    378369      #endif 
    379370    } 
    380     //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_volume, result[0]); 
     371    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    381372#endif // !USE_OPENCL 
    382373 
    … …  
    693684    // Note: weight==0 must always be excluded 
    694685    if (weight > cutoff) { 
    695       weighted_volume += weight * CALL_VOLUME(local_values.table); 
     686      pd_norm += weight * CALL_VOLUME(local_values.table); 
    696687      #if COMPUTE_F1_F2 
    697688      weight_norm += weight; 
    698689      #endif 
    699       if (effective_radius_type != 0) { 
    700         weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 
    701       } 
    702690      BUILD_ROTATION(); 
    703691 
    … …  
    757745        #ifdef USE_OPENCL 
    758746          #if COMPUTE_F1_F2 
     747            this_F1 += weight * F1; 
    759748            this_F2 += weight * F2; 
    760             this_F1 += weight * F1; 
    761749          #else 
    762750            this_result += weight * scattering; 
    … …  
    764752        #else // !USE_OPENCL 
    765753          #if COMPUTE_F1_F2 
    766             result[2*q_index+0] += weight * F2; 
    767             result[2*q_index+1] += weight * F1; 
     754            result[q_index] += weight * F2; 
     755            result[q_index+nq+1] += weight * F1; 
    768756          #else 
    769757            result[q_index] += weight * scattering; 
    … …  
    794782#ifdef USE_OPENCL 
    795783  #if COMPUTE_F1_F2 
    796     result[2*q_index+0] = this_F2; 
    797     result[2*q_index+1] = this_F1; 
     784    result[q_index] = this_F2; 
     785    result[q_index+nq+1] = this_F1; 
    798786    if (q_index == 0) { 
    799       result[2*nq+0] = weight_norm; 
    800       result[2*nq+1] = weighted_volume; 
    801       result[2*nq+2] = weighted_radius; 
     787      result[nq] = pd_norm; 
     788      result[2*nq+1] = weight_norm; 
    802789    } 
    803790  #else 
    804791    result[q_index] = this_result; 
    805     if (q_index == 0) { 
    806       result[nq+0] = weight_norm; 
    807       result[nq+1] = weighted_volume; 
    808       result[nq+2] = weighted_radius; 
    809     } 
     792    if (q_index == 0) result[nq] = pd_norm; 
    810793  #endif 
    811794 
    812 //if (q_index == 0) printf("res: %g/%g\n", result[0], weigthed_volume); 
     795//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    813796#else // !USE_OPENCL 
    814797  #if COMPUTE_F1_F2 
    815     result[2*nq] = weight_norm; 
    816     result[2*nq+1] = weighted_volume; 
    817     result[2*nq+2] = weighted_radius; 
     798    result[nq] = pd_norm; 
     799    result[2*nq+1] = weight_norm; 
    818800  #else 
    819     result[nq] = weight_norm; 
    820     result[nq+1] = weighted_volume; 
    821     result[nq+2] = weighted_radius; 
     801    result[nq] = pd_norm; 
    822802  #endif 
    823 //printf("res: %g/%g\n", result[0], weighted_volume); 
     803//printf("res: %g/%g\n", result[0], pd_norm); 
    824804#endif // !USE_OPENCL 
    825805 

  • sasmodels/kernelcl.py

    r5399809 rc036ddb  
    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 
    484483        if self.is_2d: 
    485             width = ((self.nq+15+extra_q)//16)*16 
     484            # Note: 16 rather than 15 because result is 1 longer than input. 
     485            width = ((self.nq+16)//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             width = ((self.nq+31+extra_q)//32)*32 
     490            # Note: 32 rather than 31 because result is 1 longer than input. 
     491            width = ((self.nq+32)//32)*32 
    491492            self.q = np.empty(width, dtype=dtype) 
    492493            self.q[:self.nq] = q_vectors[0] 
    … …  
    538539        self.dim = '2d' if q_input.is_2d else '1d' 
    539540        # leave room for f1/f2 results in case we need to compute beta for 1d models 
    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) 
     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) 
    543544 
    544545        # Inputs and outputs for each kernel call 
    … …  
    548549 
    549550        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    550                                   q_input.global_size[0] * nout * dtype.itemsize) 
     551                                  q_input.global_size[0] * num_returns * dtype.itemsize) 
    551552        self.q_input = q_input # allocated by GpuInput above 
    552553 
    … …  
    557558                     else np.float32)  # will never get here, so use np.float32 
    558559 
    559     def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     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): 
    560586        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    561587        context = self.queue.context 
    … …  
    571597            details_b, values_b, self.q_input.q_b, self.result_b, 
    572598            self.real(cutoff), 
    573             np.uint32(effective_radius_type), 
    574599        ] 
    575600        #print("Calling OpenCL") 

  • sasmodels/kerneldll.py

    r6e7ba14 r2f8cbb9  
    315315 
    316316        # int, int, int, int*, double*, double*, double*, double*, double 
    317         argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32] 
     317        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type] 
    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         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) 
     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) 
    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 _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    392         # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     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): 
    393417        kernel = self.kernel[1 if magnetic else 0] 
    394418        args = [ 
    … …  
    401425            self.result.ctypes.data,   # results 
    402426            self.real(cutoff), # cutoff 
    403             effective_radius_type, # cutoff 
    404427        ] 
    405428        #print("Calling DLL") 

  • sasmodels/kernelpy.py

    r5399809 r108e70e  
    1212 
    1313import numpy as np  # type: ignore 
    14 from numpy import pi 
    15 try: 
    16     from numpy import cbrt 
    17 except ImportError: 
    18     def cbrt(x): return x ** (1.0/3.0) 
    1914 
    2015from .generate import F64 
    … …  
    163158 
    164159        # Generate a closure which calls the form_volume if it exists. 
    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): 
     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): 
    177165        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    178166        if magnetic: 
    … …  
    180168        #print("Calling python kernel") 
    181169        #call_details.show(values) 
    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) 
     170        res = _loops(self._parameter_vector, self._form, self._volume, 
     171                     self.q_input.nq, call_details, values, cutoff) 
     172        return res 
    187173 
    188174    def release(self): 
    … …  
    197183           form,          # type: Callable[[], np.ndarray] 
    198184           form_volume,   # type: Callable[[], float] 
    199            form_radius,   # type: Callable[[], float] 
    200185           nq,            # type: int 
    201186           call_details,  # type: details.CallDetails 
    202187           values,        # type: np.ndarray 
    203            cutoff,        # type: float 
     188           cutoff         # type: float 
    204189          ): 
    205190    # type: (...) -> None 
    … …  
    213198    #                                                              # 
    214199    ################################################################ 
    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. 
    222200    n_pars = len(parameters) 
    223201    parameters[:] = values[2:n_pars+2] 
    224  
    225202    if call_details.num_active == 0: 
    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 
     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 
    280253 
    281254 

  • sasmodels/modelinfo.py

    r5399809 r7b9e4dd  
    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) 
    793792    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    794793    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    795794    info.source = getattr(kernel_module, 'source', []) 
    796795    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 
    800796    # TODO: check the structure of the tests 
    801797    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 
    802800    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
    803801    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    … …  
    925923    #: use those functions.  Form factors are indicated by providing 
    926924    #: 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. 
    930925    source = None           # type: List[str] 
    931926    #: The set of tests that must pass.  The format of the tests is described 
    … …  
    948943    #: each parameter set.  Multiplicity parameters will be received as 
    949944    #: 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. 
    950952    c_code = None 
    951953    #: Returns the form volume for python-based models.  Form volume is needed 

  • sasmodels/models/barbell.c

    rd277229 r71b751d  
    6262} 
    6363 
    64 static double 
    65 radius_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  
    71 static double 
    72 radius_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  
    78 static double 
    79 effective_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  
    9264static void 
    9365Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 

  • sasmodels/models/barbell.py

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

  • sasmodels/models/capped_cylinder.c

    rd277229 r71b751d  
    8484} 
    8585 
    86 static double 
    87 radius_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  
    93 static double 
    94 radius_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  
    100 static double 
    101 effective_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  
    11486static void 
    11587Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 

  • sasmodels/models/capped_cylinder.py

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

  • sasmodels/models/core_multi_shell.c

    rd277229 r71b751d  
    1919} 
    2020 
    21 static double 
    22 outer_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  
    32 static double 
    33 effective_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  
    4821static void 
    4922Fq(double q, double *F1, double *F2, double core_sld, double core_radius, 

  • sasmodels/models/core_multi_shell.py

    rd277229 r71b751d  
    100100source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 
    101101have_Fq = True 
    102 effective_radius_type = ["outer radius", "core radius"] 
    103102 
    104103def random(): 
    … …  
    145144    return np.asarray(z), np.asarray(rho) 
    146145 
    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 
     146def 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 
    151150 
    152151demo = dict(sld_core=6.4, 

  • sasmodels/models/core_shell_bicelle.c

    rd277229 r71b751d  
    3434 
    3535    return t; 
    36 } 
    37  
    38 static double 
    39 radius_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  
    45 static double 
    46 radius_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  
    53 static double 
    54 effective_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     } 
    6536} 
    6637 

  • sasmodels/models/core_shell_bicelle.py

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

  • sasmodels/models/core_shell_bicelle_elliptical.c

    rd277229 r71b751d  
    88{ 
    99    return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face); 
    10 } 
    11  
    12 static double 
    13 radius_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  
    19 static double 
    20 radius_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  
    28 static double 
    29 effective_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     } 
    4410} 
    4511 

  • sasmodels/models/core_shell_bicelle_elliptical.py

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

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    rd277229 r71b751d  
    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  
    13 static double 
    14 radius_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  
    20 static double 
    21 radius_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  
    29 static double 
    30 effective_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     } 
    4511} 
    4612 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

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

  • sasmodels/models/core_shell_cylinder.c

    rd277229 r71b751d  
    1111{ 
    1212    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
    13 } 
    14  
    15 static double 
    16 radius_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  
    22 static double 
    23 radius_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  
    30 static double 
    31 effective_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     } 
    4613} 
    4714 

  • sasmodels/models/core_shell_cylinder.py

    rd277229 r71b751d  
    126126source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 
    127127have_Fq = True 
    128 effective_radius_type = ["equivalent sphere","outer radius","half outer length","half min outer dimension", 
    129                          "half max outer dimension","half outer diagonal"] 
    130128 
    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.) 
     129def 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.) 
    139137 
    140138def VR(radius, thickness, length): 

  • sasmodels/models/core_shell_ellipsoid.c

    r3c60146 r71b751d  
    3636    double vol = M_4PI_3*equat_shell*equat_shell*polar_shell; 
    3737    return vol; 
    38 } 
    39  
    40 static double 
    41 radius_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  
    47 static double 
    48 radius_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  
    69 static double 
    70 effective_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     } 
    8438} 
    8539 

  • sasmodels/models/core_shell_ellipsoid.py

    rd277229 r71b751d  
    146146source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    147147have_Fq = True 
    148 effective_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) 
     148 
     149def 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) 
    158157 
    159158def random(): 

  • sasmodels/models/core_shell_parallelepiped.c

    rd277229 r71b751d  
    2525        2.0 * length_a * length_b * thick_rim_c; 
    2626#endif 
    27 } 
    28  
    29 static double 
    30 radius_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  
    37 static double 
    38 radius_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  
    44 static double 
    45 effective_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     } 
    6327} 
    6428 

  • sasmodels/models/core_shell_parallelepiped.py

    rd277229 r71b751d  
    227227source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] 
    228228have_Fq = True 
    229 effective_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) 
     229 
     230 
     231def 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) 
    242241 
    243242# VR defaults to 1.0 

  • sasmodels/models/core_shell_sphere.c

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

  • sasmodels/models/core_shell_sphere.py

    rd277229 r71b751d  
    7878source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 
    7979have_Fq = True 
    80 effective_radius_type = ["outer radius", "core radius"] 
    8180 
    8281demo = dict(scale=1, background=0, radius=60, thickness=10, 
    8382            sld_core=1.0, sld_shell=2.0, sld_solvent=0.0) 
    8483 
    85 #def ER(radius, thickness): 
    86 #    """ 
    87 #        Equivalent radius 
    88 #        @param radius: core radius 
    89 #        @param thickness: shell thickness 
    90 #    """ 
    91 #    return radius + thickness 
     84def ER(radius, thickness): 
     85    """ 
     86        Equivalent radius 
     87        @param radius: core radius 
     88        @param thickness: shell thickness 
     89    """ 
     90    return radius + thickness 
    9291 
    9392def VR(radius, thickness): 

  • sasmodels/models/cylinder.c

    rd277229 r71b751d  
    1111{ 
    1212    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
    13 } 
    14  
    15 static double 
    16 radius_from_volume(double radius, double length) 
    17 { 
    18     return cbrt(0.75*radius*radius*length); 
    19 } 
    20  
    21 static double 
    22 radius_from_diagonal(double radius, double length) 
    23 { 
    24     return sqrt(radius*radius + 0.25*length*length); 
    25 } 
    26  
    27 static double 
    28 effective_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     } 
    4313} 
    4414 

  • sasmodels/models/cylinder.py

    rd277229 r71b751d  
    139139source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
    140140have_Fq = True 
    141 effective_radius_type = ["equivalent sphere","radius","half length","half min dimension","half max dimension","half diagonal"] 
    142141 
    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.) 
     142def 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.) 
    149148 
    150149def random(): 

  • sasmodels/models/ellipsoid.c

    rd277229 r71b751d  
    44    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    55} 
    6  
    7 static double 
    8 radius_from_volume(double radius_polar, double radius_equatorial) 
    9 { 
    10     return cbrt(radius_polar*radius_equatorial*radius_equatorial); 
    11 } 
    12  
    13 static double 
    14 radius_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  
    33 static double 
    34 effective_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  
    476 
    487static void 

  • sasmodels/models/ellipsoid.py

    rd277229 r0168844  
    169169source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
    170170have_Fq = True 
    171 effective_radius_type = ["equivalent sphere","average curvature", "min radius", "max radius"] 
     171 
     172def 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 
    172193 
    173194def random(): 

  • sasmodels/models/elliptical_cylinder.c

    rd277229 r71b751d  
    33{ 
    44    return M_PI * radius_minor * radius_minor * r_ratio * length; 
    5 } 
    6  
    7 static double 
    8 radius_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  
    14 static double 
    15 radius_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  
    21 static double 
    22 radius_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  
    28 static double 
    29 radius_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  
    35 static double 
    36 effective_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     } 
    575} 
    586 

  • sasmodels/models/elliptical_cylinder.py

    rd277229 r71b751d  
    123123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
    124124have_Fq = True 
    125 effective_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"] 
    127125 
    128126demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 
    … …  
    130128            theta_pd=10, phi_pd=2, psi_pd=3) 
    131129 
    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.) 
     130def 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.) 
    143141 
    144142def random(): 

  • sasmodels/models/fuzzy_sphere.py

    rd277229 r71b751d  
    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],    "volume",       "std deviation of Gaussian convolution for interface (must be << radius)"], 
     80              ["fuzziness",   "Ang",        10, [0, inf],    "",       "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","fuzzy_sphere.c"] 
     84source = ["lib/sas_3j1x_x.c"] 
    8585have_Fq = True 
    86 effective_radius_type = ["radius","radius + fuzziness"] 
    8786 
    88 #def ER(radius): 
    89 #    """ 
    90 #    Return radius 
    91 #    """ 
    92 #    return radius 
     87c_code = """ 
     88static double form_volume(double radius) 
     89{ 
     90    return M_4PI_3*cube(radius); 
     91} 
     92 
     93static 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 
     106def ER(radius): 
     107    """ 
     108    Return radius 
     109    """ 
     110    return radius 
    93111 
    94112# VR defaults to 1.0 

  • sasmodels/models/hollow_cylinder.c

    rd277229 r71b751d  
    2222} 
    2323 
    24 static double 
    25 radius_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  
    31 static double 
    32 radius_from_diagonal(double radius, double thickness, double length) 
    33 { 
    34     return sqrt(square(radius + thickness) + 0.25*square(length)); 
    35 } 
    36  
    37 static double 
    38 effective_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 } 
    5424 
    5525static void 

  • sasmodels/models/hollow_cylinder.py

    rd277229 r71b751d  
    9090source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    9191have_Fq = True 
    92 effective_radius_type = ["equivalent sphere","outer radius","half length", 
    93                          "half outer min dimension","half outer max dimension","half outer diagonal"] 
    9492 
    9593# pylint: disable=W0613 
    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 
     94def 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 
    113111 
    114112def VR(radius, thickness, length): 

  • sasmodels/models/hollow_rectangular_prism.c

    rd277229 r71b751d  
    1111    double vol_shell = vol_total - vol_core; 
    1212    return vol_shell; 
    13 } 
    14  
    15 static double 
    16 effective_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     } 
    3313} 
    3414 

  • sasmodels/models/hollow_rectangular_prism.py

    rd277229 r71b751d  
    143143source = ["lib/gauss76.c", "hollow_rectangular_prism.c"] 
    144144have_Fq = True 
    145 effective_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.) 
     145 
     146def 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.) 
    161159 
    162160def VR(length_a, b2a_ratio, c2a_ratio, thickness): 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rd277229 r71b751d  
    66    double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
    77    return vol_shell; 
    8 } 
    9  
    10 static double 
    11 effective_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     } 
    288} 
    299 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.py

    rd277229 r71b751d  
    103103source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"] 
    104104have_Fq = True 
    105 effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
    106                          "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 
    107105 
    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.) 
     106def 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.) 
    120118 
    121119def VR(length_a, b2a_ratio, c2a_ratio): 

  • sasmodels/models/mono_gauss_coil.py

    rd277229 r2d81cfe  
    5454 
    5555import numpy as np 
    56 from numpy import inf 
     56from numpy import inf, exp, errstate 
    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], "volume", "Radius of gyration"], 
     71    ["rg", "Ang", 75.0, [0.0, inf], "", "Radius of gyration"], 
    7272    ] 
    73  
    74 source = ["mono_gauss_coil.c"] 
    75 have_Fq = False 
    76 effective_radius_type = ["R_g","2R_g","3R_g","(5/3)^0.5*R_g"] 
    77  
    78  
    7973# pylint: enable=bad-whitespace, line-too-long 
    8074 
    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 
     75# NB: Scale and Background are implicit parameters on every model 
     76def 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 
     84Iq.vectorized = True # Iq accepts an array of q values 
    9185 
    9286def random(): 

  • sasmodels/models/multilayer_vesicle.c

    rd277229 r71b751d  
    4747} 
    4848 
    49 static double 
    50 effective_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  
    5549static void 
    5650Fq(double q, 

  • sasmodels/models/multilayer_vesicle.py

    rd277229 r71b751d  
    146146source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
    147147have_Fq = True 
    148 effective_radius_type = ["outer radius"] 
    149148 
    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 
     149def 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 
    153152 
    154153def random(): 

  • sasmodels/models/onion.c

    rd277229 r71b751d  
    4040} 
    4141 
    42 static double 
    43 effective_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  
    5342static void 
    5443Fq(double q, double *F1, double *F2, double sld_core, double radius_core, double sld_solvent, 

  • sasmodels/models/onion.py

    rd277229 r71b751d  
    316316single = False 
    317317have_Fq = True 
    318 effective_radius_type = ["outer radius"] 
    319318 
    320319profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    … …  
    367366    return np.asarray(z), np.asarray(rho) 
    368367 
    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 
     368def 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 
    373372 
    374373demo = { 

  • sasmodels/models/parallelepiped.c

    rd277229 r71b751d  
    33{ 
    44    return length_a * length_b * length_c; 
    5 } 
    6  
    7 static double 
    8 effective_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     } 
    255} 
    266 

  • sasmodels/models/parallelepiped.py

    rd277229 r71b751d  
    231231source = ["lib/gauss76.c", "parallelepiped.c"] 
    232232have_Fq = True 
    233 effective_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.) 
     233 
     234def 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.) 
    251249 
    252250# VR defaults to 1.0 

  • sasmodels/models/pringle.c

    rd277229 r74768cb  
    104104} 
    105105 
    106 static double 
    107 effective_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  
    116106double Iq( 
    117107    double q, 

  • sasmodels/models/pringle.py

    rd277229 ref07e95  
    7474source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", \ 
    7575          "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"] 
    76 effective_radius_type = ["equivalent sphere","radius"] 
    7776 
    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.) 
     77def 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.) 
    8584 
    8685def random(): 

  • sasmodels/models/raspberry.c

    rd277229 r71b751d  
    1 double form_volume(double radius_lg, double radius_sm, double penetration); 
     1double form_volume(double radius_lg); 
    22 
    33double Iq(double q, 
    … …  
    66          double radius_lg, double radius_sm, double penetration); 
    77 
    8 double form_volume(double radius_lg, double radius_sm, double penetration) 
     8double form_volume(double radius_lg) 
    99{ 
    1010    //Because of the complex structure, volume normalization must 
    … …  
    1212    double volume=1.0; 
    1313    return volume; 
    14 } 
    15  
    16 static double 
    17 effective_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     } 
    2414} 
    2515 

  • sasmodels/models/raspberry.py

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

  • sasmodels/models/rectangular_prism.c

    rd277229 r71b751d  
    33{ 
    44    return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio); 
    5 } 
    6  
    7 static double 
    8 effective_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     } 
    255} 
    266 

  • sasmodels/models/rectangular_prism.py

    rd277229 r71b751d  
    136136source = ["lib/gauss76.c", "rectangular_prism.c"] 
    137137have_Fq = True 
    138 effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 
    139                          "equivalent circular cross-section","half ab diagonal","half diagonal"] 
    140138 
    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.) 
     139def 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.) 
    153151 
    154152def random(): 

  • sasmodels/models/sphere.py

    rd277229 r71b751d  
    6767             ] 
    6868 
    69 source = ["lib/sas_3j1x_x.c","sphere.c"] 
     69source = ["lib/sas_3j1x_x.c"] 
    7070have_Fq = True 
    71 effective_radius_type = ["radius"] 
    7271 
    73 #def ER(radius): 
    74 #    """ 
    75 #    Return equivalent radius (ER) 
    76 #    """ 
    77 #    return radius 
     72c_code = """ 
     73static double form_volume(double radius) 
     74{ 
     75    return M_4PI_3*cube(radius); 
     76} 
     77 
     78static 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 
     88def ER(radius): 
     89    """ 
     90    Return equivalent radius (ER) 
     91    """ 
     92    return radius 
    7893 
    7994# VR defaults to 1.0 

  • sasmodels/models/spherical_sld.c

    rd277229 r71b751d  
    1010    } 
    1111    return M_4PI_3*cube(r); 
    12 } 
    13  
    14 static double 
    15 effective_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; 
    2312} 
    2413 

  • sasmodels/models/spherical_sld.py

    rd277229 r71b751d  
    210210single = False  # TODO: fix low q behaviour 
    211211have_Fq = True 
    212 effective_radius_type = ["outer radius"] 
    213212 
    214213profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    … …  
    256255 
    257256 
    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 
     257def 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 
    264263 
    265264 

  • sasmodels/models/triaxial_ellipsoid.c

    rd277229 r71b751d  
    77} 
    88 
    9 static double 
    10 radius_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  
    15 static double 
    16 radius_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  
    22 static double 
    23 radius_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  
    29 static double 
    30 effective_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 } 
    409 
    4110static void 

  • sasmodels/models/triaxial_ellipsoid.py

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

  • sasmodels/models/vesicle.c

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

  • sasmodels/models/vesicle.py

    rd277229 r71b751d  
    9595source = ["lib/sas_3j1x_x.c", "vesicle.c"] 
    9696have_Fq = True 
    97 effective_radius_type = ["outer radius"] 
    9897 
    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 
     98def 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 
    107106 
    108107def VR(radius, thickness): 

  • sasmodels/product.py

    r33d7be3 rc0131d44  
    2424# pylint: disable=unused-import 
    2525try: 
    26     from typing import Tuple, Callable, Union 
     26    from typing import Tuple, Callable 
    2727except ImportError: 
    2828    pass 
    … …  
    3737#] 
    3838 
    39 RADIUS_ID = "radius_effective" 
    40 VOLFRAC_ID = "volfraction" 
     39ER_ID = "radius_effective" 
     40VF_ID = "volfraction" 
     41 
    4142def make_extra_pars(p_info): 
    4243    pars = [] 
    … …  
    5051                "Structure factor calculation") 
    5152        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) 
    6153    return pars 
    6254 
    … …  
    7365    # have any magnetic parameters 
    7466    if not len(s_info.parameters.kernel_parameters) >= 2: 
    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)) 
     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)) 
    8072    if not s_info.parameters.magnetism_index == []: 
    8173        raise TypeError("S should not have SLD parameters") 
    … …  
    8375    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    8476 
    85     # Create list of parameters for the combined model.  If there 
     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 
    8679    # are any names in P that overlap with those in S, modify the name in S 
    8780    # to distinguish it. 
    8881    p_set = set(p.id for p in p_pars.kernel_parameters) 
    8982    s_list = [(_tag_parameter(par) if par.id in p_set else par) 
    90               for par in s_pars.kernel_parameters] 
     83              for par in s_pars.kernel_parameters[1:]] 
    9184    # Check if still a collision after renaming.  This could happen if for 
    9285    # example S has volfrac and P has both volfrac and volfrac_S. 
    … …  
    9487        raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 
    9588 
    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  
    10089    translate_name = dict((old.id, new.id) for old, new 
    101                           in zip(s_pars.kernel_parameters, s_list)) 
     90                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    10291    combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info) 
    10392    parameters = ParameterTable(combined_pars) 
    … …  
    10594    def random(): 
    10695        combined_pars = p_info.random() 
    107         s_names = set(par.id for par in s_pars.kernel_parameters) 
     96        s_names = set(par.id for par in s_pars.kernel_parameters[1:]) 
    10897        combined_pars.update((translate_name[k], v) 
    10998                             for k, v in s_info.random().items() 
    … …  
    168157    return par 
    169158 
    170 def _intermediates(P, S, effective_radius): 
    171     # type: (np.ndarray, np.ndarray, float) -> OrderedDict[str, np.ndarray] 
     159def _intermediates(P, S): 
     160    # type: (np.ndarray, np.ndarray) -> OrderedDict[str, np.ndarray] 
    172161    """ 
    173162    Returns intermediate results for standard product (P(Q)*S(Q)) 
    … …  
    176165        ("P(Q)", P), 
    177166        ("S(Q)", S), 
    178         ("effective_radius", effective_radius), 
    179167    )) 
    180168 
    181 def _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. 
     169def _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 
    191173    """ 
    192174    # TODO: 1. include calculated Q vector 
    … …  
    197179        ("beta(Q)", F1**2 / F2), 
    198180        ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), 
    199         ("effective_radius", effective_radius), 
    200181        # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), 
    201182    )) 
    … …  
    281262        p_values = np.hstack(p_values).astype(self.p_kernel.dtype) 
    282263 
     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 
    283269        # Construct the calling parameters for S. 
    284         s_npars = s_info.parameters.npars 
     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 
    285277        s_length = call_details.length[p_npars:p_npars+s_npars] 
    286278        s_offset = call_details.offset[p_npars:p_npars+s_npars] 
    … …  
    288280        s_offset = np.hstack((nweights, s_offset)) 
    289281        s_details = make_details(s_info, s_length, s_offset, nweights+1) 
     282        v, w = weights[:nweights], weights[nweights:] 
    290283        s_values = [ 
    291             # scale=1, background=0, 
    292             [1., 0.], 
    293             values[2+p_npars:2+p_npars+s_npars], 
    294             weights, 
     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] 
    295294        ] 
    296295        spacer = (32 - sum(len(v) for v in s_values)%32)%32 
    … …  
    299298 
    300299        # beta mode is the first parameter after the structure factor pars 
    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 
     300        beta_index = 2+p_npars+s_npars 
     301        beta_mode = values[beta_index] 
    312302 
    313303        # Call the kernels 
     304        s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    314305        scale, background = values[0], values[1] 
    315306        if beta_mode: 
    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) 
     307            F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic) 
    322308            combined_scale = scale*volfrac/volume_avg 
    323309 
    324             self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background, effective_radius) 
     310            self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background) 
    325311            final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background 
    326312 
    327313        else: 
    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) 
     314            p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic) 
     315 
     316            self.results = lambda: _intermediates(p_result, s_result) 
    336317            final_result = scale*(p_result*s_result) + background 
    337318 

  • 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.effective_radius_type is not None 
     274    attrs['is_form_factor'] = model_info.ER 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             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] 
     383            self.details[p.id] = [p.units, p.limits[0], p.limits[1]] 
    387384            if p.polydisperse: 
    388385                self.details[p.id+".width"] = [ 
Note: See TracChangeset for help on using the changeset viewer.