Changes in / [db472a5:0159b5e] in sasmodels
- Files:
-
- 4 deleted
- 68 edited
Legend:
- Unmodified
- Added
- Removed
-
explore/beta/sasfit_compare.py
r2a12351b r7b0abf8 208 208 F1, F2 = np.zeros_like(q), np.zeros_like(q) 209 209 for k, Rpk in enumerate(Rp_val): 210 print("ellipsoid cycle", k, "of", len(Rp_val))211 210 for i, Rei in enumerate(Re_val): 212 211 theory = ellipsoid_theta(q, Rpk, Rei, sld, sld_solvent) … … 344 343 #Sq = Iq/Pq 345 344 #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) 348 347 349 348 def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): … … 402 401 radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, 403 402 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, 406 405 volfraction=0.15, 407 406 radius_effective=270.7543927018, 408 407 #radius_effective=12.59921049894873, 409 408 ) 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) 411 410 actual = ellipsoid_pe(q, norm='sasview', **pars) 412 411 title = " ".join(("sasmodels", model, pd_type)) -
explore/precision.py
ree60aa7 r2a7e20e 357 357 ) 358 358 add_function( 359 name=" gauss_coil",359 name="debye", 360 360 mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4, 361 361 np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4, 362 362 ocl_function=make_ocl(""" 363 363 const double qsq = q*q; 364 // For double: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide) 365 // For single: use O(7) Taylor with 0.8 cutoff (7 mad) 366 if (qsq < 0.0) { 364 if (qsq < 1.0) { // Pade approximation 367 365 const double x = qsq; 368 366 if (0) { // 0.36 single … … 374 372 const double B1=3./8., B2=3./56., B3=1./336.; 375 373 return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.); 376 } else if ( 0) { // 1.0 for single, 0.25 for double374 } else if (1) { // 1.0 for single, 0.25 for double 377 375 // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 378 376 const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; … … 387 385 /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 388 386 } 389 } else if (qsq < 0.8) {387 } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double 390 388 const double x = qsq; 391 389 const double C0 = +1.; -
sasmodels/core.py
ree60aa7 r71b751d 363 363 model = load_model("cylinder@hardsphere*sphere") 364 364 actual = [p.name for p in model.info.parameters.kernel_parameters] 365 target = ("sld sld_solvent radius length theta phi" 366 " radius_effective volfraction " 367 " structure_factor_mode radius_effective_mode" 368 " A_sld A_sld_solvent A_radius").split() 365 target = ("sld sld_solvent radius length theta phi volfraction" 366 " beta_mode A_sld A_sld_solvent A_radius").split() 369 367 assert target == actual, "%s != %s"%(target, actual) 370 368 -
sasmodels/generate.py
r5399809 r6e45516 802 802 for p in partable.kernel_parameters)) 803 803 # Define the function calls 804 call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, v) 0.0"805 804 if partable.form_volume_parameters: 806 805 refs = _call_pars("_v.", partable.form_volume_parameters) 807 806 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))810 807 else: 811 808 # Model doesn't have volume. We could make the kernel run a little … … 814 811 call_volume = "#define CALL_VOLUME(v) 1.0" 815 812 source.append(call_volume) 816 source.append(call_effective_radius)817 813 model_refs = _call_pars("_v.", partable.iq_parameters) 818 814 -
sasmodels/kernel.py
ree60aa7 r2d81cfe 41 41 dtype = None # type: np.dtype 42 42 43 def Iq(self, call_details, values, cutoff, magnetic):43 def __call__(self, call_details, values, cutoff, magnetic): 44 44 # 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__") 87 46 88 47 def release(self): 89 48 # type: () -> None 90 49 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.ndarray94 raise NotImplementedError() -
sasmodels/kernel_header.c
r296c52b r108e70e 68 68 #define NEED_EXPM1 69 69 #define NEED_TGAMMA 70 #define NEED_CBRT71 70 // expf missing from windows? 72 71 #define expf exp … … 86 85 # define pown(a,b) pow(a,b) 87 86 #endif // !USE_OPENCL 88 89 #if defined(NEED_CBRT)90 #define cbrt(_x) pow(_x, 0.33333333333333333333333)91 #endif92 87 93 88 #if defined(NEED_EXPM1) -
sasmodels/kernel_iq.c
r70530778 r70530778 27 27 // parameters in the parameter table. 28 28 // CALL_VOLUME(table) : call the form volume function 29 // CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function30 29 // CALL_IQ(q, table) : call the Iq function for 1D calcs. 31 30 // CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. … … 285 284 global const double *q, // nq q values, with padding to boundary 286 285 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 289 287 ) 290 288 { … … 346 344 #ifdef USE_OPENCL 347 345 #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]); 353 350 #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]); 357 352 double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 358 353 #endif 359 354 #else // !USE_OPENCL 360 355 #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]); 364 358 #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]); 368 360 #endif 369 361 if (pd_start == 0) { … … 372 364 #endif 373 365 #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; 376 367 #else 377 368 for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 378 369 #endif 379 370 } 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]); 381 372 #endif // !USE_OPENCL 382 373 … … 693 684 // Note: weight==0 must always be excluded 694 685 if (weight > cutoff) { 695 weighted_volume+= weight * CALL_VOLUME(local_values.table);686 pd_norm += weight * CALL_VOLUME(local_values.table); 696 687 #if COMPUTE_F1_F2 697 688 weight_norm += weight; 698 689 #endif 699 if (effective_radius_type != 0) {700 weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table);701 }702 690 BUILD_ROTATION(); 703 691 … … 757 745 #ifdef USE_OPENCL 758 746 #if COMPUTE_F1_F2 747 this_F1 += weight * F1; 759 748 this_F2 += weight * F2; 760 this_F1 += weight * F1;761 749 #else 762 750 this_result += weight * scattering; … … 764 752 #else // !USE_OPENCL 765 753 #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; 768 756 #else 769 757 result[q_index] += weight * scattering; … … 794 782 #ifdef USE_OPENCL 795 783 #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; 798 786 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; 802 789 } 803 790 #else 804 791 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; 810 793 #endif 811 794 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); 813 796 #else // !USE_OPENCL 814 797 #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; 818 800 #else 819 result[nq] = weight_norm; 820 result[nq+1] = weighted_volume; 821 result[nq+2] = weighted_radius; 801 result[nq] = pd_norm; 822 802 #endif 823 //printf("res: %g/%g\n", result[0], weighted_volume);803 //printf("res: %g/%g\n", result[0], pd_norm); 824 804 #endif // !USE_OPENCL 825 805 -
sasmodels/kernelcl.py
r5399809 rc036ddb 481 481 # at this point, so instead using 32, which is good on the set of 482 482 # architectures tested so far. 483 extra_q = 3 # total weight, weighted volume and weighted radius484 483 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 486 486 self.q = np.empty((width, 2), dtype=dtype) 487 487 self.q[:self.nq, 0] = q_vectors[0] 488 488 self.q[:self.nq, 1] = q_vectors[1] 489 489 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 491 492 self.q = np.empty(width, dtype=dtype) 492 493 self.q[:self.nq] = q_vectors[0] … … 538 539 self.dim = '2d' if q_input.is_2d else '1d' 539 540 # leave room for f1/f2 results in case we need to compute beta for 1d models 540 n out = 2 if self.info.have_Fq and self.dim == '1d' else 1541 # plus 3 weight, volume, radius542 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) 543 544 544 545 # Inputs and outputs for each kernel call … … 548 549 549 550 self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 550 q_input.global_size[0] * n out* dtype.itemsize)551 q_input.global_size[0] * num_returns * dtype.itemsize) 551 552 self.q_input = q_input # allocated by GpuInput above 552 553 … … 557 558 else np.float32) # will never get here, so use np.float32 558 559 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): 560 586 # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 561 587 context = self.queue.context … … 571 597 details_b, values_b, self.q_input.q_b, self.result_b, 572 598 self.real(cutoff), 573 np.uint32(effective_radius_type),574 599 ] 575 600 #print("Calling OpenCL") -
sasmodels/kerneldll.py
r6e7ba14 r2f8cbb9 315 315 316 316 # 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] 318 318 names = [generate.kernel_name(self.info, variant) 319 319 for variant in ("Iq", "Iqxy", "Imagnetic")] … … 382 382 self.dim = '2d' if q_input.is_2d else '1d' 383 383 # leave room for f1/f2 results in case we need to compute beta for 1d models 384 n out = 2 if self.info.have_Fq else 1385 # plus 3 weight, volume, radius386 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) 387 387 self.real = (np.float32 if self.q_input.dtype == generate.F32 388 388 else np.float64 if self.q_input.dtype == generate.F64 389 389 else np.float128) 390 390 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): 393 417 kernel = self.kernel[1 if magnetic else 0] 394 418 args = [ … … 401 425 self.result.ctypes.data, # results 402 426 self.real(cutoff), # cutoff 403 effective_radius_type, # cutoff404 427 ] 405 428 #print("Calling DLL") -
sasmodels/kernelpy.py
r5399809 r108e70e 12 12 13 13 import numpy as np # type: ignore 14 from numpy import pi15 try:16 from numpy import cbrt17 except ImportError:18 def cbrt(x): return x ** (1.0/3.0)19 14 20 15 from .generate import F64 … … 163 158 164 159 # 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): 177 165 # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 178 166 if magnetic: … … 180 168 #print("Calling python kernel") 181 169 #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 187 173 188 174 def release(self): … … 197 183 form, # type: Callable[[], np.ndarray] 198 184 form_volume, # type: Callable[[], float] 199 form_radius, # type: Callable[[], float]200 185 nq, # type: int 201 186 call_details, # type: details.CallDetails 202 187 values, # type: np.ndarray 203 cutoff ,# type: float188 cutoff # type: float 204 189 ): 205 190 # type: (...) -> None … … 213 198 # # 214 199 ################################################################ 215 216 # WARNING: Trickery ahead217 # The parameters[] vector is embedded in the closures for form(),218 # form_volume() and form_radius(). We set the initial vector from219 # the values for the model parameters. As we loop through the polydispesity220 # mesh, we update the components with the polydispersity values before221 # calling the respective functions.222 200 n_pars = len(parameters) 223 201 parameters[:] = values[2:n_pars+2] 224 225 202 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 280 253 281 254 -
sasmodels/modelinfo.py
r5399809 r7b9e4dd 790 790 info.structure_factor = getattr(kernel_module, 'structure_factor', False) 791 791 # TODO: find Fq by inspection 792 info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None)793 792 info.have_Fq = getattr(kernel_module, 'have_Fq', False) 794 793 info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 795 794 info.source = getattr(kernel_module, 'source', []) 796 795 info.c_code = getattr(kernel_module, 'c_code', None) 797 info.effective_radius = getattr(kernel_module, 'effective_radius', None)798 info.ER = None # CRUFT799 info.VR = None # CRUFT800 796 # TODO: check the structure of the tests 801 797 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 802 800 info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 803 801 info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore … … 925 923 #: use those functions. Form factors are indicated by providing 926 924 #: 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.930 925 source = None # type: List[str] 931 926 #: The set of tests that must pass. The format of the tests is described … … 948 943 #: each parameter set. Multiplicity parameters will be received as 949 944 #: 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. 950 952 c_code = None 951 953 #: Returns the form volume for python-based models. Form volume is needed -
sasmodels/models/barbell.c
ree60aa7 r71b751d 62 62 } 63 63 64 static double65 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 double72 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 double79 effective_radius(int mode, double radius_bell, double radius, double length)80 {81 switch (mode) {82 case 1: // equivalent sphere83 return radius_from_volume(radius_bell, radius , length);84 case 2: // radius85 return radius;86 case 3: // half length87 return 0.5*length;88 case 4: // half total length89 return radius_from_totallength(radius_bell,radius,length);90 }91 }92 93 64 static void 94 65 Fq(double q,double *F1, double *F2, double sld, double solvent_sld, -
sasmodels/models/barbell.py
ree60aa7 r71b751d 117 117 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 118 118 have_Fq = True 119 effective_radius_type = [120 "equivalent sphere", "radius", "half length", "half total length",121 ]122 119 123 120 def random(): -
sasmodels/models/capped_cylinder.c
ree60aa7 r71b751d 84 84 } 85 85 86 static double87 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 double94 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 double101 effective_radius(int mode, double radius, double radius_cap, double length)102 {103 switch (mode) {104 case 1: // equivalent sphere105 return radius_from_volume(radius, radius_cap, length);106 case 2: // radius107 return radius;108 case 3: // half length109 return 0.5*length;110 case 4: // half total length111 return radius_from_totallength(radius, radius_cap,length);112 }113 }114 115 86 static void 116 87 Fq(double q,double *F1, double *F2, double sld, double solvent_sld, -
sasmodels/models/capped_cylinder.py
ree60aa7 r71b751d 137 137 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] 138 138 have_Fq = True 139 effective_radius_type = [140 "equivalent sphere", "radius", "half length", "half total length",141 ]142 139 143 140 def random(): -
sasmodels/models/core_multi_shell.c
ree60aa7 r71b751d 9 9 10 10 static double 11 outer_radius(double core_radius, double fp_n, double thickness[])11 form_volume(double core_radius, double fp_n, double thickness[]) 12 12 { 13 13 double r = core_radius; … … 16 16 r += thickness[i]; 17 17 } 18 return r; 19 } 20 21 static double 22 form_volume(double core_radius, double fp_n, double thickness[]) 23 { 24 return M_4PI_3 * cube(outer_radius(core_radius, fp_n, thickness)); 25 } 26 27 static double 28 effective_radius(int mode, double core_radius, double fp_n, double thickness[]) 29 { 30 switch (mode) { 31 case 1: // outer radius 32 return outer_radius(core_radius, fp_n, thickness); 33 case 2: // core radius 34 return core_radius; 35 } 18 return M_4PI_3 * cube(r); 36 19 } 37 20 -
sasmodels/models/core_multi_shell.py
ree60aa7 r71b751d 100 100 source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 101 101 have_Fq = True 102 effective_radius_type = ["outer radius", "core radius"]103 102 104 103 def random(): … … 145 144 return np.asarray(z), np.asarray(rho) 146 145 146 def ER(radius, n, thickness): 147 """Effective radius""" 148 n = int(n[0]+0.5) # n is a control parameter and is not polydisperse 149 return np.sum(thickness[:n], axis=0) + radius 150 147 151 demo = dict(sld_core=6.4, 148 152 radius=60, -
sasmodels/models/core_shell_bicelle.c
ree60aa7 r71b751d 34 34 35 35 return t; 36 }37 38 static double39 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 double46 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 double54 effective_radius(int mode, double radius, double thick_rim, double thick_face, double length)55 {56 switch (mode) {57 case 1: // equivalent sphere58 return radius_from_volume(radius, thick_rim, thick_face, length);59 case 2: // outer rim radius60 return radius + thick_rim;61 case 3: // half outer thickness62 return 0.5*length + thick_face;63 case 4: // half diagonal64 return radius_from_diagonal(radius,thick_rim,thick_face,length);65 }66 36 } 67 37 -
sasmodels/models/core_shell_bicelle.py
ree60aa7 r71b751d 155 155 "core_shell_bicelle.c"] 156 156 have_Fq = True 157 effective_radius_type = [158 "equivalent sphere", "outer rim radius",159 "half outer thickness", "half diagonal",160 ]161 157 162 158 def random(): -
sasmodels/models/core_shell_bicelle_elliptical.c
ree60aa7 r71b751d 8 8 { 9 9 return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face); 10 }11 12 static double13 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 double20 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 double29 effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length)30 {31 switch (mode) {32 case 1: // equivalent sphere33 return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length);34 case 2: // outer rim average radius35 return 0.5*r_minor*(1.0 + x_core) + thick_rim;36 case 3: // outer rim min radius37 return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);38 case 4: // outer max radius39 return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);40 case 5: // half outer thickness41 return 0.5*length + thick_face;42 case 6: // half diagonal43 return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length);44 }45 10 } 46 11 -
sasmodels/models/core_shell_bicelle_elliptical.py
ree60aa7 r71b751d 147 147 "core_shell_bicelle_elliptical.c"] 148 148 have_Fq = True 149 effective_radius_type = [150 "equivalent sphere", "outer rim average radius", "outer rim min radius",151 "outer max radius", "half outer thickness", "half diagonal",152 ]153 149 154 150 def random(): -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c
rd299327 r71b751d 9 9 return M_PI*( (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 10 10 square(r_minor)*x_core*2.0*thick_face ); 11 }12 13 static double14 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 double21 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 double30 effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length)31 {32 switch (mode) {33 case 1: // equivalent sphere34 return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length);35 case 2: // outer rim average radius36 return 0.5*r_minor*(1.0 + x_core) + thick_rim;37 case 3: // outer rim min radius38 return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);39 case 4: // outer max radius40 return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim);41 case 5: // half outer thickness42 return 0.5*length + thick_face;43 case 6: // half diagonal (this ignores the missing "corners", so may give unexpected answer if thick_face44 // or thick_rim is extremely large)45 return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length);46 }47 11 } 48 12 -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py
ree60aa7 r71b751d 160 160 "core_shell_bicelle_elliptical_belt_rough.c"] 161 161 have_Fq = True 162 effective_radius_type = [163 "equivalent sphere", "outer rim average radius", "outer rim min radius",164 "outer max radius", "half outer thickness", "half diagonal",165 ]166 162 167 163 demo = dict(scale=1, background=0, -
sasmodels/models/core_shell_cylinder.c
ree60aa7 r71b751d 11 11 { 12 12 return M_PI*square(radius+thickness)*(length+2.0*thickness); 13 }14 15 static double16 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 double23 radius_from_diagonal(double radius, double thickness, double length)24 {25 const double radius_outer = radius + thickness;26 const double length_outer = length + 2.0*thickness;27 return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer);28 }29 30 static double31 effective_radius(int mode, double radius, double thickness, double length)32 {33 switch (mode) {34 case 1: // equivalent sphere35 return radius_from_volume(radius, thickness, length);36 case 2: // outer radius37 return radius + thickness;38 case 3: // half outer length39 return 0.5*length + thickness;40 case 4: // half min outer length41 return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness);42 case 5: // half max outer length43 return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness);44 case 6: // half outer diagonal45 return radius_from_diagonal(radius,thickness,length);46 }47 13 } 48 14 -
sasmodels/models/core_shell_cylinder.py
ree60aa7 re31b19a 132 132 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 133 133 have_Fq = True 134 effective_radius_type = [ 135 "equivalent sphere", "outer radius", "half outer length", 136 "half min outer dimension", "half max outer dimension", 137 "half outer diagonal", 138 ] 134 135 def ER(radius, thickness, length): 136 """ 137 Returns the effective radius used in the S*P calculation 138 """ 139 radius = radius + thickness 140 length = length + 2 * thickness 141 ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 142 return 0.5 * (ddd) ** (1. / 3.) 139 143 140 144 def random(): -
sasmodels/models/core_shell_ellipsoid.c
ree60aa7 r71b751d 36 36 double vol = M_4PI_3*equat_shell*equat_shell*polar_shell; 37 37 return vol; 38 }39 40 static double41 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 double48 radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell)49 {50 // Trivial cases51 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-144955 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_tot58 ? radius_polar_tot / radius_equat_tot59 : 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 double70 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 switch (mode) {76 case 1: // equivalent sphere77 return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell);78 case 2: // average outer curvature79 return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell);80 case 3: // min outer radius81 return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot);82 case 4: // max outer radius83 return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot);84 }85 38 } 86 39 -
sasmodels/models/core_shell_ellipsoid.py
ree60aa7 r71b751d 146 146 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 147 147 have_Fq = True 148 effective_radius_type = [ 149 "equivalent sphere", "average outer curvature", 150 "min outer radius", "max outer radius", 151 ] 148 149 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 150 """ 151 Returns the effective radius used in the S*P calculation 152 """ 153 from .ellipsoid import ER as ellipsoid_ER 154 polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell 155 equat_outer = radius_equat_core + thick_shell 156 return ellipsoid_ER(polar_outer, equat_outer) 152 157 153 158 def random(): -
sasmodels/models/core_shell_parallelepiped.c
ree60aa7 r71b751d 25 25 2.0 * length_a * length_b * thick_rim_c; 26 26 #endif 27 }28 29 static double30 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 = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);34 return cbrt(volume/M_4PI_3);35 }36 37 static double38 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 double45 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 switch (mode) {49 case 1: // equivalent sphere50 return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);51 case 2: // half outer length a52 return 0.5 * length_a + thick_rim_a;53 case 3: // half outer length b54 return 0.5 * length_b + thick_rim_b;55 case 4: // half outer length c56 return 0.5 * length_c + thick_rim_c;57 case 5: // equivalent circular cross-section58 return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b);59 case 6: // half outer ab diagonal60 return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b));61 case 7: // half outer diagonal62 return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*thick_rim_c));63 }64 27 } 65 28 -
sasmodels/models/core_shell_parallelepiped.py
ree60aa7 r71b751d 98 98 is the scattering length of the solvent. 99 99 100 .. note:: 100 .. note:: 101 101 102 102 the code actually implements two substitutions: $d(cos\alpha)$ is … … 227 227 source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] 228 228 have_Fq = True 229 effective_radius_type = [ 230 "equivalent sphere", 231 "half outer length_a", "half outer length_b", "half outer length_c", 232 "equivalent circular cross-section", 233 "half outer ab diagonal", "half outer diagonal", 234 ] 229 230 231 def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c): 232 """ 233 Return equivalent radius (ER) 234 """ 235 from .parallelepiped import ER as ER_p 236 237 a = length_a + 2*thick_rim_a 238 b = length_b + 2*thick_rim_b 239 c = length_c + 2*thick_rim_c 240 return ER_p(a, b, c) 241 242 # VR defaults to 1.0 235 243 236 244 def random(): -
sasmodels/models/core_shell_sphere.c
ree60aa7 r71b751d 3 3 { 4 4 return M_4PI_3 * cube(radius + thickness); 5 }6 7 static double8 effective_radius(int mode, double radius, double thickness)9 {10 switch (mode) {11 case 1: // outer radius12 return radius + thickness;13 case 2: // core radius14 return radius;15 }16 5 } 17 6 -
sasmodels/models/core_shell_sphere.py
ree60aa7 rda1c8d1 78 78 source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 79 79 have_Fq = True 80 effective_radius_type = ["outer radius", "core radius"]81 80 82 81 demo = dict(scale=1, background=0, radius=60, thickness=10, 83 82 sld_core=1.0, sld_shell=2.0, sld_solvent=0.0) 83 84 def ER(radius, thickness): 85 """ 86 Equivalent radius 87 @param radius: core radius 88 @param thickness: shell thickness 89 """ 90 return radius + thickness 84 91 85 92 def random(): … … 96 103 97 104 tests = [ 98 #[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],105 [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 99 106 100 107 # The SasView test result was 0.00169, with a background of 0.001 -
sasmodels/models/cylinder.c
ree60aa7 r71b751d 11 11 { 12 12 return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 13 }14 15 static double16 radius_from_volume(double radius, double length)17 {18 return cbrt(0.75*radius*radius*length);19 }20 21 static double22 radius_from_diagonal(double radius, double length)23 {24 return sqrt(radius*radius + 0.25*length*length);25 }26 27 static double28 effective_radius(int mode, double radius, double length)29 {30 switch (mode) {31 case 1:32 return radius_from_volume(radius, length);33 case 2:34 return radius;35 case 3:36 return 0.5*length;37 case 4:38 return (radius < 0.5*length ? radius : 0.5*length);39 case 5:40 return (radius > 0.5*length ? radius : 0.5*length);41 case 6:42 return radius_from_diagonal(radius,length);43 }44 13 } 45 14 -
sasmodels/models/cylinder.py
ree60aa7 r71b751d 139 139 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 140 140 have_Fq = True 141 effective_radius_type = [ 142 "equivalent sphere", "radius", 143 "half length", "half min dimension", "half max dimension", "half diagonal", 144 ] 141 142 def ER(radius, length): 143 """ 144 Return equivalent radius (ER) 145 """ 146 ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 147 return 0.5 * (ddd) ** (1. / 3.) 145 148 146 149 def random(): -
sasmodels/models/ellipsoid.c
ree60aa7 r71b751d 4 4 return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 5 5 } 6 7 static double8 radius_from_volume(double radius_polar, double radius_equatorial)9 {10 return cbrt(radius_polar*radius_equatorial*radius_equatorial);11 }12 13 static double14 radius_from_curvature(double radius_polar, double radius_equatorial)15 {16 // Trivial cases17 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-144921 const double ratio = (radius_polar < radius_equatorial22 ? radius_polar / radius_equatorial23 : 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 double34 effective_radius(int mode, double radius_polar, double radius_equatorial)35 {36 switch (mode) {37 case 1: // equivalent sphere38 return radius_from_volume(radius_polar, radius_equatorial);39 case 2: // average curvature40 return radius_from_curvature(radius_polar, radius_equatorial);41 case 3: // min radius42 return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial);43 case 4: // max radius44 return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial);45 }46 }47 48 6 49 7 static void -
sasmodels/models/ellipsoid.py
ree60aa7 r0168844 169 169 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 170 170 have_Fq = True 171 effective_radius_type = [ 172 "equivalent sphere", "average curvature", "min radius", "max radius", 173 ] 171 172 def ER(radius_polar, radius_equatorial): 173 # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 174 ee = np.empty_like(radius_polar) 175 idx = radius_polar > radius_equatorial 176 ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2 177 idx = radius_polar < radius_equatorial 178 ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2 179 valid = (radius_polar * radius_equatorial != 0) & (radius_polar != radius_equatorial) 180 bd = 1.0 - ee[valid] 181 e1 = np.sqrt(ee[valid]) 182 b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd)) 183 bL = (1.0 + e1) / (1.0 - e1) 184 b2 = 1.0 + bd / 2 / e1 * np.log(bL) 185 delta = 0.75 * b1 * b2 186 ddd = 2.0 * (delta + 1.0) * (radius_polar * radius_equatorial**2)[valid] 187 188 r = np.zeros_like(radius_polar) 189 r[valid] = 0.5 * cbrt(ddd) 190 idx = radius_polar == radius_equatorial 191 r[idx] = radius_polar[idx] 192 return r 174 193 175 194 def random(): -
sasmodels/models/elliptical_cylinder.c
ree60aa7 r71b751d 3 3 { 4 4 return M_PI * radius_minor * radius_minor * r_ratio * length; 5 }6 7 static double8 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 double15 radius_from_min_dimension(double radius_minor, double r_ratio, double hlength)16 {17 const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);18 return (rad_min < hlength ? rad_min : hlength);19 }20 21 static double22 radius_from_max_dimension(double radius_minor, double r_ratio, double hlength)23 {24 const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);25 return (rad_max > hlength ? rad_max : hlength);26 }27 28 static double29 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 double36 effective_radius(int mode, double radius_minor, double r_ratio, double length)37 {38 switch (mode) {39 case 1: // equivalent sphere40 return radius_from_volume(radius_minor, r_ratio, length);41 case 2: // average radius42 return 0.5*radius_minor*(1.0 + r_ratio);43 case 3: // min radius44 return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);45 case 4: // max radius46 return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);47 case 5: // equivalent circular cross-section48 return sqrt(radius_minor*radius_minor*r_ratio);49 case 6: // half length50 return 0.5*length;51 case 7: // half min dimension52 return radius_from_min_dimension(radius_minor,r_ratio,0.5*length);53 case 8: // half max dimension54 return radius_from_max_dimension(radius_minor,r_ratio,0.5*length);55 case 9: // half diagonal56 return radius_from_diagonal(radius_minor,r_ratio,length);57 }58 5 } 59 6 -
sasmodels/models/elliptical_cylinder.py
ree60aa7 r71b751d 123 123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 124 124 have_Fq = True 125 effective_radius_type = [126 "equivalent sphere", "average radius", "min radius", "max radius",127 "equivalent circular cross-section",128 "half length", "half min dimension", "half max dimension", "half diagonal",129 ]130 125 131 126 demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 132 127 sld=4.0, sld_solvent=1.0, theta=10.0, phi=20, psi=30, 133 128 theta_pd=10, phi_pd=2, psi_pd=3) 129 130 def ER(radius_minor, axis_ratio, length): 131 """ 132 Equivalent radius 133 @param radius_minor: Ellipse minor radius 134 @param axis_ratio: Ratio of major radius over minor radius 135 @param length: Length of the cylinder 136 """ 137 radius = sqrt(radius_minor * radius_minor * axis_ratio) 138 ddd = 0.75 * radius * (2 * radius * length 139 + (length + radius) * (length + pi * radius)) 140 return 0.5 * (ddd) ** (1. / 3.) 134 141 135 142 def random(): … … 155 162 156 163 tests = [ 157 #[{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024],158 #[{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1],164 [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], 165 [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], 159 166 160 167 # The SasView test result was 0.00169, with a background of 0.001 -
sasmodels/models/fractal_core_shell.py
r2cc8aa2 reb3eb38 134 134 return radius + thickness 135 135 136 #tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0],137 tests = [ 136 tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 137 138 138 # # At some point the SasView 3.x test result was deemed incorrect. The 139 139 #following tests were verified against NIST IGOR macros ver 7.850. -
sasmodels/models/fuzzy_sphere.py
ree60aa7 r71b751d 78 78 ["sld_solvent", "1e-6/Ang^2", 3, [-inf, inf], "sld", "Solvent scattering length density"], 79 79 ["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)"], 81 81 ] 82 82 # pylint: enable=bad-whitespace,line-too-long 83 83 84 source = ["lib/sas_3j1x_x.c" ,"fuzzy_sphere.c"]84 source = ["lib/sas_3j1x_x.c"] 85 85 have_Fq = True 86 effective_radius_type = ["radius", "radius + fuzziness"] 86 87 c_code = """ 88 static double form_volume(double radius) 89 { 90 return M_4PI_3*cube(radius); 91 } 92 93 static void Fq(double q, double *F1, double *F2, double sld, double sld_solvent, 94 double radius, double fuzziness) 95 { 96 const double qr = q*radius; 97 const double bes = sas_3j1x_x(qr); 98 const double qf = exp(-0.5*square(q*fuzziness)); 99 const double contrast = (sld - sld_solvent); 100 const double form = contrast * form_volume(radius) * bes * qf; 101 *F1 = 1.0e-2*form; 102 *F2 = 1.0e-4*form*form; 103 } 104 """ 105 106 def ER(radius): 107 """ 108 Return radius 109 """ 110 return radius 111 112 # VR defaults to 1.0 87 113 88 114 def random(): -
sasmodels/models/hollow_cylinder.c
ree60aa7 r71b751d 15 15 } 16 16 17 // TODO: interface to form_volume/shell_volume not yet settled18 static double19 shell_volume(double *total, double radius, double thickness, double length)20 {21 *total = M_PI*length*square(radius+thickness);22 return *total - M_PI*length*radius*radius;23 }24 25 17 static double 26 18 form_volume(double radius, double thickness, double length) 27 19 { 28 double total;29 return shell_volume(&total, radius, thickness, length);20 double v_shell = M_PI*length*(square(radius+thickness) - radius*radius); 21 return v_shell; 30 22 } 31 23 32 static double33 radius_from_volume(double radius, double thickness, double length)34 {35 const double volume_outer_cyl = M_PI*square(radius + thickness)*length;36 return cbrt(0.75*volume_outer_cyl/M_PI);37 }38 39 static double40 radius_from_diagonal(double radius, double thickness, double length)41 {42 return sqrt(square(radius + thickness) + 0.25*square(length));43 }44 45 static double46 effective_radius(int mode, double radius, double thickness, double length)47 {48 switch (mode) {49 case 1: // equivalent sphere50 return radius_from_volume(radius, thickness, length);51 case 2: // outer radius52 return radius + thickness;53 case 3: // half length54 return 0.5*length;55 case 4: // half outer min dimension56 return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length);57 case 5: // half outer max dimension58 return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length);59 case 6: // half outer diagonal60 return radius_from_diagonal(radius,thickness,length);61 }62 }63 24 64 25 static void -
sasmodels/models/hollow_cylinder.py
ree60aa7 r455aaa1 100 100 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 101 101 have_Fq = True 102 effective_radius_type = [ 103 "equivalent sphere", "outer radius", "half length", 104 "half outer min dimension", "half outer max dimension", 105 "half outer diagonal", 106 ] 102 103 # pylint: disable=W0613 104 def ER(radius, thickness, length): 105 """ 106 :param radius: Cylinder core radius 107 :param thickness: Cylinder wall thickness 108 :param length: Cylinder length 109 :return: Effective radius 110 """ 111 router = radius + thickness 112 if router == 0 or length == 0: 113 return 0.0 114 len1 = router 115 len2 = length/2.0 116 term1 = len1*len1*2.0*len2/2.0 117 term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0) 118 ddd = 3.0*term1*term2 119 diam = pow(ddd, (1.0/3.0)) 120 return diam 121 122 def VR(radius, thickness, length): 123 """ 124 :param radius: Cylinder radius 125 :param thickness: Cylinder wall thickness 126 :param length: Cylinder length 127 :return: Volf ratio for P(q)*S(q) 128 """ 129 router = radius + thickness 130 vol_core = pi*radius*radius*length 131 vol_total = pi*router*router*length 132 vol_shell = vol_total - vol_core 133 return vol_total, vol_shell 107 134 108 135 def random(): … … 135 162 tests = [ 136 163 [{}, 0.00005, 1764.926], 137 #[{}, 'VR', 0.55555556],164 [{}, 'VR', 0.55555556], 138 165 [{}, 0.001, 1756.76], 139 166 [{}, (qx, qy), 2.36885476192], -
sasmodels/models/hollow_rectangular_prism.c
ree60aa7 r71b751d 1 // TODO: interface to form_volume/shell_volume not yet settled2 1 static double 3 shell_volume(double *total,double length_a, double b2a_ratio, double c2a_ratio, double thickness)2 form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 4 3 { 5 4 double length_b = length_a * b2a_ratio; … … 9 8 double c_core = length_c - 2.0*thickness; 10 9 double vol_core = a_core * b_core * c_core; 11 *total = length_a * length_b * length_c; 12 return *total - vol_core; 13 } 14 15 static double 16 form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 17 { 18 double total; 19 return shell_volume(&total, length_a, b2a_ratio, c2a_ratio, thickness); 20 } 21 22 static double 23 effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 24 //effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 25 // "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 26 // NOTE length_a is external dimension! 27 { 28 switch (mode) { 29 case 1: // equivalent sphere 30 return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI); 31 case 2: // half length_a 32 return 0.5 * length_a; 33 case 3: // half length_b 34 return 0.5 * length_a*b2a_ratio; 35 case 4: // half length_c 36 return 0.5 * length_a*c2a_ratio; 37 case 5: // equivalent outer circular cross-section 38 return length_a*sqrt(b2a_ratio/M_PI); 39 case 6: // half ab diagonal 40 return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 41 case 7: // half diagonal 42 return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 43 } 10 double vol_total = length_a * length_b * length_c; 11 double vol_shell = vol_total - vol_core; 12 return vol_shell; 44 13 } 45 14 -
sasmodels/models/hollow_rectangular_prism.py
ree60aa7 r455aaa1 132 132 "Solvent scattering length density"], 133 133 ["length_a", "Ang", 35, [0, inf], "volume", 134 "Shorte st, external, size of the parallelepiped"],134 "Shorter side of the parallelepiped"], 135 135 ["b2a_ratio", "Ang", 1, [0, inf], "volume", 136 136 "Ratio sides b/a"], … … 149 149 source = ["lib/gauss76.c", "hollow_rectangular_prism.c"] 150 150 have_Fq = True 151 effective_radius_type = [ 152 "equivalent sphere", "half length_a", "half length_b", "half length_c", 153 "equivalent outer circular cross-section", 154 "half ab diagonal", "half diagonal", 155 ] 151 152 def ER(length_a, b2a_ratio, c2a_ratio, thickness): 153 """ 154 Return equivalent radius (ER) 155 thickness parameter not used 156 """ 157 b_side = length_a * b2a_ratio 158 c_side = length_a * c2a_ratio 159 160 # surface average radius (rough approximation) 161 surf_rad = sqrt(length_a * b_side / pi) 162 163 ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 164 return 0.5 * (ddd) ** (1. / 3.) 165 166 def VR(length_a, b2a_ratio, c2a_ratio, thickness): 167 """ 168 Return shell volume and total volume 169 """ 170 b_side = length_a * b2a_ratio 171 c_side = length_a * c2a_ratio 172 a_core = length_a - 2.0*thickness 173 b_core = b_side - 2.0*thickness 174 c_core = c_side - 2.0*thickness 175 vol_core = a_core * b_core * c_core 176 vol_total = length_a * b_side * c_side 177 vol_shell = vol_total - vol_core 178 return vol_total, vol_shell 179 156 180 157 181 def random(): -
sasmodels/models/hollow_rectangular_prism_thin_walls.c
ree60aa7 r71b751d 1 // TODO: interface to form_volume/shell_volume not yet settled2 1 static double 3 shell_volume(double *total,double length_a, double b2a_ratio, double c2a_ratio)2 form_volume(double length_a, double b2a_ratio, double c2a_ratio) 4 3 { 5 4 double length_b = length_a * b2a_ratio; 6 5 double length_c = length_a * c2a_ratio; 7 6 double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 8 *total = length_a * length_b * length_c;9 7 return vol_shell; 10 }11 12 static double13 form_volume(double length_a, double b2a_ratio, double c2a_ratio)14 {15 double total;16 return shell_volume(&total, length_a, b2a_ratio, c2a_ratio);17 }18 19 20 static double21 effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio)22 {23 switch (mode) {24 case 1: // equivalent sphere25 return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI);26 case 2: // half length_a27 return 0.5 * length_a;28 case 3: // half length_b29 return 0.5 * length_a*b2a_ratio;30 case 4: // half length_c31 return 0.5 * length_a*c2a_ratio;32 case 5: // equivalent outer circular cross-section33 return length_a*sqrt(b2a_ratio/M_PI);34 case 6: // half ab diagonal35 return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio)));36 case 7: // half diagonal37 return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio)));38 }39 8 } 40 9 -
sasmodels/models/hollow_rectangular_prism_thin_walls.py
ree60aa7 r6e7d7b6 109 109 source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"] 110 110 have_Fq = True 111 effective_radius_type = [ 112 "equivalent sphere", "half length_a", "half length_b", "half length_c", 113 "equivalent outer circular cross-section", 114 "half ab diagonal", "half diagonal", 115 ] 111 112 def ER(length_a, b2a_ratio, c2a_ratio): 113 """ 114 Return equivalent radius (ER) 115 """ 116 b_side = length_a * b2a_ratio 117 c_side = length_a * c2a_ratio 118 119 # surface average radius (rough approximation) 120 surf_rad = sqrt(length_a * b_side / pi) 121 122 ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 123 return 0.5 * (ddd) ** (1. / 3.) 124 125 def VR(length_a, b2a_ratio, c2a_ratio): 126 """ 127 Return shell volume and total volume 128 """ 129 b_side = length_a * b2a_ratio 130 c_side = length_a * c2a_ratio 131 vol_total = length_a * b_side * c_side 132 vol_shell = 2.0 * (length_a*b_side + length_a*c_side + b_side*c_side) 133 return vol_shell, vol_total 116 134 117 135 -
sasmodels/models/mono_gauss_coil.py
ree60aa7 r2d81cfe 54 54 55 55 import numpy as np 56 from numpy import inf 56 from numpy import inf, exp, errstate 57 57 58 58 name = "mono_gauss_coil" … … 69 69 parameters = [ 70 70 ["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"], 72 72 ] 73 73 # pylint: enable=bad-whitespace, line-too-long 74 74 75 source = ["mono_gauss_coil.c"] 76 have_Fq = False 77 effective_radius_type = ["R_g", "2R_g", "3R_g", "(5/3)^0.5*R_g"] 75 # NB: Scale and Background are implicit parameters on every model 76 def Iq(q, i_zero, rg): 77 # pylint: disable = missing-docstring 78 z = (q * rg)**2 78 79 80 with errstate(invalid='ignore'): 81 inten = (i_zero * 2.0) * (exp(-z) + z - 1.0)/z**2 82 inten[q == 0] = i_zero 83 return inten 84 Iq.vectorized = True # Iq accepts an array of q values 79 85 80 86 def random(): -
sasmodels/models/multilayer_vesicle.c
ree60aa7 r71b751d 47 47 } 48 48 49 static double50 effective_radius(int mode, double radius, double thick_shell, double thick_solvent, double fp_n_shells)51 {52 // case 1: outer radius53 return radius + fp_n_shells*thick_shell + (fp_n_shells - 1.0)*thick_solvent;54 }55 56 49 static void 57 50 Fq(double q, -
sasmodels/models/multilayer_vesicle.py
ree60aa7 r71b751d 146 146 source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 147 147 have_Fq = True 148 effective_radius_type = ["outer radius"] 148 149 def ER(radius, thick_shell, thick_solvent, n_shells): 150 n_shells = int(n_shells+0.5) 151 return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent 149 152 150 153 def random(): -
sasmodels/models/onion.c
ree60aa7 r71b751d 30 30 31 31 static double 32 outer_radius(double radius_core, double n_shells, double thickness[])32 form_volume(double radius_core, double n_shells, double thickness[]) 33 33 { 34 34 int n = (int)(n_shells+0.5); … … 37 37 r += thickness[i]; 38 38 } 39 return r; 40 } 41 42 static double 43 form_volume(double radius_core, double n_shells, double thickness[]) 44 { 45 return M_4PI_3*cube(outer_radius(radius_core, n_shells, thickness)); 46 } 47 48 static double 49 effective_radius(int mode, double radius_core, double n_shells, double thickness[]) 50 { 51 // case 1: outer radius 52 return outer_radius(radius_core, n_shells, thickness); 39 return M_4PI_3*cube(r); 53 40 } 54 41 -
sasmodels/models/onion.py
ree60aa7 r71b751d 316 316 single = False 317 317 have_Fq = True 318 effective_radius_type = ["outer radius"]319 318 320 319 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] … … 366 365 367 366 return np.asarray(z), np.asarray(rho) 367 368 def ER(radius_core, n_shells, thickness): 369 """Effective radius""" 370 n = int(n_shells[0]+0.5) 371 return np.sum(thickness[:n], axis=0) + radius_core 368 372 369 373 demo = { -
sasmodels/models/parallelepiped.c
ree60aa7 r71b751d 5 5 } 6 6 7 static double8 effective_radius(int mode, double length_a, double length_b, double length_c)9 {10 switch (mode) {11 case 1: // equivalent sphere12 return cbrt(0.75*length_a*length_b*length_c/M_PI);13 case 2: // half length_a14 return 0.5 * length_a;15 case 3: // half length_b16 return 0.5 * length_b;17 case 4: // half length_c18 return 0.5 * length_c;19 case 5: // equivalent circular cross-section20 return sqrt(length_a*length_b/M_PI);21 case 6: // half ab diagonal22 return 0.5*sqrt(length_a*length_a + length_b*length_b);23 case 7: // half diagonal24 return 0.5*sqrt(length_a*length_a + length_b*length_b + length_c*length_c);25 }26 }27 7 28 8 static void -
sasmodels/models/parallelepiped.py
ree60aa7 r71b751d 140 140 ensuring that the inequality $A < B < C$ is not violated, The calculation 141 141 will not report an error, but the results may be not correct. 142 142 143 143 .. _parallelepiped-orientation: 144 144 … … 231 231 source = ["lib/gauss76.c", "parallelepiped.c"] 232 232 have_Fq = True 233 effective_radius_type = [ 234 "equivalent sphere", "half length_a", "half length_b", "half length_c", 235 "equivalent circular cross-section", "half ab diagonal", "half diagonal", 236 ] 233 234 def ER(length_a, length_b, length_c): 235 """ 236 Return effective radius (ER) for P(q)*S(q) 237 """ 238 # now that axes can be in any size order, need to sort a,b,c 239 # where a~b and c is either much smaller or much larger 240 abc = np.vstack((length_a, length_b, length_c)) 241 abc = np.sort(abc, axis=0) 242 selector = (abc[1] - abc[0]) > (abc[2] - abc[1]) 243 length = np.where(selector, abc[0], abc[2]) 244 # surface average radius (rough approximation) 245 radius = sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi) 246 247 ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius)) 248 return 0.5 * (ddd) ** (1. / 3.) 249 250 # VR defaults to 1.0 251 237 252 238 253 def random(): -
sasmodels/models/pearl_necklace.py
r2cc8aa2 r2d81cfe 140 140 thick_string_pd=0.2, thick_string_pd_n=5, 141 141 ) 142 # ER function is not being used here, not that it is likely very sensible to 143 # include an S(Q) with this model, the default in sasview 5.0 would be to the 144 # "unconstrained" radius_effective. 145 #tests = [[{}, 0.001, 17380.245], [{}, 'ER', 115.39502]] 146 tests = [[{}, 0.001, 17380.245]] 142 143 tests = [[{}, 0.001, 17380.245], [{}, 'ER', 115.39502]] -
sasmodels/models/pringle.c
ree60aa7 r74768cb 104 104 } 105 105 106 static double107 effective_radius(int mode, double radius, double thickness, double alpha, double beta)108 {109 switch (mode) {110 case 1: // equivalent sphere111 return cbrt(0.75*radius*radius*thickness);112 case 2: // radius113 return radius;114 }115 }116 117 106 double Iq( 118 107 double q, -
sasmodels/models/pringle.py
ree60aa7 ref07e95 72 72 73 73 74 source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", 74 source = ["lib/polevl.c", "lib/sas_J0.c", "lib/sas_J1.c", \ 75 75 "lib/sas_JN.c", "lib/gauss76.c", "pringle.c"] 76 effective_radius_type = ["equivalent sphere", "radius"] 76 77 def ER(radius, thickness, alpha, beta): 78 """ 79 Return equivalent radius (ER) 80 """ 81 ddd = 0.75 * radius * (2 * radius * thickness + (thickness + radius) \ 82 * (thickness + pi * radius)) 83 return 0.5 * (ddd) ** (1. / 3.) 77 84 78 85 def random(): -
sasmodels/models/raspberry.c
ree60aa7 r71b751d 1 double form_volume(double radius_lg , double radius_sm, double penetration);1 double form_volume(double radius_lg); 2 2 3 3 double Iq(double q, … … 6 6 double radius_lg, double radius_sm, double penetration); 7 7 8 double form_volume(double radius_lg , double radius_sm, double penetration)8 double form_volume(double radius_lg) 9 9 { 10 10 //Because of the complex structure, volume normalization must … … 12 12 double volume=1.0; 13 13 return volume; 14 }15 16 static double17 effective_radius(int mode, double radius_lg, double radius_sm, double penetration)18 {19 switch (mode) {20 case 1: // radius_large21 return radius_lg;22 case 2: // radius_outer23 return radius_lg + 2.0*radius_sm - penetration;24 }25 14 } 26 15 -
sasmodels/models/raspberry.py
ree60aa7 ref07e95 145 145 ["radius_lg", "Ang", 5000, [0, inf], "volume", 146 146 "radius of large spheres"], 147 ["radius_sm", "Ang", 100, [0, inf], " volume",147 ["radius_sm", "Ang", 100, [0, inf], "", 148 148 "radius of small spheres"], 149 ["penetration", "Ang", 0, [-1, 1], " volume",149 ["penetration", "Ang", 0, [-1, 1], "", 150 150 "fractional penetration depth of small spheres into large sphere"], 151 151 ] 152 152 153 153 source = ["lib/sas_3j1x_x.c", "raspberry.c"] 154 effective_radius_type = ["radius_large", "radius_outer"]155 154 156 155 def random(): -
sasmodels/models/rectangular_prism.c
ree60aa7 r71b751d 3 3 { 4 4 return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio); 5 }6 7 static double8 effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio)9 {10 switch (mode) {11 case 1: // equivalent sphere12 return cbrt(0.75*cube(length_a)*b2a_ratio*c2a_ratio/M_PI);13 case 2: // half length_a14 return 0.5 * length_a;15 case 3: // half length_b16 return 0.5 * length_a*b2a_ratio;17 case 4: // half length_c18 return 0.5 * length_a*c2a_ratio;19 case 5: // equivalent circular cross-section20 return length_a*sqrt(b2a_ratio/M_PI);21 case 6: // half ab diagonal22 return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio)));23 case 7: // half diagonal24 return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio)));25 }26 5 } 27 6 -
sasmodels/models/rectangular_prism.py
ree60aa7 r71b751d 136 136 source = ["lib/gauss76.c", "rectangular_prism.c"] 137 137 have_Fq = True 138 effective_radius_type = [ 139 "equivalent sphere", "half length_a", "half length_b", "half length_c", 140 "equivalent circular cross-section", "half ab diagonal", "half diagonal", 141 ] 138 139 def ER(length_a, b2a_ratio, c2a_ratio): 140 """ 141 Return equivalent radius (ER) 142 """ 143 b_side = length_a * b2a_ratio 144 c_side = length_a * c2a_ratio 145 146 # surface average radius (rough approximation) 147 surf_rad = sqrt(length_a * b_side / pi) 148 149 ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 150 return 0.5 * (ddd) ** (1. / 3.) 142 151 143 152 def random(): -
sasmodels/models/sphere.py
ree60aa7 r71b751d 67 67 ] 68 68 69 source = ["lib/sas_3j1x_x.c" ,"sphere.c"]69 source = ["lib/sas_3j1x_x.c"] 70 70 have_Fq = True 71 effective_radius_type = ["radius"] 71 72 c_code = """ 73 static double form_volume(double radius) 74 { 75 return M_4PI_3*cube(radius); 76 } 77 78 static void Fq(double q, double *f1, double *f2, double sld, double sld_solvent, double radius) 79 { 80 const double bes = sas_3j1x_x(q*radius); 81 const double contrast = (sld - sld_solvent); 82 const double form = contrast * form_volume(radius) * bes; 83 *f1 = 1.0e-2*form; 84 *f2 = 1.0e-4*form*form; 85 } 86 """ 87 88 def ER(radius): 89 """ 90 Return equivalent radius (ER) 91 """ 92 return radius 93 94 # VR defaults to 1.0 72 95 73 96 def random(): … … 83 106 "radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, 84 107 0.2, 0.228843], 85 #[{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.],86 #[{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "VR", 1.],108 [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.], 109 [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "VR", 1.], 87 110 ] -
sasmodels/models/spherical_sld.c
ree60aa7 r71b751d 1 static double 2 outer_radius(double fp_n_shells, double thickness[], double interface[]) 1 static double form_volume( 2 double fp_n_shells, 3 double thickness[], 4 double interface[]) 3 5 { 4 6 int n_shells= (int)(fp_n_shells + 0.5); … … 7 9 r += thickness[i] + interface[i]; 8 10 } 9 return r; 10 } 11 12 static double form_volume( 13 double fp_n_shells, 14 double thickness[], 15 double interface[]) 16 { 17 return M_4PI_3*cube(outer_radius(fp_n_shells, thickness, interface)); 18 } 19 20 static double 21 effective_radius(int mode, double fp_n_shells, double thickness[], double interface[]) 22 { 23 // case 1: outer radius 24 return outer_radius(fp_n_shells, thickness, interface); 11 return M_4PI_3*cube(r); 25 12 } 26 13 -
sasmodels/models/spherical_sld.py
ree60aa7 r5601947 22 22 23 23 0: erf($\nu z$) 24 24 25 25 1: Rpow($z^\nu$) 26 26 27 27 2: Lpow($z^\nu$) 28 28 29 29 3: Rexp($-\nu z$) 30 30 31 31 4: Lexp($-\nu z$) 32 32 … … 217 217 ["interface[n_shells]", "Ang", 50.0, [0, inf], "volume", "thickness of the interface"], 218 218 ["shape[n_shells]", "", 0, [SHAPES], "", "interface shape"], 219 ["nu[n_shells]", "", 2.5, [ 1, inf], "", "interface shape exponent"],219 ["nu[n_shells]", "", 2.5, [0, inf], "", "interface shape exponent"], 220 220 ["n_steps", "", 35, [0, inf], "", "number of steps in each interface (must be an odd integer)"], 221 221 ] … … 224 224 single = False # TODO: fix low q behaviour 225 225 have_Fq = True 226 effective_radius_type = ["outer radius"]227 226 228 227 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] … … 270 269 271 270 271 def ER(n_shells, thickness, interface): 272 """Effective radius""" 273 n_shells = int(n_shells + 0.5) 274 total = (np.sum(thickness[:n_shells], axis=1) 275 + np.sum(interface[:n_shells], axis=1)) 276 return total 277 278 272 279 demo = { 273 280 "n_shells": 5, -
sasmodels/models/triaxial_ellipsoid.c
ree60aa7 r71b751d 7 7 } 8 8 9 static double10 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 double16 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 double23 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 double30 effective_radius(int mode, double radius_equat_minor, double radius_equat_major, double radius_polar)31 {32 switch (mode) {33 case 1: // equivalent sphere34 return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar);35 case 2: // min radius36 return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar);37 case 3: // max radius38 return radius_from_max_dimension(radius_equat_minor,radius_equat_major, radius_polar);39 }40 }41 9 42 10 static void -
sasmodels/models/triaxial_ellipsoid.py
ree60aa7 r71b751d 158 158 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 159 159 have_Fq = True 160 effective_radius_type = ["equivalent sphere", "min radius", "max radius"] 160 161 def ER(radius_equat_minor, radius_equat_major, radius_polar): 162 """ 163 Returns the effective radius used in the S*P calculation 164 """ 165 from .ellipsoid import ER as ellipsoid_ER 166 167 # now that radii can be in any size order, radii need sorting a,b,c 168 # where a~b and c is either much smaller or much larger 169 radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar)) 170 radii = np.sort(radii, axis=0) 171 selector = (radii[1] - radii[0]) > (radii[2] - radii[1]) 172 polar = np.where(selector, radii[0], radii[2]) 173 equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2])) 174 return ellipsoid_ER(polar, equatorial) 161 175 162 176 def random(): -
sasmodels/models/vesicle.c
ree60aa7 r71b751d 1 // TODO: interface to form_volume/shell_volume not yet settled2 static double3 shell_volume(double *total, double radius, double thickness)4 {5 //note that for the vesicle model, the volume is ONLY the shell volume6 *total = M_4PI_3 * cube(radius+thickness);7 return *total - M_4PI_3*cube(radius);8 }9 10 1 static double 11 2 form_volume(double radius, double thickness) 12 3 { 13 4 //note that for the vesicle model, the volume is ONLY the shell volume 14 double total; 15 return shell_volume(&total, radius, thickness); 5 return M_4PI_3*(cube(radius+thickness) - cube(radius)); 16 6 } 17 7 18 static double19 effective_radius(int mode, double radius, double thickness)20 {21 // case 1: outer radius22 return radius + thickness;23 }24 8 25 9 static void -
sasmodels/models/vesicle.py
ree60aa7 rb477605 100 100 source = ["lib/sas_3j1x_x.c", "vesicle.c"] 101 101 have_Fq = True 102 effective_radius_type = ["outer radius"] 102 103 def ER(radius, thickness): 104 ''' 105 returns the effective radius used in the S*P calculation 106 107 :param radius: core radius 108 :param thickness: shell thickness 109 ''' 110 return radius + thickness 111 112 def VR(radius, thickness): 113 ''' 114 returns the volumes of the shell and of the whole sphere including the 115 core plus shell - is used to normalize when including polydispersity. 116 117 :param radius: core radius 118 :param thickness: shell thickness 119 :return whole: volume of core and shell 120 :return whole-core: volume of the shell 121 ''' 122 123 whole = 4./3. * pi * (radius + thickness)**3 124 core = 4./3. * pi * radius**3 125 return whole, whole - core 103 126 104 127 def random(): … … 128 151 [{}, 0.100600200401, 1.77063682331], 129 152 [{}, 0.5, 0.00355351388906], 130 #[{}, 'ER', 130.],131 #[{}, 'VR', 0.54483386436],153 [{}, 'ER', 130.], 154 [{}, 'VR', 0.54483386436], 132 155 ] -
sasmodels/product.py
r33d7be3 rc0131d44 24 24 # pylint: disable=unused-import 25 25 try: 26 from typing import Tuple, Callable , Union26 from typing import Tuple, Callable 27 27 except ImportError: 28 28 pass … … 37 37 #] 38 38 39 RADIUS_ID = "radius_effective" 40 VOLFRAC_ID = "volfraction" 39 ER_ID = "radius_effective" 40 VF_ID = "volfraction" 41 41 42 def make_extra_pars(p_info): 42 43 pars = [] … … 50 51 "Structure factor calculation") 51 52 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)61 53 return pars 62 54 … … 73 65 # have any magnetic parameters 74 66 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 == V OLFRAC_ID:79 raise TypeError("S needs {} as second parameter".format(V OLFRAC_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)) 80 72 if not s_info.parameters.magnetism_index == []: 81 73 raise TypeError("S should not have SLD parameters") … … 83 75 s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 84 76 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 86 79 # are any names in P that overlap with those in S, modify the name in S 87 80 # to distinguish it. 88 81 p_set = set(p.id for p in p_pars.kernel_parameters) 89 82 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:]] 91 84 # Check if still a collision after renaming. This could happen if for 92 85 # example S has volfrac and P has both volfrac and volfrac_S. … … 94 87 raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 95 88 96 # make sure effective radius is not a polydisperse parameter in product97 s_list[0] = copy(s_list[0])98 s_list[0].polydisperse = False99 100 89 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)) 102 91 combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info) 103 92 parameters = ParameterTable(combined_pars) … … 105 94 def random(): 106 95 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:]) 108 97 combined_pars.update((translate_name[k], v) 109 98 for k, v in s_info.random().items() … … 168 157 return par 169 158 170 def _intermediates(P, S , effective_radius):171 # type: (np.ndarray, np.ndarray , float) -> OrderedDict[str, np.ndarray]159 def _intermediates(P, S): 160 # type: (np.ndarray, np.ndarray) -> OrderedDict[str, np.ndarray] 172 161 """ 173 162 Returns intermediate results for standard product (P(Q)*S(Q)) … … 176 165 ("P(Q)", P), 177 166 ("S(Q)", S), 178 ("effective_radius", effective_radius),179 167 )) 180 168 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. 169 def _intermediates_beta(F1, F2, S, scale, bg): 170 # type: (np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray) -> OrderedDict[str, np.ndarray] 171 """ 172 Returns intermediate results for beta approximation-enabled product 191 173 """ 192 174 # TODO: 1. include calculated Q vector … … 197 179 ("beta(Q)", F1**2 / F2), 198 180 ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), 199 ("effective_radius", effective_radius),200 181 # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), 201 182 )) … … 281 262 p_values = np.hstack(p_values).astype(self.p_kernel.dtype) 282 263 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 283 269 # 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 285 277 s_length = call_details.length[p_npars:p_npars+s_npars] 286 278 s_offset = call_details.offset[p_npars:p_npars+s_npars] … … 288 280 s_offset = np.hstack((nweights, s_offset)) 289 281 s_details = make_details(s_info, s_length, s_offset, nweights+1) 282 v, w = weights[:nweights], weights[nweights:] 290 283 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] 295 294 ] 296 295 spacer = (32 - sum(len(v) for v in s_values)%32)%32 … … 299 298 300 299 # 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] 312 302 313 303 # Call the kernels 304 s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 314 305 scale, background = values[0], values[1] 315 306 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) 322 308 combined_scale = scale*volfrac/volume_avg 323 309 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) 325 311 final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background 326 312 327 313 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) 336 317 final_result = scale*(p_result*s_result) + background 337 318 -
sasmodels/sasview_model.py
rc952e59 r84f2962 272 272 attrs['category'] = model_info.category 273 273 attrs['is_structure_factor'] = model_info.structure_factor 274 attrs['is_form_factor'] = model_info. effective_radius_typeis not None274 attrs['is_form_factor'] = model_info.ER is not None 275 275 attrs['is_multiplicity_model'] = variants[0] > 1 276 276 attrs['multiplicity_info'] = variants
Note: See TracChangeset
for help on using the changeset viewer.