Changeset 975ec8e in sasview for sansmodels/src/sans/models/c_models


Ignore:
Timestamp:
Aug 14, 2009 5:58:58 PM (15 years ago)
Author:
Jae Cho <jhjcho@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
d5bd424
Parents:
cad821b
Message:

working on 2D models. Still need smore corrections and unit tests.

Location:
sansmodels/src/sans/models/c_models
Files:
2 added
18 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/sans/models/c_models/CBinaryHSModel.cpp

    r9bd69098 r975ec8e  
    8686         
    8787        // Initialize parameter dictionary 
    88         PyDict_SetItemString(self->params,"vol_frac_ls",Py_BuildValue("d",0.200000)); 
     88        PyDict_SetItemString(self->params,"vol_frac_ls",Py_BuildValue("d",0.100000)); 
    8989        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.001000)); 
    9090        PyDict_SetItemString(self->params,"vol_frac_ss",Py_BuildValue("d",0.200000)); 
    … …  
    9393        PyDict_SetItemString(self->params,"ss_sld",Py_BuildValue("d",0.000000)); 
    9494        PyDict_SetItemString(self->params,"s_radius",Py_BuildValue("d",25.000000)); 
    95         PyDict_SetItemString(self->params,"l_radius",Py_BuildValue("d",160.000000)); 
     95        PyDict_SetItemString(self->params,"l_radius",Py_BuildValue("d",100.000000)); 
    9696        // Initialize dispersion / averaging parameter dict 
    9797        DispersionVisitor* visitor = new DispersionVisitor(); 
    … …  
    170170    return PyArray_Return(result);  
    171171 } 
    172 /** 
    173  * Function to call to evaluate model 
    174  * @param args: input numpy array  [q[],phi[]] 
    175  * @return: numpy array object  
    176  */ 
    177 static PyObject * evaluateTwoDim( BinaryHSModel* model,  
    178                               PyArrayObject *q, PyArrayObject *phi) 
    179  { 
    180     PyArrayObject *result; 
    181     //check validity of input vectors 
    182     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    183         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    184         || phi->dimensions[0] != q->dimensions[0]){ 
    185       
    186         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    187         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    188         return NULL; 
    189     } 
    190         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    191  
    192         if (result == NULL){ 
    193             const char * message= "Could not create result "; 
    194         PyErr_SetString(PyExc_RuntimeError , message); 
    195             return NULL; 
    196         } 
    197          
    198     for (int i = 0; i < q->dimensions[0]; i++) { 
    199       double q_value = *(double *)(q->data + i*q->strides[0]); 
    200       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    201       double *result_value = (double *)(result->data + i*result->strides[0]); 
    202       if (q_value == 0) 
    203           *result_value = 0.0; 
    204       else 
    205           *result_value = model->evaluate_rphi(q_value, phi_value); 
    206     } 
    207     return PyArray_Return(result);  
    208  } 
     172 
    209173 /** 
    210174 * Function to call to evaluate model 

  • sansmodels/src/sans/models/c_models/CCylinderModel.cpp

    r9bd69098 r975ec8e  
    175175    return PyArray_Return(result);  
    176176 } 
    177 /** 
    178  * Function to call to evaluate model 
    179  * @param args: input numpy array  [q[],phi[]] 
    180  * @return: numpy array object  
    181  */ 
    182 static PyObject * evaluateTwoDim( CylinderModel* model,  
    183                               PyArrayObject *q, PyArrayObject *phi) 
    184  { 
    185     PyArrayObject *result; 
    186     //check validity of input vectors 
    187     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    188         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    189         || phi->dimensions[0] != q->dimensions[0]){ 
    190       
    191         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    192         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    193         return NULL; 
    194     } 
    195         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    196  
    197         if (result == NULL){ 
    198             const char * message= "Could not create result "; 
    199         PyErr_SetString(PyExc_RuntimeError , message); 
    200             return NULL; 
    201         } 
    202          
    203     for (int i = 0; i < q->dimensions[0]; i++) { 
    204       double q_value = *(double *)(q->data + i*q->strides[0]); 
    205       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    206       double *result_value = (double *)(result->data + i*result->strides[0]); 
    207       if (q_value == 0) 
    208           *result_value = 0.0; 
    209       else 
    210           *result_value = model->evaluate_rphi(q_value, phi_value); 
    211     } 
    212     return PyArray_Return(result);  
    213  } 
     177 
    214178 /** 
    215179 * Function to call to evaluate model 

  • sansmodels/src/sans/models/c_models/CEllipsoidModel.cpp

    r9bd69098 r975ec8e  
    175175    return PyArray_Return(result);  
    176176 } 
    177 /** 
    178  * Function to call to evaluate model 
    179  * @param args: input numpy array  [q[],phi[]] 
    180  * @return: numpy array object  
    181  */ 
    182 static PyObject * evaluateTwoDim( EllipsoidModel* model,  
    183                               PyArrayObject *q, PyArrayObject *phi) 
    184  { 
    185     PyArrayObject *result; 
    186     //check validity of input vectors 
    187     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    188         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    189         || phi->dimensions[0] != q->dimensions[0]){ 
    190       
    191         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    192         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    193         return NULL; 
    194     } 
    195         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    196  
    197         if (result == NULL){ 
    198             const char * message= "Could not create result "; 
    199         PyErr_SetString(PyExc_RuntimeError , message); 
    200             return NULL; 
    201         } 
    202          
    203     for (int i = 0; i < q->dimensions[0]; i++) { 
    204       double q_value = *(double *)(q->data + i*q->strides[0]); 
    205       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    206       double *result_value = (double *)(result->data + i*result->strides[0]); 
    207       if (q_value == 0) 
    208           *result_value = 0.0; 
    209       else 
    210           *result_value = model->evaluate_rphi(q_value, phi_value); 
    211     } 
    212     return PyArray_Return(result);  
    213  } 
     177 
    214178 /** 
    215179 * Function to call to evaluate model 

  • sansmodels/src/sans/models/c_models/CEllipticalCylinderModel.cpp

    r9bd69098 r975ec8e  
    183183    return PyArray_Return(result);  
    184184 } 
    185 /** 
    186  * Function to call to evaluate model 
    187  * @param args: input numpy array  [q[],phi[]] 
    188  * @return: numpy array object  
    189  */ 
    190 static PyObject * evaluateTwoDim( EllipticalCylinderModel* model,  
    191                               PyArrayObject *q, PyArrayObject *phi) 
    192  { 
    193     PyArrayObject *result; 
    194     //check validity of input vectors 
    195     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    196         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    197         || phi->dimensions[0] != q->dimensions[0]){ 
    198       
    199         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    200         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    201         return NULL; 
    202     } 
    203         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    204  
    205         if (result == NULL){ 
    206             const char * message= "Could not create result "; 
    207         PyErr_SetString(PyExc_RuntimeError , message); 
    208             return NULL; 
    209         } 
    210          
    211     for (int i = 0; i < q->dimensions[0]; i++) { 
    212       double q_value = *(double *)(q->data + i*q->strides[0]); 
    213       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    214       double *result_value = (double *)(result->data + i*result->strides[0]); 
    215       if (q_value == 0) 
    216           *result_value = 0.0; 
    217       else 
    218           *result_value = model->evaluate_rphi(q_value, phi_value); 
    219     } 
    220     return PyArray_Return(result);  
    221  } 
     185 
    222186 /** 
    223187 * Function to call to evaluate model 

  • sansmodels/src/sans/models/c_models/CLamellarModel.cpp

    r9bd69098 r975ec8e  
    8686         
    8787        // Initialize parameter dictionary 
     88        PyDict_SetItemString(self->params,"sld_sol",Py_BuildValue("d",0.000006)); 
    8889        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000)); 
    89         PyDict_SetItemString(self->params,"sigma",Py_BuildValue("d",0.150000)); 
     90        PyDict_SetItemString(self->params,"bi_thick",Py_BuildValue("d",50.000000)); 
    9091        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000)); 
    91         PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000005)); 
    92         PyDict_SetItemString(self->params,"delta",Py_BuildValue("d",50.000000)); 
     92        PyDict_SetItemString(self->params,"sld_bi",Py_BuildValue("d",0.000001)); 
    9393        // Initialize dispersion / averaging parameter dict 
    9494        DispersionVisitor* visitor = new DispersionVisitor(); 
    9595        PyObject * disp_dict; 
     96        disp_dict = PyDict_New(); 
     97        self->model->bi_thick.dispersion->accept_as_source(visitor, self->model->bi_thick.dispersion, disp_dict); 
     98        PyDict_SetItemString(self->dispersion, "bi_thick", disp_dict); 
    9699 
    97100 
    … …  
    161164    return PyArray_Return(result);  
    162165 } 
    163 /** 
    164  * Function to call to evaluate model 
    165  * @param args: input numpy array  [q[],phi[]] 
    166  * @return: numpy array object  
    167  */ 
    168 static PyObject * evaluateTwoDim( LamellarModel* model,  
    169                               PyArrayObject *q, PyArrayObject *phi) 
    170  { 
    171     PyArrayObject *result; 
    172     //check validity of input vectors 
    173     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    174         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    175         || phi->dimensions[0] != q->dimensions[0]){ 
    176       
    177         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    178         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    179         return NULL; 
    180     } 
    181         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    182  
    183         if (result == NULL){ 
    184             const char * message= "Could not create result "; 
    185         PyErr_SetString(PyExc_RuntimeError , message); 
    186             return NULL; 
    187         } 
    188          
    189     for (int i = 0; i < q->dimensions[0]; i++) { 
    190       double q_value = *(double *)(q->data + i*q->strides[0]); 
    191       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    192       double *result_value = (double *)(result->data + i*result->strides[0]); 
    193       if (q_value == 0) 
    194           *result_value = 0.0; 
    195       else 
    196           *result_value = model->evaluate_rphi(q_value, phi_value); 
    197     } 
    198     return PyArray_Return(result);  
    199  } 
     166 
    200167 /** 
    201168 * Function to call to evaluate model 
    … …  
    260227         
    261228            // Reader parameter dictionary 
     229    self->model->sld_sol = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_sol") ); 
    262230    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
    263     self->model->sigma = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma") ); 
     231    self->model->bi_thick = PyFloat_AsDouble( PyDict_GetItemString(self->params, "bi_thick") ); 
    264232    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") ); 
    265     self->model->contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") ); 
    266     self->model->delta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "delta") ); 
     233    self->model->sld_bi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_bi") ); 
    267234    // Read in dispersion parameters 
    268235    PyObject* disp_dict; 
    269236    DispersionVisitor* visitor = new DispersionVisitor(); 
     237    disp_dict = PyDict_GetItemString(self->dispersion, "bi_thick"); 
     238    self->model->bi_thick.dispersion->accept_as_destination(visitor, self->model->bi_thick.dispersion, disp_dict); 
    270239 
    271240         
    … …  
    330299         
    331300            // Reader parameter dictionary 
     301    self->model->sld_sol = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_sol") ); 
    332302    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
    333     self->model->sigma = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma") ); 
     303    self->model->bi_thick = PyFloat_AsDouble( PyDict_GetItemString(self->params, "bi_thick") ); 
    334304    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") ); 
    335     self->model->contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") ); 
    336     self->model->delta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "delta") ); 
     305    self->model->sld_bi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_bi") ); 
    337306    // Read in dispersion parameters 
    338307    PyObject* disp_dict; 
    339308    DispersionVisitor* visitor = new DispersionVisitor(); 
     309    disp_dict = PyDict_GetItemString(self->dispersion, "bi_thick"); 
     310    self->model->bi_thick.dispersion->accept_as_destination(visitor, self->model->bi_thick.dispersion, disp_dict); 
    340311 
    341312         
    … …  
    389360         
    390361            // Reader parameter dictionary 
     362    self->model->sld_sol = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_sol") ); 
    391363    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
    392     self->model->sigma = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma") ); 
     364    self->model->bi_thick = PyFloat_AsDouble( PyDict_GetItemString(self->params, "bi_thick") ); 
    393365    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") ); 
    394     self->model->contrast = PyFloat_AsDouble( PyDict_GetItemString(self->params, "contrast") ); 
    395     self->model->delta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "delta") ); 
     366    self->model->sld_bi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sld_bi") ); 
    396367    // Read in dispersion parameters 
    397368    PyObject* disp_dict; 
    398369    DispersionVisitor* visitor = new DispersionVisitor(); 
     370    disp_dict = PyDict_GetItemString(self->dispersion, "bi_thick"); 
     371    self->model->bi_thick.dispersion->accept_as_destination(visitor, self->model->bi_thick.dispersion, disp_dict); 
    399372 
    400373         
    … …  
    452425        // Ugliness necessary to go from python to C 
    453426            // TODO: refactor this 
    454  { 
     427    if (!strcmp(par_name, "bi_thick")) { 
     428        self->model->bi_thick.dispersion = dispersion; 
     429    } else { 
    455430            PyErr_SetString(CLamellarModelError, 
    456431                "CLamellarModel.set_dispersion expects a valid parameter name."); 

  • sansmodels/src/sans/models/c_models/COblateModel.cpp

    r9bd69098 r975ec8e  
    178178    return PyArray_Return(result);  
    179179 } 
    180 /** 
    181  * Function to call to evaluate model 
    182  * @param args: input numpy array  [q[],phi[]] 
    183  * @return: numpy array object  
    184  */ 
    185 static PyObject * evaluateTwoDim( OblateModel* model,  
    186                               PyArrayObject *q, PyArrayObject *phi) 
    187  { 
    188     PyArrayObject *result; 
    189     //check validity of input vectors 
    190     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    191         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    192         || phi->dimensions[0] != q->dimensions[0]){ 
    193       
    194         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    195         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    196         return NULL; 
    197     } 
    198         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    199  
    200         if (result == NULL){ 
    201             const char * message= "Could not create result "; 
    202         PyErr_SetString(PyExc_RuntimeError , message); 
    203             return NULL; 
    204         } 
    205          
    206     for (int i = 0; i < q->dimensions[0]; i++) { 
    207       double q_value = *(double *)(q->data + i*q->strides[0]); 
    208       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    209       double *result_value = (double *)(result->data + i*result->strides[0]); 
    210       if (q_value == 0) 
    211           *result_value = 0.0; 
    212       else 
    213           *result_value = model->evaluate_rphi(q_value, phi_value); 
    214     } 
    215     return PyArray_Return(result);  
    216  } 
     180 
    217181 /** 
    218182 * Function to call to evaluate model 

  • sansmodels/src/sans/models/c_models/CParallelepipedModel.cpp

    r9bd69098 r975ec8e  
    8888        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000)); 
    8989        PyDict_SetItemString(self->params,"longer_edgeB",Py_BuildValue("d",75.000000)); 
     90        PyDict_SetItemString(self->params,"parallel_psi",Py_BuildValue("d",0.000000)); 
    9091        PyDict_SetItemString(self->params,"longuest_edgeC",Py_BuildValue("d",400.000000)); 
    91         PyDict_SetItemString(self->params,"parallel_phi",Py_BuildValue("d",1.000000)); 
    92         PyDict_SetItemString(self->params,"parallel_theta",Py_BuildValue("d",1.000000)); 
     92        PyDict_SetItemString(self->params,"parallel_phi",Py_BuildValue("d",0.000000)); 
     93        PyDict_SetItemString(self->params,"parallel_theta",Py_BuildValue("d",0.000000)); 
    9394        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000)); 
    9495        PyDict_SetItemString(self->params,"short_edgeA",Py_BuildValue("d",35.000000)); 
    … …  
    109110        self->model->parallel_phi.dispersion->accept_as_source(visitor, self->model->parallel_phi.dispersion, disp_dict); 
    110111        PyDict_SetItemString(self->dispersion, "parallel_phi", disp_dict); 
     112        disp_dict = PyDict_New(); 
     113        self->model->parallel_psi.dispersion->accept_as_source(visitor, self->model->parallel_psi.dispersion, disp_dict); 
     114        PyDict_SetItemString(self->dispersion, "parallel_psi", disp_dict); 
    111115        disp_dict = PyDict_New(); 
    112116        self->model->parallel_theta.dispersion->accept_as_source(visitor, self->model->parallel_theta.dispersion, disp_dict); 
    … …  
    179183    return PyArray_Return(result);  
    180184 } 
    181 /** 
    182  * Function to call to evaluate model 
    183  * @param args: input numpy array  [q[],phi[]] 
    184  * @return: numpy array object  
    185  */ 
    186 static PyObject * evaluateTwoDim( ParallelepipedModel* model,  
    187                               PyArrayObject *q, PyArrayObject *phi) 
    188  { 
    189     PyArrayObject *result; 
    190     //check validity of input vectors 
    191     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    192         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    193         || phi->dimensions[0] != q->dimensions[0]){ 
    194       
    195         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    196         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    197         return NULL; 
    198     } 
    199         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    200  
    201         if (result == NULL){ 
    202             const char * message= "Could not create result "; 
    203         PyErr_SetString(PyExc_RuntimeError , message); 
    204             return NULL; 
    205         } 
    206          
    207     for (int i = 0; i < q->dimensions[0]; i++) { 
    208       double q_value = *(double *)(q->data + i*q->strides[0]); 
    209       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    210       double *result_value = (double *)(result->data + i*result->strides[0]); 
    211       if (q_value == 0) 
    212           *result_value = 0.0; 
    213       else 
    214           *result_value = model->evaluate_rphi(q_value, phi_value); 
    215     } 
    216     return PyArray_Return(result);  
    217  } 
     185 
    218186 /** 
    219187 * Function to call to evaluate model 
    … …  
    280248    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
    281249    self->model->longer_edgeB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longer_edgeB") ); 
     250    self->model->parallel_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_psi") ); 
    282251    self->model->longuest_edgeC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longuest_edgeC") ); 
    283252    self->model->parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") ); 
    … …  
    297266    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_phi"); 
    298267    self->model->parallel_phi.dispersion->accept_as_destination(visitor, self->model->parallel_phi.dispersion, disp_dict); 
     268    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_psi"); 
     269    self->model->parallel_psi.dispersion->accept_as_destination(visitor, self->model->parallel_psi.dispersion, disp_dict); 
    299270    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_theta"); 
    300271    self->model->parallel_theta.dispersion->accept_as_destination(visitor, self->model->parallel_theta.dispersion, disp_dict); 
    … …  
    363334    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
    364335    self->model->longer_edgeB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longer_edgeB") ); 
     336    self->model->parallel_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_psi") ); 
    365337    self->model->longuest_edgeC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longuest_edgeC") ); 
    366338    self->model->parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") ); 
    … …  
    380352    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_phi"); 
    381353    self->model->parallel_phi.dispersion->accept_as_destination(visitor, self->model->parallel_phi.dispersion, disp_dict); 
     354    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_psi"); 
     355    self->model->parallel_psi.dispersion->accept_as_destination(visitor, self->model->parallel_psi.dispersion, disp_dict); 
    382356    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_theta"); 
    383357    self->model->parallel_theta.dispersion->accept_as_destination(visitor, self->model->parallel_theta.dispersion, disp_dict); 
    … …  
    435409    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
    436410    self->model->longer_edgeB = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longer_edgeB") ); 
     411    self->model->parallel_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_psi") ); 
    437412    self->model->longuest_edgeC = PyFloat_AsDouble( PyDict_GetItemString(self->params, "longuest_edgeC") ); 
    438413    self->model->parallel_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "parallel_phi") ); 
    … …  
    452427    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_phi"); 
    453428    self->model->parallel_phi.dispersion->accept_as_destination(visitor, self->model->parallel_phi.dispersion, disp_dict); 
     429    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_psi"); 
     430    self->model->parallel_psi.dispersion->accept_as_destination(visitor, self->model->parallel_psi.dispersion, disp_dict); 
    454431    disp_dict = PyDict_GetItemString(self->dispersion, "parallel_theta"); 
    455432    self->model->parallel_theta.dispersion->accept_as_destination(visitor, self->model->parallel_theta.dispersion, disp_dict); 
    … …  
    517494    } else    if (!strcmp(par_name, "parallel_phi")) { 
    518495        self->model->parallel_phi.dispersion = dispersion; 
     496    } else    if (!strcmp(par_name, "parallel_psi")) { 
     497        self->model->parallel_psi.dispersion = dispersion; 
    519498    } else    if (!strcmp(par_name, "parallel_theta")) { 
    520499        self->model->parallel_theta.dispersion = dispersion; 

  • sansmodels/src/sans/models/c_models/CSphereModel.cpp

    r9bd69098 r975ec8e  
    163163    return PyArray_Return(result);  
    164164 } 
    165 /** 
    166  * Function to call to evaluate model 
    167  * @param args: input numpy array  [q[],phi[]] 
    168  * @return: numpy array object  
    169  */ 
    170 static PyObject * evaluateTwoDim( SphereModel* model,  
    171                               PyArrayObject *q, PyArrayObject *phi) 
    172  { 
    173     PyArrayObject *result; 
    174     //check validity of input vectors 
    175     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    176         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    177         || phi->dimensions[0] != q->dimensions[0]){ 
    178       
    179         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    180         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    181         return NULL; 
    182     } 
    183         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    184  
    185         if (result == NULL){ 
    186             const char * message= "Could not create result "; 
    187         PyErr_SetString(PyExc_RuntimeError , message); 
    188             return NULL; 
    189         } 
    190          
    191     for (int i = 0; i < q->dimensions[0]; i++) { 
    192       double q_value = *(double *)(q->data + i*q->strides[0]); 
    193       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    194       double *result_value = (double *)(result->data + i*result->strides[0]); 
    195       if (q_value == 0) 
    196           *result_value = 0.0; 
    197       else 
    198           *result_value = model->evaluate_rphi(q_value, phi_value); 
    199     } 
    200     return PyArray_Return(result);  
    201  } 
     165 
    202166 /** 
    203167 * Function to call to evaluate model 

  • sansmodels/src/sans/models/c_models/CStackedDisksModel.cpp

    r9bd69098 r975ec8e  
    8686         
    8787        // Initialize parameter dictionary 
     88        PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.000004)); 
     89        PyDict_SetItemString(self->params,"core_thick",Py_BuildValue("d",10.000000)); 
     90        PyDict_SetItemString(self->params,"layer_thick",Py_BuildValue("d",15.000000)); 
     91        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",0.000000)); 
     92        PyDict_SetItemString(self->params,"layer_sld",Py_BuildValue("d",-0.000000)); 
     93        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.000000)); 
     94        PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.000005)); 
    8895        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",0.010000)); 
    89         PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.000000)); 
    90         PyDict_SetItemString(self->params,"spacing",Py_BuildValue("d",0.000000)); 
    91         PyDict_SetItemString(self->params,"nlayers",Py_BuildValue("d",1.000000)); 
    92         PyDict_SetItemString(self->params,"layer_sld",Py_BuildValue("d",-0.000000)); 
    93         PyDict_SetItemString(self->params,"solvent_sld",Py_BuildValue("d",0.000005)); 
    94         PyDict_SetItemString(self->params,"thickness",Py_BuildValue("d",15.000000)); 
    95         PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.000000)); 
    96         PyDict_SetItemString(self->params,"length",Py_BuildValue("d",10.000000)); 
    97         PyDict_SetItemString(self->params,"core_sld",Py_BuildValue("d",0.000004)); 
    9896        PyDict_SetItemString(self->params,"radius",Py_BuildValue("d",3000.000000)); 
    9997        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.001000)); 
     98        PyDict_SetItemString(self->params,"sigma_d",Py_BuildValue("d",0.000000)); 
     99        PyDict_SetItemString(self->params,"n_stacking",Py_BuildValue("d",1.000000)); 
    100100        // Initialize dispersion / averaging parameter dict 
    101101        DispersionVisitor* visitor = new DispersionVisitor(); 
    102102        PyObject * disp_dict; 
    103103        disp_dict = PyDict_New(); 
    104         self->model->length.dispersion->accept_as_source(visitor, self->model->length.dispersion, disp_dict); 
    105         PyDict_SetItemString(self->dispersion, "length", disp_dict); 
     104        self->model->core_thick.dispersion->accept_as_source(visitor, self->model->core_thick.dispersion, disp_dict); 
     105        PyDict_SetItemString(self->dispersion, "core_thick", disp_dict); 
     106        disp_dict = PyDict_New(); 
     107        self->model->layer_thick.dispersion->accept_as_source(visitor, self->model->layer_thick.dispersion, disp_dict); 
     108        PyDict_SetItemString(self->dispersion, "layer_thick", disp_dict); 
    106109        disp_dict = PyDict_New(); 
    107110        self->model->radius.dispersion->accept_as_source(visitor, self->model->radius.dispersion, disp_dict); 
    … …  
    180183    return PyArray_Return(result);  
    181184 } 
    182 /** 
    183  * Function to call to evaluate model 
    184  * @param args: input numpy array  [q[],phi[]] 
    185  * @return: numpy array object  
    186  */ 
    187 static PyObject * evaluateTwoDim( StackedDisksModel* model,  
    188                               PyArrayObject *q, PyArrayObject *phi) 
    189  { 
    190     PyArrayObject *result; 
    191     //check validity of input vectors 
    192     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    193         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    194         || phi->dimensions[0] != q->dimensions[0]){ 
    195       
    196         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    197         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    198         return NULL; 
    199     } 
    200         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    201  
    202         if (result == NULL){ 
    203             const char * message= "Could not create result "; 
    204         PyErr_SetString(PyExc_RuntimeError , message); 
    205             return NULL; 
    206         } 
    207          
    208     for (int i = 0; i < q->dimensions[0]; i++) { 
    209       double q_value = *(double *)(q->data + i*q->strides[0]); 
    210       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    211       double *result_value = (double *)(result->data + i*result->strides[0]); 
    212       if (q_value == 0) 
    213           *result_value = 0.0; 
    214       else 
    215           *result_value = model->evaluate_rphi(q_value, phi_value); 
    216     } 
    217     return PyArray_Return(result);  
    218  } 
     185 
    219186 /** 
    220187 * Function to call to evaluate model 
    … …  
    279246         
    280247            // Reader parameter dictionary 
     248    self->model->core_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "core_sld") ); 
     249    self->model->core_thick = PyFloat_AsDouble( PyDict_GetItemString(self->params, "core_thick") ); 
     250    self->model->layer_thick = PyFloat_AsDouble( PyDict_GetItemString(self->params, "layer_thick") ); 
     251    self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") ); 
     252    self->model->layer_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "layer_sld") ); 
     253    self->model->axis_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_phi") ); 
     254    self->model->solvent_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "solvent_sld") ); 
    281255    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
    282     self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") ); 
    283     self->model->spacing = PyFloat_AsDouble( PyDict_GetItemString(self->params, "spacing") ); 
    284     self->model->nlayers = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nlayers") ); 
    285     self->model->layer_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "layer_sld") ); 
    286     self->model->solvent_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "solvent_sld") ); 
    287     self->model->thickness = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thickness") ); 
    288     self->model->axis_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_phi") ); 
    289     self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") ); 
    290     self->model->core_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "core_sld") ); 
    291256    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") ); 
    292257    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") ); 
     258    self->model->sigma_d = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma_d") ); 
     259    self->model->n_stacking = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_stacking") ); 
    293260    // Read in dispersion parameters 
    294261    PyObject* disp_dict; 
    295262    DispersionVisitor* visitor = new DispersionVisitor(); 
    296     disp_dict = PyDict_GetItemString(self->dispersion, "length"); 
    297     self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict); 
     263    disp_dict = PyDict_GetItemString(self->dispersion, "core_thick"); 
     264    self->model->core_thick.dispersion->accept_as_destination(visitor, self->model->core_thick.dispersion, disp_dict); 
     265    disp_dict = PyDict_GetItemString(self->dispersion, "layer_thick"); 
     266    self->model->layer_thick.dispersion->accept_as_destination(visitor, self->model->layer_thick.dispersion, disp_dict); 
    298267    disp_dict = PyDict_GetItemString(self->dispersion, "radius"); 
    299268    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict); 
    … …  
    364333         
    365334            // Reader parameter dictionary 
     335    self->model->core_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "core_sld") ); 
     336    self->model->core_thick = PyFloat_AsDouble( PyDict_GetItemString(self->params, "core_thick") ); 
     337    self->model->layer_thick = PyFloat_AsDouble( PyDict_GetItemString(self->params, "layer_thick") ); 
     338    self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") ); 
     339    self->model->layer_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "layer_sld") ); 
     340    self->model->axis_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_phi") ); 
     341    self->model->solvent_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "solvent_sld") ); 
    366342    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
    367     self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") ); 
    368     self->model->spacing = PyFloat_AsDouble( PyDict_GetItemString(self->params, "spacing") ); 
    369     self->model->nlayers = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nlayers") ); 
    370     self->model->layer_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "layer_sld") ); 
    371     self->model->solvent_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "solvent_sld") ); 
    372     self->model->thickness = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thickness") ); 
    373     self->model->axis_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_phi") ); 
    374     self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") ); 
    375     self->model->core_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "core_sld") ); 
    376343    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") ); 
    377344    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") ); 
     345    self->model->sigma_d = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma_d") ); 
     346    self->model->n_stacking = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_stacking") ); 
    378347    // Read in dispersion parameters 
    379348    PyObject* disp_dict; 
    380349    DispersionVisitor* visitor = new DispersionVisitor(); 
    381     disp_dict = PyDict_GetItemString(self->dispersion, "length"); 
    382     self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict); 
     350    disp_dict = PyDict_GetItemString(self->dispersion, "core_thick"); 
     351    self->model->core_thick.dispersion->accept_as_destination(visitor, self->model->core_thick.dispersion, disp_dict); 
     352    disp_dict = PyDict_GetItemString(self->dispersion, "layer_thick"); 
     353    self->model->layer_thick.dispersion->accept_as_destination(visitor, self->model->layer_thick.dispersion, disp_dict); 
    383354    disp_dict = PyDict_GetItemString(self->dispersion, "radius"); 
    384355    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict); 
    … …  
    438409         
    439410            // Reader parameter dictionary 
     411    self->model->core_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "core_sld") ); 
     412    self->model->core_thick = PyFloat_AsDouble( PyDict_GetItemString(self->params, "core_thick") ); 
     413    self->model->layer_thick = PyFloat_AsDouble( PyDict_GetItemString(self->params, "layer_thick") ); 
     414    self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") ); 
     415    self->model->layer_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "layer_sld") ); 
     416    self->model->axis_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_phi") ); 
     417    self->model->solvent_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "solvent_sld") ); 
    440418    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
    441     self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") ); 
    442     self->model->spacing = PyFloat_AsDouble( PyDict_GetItemString(self->params, "spacing") ); 
    443     self->model->nlayers = PyFloat_AsDouble( PyDict_GetItemString(self->params, "nlayers") ); 
    444     self->model->layer_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "layer_sld") ); 
    445     self->model->solvent_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "solvent_sld") ); 
    446     self->model->thickness = PyFloat_AsDouble( PyDict_GetItemString(self->params, "thickness") ); 
    447     self->model->axis_phi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_phi") ); 
    448     self->model->length = PyFloat_AsDouble( PyDict_GetItemString(self->params, "length") ); 
    449     self->model->core_sld = PyFloat_AsDouble( PyDict_GetItemString(self->params, "core_sld") ); 
    450419    self->model->radius = PyFloat_AsDouble( PyDict_GetItemString(self->params, "radius") ); 
    451420    self->model->background = PyFloat_AsDouble( PyDict_GetItemString(self->params, "background") ); 
     421    self->model->sigma_d = PyFloat_AsDouble( PyDict_GetItemString(self->params, "sigma_d") ); 
     422    self->model->n_stacking = PyFloat_AsDouble( PyDict_GetItemString(self->params, "n_stacking") ); 
    452423    // Read in dispersion parameters 
    453424    PyObject* disp_dict; 
    454425    DispersionVisitor* visitor = new DispersionVisitor(); 
    455     disp_dict = PyDict_GetItemString(self->dispersion, "length"); 
    456     self->model->length.dispersion->accept_as_destination(visitor, self->model->length.dispersion, disp_dict); 
     426    disp_dict = PyDict_GetItemString(self->dispersion, "core_thick"); 
     427    self->model->core_thick.dispersion->accept_as_destination(visitor, self->model->core_thick.dispersion, disp_dict); 
     428    disp_dict = PyDict_GetItemString(self->dispersion, "layer_thick"); 
     429    self->model->layer_thick.dispersion->accept_as_destination(visitor, self->model->layer_thick.dispersion, disp_dict); 
    457430    disp_dict = PyDict_GetItemString(self->dispersion, "radius"); 
    458431    self->model->radius.dispersion->accept_as_destination(visitor, self->model->radius.dispersion, disp_dict); 
    … …  
    516489        // Ugliness necessary to go from python to C 
    517490            // TODO: refactor this 
    518     if (!strcmp(par_name, "length")) { 
    519         self->model->length.dispersion = dispersion; 
     491    if (!strcmp(par_name, "core_thick")) { 
     492        self->model->core_thick.dispersion = dispersion; 
     493    } else    if (!strcmp(par_name, "layer_thick")) { 
     494        self->model->layer_thick.dispersion = dispersion; 
    520495    } else    if (!strcmp(par_name, "radius")) { 
    521496        self->model->radius.dispersion = dispersion; 

  • sansmodels/src/sans/models/c_models/CTriaxialEllipsoidModel.cpp

    r9bd69098 r975ec8e  
    8787        // Initialize parameter dictionary 
    8888        PyDict_SetItemString(self->params,"scale",Py_BuildValue("d",1.000000)); 
    89         PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",1.000000)); 
     89        PyDict_SetItemString(self->params,"axis_psi",Py_BuildValue("d",0.000000)); 
     90        PyDict_SetItemString(self->params,"axis_theta",Py_BuildValue("d",0.000000)); 
    9091        PyDict_SetItemString(self->params,"semi_axisA",Py_BuildValue("d",100.000000)); 
    9192        PyDict_SetItemString(self->params,"semi_axisB",Py_BuildValue("d",35.000000)); 
    9293        PyDict_SetItemString(self->params,"semi_axisC",Py_BuildValue("d",400.000000)); 
    93         PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",1.000000)); 
     94        PyDict_SetItemString(self->params,"axis_phi",Py_BuildValue("d",0.000000)); 
    9495        PyDict_SetItemString(self->params,"background",Py_BuildValue("d",0.000000)); 
    9596        PyDict_SetItemString(self->params,"contrast",Py_BuildValue("d",0.000005)); 
    … …  
    103104        self->model->axis_phi.dispersion->accept_as_source(visitor, self->model->axis_phi.dispersion, disp_dict); 
    104105        PyDict_SetItemString(self->dispersion, "axis_phi", disp_dict); 
     106        disp_dict = PyDict_New(); 
     107        self->model->axis_psi.dispersion->accept_as_source(visitor, self->model->axis_psi.dispersion, disp_dict); 
     108        PyDict_SetItemString(self->dispersion, "axis_psi", disp_dict); 
    105109 
    106110 
    … …  
    170174    return PyArray_Return(result);  
    171175 } 
    172 /** 
    173  * Function to call to evaluate model 
    174  * @param args: input numpy array  [q[],phi[]] 
    175  * @return: numpy array object  
    176  */ 
    177 static PyObject * evaluateTwoDim( TriaxialEllipsoidModel* model,  
    178                               PyArrayObject *q, PyArrayObject *phi) 
    179  { 
    180     PyArrayObject *result; 
    181     //check validity of input vectors 
    182     if (q->nd != 1 || q->descr->type_num != PyArray_DOUBLE 
    183         || phi->nd != 1 || phi->descr->type_num != PyArray_DOUBLE 
    184         || phi->dimensions[0] != q->dimensions[0]){ 
    185       
    186         //const char * message= "Invalid array: q->nd=%d,type_num=%d\n",q->nd,q->descr->type_num; 
    187         PyErr_SetString(PyExc_ValueError ,"wrong input");  
    188         return NULL; 
    189     } 
    190         result= (PyArrayObject *)PyArray_FromDims(q->nd,(int*)(q->dimensions), PyArray_DOUBLE); 
    191  
    192         if (result == NULL){ 
    193             const char * message= "Could not create result "; 
    194         PyErr_SetString(PyExc_RuntimeError , message); 
    195             return NULL; 
    196         } 
    197          
    198     for (int i = 0; i < q->dimensions[0]; i++) { 
    199       double q_value = *(double *)(q->data + i*q->strides[0]); 
    200       double phi_value = *(double *)(phi->data + i*phi->strides[0]); 
    201       double *result_value = (double *)(result->data + i*result->strides[0]); 
    202       if (q_value == 0) 
    203           *result_value = 0.0; 
    204       else 
    205           *result_value = model->evaluate_rphi(q_value, phi_value); 
    206     } 
    207     return PyArray_Return(result);  
    208  } 
     176 
    209177 /** 
    210178 * Function to call to evaluate model 
    … …  
    270238            // Reader parameter dictionary 
    271239    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
     240    self->model->axis_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_psi") ); 
    272241    self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") ); 
    273242    self->model->semi_axisA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "semi_axisA") ); 
    … …  
    284253    disp_dict = PyDict_GetItemString(self->dispersion, "axis_phi"); 
    285254    self->model->axis_phi.dispersion->accept_as_destination(visitor, self->model->axis_phi.dispersion, disp_dict); 
     255    disp_dict = PyDict_GetItemString(self->dispersion, "axis_psi"); 
     256    self->model->axis_psi.dispersion->accept_as_destination(visitor, self->model->axis_psi.dispersion, disp_dict); 
    286257 
    287258         
    … …  
    347318            // Reader parameter dictionary 
    348319    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
     320    self->model->axis_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_psi") ); 
    349321    self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") ); 
    350322    self->model->semi_axisA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "semi_axisA") ); 
    … …  
    361333    disp_dict = PyDict_GetItemString(self->dispersion, "axis_phi"); 
    362334    self->model->axis_phi.dispersion->accept_as_destination(visitor, self->model->axis_phi.dispersion, disp_dict); 
     335    disp_dict = PyDict_GetItemString(self->dispersion, "axis_psi"); 
     336    self->model->axis_psi.dispersion->accept_as_destination(visitor, self->model->axis_psi.dispersion, disp_dict); 
    363337 
    364338         
    … …  
    413387            // Reader parameter dictionary 
    414388    self->model->scale = PyFloat_AsDouble( PyDict_GetItemString(self->params, "scale") ); 
     389    self->model->axis_psi = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_psi") ); 
    415390    self->model->axis_theta = PyFloat_AsDouble( PyDict_GetItemString(self->params, "axis_theta") ); 
    416391    self->model->semi_axisA = PyFloat_AsDouble( PyDict_GetItemString(self->params, "semi_axisA") ); 
    … …  
    427402    disp_dict = PyDict_GetItemString(self->dispersion, "axis_phi"); 
    428403    self->model->axis_phi.dispersion->accept_as_destination(visitor, self->model->axis_phi.dispersion, disp_dict); 
     404    disp_dict = PyDict_GetItemString(self->dispersion, "axis_psi"); 
     405    self->model->axis_psi.dispersion->accept_as_destination(visitor, self->model->axis_psi.dispersion, disp_dict); 
    429406 
    430407         
    … …  
    486463    } else    if (!strcmp(par_name, "axis_phi")) { 
    487464        self->model->axis_phi.dispersion = dispersion; 
     465    } else    if (!strcmp(par_name, "axis_psi")) { 
     466        self->model->axis_psi.dispersion = dispersion; 
    488467    } else { 
    489468            PyErr_SetString(CTriaxialEllipsoidModelError, 

  • sansmodels/src/sans/models/c_models/c_models.cpp

    reba9885 r975ec8e  
    3030void addCLamellarPSModel(PyObject *module); 
    3131void addCLamellarPSHGModel(PyObject *module); 
     32void addCCoreShellSpheroidModel(PyObject *module); 
    3233void addCOblateModel(PyObject *module); 
    3334void addCProlateModel(PyObject *module); 
    … …  
    225226        addCLamellarPSModel(m); 
    226227        addCLamellarPSHGModel(m); 
     228        addCCoreShellSpheroidModel(m); 
    227229        addCOblateModel(m); 
    228230        addCProlateModel(m); 

  • sansmodels/src/sans/models/c_models/lamellar.cpp

    r42f193a r975ec8e  
    3434LamellarModel :: LamellarModel() { 
    3535        scale      = Parameter(1.0); 
    36         delta     = Parameter(50.0); 
    37         delta.set_min(0.0); 
    38         sigma    = Parameter(0.15); 
    39         sigma.set_min(0.0); 
    40         contrast   = Parameter(5.3e-6); 
     36        bi_thick     = Parameter(50.0, true); 
     37        bi_thick.set_min(0.0); 
     38        sld_bi    = Parameter(1.0e-6); 
     39        sld_sol    = Parameter(6.3e-6); 
    4140        background = Parameter(0.0); 
    4241 
    … …  
    5554        // Add the background after averaging 
    5655        dp[0] = scale(); 
    57         dp[1] = delta(); 
    58         dp[2] = sigma(); 
    59         dp[3] = contrast(); 
    60         dp[4] = background(); 
     56        dp[1] = bi_thick(); 
     57        dp[2] = sld_bi(); 
     58        dp[3] = sld_sol(); 
     59        dp[4] = 0.0; 
    6160 
    62         return LamellarFF(dp, q); 
     61        // Get the dispersion points for the bi_thick 
     62        vector<WeightPoint> weights_bi_thick; 
     63        bi_thick.get_weights(weights_bi_thick); 
     64        // Perform the computation, with all weight points 
     65        double sum = 0.0; 
     66        double norm = 0.0; 
     67 
     68        // Loop over short_edgeA weight points 
     69        for(int i=0; i< (int)weights_bi_thick.size(); i++) { 
     70                dp[1] = weights_bi_thick[i].value; 
     71 
     72                sum += weights_bi_thick[i].weight * lamellar_kernel(dp, q); 
     73                norm += weights_bi_thick[i].weight; 
     74 
     75        } 
     76 
     77        return sum/norm + background(); 
    6378} 
    6479 

  • sansmodels/src/sans/models/c_models/models.hh

    r2c4b289 r975ec8e  
    4646class ParallelepipedModel{ 
    4747public: 
    48         // TODO: add 2D  
     48        // TODO: add 2D 
    4949        // Model parameters 
    5050        Parameter scale; 
    … …  
    5656        Parameter parallel_theta; 
    5757        Parameter parallel_phi; 
     58        Parameter parallel_psi; 
    5859 
    5960        // Constructor 
    … …  
    277278        Parameter axis_theta; 
    278279        Parameter axis_phi; 
    279          
     280        Parameter axis_psi; 
     281 
    280282        // Constructor 
    281283        TriaxialEllipsoidModel(); 
    … …  
    298300        Parameter axis_theta; 
    299301        Parameter axis_phi; 
    300          
     302 
    301303        // Constructor 
    302304        FlexibleCylinderModel(); 
    … …  
    312314        // Model parameters 
    313315        Parameter scale; 
    314         Parameter length; 
    315         Parameter radius; 
    316         Parameter thickness; 
     316        Parameter core_thick; 
     317        Parameter radius; 
     318        Parameter layer_thick; 
    317319        Parameter core_sld; 
    318320        Parameter layer_sld; 
    319321        Parameter solvent_sld; 
    320         Parameter nlayers; 
    321         Parameter spacing; 
    322         Parameter background; 
    323         Parameter axis_theta; 
    324         Parameter axis_phi; 
    325          
     322        Parameter n_stacking; 
     323        Parameter sigma_d; 
     324        Parameter background; 
     325        Parameter axis_theta; 
     326        Parameter axis_phi; 
     327 
    326328        // Constructor 
    327329        StackedDisksModel(); 
    … …  
    337339        // Model parameters 
    338340        Parameter scale; 
    339         Parameter delta; 
    340         Parameter sigma; 
    341         Parameter contrast; 
    342         Parameter background; 
    343          
     341        Parameter bi_thick; 
     342        Parameter sld_bi; 
     343        Parameter sld_sol; 
     344        Parameter background; 
     345 
    344346        // Constructor 
    345347        LamellarModel(); 
    … …  
    362364        Parameter sld_solvent; 
    363365        Parameter background; 
    364          
     366 
    365367        // Constructor 
    366368        LamellarFFHGModel(); 
    … …  
    386388        Parameter caille; 
    387389        Parameter background; 
    388          
     390 
    389391        // Constructor 
    390392        LamellarPSModel(); 
    … …  
    409411        Parameter caille; 
    410412        Parameter background; 
    411          
     413 
    412414        // Constructor 
    413415        LamellarPSHGModel(); 
     416 
     417        // Operators to get I(Q) 
     418        double operator()(double q); 
     419        double operator()(double qx, double qy); 
     420        double evaluate_rphi(double q, double phi); 
     421}; 
     422 
     423class CoreShellSpheroidModel{ 
     424public: 
     425        // Model parameters 
     426        Parameter scale; 
     427        Parameter equat_core; 
     428        Parameter polar_core; 
     429        Parameter equat_shell; 
     430        Parameter polar_shell; 
     431        Parameter contrast; 
     432        Parameter sld_solvent; 
     433        Parameter background; 
     434        Parameter axis_theta; 
     435        Parameter axis_phi; 
     436 
     437        // Constructor 
     438        CoreShellSpheroidModel(); 
    414439 
    415440        // Operators to get I(Q) 
    … …  
    432457        Parameter axis_theta; 
    433458        Parameter axis_phi; 
    434          
     459 
    435460        // Constructor 
    436461        OblateModel(); 
    … …  
    454479        Parameter axis_theta; 
    455480        Parameter axis_phi; 
    456          
     481 
    457482        // Constructor 
    458483        ProlateModel(); 
    … …  
    474499        Parameter axis_theta; 
    475500        Parameter axis_phi; 
    476          
     501 
    477502        //Constructor 
    478503        HollowCylinderModel(); 
    479          
     504 
    480505        //Operators to get I(Q) 
    481506        double operator()(double q); 
    … …  
    495520        Parameter n_pairs; 
    496521        Parameter background; 
    497          
     522 
    498523        //Constructor 
    499524        MultiShellModel(); 
    500          
     525 
    501526        //Operators to get I(Q) 
    502527        double operator()(double q); 
    … …  
    514539        Parameter shell_sld; 
    515540        Parameter background; 
    516          
     541 
    517542        //Constructor 
    518543        VesicleModel(); 
    519          
     544 
    520545        //Operators to get I(Q) 
    521546        double operator()(double q); 
    … …  
    535560        Parameter solvent_sld; 
    536561        Parameter background; 
    537          
     562 
    538563        //Constructor 
    539564        BinaryHSModel(); 
    540          
     565 
    541566        //Operators to get I(Q) 
    542567        double operator()(double q); 
    … …  
    556581        Parameter solvent_sld; 
    557582        Parameter background; 
    558          
     583 
    559584        //Constructor 
    560585        BinaryHSPSF11Model(); 
    561          
     586 
    562587        //Operators to get I(Q) 
    563588        double operator()(double q); 

  • sansmodels/src/sans/models/c_models/oblate.cpp

    r96b59384 r975ec8e  
    121121 * @return: function value 
    122122 */ 
    123  
     123/* 
    124124double OblateModel :: operator()(double qx, double qy) { 
    125125        double q = sqrt(qx*qx + qy*qy); 
    … …  
    127127        return (*this).operator()(q); 
    128128} 
    129  
     129*/ 
    130130 
    131131/** 
    … …  
    140140} 
    141141 
    142 /* disable below for now 
    143 double OblateShellModel :: operator()(double qx, double qy) { 
     142 
     143double OblateModel :: operator()(double qx, double qy) { 
    144144        OblateParameters dp; 
    145145        // Fill parameter array 
    … …  
    233233        return sum/norm + background(); 
    234234} 
    235 */// 
     235 

  • sansmodels/src/sans/models/c_models/parallelepiped.cpp

    r8a48713 r975ec8e  
    4545        parallel_theta  = Parameter(0.0, true); 
    4646        parallel_phi    = Parameter(0.0, true); 
     47        parallel_psi    = Parameter(0.0, true); 
    4748} 
    4849 
    … …  
    5455 */ 
    5556double ParallelepipedModel :: operator()(double q) { 
    56         double dp[5]; 
     57        double dp[6]; 
    5758 
    5859        // Fill parameter array for IGOR library 
    … …  
    6364        dp[3] = longuest_edgeC(); 
    6465        dp[4] = contrast(); 
    65         //dp[5] = background(); 
    6666        dp[5] = 0.0; 
    6767 
    … …  
    6969        vector<WeightPoint> weights_short_edgeA; 
    7070        short_edgeA.get_weights(weights_short_edgeA); 
    71          
     71 
    7272        // Get the dispersion points for the longer_edgeB 
    7373        vector<WeightPoint> weights_longer_edgeB; 
    … …  
    8383        double sum = 0.0; 
    8484        double norm = 0.0; 
    85          
     85 
    8686        // Loop over short_edgeA weight points 
    8787        for(int i=0; i< (int)weights_short_edgeA.size(); i++) { 
    … …  
    9696                                dp[3] = weights_longuest_edgeC[j].value; 
    9797 
    98                                 sum += weights_short_edgeA[i].weight * weights_longer_edgeB[j].weight  
     98                                sum += weights_short_edgeA[i].weight * weights_longer_edgeB[j].weight 
    9999                                        * weights_longuest_edgeC[k].weight * Parallelepiped(dp, q); 
    100100 
    … …  
    124124        dp.parallel_theta  = parallel_theta(); 
    125125        dp.parallel_phi    = parallel_phi(); 
    126          
     126        dp.parallel_psi    = parallel_psi(); 
     127 
    127128 
    128129        // Get the dispersion points for the short_edgeA 
    … …  
    146147        parallel_phi.get_weights(weights_parallel_phi); 
    147148 
     149        // Get angular averaging for psi 
     150        vector<WeightPoint> weights_parallel_psi; 
     151        parallel_psi.get_weights(weights_parallel_psi); 
    148152 
    149153        // Perform the computation, with all weight points 
    … …  
    171175                                                dp.parallel_phi = weights_parallel_phi[m].value; 
    172176 
    173                                                 double _ptvalue = weights_short_edgeA[i].weight 
    174                                                         * weights_longer_edgeB[j].weight 
    175                                                         * weights_longuest_edgeC[k].weight 
    176                                                         * weights_parallel_theta[l].weight 
    177                                                         * weights_parallel_phi[m].weight 
    178                                                         * parallelepiped_analytical_2DXY(&dp, qx, qy); 
    179                                                 if (weights_parallel_theta.size()>1) { 
    180                                                         _ptvalue *= sin(weights_parallel_theta[l].value); 
     177                                                // Average over phi distribution 
     178                                                for(int n=0; n< (int)weights_parallel_psi.size(); n++) { 
     179                                                        dp.parallel_psi = weights_parallel_psi[n].value; 
     180 
     181                                                        double _ptvalue = weights_short_edgeA[i].weight 
     182                                                                * weights_longer_edgeB[j].weight 
     183                                                                * weights_longuest_edgeC[k].weight 
     184                                                                * weights_parallel_theta[l].weight 
     185                                                                * weights_parallel_phi[m].weight 
     186                                                                * weights_parallel_psi[n].weight 
     187                                                                * parallelepiped_analytical_2DXY(&dp, qx, qy); 
     188                                                        if (weights_parallel_theta.size()>1) { 
     189                                                                _ptvalue *= sin(weights_parallel_theta[l].value); 
     190                                                        } 
     191                                                        sum += _ptvalue; 
     192 
     193                                                        norm += weights_short_edgeA[i].weight 
     194                                                                * weights_longer_edgeB[j].weight 
     195                                                                * weights_longuest_edgeC[k].weight 
     196                                                                * weights_parallel_theta[l].weight 
     197                                                                * weights_parallel_phi[m].weight 
     198                                                                * weights_parallel_psi[n].weight; 
    181199                                                } 
    182                                                 sum += _ptvalue; 
    183  
    184                                                 norm += weights_short_edgeA[i].weight 
    185                                                         * weights_longer_edgeB[j].weight 
    186                                                         * weights_longuest_edgeC[k].weight 
    187                                                         * weights_parallel_theta[l].weight 
    188                                                         * weights_parallel_phi[m].weight; 
    189200                                        } 
    190201 

  • sansmodels/src/sans/models/c_models/stackeddisks.cpp

    r9188cc1 r975ec8e  
    3737        radius     = Parameter(3000.0, true); 
    3838        radius.set_min(0.0); 
    39         length     = Parameter(10.0, true); 
    40         length.set_min(0.0); 
    41         thickness     = Parameter(15.0); 
    42         thickness.set_min(0.0); 
     39        core_thick  = Parameter(10.0, true); 
     40        core_thick.set_min(0.0); 
     41        layer_thick     = Parameter(15.0); 
     42        layer_thick.set_min(0.0); 
    4343        core_sld = Parameter(4.0e-6); 
    4444        layer_sld  = Parameter(-4.0e-7); 
    4545        solvent_sld  = Parameter(5.0e-6); 
    46         nlayers   = Parameter(1); 
    47         spacing   = Parameter(0); 
     46        n_stacking   = Parameter(1); 
     47        sigma_d   = Parameter(0); 
    4848        background = Parameter(0.001); 
    4949        axis_theta  = Parameter(0.0, true); 
    … …  
    6464        dp[0] = scale(); 
    6565        dp[1] = radius(); 
    66         dp[2] = length(); 
    67         dp[3] = thickness(); 
     66        dp[2] = core_thick(); 
     67        dp[3] = layer_thick(); 
    6868        dp[4] = core_sld(); 
    6969        dp[5] = layer_sld(); 
    7070        dp[6] = solvent_sld(); 
    71         dp[7] = nlayers(); 
    72         dp[8] = spacing(); 
     71        dp[7] = n_stacking(); 
     72        dp[8] = sigma_d(); 
    7373        dp[9] = 0.0; 
    74  
    75         // Get the dispersion points for the length 
    76         vector<WeightPoint> weights_length; 
    77         length.get_weights(weights_length); 
    7874 
    7975        // Get the dispersion points for the radius 
    … …  
    8177        radius.get_weights(weights_radius); 
    8278 
    83         // Get the dispersion points for the thickness 
    84         vector<WeightPoint> weights_thickness; 
    85         thickness.get_weights(weights_thickness); 
     79        // Get the dispersion points for the core_thick 
     80        vector<WeightPoint> weights_core_thick; 
     81        core_thick.get_weights(weights_core_thick); 
     82 
     83        // Get the dispersion points for the layer_thick 
     84        vector<WeightPoint> weights_layer_thick; 
     85        layer_thick.get_weights(weights_layer_thick); 
    8686 
    8787        // Perform the computation, with all weight points 
    … …  
    9090 
    9191        // Loop over length weight points 
    92         for(int i=0; i< (int)weights_length.size(); i++) { 
    93                 dp[1] = weights_length[i].value; 
     92        for(int i=0; i< (int)weights_radius.size(); i++) { 
     93                dp[1] = weights_radius[i].value; 
    9494 
    9595                // Loop over radius weight points 
    96                 for(int j=0; j< (int)weights_radius.size(); j++) { 
    97                         dp[2] = weights_radius[j].value; 
     96                for(int j=0; j< (int)weights_core_thick.size(); j++) { 
     97                        dp[2] = weights_core_thick[j].value; 
    9898 
    9999                        // Loop over thickness weight points 
    100                         for(int k=0; k< (int)weights_radius.size(); k++) { 
    101                                 dp[3] = weights_radius[k].value; 
    102  
    103                                 sum += weights_length[i].weight 
    104                                         * weights_radius[j].weight * weights_thickness[k].weight* StackedDiscs(dp, q); 
    105                                 norm += weights_length[i].weight 
    106                                         * weights_radius[j].weight* weights_thickness[k].weight; 
     100                        for(int k=0; k< (int)weights_layer_thick.size(); k++) { 
     101                                dp[3] = weights_layer_thick[k].value; 
     102 
     103                                sum += weights_radius[i].weight 
     104                                        * weights_core_thick[j].weight * weights_layer_thick[k].weight* StackedDiscs(dp, q); 
     105                                norm += weights_radius[i].weight 
     106                                        * weights_core_thick[j].weight* weights_layer_thick[k].weight; 
    107107                        } 
    108108                } 
    … …  
    121121        // Fill parameter array 
    122122        dp.scale      = scale(); 
    123         dp.length     = length(); 
     123        dp.core_thick    = core_thick(); 
    124124        dp.radius         = radius(); 
    125         dp.thickness  = thickness(); 
     125        dp.core_thick  = core_thick(); 
    126126        dp.core_sld   = core_sld(); 
    127127        dp.layer_sld  = layer_sld(); 
    128128        dp.solvent_sld= solvent_sld(); 
    129         dp.nlayers        = nlayers(); 
    130         dp.spacing    = spacing(); 
     129        dp.n_stacking     = n_stacking(); 
     130        dp.sigma_d   = sigma_d(); 
    131131        dp.background = 0.0; 
    132132        dp.axis_theta = axis_theta(); 
    … …  
    134134 
    135135        // Get the dispersion points for the length 
    136         vector<WeightPoint> weights_length; 
    137         length.get_weights(weights_length); 
     136        vector<WeightPoint> weights_core_thick; 
     137        core_thick.get_weights(weights_core_thick); 
    138138 
    139139        // Get the dispersion points for the radius 
    … …  
    142142 
    143143        // Get the dispersion points for the thickness 
    144         vector<WeightPoint> weights_thickness; 
    145         thickness.get_weights(weights_thickness); 
     144        vector<WeightPoint> weights_layer_thick; 
     145        layer_thick.get_weights(weights_layer_thick); 
    146146 
    147147        // Get angular averaging for theta 
    … …  
    158158 
    159159        // Loop over length weight points 
    160         for(int i=0; i< (int)weights_length.size(); i++) { 
    161                 dp.length = weights_length[i].value; 
     160        for(int i=0; i< (int)weights_core_thick.size(); i++) { 
     161                dp.core_thick = weights_core_thick[i].value; 
    162162 
    163163                // Loop over radius weight points 
    … …  
    166166 
    167167                                // Loop over thickness weight points 
    168                                 for(int k=0; k< (int)weights_thickness.size(); k++) { 
    169                                 dp.thickness = weights_thickness[k].value; 
     168                                for(int k=0; k< (int)weights_layer_thick.size(); k++) { 
     169                                dp.layer_thick = weights_layer_thick[k].value; 
    170170 
    171171                                        for(int l=0; l< (int)weights_theta.size(); l++) { 
    … …  
    176176                                                        dp.axis_phi = weights_phi[m].value; 
    177177 
    178                                                         double _ptvalue = weights_length[i].weight 
     178                                                        double _ptvalue = weights_core_thick[i].weight 
    179179                                                                * weights_radius[j].weight 
    180                                                                 * weights_thickness[k].weight 
     180                                                                * weights_layer_thick[k].weight 
    181181                                                                * weights_theta[l].weight 
    182182                                                                * weights_phi[m].weight 
    … …  
    187187                                                        sum += _ptvalue; 
    188188 
    189                                                         norm += weights_length[i].weight 
     189                                                        norm += weights_core_thick[i].weight 
    190190                                                                * weights_radius[j].weight 
    191                                                                 * weights_thickness[k].weight 
     191                                                                * weights_layer_thick[k].weight 
    192192                                                                * weights_theta[l].weight 
    193193                                                                * weights_phi[m].weight; 

  • sansmodels/src/sans/models/c_models/triaxialellipsoid.cpp

    r9188cc1 r975ec8e  
    4444        axis_theta  = Parameter(0.0, true); 
    4545        axis_phi    = Parameter(0.0, true); 
     46        axis_psi    = Parameter(0.0, true); 
    4647} 
    4748 
    … …  
    5354 */ 
    5455double TriaxialEllipsoidModel :: operator()(double q) { 
    55         double dp[5]; 
     56        double dp[6]; 
    5657 
    5758        // Fill parameter array for IGOR library 
    … …  
    119120        dp.axis_theta  = axis_theta(); 
    120121        dp.axis_phi    = axis_phi(); 
     122        dp.axis_psi    = axis_psi(); 
    121123 
    122124        // Get the dispersion points for the semi_axis A 
    … …  
    140142        axis_phi.get_weights(weights_phi); 
    141143 
     144        // Get angular averaging for psi 
     145        vector<WeightPoint> weights_psi; 
     146        axis_psi.get_weights(weights_psi); 
     147 
    142148        // Perform the computation, with all weight points 
    143149        double sum = 0.0; 
    … …  
    163169                                        for(int m=0; m <(int)weights_phi.size(); m++) { 
    164170                                                dp.axis_phi = weights_phi[m].value; 
    165  
    166                                                 double _ptvalue = weights_semi_axisA[i].weight 
    167                                                         * weights_semi_axisB[j].weight 
    168                                                         * weights_semi_axisC[k].weight 
    169                                                         * weights_theta[l].weight 
    170                                                         * weights_phi[m].weight 
    171                                                         * triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy); 
    172                                                 if (weights_theta.size()>1) { 
    173                                                         _ptvalue *= sin(weights_theta[k].value); 
     171                                                // Average over psi distribution 
     172                                                for(int n=0; n <(int)weights_psi.size(); n++) { 
     173                                                        dp.axis_psi = weights_psi[n].value; 
     174 
     175                                                        double _ptvalue = weights_semi_axisA[i].weight 
     176                                                                * weights_semi_axisB[j].weight 
     177                                                                * weights_semi_axisC[k].weight 
     178                                                                * weights_theta[l].weight 
     179                                                                * weights_phi[m].weight 
     180                                                                * weights_psi[n].weight 
     181                                                                * triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy); 
     182                                                        if (weights_theta.size()>1) { 
     183                                                                _ptvalue *= sin(weights_theta[k].value); 
     184                                                        } 
     185                                                        sum += _ptvalue; 
     186 
     187                                                        norm += weights_semi_axisA[i].weight 
     188                                                                * weights_semi_axisB[j].weight 
     189                                                                * weights_semi_axisC[k].weight 
     190                                                                * weights_theta[l].weight 
     191                                                                * weights_phi[m].weight 
     192                                                                * weights_psi[m].weight; 
    174193                                                } 
    175                                                 sum += _ptvalue; 
    176  
    177                                                 norm += weights_semi_axisA[i].weight 
    178                                                         * weights_semi_axisB[j].weight 
    179                                                         * weights_semi_axisC[k].weight 
    180                                                         * weights_theta[l].weight 
    181                                                         * weights_phi[m].weight; 
    182194                                        } 
    183195 
Note: See TracChangeset for help on using the changeset viewer.