Changes in / [502c7b8:cdd676e] in sasmodels


Ignore:
Files:
1 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • explore/beta/sasfit_compare.py

    r01c8d9e rcdd676e  
    33import sys, os 
    44BETA_DIR = os.path.dirname(os.path.realpath(__file__)) 
    5 #SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 
    6 SASMODELS_DIR = r"C:\Source\sasmodels" 
     5SASMODELS_DIR = os.path.dirname(os.path.dirname(BETA_DIR)) 
    76sys.path.insert(0, SASMODELS_DIR) 
    87 
    … …  
    6261    calculator = DirectModel(data, model,cutoff=0) 
    6362    calculator.pars = pars.copy() 
    64     calculator.pars.setdefault('background', 0) 
     63    calculator.pars.setdefault('background', 1e-3) 
    6564    return calculator 
    6665 
    … …  
    232231        #     = F2/total_volume 
    233232        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 
    236233    elif norm == 'yun': 
    237234        F1 *= 1e-6  # Yun is using sld in 1/A^2, not 1e-6/A^2 
    … …  
    339336    I = build_model(Pname+"@hardsphere", q) 
    340337    Pq = P(**Ppars)*pars.get('volfraction', 1) 
    341     Sq = S(**Spars) 
     338    #Sq = S(**Spars) 
    342339    Iq = I(**Ipars) 
    343340    #Iq = Pq*Sq*pars.get('volfraction', 1) 
    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) 
     341    Sq = Iq/Pq 
     342    return Theory(Q=q, F1=None, F2=None, P=Pq, S=Sq, I=Iq, Seff=None, Ibeta=None) 
    348343 
    349344def compare(title, target, actual, fields='F1 F2 P S I Seff Ibeta'): 
    … …  
    405400        radius_equatorial_pd=.1, radius_equatorial_pd_type=pd_type, 
    406401        volfraction=0.15, 
    407         radius_effective=270.7543927018, 
    408402        #radius_effective=12.59921049894873, 
    409403        ) 
    410     target = sasmodels_theory(q, model, beta_mode=1, **pars) 
     404    target = sasmodels_theory(q, model, **pars) 
    411405    actual = ellipsoid_pe(q, norm='sasview', **pars) 
    412406    title = " ".join(("sasmodels", model, pd_type)) 
    … …  
    456450    S = data[5] 
    457451    Seff = data[6] 
    458     target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 
     452    target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, Seff=Seff) 
    459453    actual = sphere_r(Q, norm='yun', **pars) 
    460454    title = " ".join(("yun", "sphere", "10% dispersion 10% Vf")) 
    … …  
    468462    S = data[5] 
    469463    Seff = data[6] 
    470     target = Theory(Q=Q, F1=F1, P=P, S=S, Seff=Seff) 
     464    target = Theory(Q=Q, F1=F1, F2=F2, P=P, S=S, Seff=Seff) 
    471465    actual = sphere_r(Q, norm='yun', **pars) 
    472466    title = " ".join(("yun", "sphere", "15% dispersion 10% Vf")) 

  • sasmodels/core.py

    r01c8d9e r3221de0  
    140140    used with functions in generate to build the docs or extract model info. 
    141141    """ 
    142  
    143142    if "+" in model_string: 
    144143        parts = [load_model_info(part) 
    … …  
    206205 
    207206    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 
    267268    if platform is None: 
    268269        platform = "ocl" 
    … …  
    290291    else: 
    291292        numpy_dtype = np.dtype(dtype) 
     293 
    292294    # Make sure that the type is supported by opencl, otherwise use dll 
    293295    if platform == "ocl": 

  • sasmodels/data.py

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

  • sasmodels/generate.py

    r01c8d9e rd86f0fc  
    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(ac|abc|xy))\s*[(]", 
     674_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|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) 
    705 def 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 
    711703 
    712704def _add_source(source, code, path, lineno=1): 
    … …  
    738730    # dispersion.  Need to be careful that necessary parameters are available 
    739731    # for computing volume even if we allow non-disperse volume parameters. 
     732 
    740733    partable = model_info.parameters 
     734 
    741735    # Load templates and user code 
    742736    kernel_header = load_template('kernel_header.c') 
    743737    kernel_code = load_template('kernel_iq.c') 
    744738    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
     739 
    745740    # Build initial sources 
    746741    source = [] 
    … …  
    748743    for path, code in user_code: 
    749744        _add_source(source, code, path) 
     745 
    750746    if model_info.c_code: 
    751747        _add_source(source, model_info.c_code, model_info.filename, 
    … …  
    793789    source.append("\\\n".join(p.as_definition() 
    794790                              for p in partable.kernel_parameters)) 
     791 
    795792    # Define the function calls 
    796793    if partable.form_volume_parameters: 
    … …  
    803800        call_volume = "#define CALL_VOLUME(v) 1.0" 
    804801    source.append(call_volume) 
     802 
    805803    model_refs = _call_pars("_v.", partable.iq_parameters) 
    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 
     804    pars = ",".join(["_q"] + model_refs) 
     805    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
    814806    if xy_mode == 'qabc': 
    815807        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    … …  
    839831    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    840832               if p.type == 'sld'] 
     833 
    841834    # Fill in definitions for numbers of parameters 
    842     source.append("#define BETA %d" %(1 if BETA else 0)) 
    843835    source.append("#define MAX_PD %s"%partable.max_pd) 
    844836    source.append("#define NUM_PARS %d"%partable.npars) 
    … …  
    847839    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    848840    source.append("#define PROJECTION %d"%PROJECTION) 
     841 
    849842    # TODO: allow mixed python/opencl kernels? 
     843 
    850844    ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    851845    dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    852  
    853846    result = { 
    854847        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    855848        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    856849    } 
    857     #print(result['dll']) 
     850 
    858851    return result 
    859852 
    … …  
    10751068        model_info = make_model_info(kernel_module) 
    10761069        source = make_source(model_info) 
    1077         #print(source['dll']) 
     1070        print(source['dll']) 
    10781071 
    10791072 

  • sasmodels/kernel_iq.c

    r01c8d9e r7c35fda  
    3636//  PROJECTION : equirectangular=1, sinusoidal=2 
    3737//      see explore/jitter.py for definitions. 
    38  
    3938 
    4039#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
    … …  
    271270#endif // _QABC_SECTION 
    272271 
     272 
    273273// ==================== KERNEL CODE ======================== 
     274 
    274275kernel 
    275276void KERNEL_NAME( 
    … …  
    338339  // The code differs slightly between opencl and dll since opencl is only 
    339340  // seeing one q value (stored in the variable "this_result") while the dll 
    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 
     341  // version must loop over all q. 
    345342  #ifdef USE_OPENCL 
    346343    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    347344    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]; 
    350345  #else // !USE_OPENCL 
    351     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]);  
     346    double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    352347    if (pd_start == 0) { 
    353348      #ifdef USE_OPENMP 
    … …  
    355350      #endif 
    356351      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 
    360352    } 
    361353    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    362354#endif // !USE_OPENCL 
    363  
    364  
    365  
    366355 
    367356 
    … …  
    453442// inner loop and defines the macros that use them. 
    454443 
    455  
    456444#if defined(CALL_IQ) 
    457445  // unoriented 1D 
    458446  double qk; 
    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 
     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) 
    473451 
    474452#elif defined(CALL_IQ_A) 
    … …  
    669647      pd_norm += weight * CALL_VOLUME(local_values.table); 
    670648      BUILD_ROTATION(); 
    671 #if BETA 
    672     if (weight > cutoff) { 
    673       weight_norm += weight;} 
    674 #endif 
     649 
    675650#ifndef USE_OPENCL 
    676651      // DLL needs to explicitly loop over the q values. 
    … …  
    718693          } 
    719694        #else  // !MAGNETIC 
    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 
     695          const double scattering = CALL_KERNEL(); 
    727696        #endif // !MAGNETIC 
    728697//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 
    729698 
    730699        #ifdef USE_OPENCL 
    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 
     700          this_result += weight * scattering; 
    737701        #else // !USE_OPENCL 
    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 
     702          result[q_index] += weight * scattering; 
    745703        #endif // !USE_OPENCL 
    746704      } 
    747705    } 
    748706  } 
     707 
    749708// close nested loops 
    750709++step; 
    … …  
    769728  result[q_index] = this_result; 
    770729  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  
    776730//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    777731#else // !USE_OPENCL 
    778732  result[nq] = pd_norm; 
    779   #if BETA 
    780   result[2*nq+1] = weight_norm; 
    781   #endif 
    782733//printf("res: %g/%g\n", result[0], pd_norm); 
    783734#endif // !USE_OPENCL 

  • sasmodels/kernelcl.py

    r01c8d9e rd86f0fc  
    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) 
    424422            timestamp = generate.ocl_timestamp(self.info) 
    425423            self.program = compile_program( 
    … …  
    541539        self.dim = '2d' if q_input.is_2d else '1d' 
    542540        # plus three for the normalization values 
    543         self.result = np.empty(2*q_input.nq+2,dtype) 
     541        self.result = np.empty(q_input.nq+1, dtype) 
    544542 
    545543        # Inputs and outputs for each kernel call 
    … …  
    557555                     else np.float16 if dtype == generate.F16 
    558556                     else np.float32)  # will never get here, so use np.float32 
    559     __call__= Iq 
    560  
    561     def Iq(self, call_details, values, cutoff, magnetic): 
     557 
     558    def __call__(self, call_details, values, cutoff, magnetic): 
    562559        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    563560        context = self.queue.context 
    … …  
    575572        ] 
    576573        #print("Calling OpenCL") 
    577         call_details.show(values) 
    578         #Call kernel and retrieve results 
     574        #call_details.show(values) 
     575        # Call kernel and retrieve results 
    579576        wait_for = None 
    580577        last_nap = time.clock() 
    … …  
    607604        return scale*self.result[:self.q_input.nq] + background 
    608605        # return self.result[:self.q_input.nq] 
    609      #NEEDS TO BE FINISHED FOR OPENCL 
    610      def beta(): 
    611          return 0 
    612606 
    613607    def release(self): 

  • sasmodels/kerneldll.py

    r01c8d9e r33969b6  
    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) 
    374373        self.kernel = kernel 
    375374        self.info = model_info 
    … …  
    377376        self.dtype = q_input.dtype 
    378377        self.dim = '2d' if q_input.is_2d else '1d' 
    379         self.result = np.empty(2*q_input.nq+2, q_input.dtype) 
     378        self.result = np.empty(q_input.nq+1, q_input.dtype) 
    380379        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    381380                     else np.float64 if self.q_input.dtype == generate.F64 
    382381                     else np.float128) 
    383382 
    384     def Iq(self, call_details, values, cutoff, magnetic): 
     383    def __call__(self, call_details, values, cutoff, magnetic): 
    385384        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    386         self._call_kernel(call_details, values, cutoff, magnetic) 
     385 
     386        kernel = self.kernel[1 if magnetic else 0] 
     387        args = [ 
     388            self.q_input.nq, # nq 
     389            None, # pd_start 
     390            None, # pd_stop pd_stride[MAX_PD] 
     391            call_details.buffer.ctypes.data, # problem 
     392            values.ctypes.data,  #pars 
     393            self.q_input.q.ctypes.data, #q 
     394            self.result.ctypes.data,   # results 
     395            self.real(cutoff), # cutoff 
     396        ] 
     397        #print("Calling DLL") 
     398        #call_details.show(values) 
     399        step = 100 
     400        for start in range(0, call_details.num_eval, step): 
     401            stop = min(start + step, call_details.num_eval) 
     402            args[1:3] = [start, stop] 
     403            kernel(*args) # type: ignore 
     404 
    387405        #print("returned",self.q_input.q, self.result) 
    388406        pd_norm = self.result[self.q_input.nq] 
    … …  
    391409        #print("scale",scale,background) 
    392410        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): 
    408         kernel = self.kernel[1 if magnetic else 0] 
    409         args = [ 
    410             self.q_input.nq, # nq 
    411             None, # pd_start 
    412             None, # pd_stop pd_stride[MAX_PD] 
    413             call_details.buffer.ctypes.data, # problem 
    414             values.ctypes.data,  # pars 
    415             self.q_input.q.ctypes.data, # q 
    416             self.result.ctypes.data,   # results 
    417             self.real(cutoff), # cutoff 
    418         ] 
    419         #print(self.beta_result.ctypes.data) 
    420 #        print("Calling DLL line 397") 
    421 #        print("values", values) 
    422         #call_details.show(values) 
    423         step = 100 
    424         for start in range(0, call_details.num_eval, step): 
    425             stop = min(start + step, call_details.num_eval) 
    426             args[1:3] = [start, stop] 
    427             kernel(*args) # type: ignore 
    428411 
    429412    def release(self): 

  • sasmodels/modelinfo.py

    r01c8d9e r95498a3  
    163163    parameter.length = length 
    164164    parameter.length_control = control 
     165 
    165166    return parameter 
    166167 
    … …  
    417418    # scale and background are implicit parameters 
    418419    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
     420 
    419421    def __init__(self, parameters): 
    420422        # type: (List[Parameter]) -> None 
    421423        self.kernel_parameters = parameters 
    422424        self._set_vector_lengths() 
     425 
    423426        self.npars = sum(p.length for p in self.kernel_parameters) 
    424427        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
    … …  
    427430        if self.nmagnetic: 
    428431            self.nvalues += 3 + 3*self.nmagnetic 
     432 
    429433        self.call_parameters = self._get_call_parameters() 
    430434        self.defaults = self._get_defaults() 
    … …  
    440444        self.form_volume_parameters = [p for p in self.kernel_parameters 
    441445                                       if p.type == 'volume'] 
     446 
    442447        # Theta offset 
    443448        offset = 0 
    … …  
    461466        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    462467                                if p.id.startswith('M0:')] 
     468 
    463469        self.pd_1d = set(p.name for p in self.call_parameters 
    464470                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    … …  
    764770        # Custom sum/multi models 
    765771        return kernel_module.model_info 
    766  
    767772    info = ModelInfo() 
    768773    #print("make parameter table", kernel_module.parameters) 
    … …  
    817822    info.lineno = {} 
    818823    _find_source_lines(info, kernel_module) 
     824 
    819825    return info 
    820826 

  • sasmodels/models/ellipsoid.c

    r01c8d9e r108e70e  
    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  
    40 static void 
    41 Fq(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; 
    7638} 
    77  
    78  
    79  
    80  
    81  
    8239 
    8340static double 

  • sasmodels/models/lib/sphere_form.c

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

  • sasmodels/models/sphere.py

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

  • sasmodels/product.py

    r01c8d9e r2d81cfe  
    1616import numpy as np  # type: ignore 
    1717 
    18 from .modelinfo import ParameterTable, ModelInfo, Parameter 
     18from .modelinfo import ParameterTable, ModelInfo 
    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     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] 
     76    combined_pars = p_pars.kernel_parameters + s_list 
    7877    parameters = ParameterTable(combined_pars) 
    7978    parameters.max_pd = p_pars.max_pd + s_pars.max_pd 
    … …  
    152151        #: Structure factor modelling interaction between particles. 
    153152        self.S = S 
    154          
    155153        #: Model precision. This is not really relevant, since it is the 
    156154        #: individual P and S models that control the effective dtype, 
    … …  
    170168        # in opencl; or both in opencl, but one in single precision and the 
    171169        # other in double precision). 
    172          
    173170        p_kernel = self.P.make_kernel(q_vectors) 
    174171        s_kernel = self.S.make_kernel(q_vectors) 
    … …  
    196193        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray 
    197194        p_info, s_info = self.info.composition[1] 
     195 
    198196        # if there are magnetic parameters, they will only be on the 
    199197        # form factor P, not the structure factor S. 
    … …  
    207205        nweights = call_details.num_weights 
    208206        weights = values[nvalues:nvalues + 2*nweights] 
     207 
    209208        # Construct the calling parameters for P. 
    210209        p_npars = p_info.parameters.npars 
    … …  
    219218        p_values.append([0.]*spacer) 
    220219        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 
    225226        # Construct the calling parameters for S. 
    226227        # The  effective radius is not in the combined parameter list, so 
    … …  
    252253        s_values.append([0.]*spacer) 
    253254        s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 
    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] 
     255 
    257256        # Call the kernels 
    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) 
     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 
    263260        #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars]) 
    264261        #call_details.show(values) 
    … …  
    268265        #s_details.show(s_values) 
    269266        #print("=>", s_result) 
     267 
     268        # remember the parts for plotting later 
     269        self.results = [p_result, s_result] 
     270 
    270271        #import pylab as plt 
    271272        #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') 
    272273        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') 
    273274        #plt.figure() 
    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 
     275 
     276        return values[0]*(p_result*s_result) + values[1] 
    285277 
    286278    def release(self): 
Note: See TracChangeset for help on using the changeset viewer.