Changeset 502c7b8 in sasmodels


Ignore:
Timestamp:
Aug 7, 2018 5:14:59 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
7b0abf8
Parents:
cdd676e (diff), 01c8d9e (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge remote-tracking branch 'greg/beta_approx' into beta_approx

Files:
6 added
14 edited

Legend:

Unmodified
Added
Removed

  • explore/beta/sasfit_compare.py

    rcdd676e r502c7b8  
    33import sys, os 
    44BETA_DIR = os.path.dirname(os.path.realpath(__file__)) 
    5 SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 
     5#SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 
     6SASMODELS_DIR = r"C:\Source\sasmodels" 
    67sys.path.insert(0, SASMODELS_DIR) 
    78 
    … …  
    6162    calculator = DirectModel(data, model,cutoff=0) 
    6263    calculator.pars = pars.copy() 
    63     calculator.pars.setdefault('background', 1e-3) 
     64    calculator.pars.setdefault('background', 0) 
    6465    return calculator 
    6566 
    … …  
    231232        #     = F2/total_volume 
    232233        IQD = F2/average_volume*1e-4*volfraction 
     234        F1 *= 1e-2  # Yun is using sld in 1/A^2, not 1e-6/A^2 
     235        F2 *= 1e-4 
    233236    elif norm == 'yun': 
    234237        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2 
    … …  
    336339    I = build_model(Pname+"@hardsphere", q) 
    337340    Pq = P(**Ppars)*pars.get('volfraction', 1) 
    338     #Sq = S(**Spars) 
     341    Sq = S(**Spars) 
    339342    Iq = I(**Ipars) 
    340343    #Iq = Pq*Sq*pars.get('volfraction', 1) 
    341     Sq = Iq/Pq 
    342     return Theory(Q=q, F1=None, F2=None, P=Pq, S=Sq, I=Iq, Seff=None, Ibeta=None) 
     344    #Sq = Iq/Pq 
     345    #Iq = None#= Sq = None 
     346    r=I._kernel.results 
     347    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r[1], Ibeta=Iq) 
    343348 
    344349def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): 
    … …  
    400405        radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type, 
    401406        volfraction=0.15, 
     407        radius_effective=270.7543927018, 
    402408        #radius_effective=12.59921049894873, 
    403409        ) 
    404     target = sasmodels_theory(q, model, **pars) 
     410    target = sasmodels_theory(q, model, beta_mode=1, **pars) 
    405411    actual = ellipsoid_pe(q, norm='sasview', **pars) 
    406412    title = " ".join(("sasmodels", model, pd_type)) 
    … …  
    450456    S = data[5] 
    451457    Seff = data[6] 
    452     target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, Seff=Seff) 
     458    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 
    453459    actual = sphere_r(Q, norm='yun', **pars) 
    454460    title = " ".join(("yun", "sphere", "10% dispersion 10% Vf")) 
    … …  
    462468    S = data[5] 
    463469    Seff = data[6] 
    464     target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, Seff=Seff) 
     470    target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 
    465471    actual = sphere_r(Q, norm='yun', **pars) 
    466472    title = " ".join(("yun", "sphere", "15% dispersion 10% Vf")) 

  • sasmodels/core.py

    r3221de0 r01c8d9e  
    140140    used with functions in generate to build the docs or extract model info. 
    141141    """ 
     142 
    142143    if "+" in model_string: 
    143144        parts = [load_model_info(part) 
    … …  
    205206 
    206207    numpy_dtype, fast, platform = parse_dtype(model_info, dtype, platform) 
    207  
    208208    source = generate.make_source(model_info) 
    209209    if platform == "dll": 
    … …  
    265265    # Assign default platform, overriding ocl with dll if OpenCL is unavailable 
    266266    # If opencl=False OpenCL is switched off 
    267  
    268267    if platform is None: 
    269268        platform = "ocl" 
    … …  
    291290    else: 
    292291        numpy_dtype = np.dtype(dtype) 
    293  
    294292    # Make sure that the type is supported by opencl, otherwise use dll 
    295293    if platform == "ocl": 

  • sasmodels/data.py

    r65fbf7c r01c8d9e  
    329329    else: 
    330330        dq = resolution * q 
    331  
    332331    data = Data1D(q, Iq, dx=dq, dy=dIq) 
    333332    data.filename = "fake data" 
    … …  
    523522                  else 'log') 
    524523        plt.xscale(xscale) 
     524         
    525525        plt.xlabel("$q$/A$^{-1}$") 
    526526        plt.yscale(yscale) 

  • sasmodels/generate.py

    rd86f0fc r01c8d9e  
    672672 
    673673# type in IQXY pattern could be single, float, double, long double, ... 
    674 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
     674_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]", 
    675675                           flags=re.MULTILINE) 
    676676def find_xy_mode(source): 
    … …  
    701701    return 'qa' 
    702702 
     703# type in IQXY pattern could be single, float, double, long double, ... 
     704_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE) 
     705def has_Fq(source): 
     706    for code in source: 
     707        m = _FQ_PATTERN.search(code) 
     708        if m is not None: 
     709            return True 
     710    return False 
    703711 
    704712def _add_source(source, code, path, lineno=1): 
    … …  
    730738    # dispersion.  Need to be careful that necessary parameters are available 
    731739    # for computing volume even if we allow non-disperse volume parameters. 
    732  
    733740    partable = model_info.parameters 
    734  
    735741    # Load templates and user code 
    736742    kernel_header = load_template('kernel_header.c') 
    737743    kernel_code = load_template('kernel_iq.c') 
    738744    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
    739  
    740745    # Build initial sources 
    741746    source = [] 
    … …  
    743748    for path, code in user_code: 
    744749        _add_source(source, code, path) 
    745  
    746750    if model_info.c_code: 
    747751        _add_source(source, model_info.c_code, model_info.filename, 
    … …  
    789793    source.append("\\\n".join(p.as_definition() 
    790794                              for p in partable.kernel_parameters)) 
    791  
    792795    # Define the function calls 
    793796    if partable.form_volume_parameters: 
    … …  
    800803        call_volume = "#define CALL_VOLUME(v) 1.0" 
    801804    source.append(call_volume) 
    802  
    803805    model_refs = _call_pars("_v.", partable.iq_parameters) 
    804     pars = ",".join(["_q"] + model_refs) 
    805     call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     806    #create varaible BETA to turn on and off beta approximation 
     807    BETA = has_Fq(source) 
     808    if not BETA: 
     809        pars = ",".join(["_q"] + model_refs) 
     810        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     811    else: 
     812        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs) 
     813        call_iq = "#define CALL_IQ(_q, _F1, _F2, _v) Fq(%s)" % pars 
    806814    if xy_mode == 'qabc': 
    807815        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    … …  
    831839    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    832840               if p.type == 'sld'] 
    833  
    834841    # Fill in definitions for numbers of parameters 
     842    source.append("#define BETA %d" %(1 if BETA else 0)) 
    835843    source.append("#define MAX_PD %s"%partable.max_pd) 
    836844    source.append("#define NUM_PARS %d"%partable.npars) 
    … …  
    839847    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    840848    source.append("#define PROJECTION %d"%PROJECTION) 
    841  
    842849    # TODO: allow mixed python/opencl kernels? 
    843  
    844850    ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    845851    dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
     852 
    846853    result = { 
    847854        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    848855        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    849856    } 
    850  
     857    #print(result['dll']) 
    851858    return result 
    852859 
    … …  
    10681075        model_info = make_model_info(kernel_module) 
    10691076        source = make_source(model_info) 
    1070         print(source['dll']) 
     1077        #print(source['dll']) 
    10711078 
    10721079 

  • sasmodels/kernel_iq.c

    r7c35fda r01c8d9e  
    3636//  PROJECTION : equirectangular=1, sinusoidal=2 
    3737//      see explore/jitter.py for definitions. 
     38 
    3839 
    3940#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
    … …  
    270271#endif // _QABC_SECTION 
    271272 
    272  
    273273// ==================== KERNEL CODE ======================== 
    274  
    275274kernel 
    276275void KERNEL_NAME( 
    … …  
    339338  // The code differs slightly between opencl and dll since opencl is only 
    340339  // seeing one q value (stored in the variable "this_result") while the dll 
    341   // version must loop over all q. 
     340  // version must loop over all q 
     341  #if BETA 
     342  double *beta_result = &(result[nq+1]); // = result + nq + 1 
     343  double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     344  #endif 
    342345  #ifdef USE_OPENCL 
    343346    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    344347    double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
     348    #if BETA 
     349      double this_beta_result = (pd_start == 0 ? 0.0 : result[nq + q_index]; 
    345350  #else // !USE_OPENCL 
    346     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     351    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);  
    347352    if (pd_start == 0) { 
    348353      #ifdef USE_OPENMP 
    … …  
    350355      #endif 
    351356      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     357      #if BETA 
     358      for (int q_index=0; q_index < nq; q_index++) beta_result[q_index] = 0.0; 
     359      #endif 
    352360    } 
    353361    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    354362#endif // !USE_OPENCL 
     363 
     364 
     365 
    355366 
    356367 
    … …  
    442453// inner loop and defines the macros that use them. 
    443454 
     455 
    444456#if defined(CALL_IQ) 
    445457  // unoriented 1D 
    446458  double qk; 
    447   #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    448   #define BUILD_ROTATION() do {} while(0) 
    449   #define APPLY_ROTATION() do {} while(0) 
    450   #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
     459  #if BETA == 0 
     460    #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     461    #define BUILD_ROTATION() do {} while(0) 
     462    #define APPLY_ROTATION() do {} while(0) 
     463    #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 
     464 
     465  // unoriented 1D Beta 
     466  #elif BETA == 1 
     467    double F1, F2; 
     468    #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     469    #define BUILD_ROTATION() do {} while(0) 
     470    #define APPLY_ROTATION() do {} while(0) 
     471    #define CALL_KERNEL() CALL_IQ(qk,F1,F2,local_values.table) 
     472  #endif 
    451473 
    452474#elif defined(CALL_IQ_A) 
    … …  
    647669      pd_norm += weight * CALL_VOLUME(local_values.table); 
    648670      BUILD_ROTATION(); 
    649  
     671#if BETA 
     672    if (weight > cutoff) { 
     673      weight_norm += weight;} 
     674#endif 
    650675#ifndef USE_OPENCL 
    651676      // DLL needs to explicitly loop over the q values. 
    … …  
    693718          } 
    694719        #else  // !MAGNETIC 
    695           const double scattering = CALL_KERNEL(); 
     720          #if defined(CALL_IQ) && BETA == 1 
     721            CALL_KERNEL(); 
     722            const double scatteringF1 = F1; 
     723            const double scatteringF2 = F2; 
     724          #else 
     725            const double scattering = CALL_KERNEL(); 
     726          #endif 
    696727        #endif // !MAGNETIC 
    697728//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 
    698729 
    699730        #ifdef USE_OPENCL 
    700           this_result += weight * scattering; 
     731          #if defined(CALL_IQ)&& BETA == 1 
     732             this_result += weight * scatteringF2; 
     733             this_beta_result += weight * scatteringF1; 
     734            #else 
     735              this_result += weight * scattering; 
     736          #endif 
    701737        #else // !USE_OPENCL 
    702           result[q_index] += weight * scattering; 
     738          #if defined(CALL_IQ)&& BETA == 1 
     739            result[q_index] += weight * scatteringF2; 
     740            beta_result[q_index] += weight * scatteringF1; 
     741            #endif 
     742            #else 
     743            result[q_index] += weight * scattering; 
     744          #endif 
    703745        #endif // !USE_OPENCL 
    704746      } 
    705747    } 
    706748  } 
    707  
    708749// close nested loops 
    709750++step; 
    … …  
    728769  result[q_index] = this_result; 
    729770  if (q_index == 0) result[nq] = pd_norm; 
     771  #if BETA 
     772  beta_result[q_index] = this_beta_result; 
     773  #endif 
     774  if (q_index == 0) result[nq] = pd_norm; 
     775 
    730776//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    731777#else // !USE_OPENCL 
    732778  result[nq] = pd_norm; 
     779  #if BETA 
     780  result[2*nq+1] = weight_norm; 
     781  #endif 
    733782//printf("res: %g/%g\n", result[0], pd_norm); 
    734783#endif // !USE_OPENCL 

  • sasmodels/kernelcl.py

    rd86f0fc r01c8d9e  
    420420        if self.program is None: 
    421421            compile_program = environment().compile_program 
     422            with open('model.c','w') as fid: 
     423                print(self.source['opencl'], file=fid) 
    422424            timestamp = generate.ocl_timestamp(self.info) 
    423425            self.program = compile_program( 
    … …  
    539541        self.dim = '2d' if q_input.is_2d else '1d' 
    540542        # plus three for the normalization values 
    541         self.result = np.empty(q_input.nq+1, dtype) 
     543        self.result = np.empty(2*q_input.nq+2,dtype) 
    542544 
    543545        # Inputs and outputs for each kernel call 
    … …  
    555557                     else np.float16 if dtype == generate.F16 
    556558                     else np.float32)  # will never get here, so use np.float32 
    557  
    558     def __call__(self, call_details, values, cutoff, magnetic): 
     559    __call__= Iq 
     560 
     561    def Iq(self, call_details, values, cutoff, magnetic): 
    559562        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    560563        context = self.queue.context 
    … …  
    572575        ] 
    573576        #print("Calling OpenCL") 
    574         #call_details.show(values) 
    575         # Call kernel and retrieve results 
     577        call_details.show(values) 
     578        #Call kernel and retrieve results 
    576579        wait_for = None 
    577580        last_nap = time.clock() 
    … …  
    604607        return scale*self.result[:self.q_input.nq] + background 
    605608        # return self.result[:self.q_input.nq] 
     609     #NEEDS TO BE FINISHED FOR OPENCL 
     610     def beta(): 
     611         return 0 
    606612 
    607613    def release(self): 

  • sasmodels/kerneldll.py

    r33969b6 r01c8d9e  
    258258        # Note: if there is a syntax error then compile raises an error 
    259259        # and the source file will not be deleted. 
    260         os.unlink(filename) 
    261         #print("saving compiled file in %r"%filename) 
     260        #os.unlink(filename) 
     261        print("saving compiled file in %r"%filename) 
    262262    return dll 
    263263 
    … …  
    371371    def __init__(self, kernel, model_info, q_input): 
    372372        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
     373        #,model_info,q_input) 
    373374        self.kernel = kernel 
    374375        self.info = model_info 
    … …  
    376377        self.dtype = q_input.dtype 
    377378        self.dim = '2d' if q_input.is_2d else '1d' 
    378         self.result = np.empty(q_input.nq+1, q_input.dtype) 
     379        self.result = np.empty(2*q_input.nq+2, q_input.dtype) 
    379380        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    380381                     else np.float64 if self.q_input.dtype == generate.F64 
    381382                     else np.float128) 
    382383 
    383     def __call__(self, call_details, values, cutoff, magnetic): 
     384    def Iq(self, call_details, values, cutoff, magnetic): 
    384385        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    385  
     386        self._call_kernel(call_details, values, cutoff, magnetic) 
     387        #print("returned",self.q_input.q, self.result) 
     388        pd_norm = self.result[self.q_input.nq] 
     389        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
     390        background = values[1] 
     391        #print("scale",scale,background) 
     392        return scale*self.result[:self.q_input.nq] + background 
     393    __call__ = Iq 
     394 
     395    def beta(self, call_details, values, cutoff, magnetic): 
     396        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
     397        self._call_kernel(call_details, values, cutoff, magnetic) 
     398        w_norm = self.result[2*self.q_input.nq + 1] 
     399        pd_norm = self.result[self.q_input.nq] 
     400        if w_norm == 0.: 
     401            w_norm = 1. 
     402        F2 = self.result[:self.q_input.nq]/w_norm 
     403        F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
     404        volume_avg = pd_norm/w_norm 
     405        return F1, F2, volume_avg 
     406 
     407    def _call_kernel(self, call_details, values, cutoff, magnetic): 
    386408        kernel = self.kernel[1 if magnetic else 0] 
    387409        args = [ 
    … …  
    390412            None, # pd_stop pd_stride[MAX_PD] 
    391413            call_details.buffer.ctypes.data, # problem 
    392             values.ctypes.data,  #pars 
    393             self.q_input.q.ctypes.data, #q 
     414            values.ctypes.data,  # pars 
     415            self.q_input.q.ctypes.data, # q 
    394416            self.result.ctypes.data,   # results 
    395417            self.real(cutoff), # cutoff 
    396418        ] 
    397         #print("Calling DLL") 
     419        #print(self.beta_result.ctypes.data) 
     420#        print("Calling DLL line 397") 
     421#        print("values", values) 
    398422        #call_details.show(values) 
    399423        step = 100 
    … …  
    403427            kernel(*args) # type: ignore 
    404428 
    405         #print("returned",self.q_input.q, self.result) 
    406         pd_norm = self.result[self.q_input.nq] 
    407         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    408         background = values[1] 
    409         #print("scale",scale,background) 
    410         return scale*self.result[:self.q_input.nq] + background 
    411  
    412429    def release(self): 
    413430        # type: () -> None 

  • sasmodels/modelinfo.py

    r95498a3 r01c8d9e  
    163163    parameter.length = length 
    164164    parameter.length_control = control 
    165  
    166165    return parameter 
    167166 
    … …  
    418417    # scale and background are implicit parameters 
    419418    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
    420  
    421419    def __init__(self, parameters): 
    422420        # type: (List[Parameter]) -> None 
    423421        self.kernel_parameters = parameters 
    424422        self._set_vector_lengths() 
    425  
    426423        self.npars = sum(p.length for p in self.kernel_parameters) 
    427424        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
    … …  
    430427        if self.nmagnetic: 
    431428            self.nvalues += 3 + 3*self.nmagnetic 
    432  
    433429        self.call_parameters = self._get_call_parameters() 
    434430        self.defaults = self._get_defaults() 
    … …  
    444440        self.form_volume_parameters = [p for p in self.kernel_parameters 
    445441                                       if p.type == 'volume'] 
    446  
    447442        # Theta offset 
    448443        offset = 0 
    … …  
    466461        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    467462                                if p.id.startswith('M0:')] 
    468  
    469463        self.pd_1d = set(p.name for p in self.call_parameters 
    470464                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    … …  
    770764        # Custom sum/multi models 
    771765        return kernel_module.model_info 
     766 
    772767    info = ModelInfo() 
    773768    #print("make parameter table", kernel_module.parameters) 
    … …  
    822817    info.lineno = {} 
    823818    _find_source_lines(info, kernel_module) 
    824  
    825819    return info 
    826820 

  • sasmodels/models/ellipsoid.c

    r108e70e r01c8d9e  
    2020    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
    2121    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
    22  
     22   
    2323    // translate a point in [-1,1] to a point in [0, 1] 
    2424    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    … …  
    3636    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    3737    return 1.0e-4 * s * s * form; 
     38}  
     39 
     40static void 
     41Fq(double q, 
     42    double *F1, 
     43    double *F2, 
     44    double sld, 
     45    double sld_solvent, 
     46    double radius_polar, 
     47    double radius_equatorial) 
     48{ 
     49    // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955) 
     50    //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT 
     51    //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT 
     52    //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT 
     53    // u-substitution of 
     54    //     u = sin, du = cos dT 
     55    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
     56    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
     57    // translate a point in [-1,1] to a point in [0, 1] 
     58    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     59    const double zm = 0.5; 
     60    const double zb = 0.5; 
     61    double total_F2 = 0.0; 
     62    double total_F1 = 0.0; 
     63    for (int i=0;i<GAUSS_N;i++) { 
     64        const double u = GAUSS_Z[i]*zm + zb; 
     65        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
     66        const double f = sas_3j1x_x(q*r); 
     67        total_F2 += GAUSS_W[i] * f * f; 
     68        total_F1 += GAUSS_W[i] * f; 
     69    } 
     70    // translate dx in [-1,1] to dx in [lower,upper] 
     71    const double form_squared_avg = total_F2*zm; 
     72    const double form_avg = total_F1*zm; 
     73    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
     74    *F2 = 1e-4 * s * s * form_squared_avg; 
     75    *F1 = 1e-2 * s * form_avg; 
    3876} 
     77 
     78 
     79 
     80 
     81 
    3982 
    4083static double 

  • sasmodels/models/lib/sphere_form.c

    r925ad6e r01c8d9e  
    1313    return 1.0e-4*square(contrast * fq); 
    1414} 
     15 

  • sasmodels/models/sphere.py

    ref07e95 r01c8d9e  
    6767             ] 
    6868 
    69 source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"] 
     69source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c", "sphere.c"] 
    7070 
    71 # No volume normalization despite having a volume parameter 
    72 # This should perhaps be volume normalized? 
    73 form_volume = """ 
    74     return sphere_volume(radius); 
    75     """ 
    76  
    77 Iq = """ 
    78     return sphere_form(q, radius, sld, sld_solvent); 
    79     """ 
    8071 
    8172def ER(radius): 

  • sasmodels/product.py

    r2d81cfe r01c8d9e  
    1616import numpy as np  # type: ignore 
    1717 
    18 from .modelinfo import ParameterTable, ModelInfo 
     18from .modelinfo import ParameterTable, ModelInfo, Parameter 
    1919from .kernel import KernelModel, Kernel 
    2020from .details import make_details, dispersion_mesh 
    … …  
    7474    translate_name = dict((old.id, new.id) for old, new 
    7575                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    76     combined_pars = p_pars.kernel_parameters + s_list 
     76    beta_parameter = Parameter("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"], "", "Structure factor dispersion calculation mode"]) 
     77    combined_pars = p_pars.kernel_parameters + s_list + [beta_parameter] 
    7778    parameters = ParameterTable(combined_pars) 
    7879    parameters.max_pd = p_pars.max_pd + s_pars.max_pd 
    … …  
    151152        #: Structure factor modelling interaction between particles. 
    152153        self.S = S 
     154         
    153155        #: Model precision. This is not really relevant, since it is the 
    154156        #: individual P and S models that control the effective dtype, 
    … …  
    168170        # in opencl; or both in opencl, but one in single precision and the 
    169171        # other in double precision). 
     172         
    170173        p_kernel = self.P.make_kernel(q_vectors) 
    171174        s_kernel = self.S.make_kernel(q_vectors) 
    … …  
    193196        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray 
    194197        p_info, s_info = self.info.composition[1] 
    195  
    196198        # if there are magnetic parameters, they will only be on the 
    197199        # form factor P, not the structure factor S. 
    … …  
    205207        nweights = call_details.num_weights 
    206208        weights = values[nvalues:nvalues + 2*nweights] 
    207  
    208209        # Construct the calling parameters for P. 
    209210        p_npars = p_info.parameters.npars 
    … …  
    218219        p_values.append([0.]*spacer) 
    219220        p_values = np.hstack(p_values).astype(self.p_kernel.dtype) 
    220  
    221221        # Call ER and VR for P since these are needed for S. 
    222222        p_er, p_vr = calc_er_vr(p_info, p_details, p_values) 
    223223        s_vr = (volfrac/p_vr if p_vr != 0. else volfrac) 
    224224        #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr)) 
    225  
    226225        # Construct the calling parameters for S. 
    227226        # The  effective radius is not in the combined parameter list, so 
    … …  
    253252        s_values.append([0.]*spacer) 
    254253        s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 
    255  
     254        # beta mode is the first parameter after the structure factor pars 
     255        beta_index = 2+p_npars+s_npars 
     256        beta_mode = values[beta_index] 
    256257        # Call the kernels 
    257         p_result = self.p_kernel(p_details, p_values, cutoff, magnetic) 
    258         s_result = self.s_kernel(s_details, s_values, cutoff, False) 
    259  
     258        if beta_mode: # beta: 
     259            F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic) 
     260        else: 
     261            p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic) 
     262        s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
    260263        #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars]) 
    261264        #call_details.show(values) 
    … …  
    265268        #s_details.show(s_values) 
    266269        #print("=>", s_result) 
    267  
    268         # remember the parts for plotting later 
    269         self.results = [p_result, s_result] 
    270  
    271270        #import pylab as plt 
    272271        #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') 
    273272        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') 
    274273        #plt.figure() 
    275  
    276         return values[0]*(p_result*s_result) + values[1] 
     274        if beta_mode:#beta 
     275            beta_factor = F1**2/F2 
     276            Sq_eff = 1+beta_factor*(s_result - 1) 
     277            self.results = [F2, Sq_eff,F1,s_result] 
     278            final_result = volfrac*values[0]*(F2 + (F1**2)*(s_result - 1))/volume_avg+values[1] 
     279            #final_result =  volfrac * values[0] * F2 * Sq_eff / volume_avg + values[1] 
     280        else: 
     281            # remember the parts for plotting later 
     282            self.results = [p_result, s_result] 
     283            final_result = values[0]*(p_result*s_result) + values[1] 
     284        return final_result 
    277285 
    278286    def release(self): 

  • doc/guide/pd/polydispersity.rst

    r29afc50 rd712a0f  
    2020  P(q) = \text{scale} \langle F^* F \rangle / V + \text{background} 
    2121 
    22 where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an 
    23 average over the size distribution. 
     22where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an  
     23average over the size distribution $f(x; \bar x, \sigma)$, giving 
     24 
     25.. math:: 
     26 
     27  P(q) = \frac{\text{scale}}{V} \int_\mathbb{R}  
     28  f(x; \bar x, \sigma) F^2(q, x)\, dx + \text{background} 
    2429 
    2530Each distribution is characterized by a center value $\bar x$ or 
    … …  
    4146with larger values of $N_\sigma$ required for heavier tailed distributions. 
    4247The scattering in general falls rapidly with $qr$ so the usual assumption 
    43 that $G(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)G(r - 3\sigma_r)$ 
     48that $f(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)f(r - 3\sigma_r)$ 
    4449will not contribute much to the average may not hold when particles are large. 
    4550This, too, will require increasing $N_\sigma$. 
    … …  
    6368 
    6469Additional distributions are under consideration. 
     70 
     71.. note:: In 2009 IUPAC decided to introduce the new term 'dispersity' to replace  
     72           the term 'polydispersity' (see `Pure Appl. Chem., (2009), 81(2),  
     73           351-353 <http://media.iupac.org/publications/pac/2009/pdf/8102x0351.pdf>`_  
     74           in order to make the terminology describing distributions of properties  
     75           unambiguous. Throughout the SasView documentation we continue to use the  
     76           term polydispersity because one of the consequences of the IUPAC change is  
     77           that orientational polydispersity would not meet their new criteria (which  
     78           requires dispersity to be dimensionless). 
    6579 
    6680Suggested Applications 

  • sasmodels/models/core_shell_sphere.py

    r2d81cfe rdc76240  
    2121.. math:: 
    2222 
    23     F^2(q) = \frac{3}{V_s}\left[ 
     23    F(q) = \frac{3}{V_s}\left[ 
    2424       V_c(\rho_c-\rho_s)\frac{\sin(qr_c)-qr_c\cos(qr_c)}{(qr_c)^3} + 
    2525       V_s(\rho_s-\rho_\text{solv})\frac{\sin(qr_s)-qr_s\cos(qr_s)}{(qr_s)^3} 
Note: See TracChangeset for help on using the changeset viewer.