Changes in / [cdd676e:502c7b8] in sasmodels


Ignore:
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed

  • explore/beta/sasfit_compare.py

    rcdd676e r01c8d9e  
    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): 
Note: See TracChangeset for help on using the changeset viewer.