Changeset 0957bb3a in sasview for src/sas/sascalc/calculator/c_extensions
- Timestamp:
- Nov 7, 2017 10:00:25 AM (7 years ago)
- Branches:
- master, magnetic_scatt, release-4.2.2, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
- Children:
- 54b0650
- Parents:
- 0fc5a03
- Location:
- src/sas/sascalc/calculator/c_extensions
- Files:
-
- 1 added
- 1 deleted
- 2 moved
-
sld2i.c (moved) (moved from src/sas/sascalc/calculator/c_extensions/sld2i.cpp) (13 diffs)
-
sld2i.h (added)
-
sld2i.hh (deleted)
-
sld2i_module.c (moved) (moved from src/sas/sascalc/calculator/c_extensions/sld2i_module.cpp) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/sas/sascalc/calculator/c_extensions/sld2i.c
rb523c0e r0957bb3a 2 2 Computes the (magnetic) scattering form sld (n and m) profile 3 3 */ 4 #include "sld2i.hh"5 4 #include <stdio.h> 6 5 #include <math.h> 7 using namespace std; 8 extern "C" { 9 #include "libfunc.h" 10 #include "librefl.h" 11 } 6 #include "sld2i.h" 7 #include "libfunc.h" 8 #include "librefl.h" 12 9 /** 13 10 * Constructor for GenI … … 28 25 * @param s_theta: angle (from x-axis) of the up spin in degree 29 26 */ 30 GenI :: GenI(int npix, double* x, double* y, double* z, double* sldn,27 void initGenI(GenI* this, int npix, double* x, double* y, double* z, double* sldn, 31 28 double* mx, double* my, double* mz, double* voli, 32 29 double in_spin, double out_spin, 33 30 double s_theta) { 34 //this->qx_val = qx;35 //this->qy_val = qy;36 //this->qz_val = qz;37 31 this->n_pix = npix; 38 32 this->x_val = x; … … 47 41 this->outspin = out_spin; 48 42 this->stheta = s_theta; 49 } ;43 } 50 44 51 45 /** 52 46 * Compute 2D anisotropic 53 47 */ 54 void GenI :: genicomXY(int npoints, double *qx, double *qy, double *I_out){55 //npoints is given negative for angular averaging 48 void genicomXY(GenI* this, int npoints, double *qx, double *qy, double *I_out){ 49 //npoints is given negative for angular averaging 56 50 // Assumes that q doesn't have qz component and sld_n is all real 57 51 //double q = 0.0; … … 76 70 //int y_size = 0; //in Ang 77 71 //int z_size = 0; //in Ang 78 72 79 73 // Loop over q-values and multiply apply matrix 80 74 81 75 for(int i=0; i<npoints; i++){ 82 76 //I_out[i] = 0.0; … … 84 78 sumj_ud = cassign(0.0, 0.0); 85 79 sumj_du = cassign(0.0, 0.0); 86 sumj_dd = cassign(0.0, 0.0); 80 sumj_dd = cassign(0.0, 0.0); 87 81 //printf ("%d ", i); 88 82 //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); 89 83 90 for(int j=0; j<n_pix; j++){ 91 if (sldn_val[j]!=0.0||mx_val[j]!=0.0||my_val[j]!=0.0||mz_val[j]!=0.0) 92 { 84 for(int j=0; j<this->n_pix; j++){ 85 if (this->sldn_val[j]!=0.0 86 ||this->mx_val[j]!=0.0 87 ||this->my_val[j]!=0.0 88 ||this->mz_val[j]!=0.0) 89 { 93 90 //anisotropic 94 91 temp_fi = cassign(0.0, 0.0); 95 b_sld = cal_msld(0, qx[i], qy[i], sldn_val[j],96 mx_val[j], my_val[j],mz_val[j],97 inspin, outspin,stheta);98 qr = (qx[i]* x_val[j] + qy[i]*y_val[j]);92 b_sld = cal_msld(0, qx[i], qy[i], this->sldn_val[j], 93 this->mx_val[j], this->my_val[j], this->mz_val[j], 94 this->inspin, this->outspin, this->stheta); 95 qr = (qx[i]*this->x_val[j] + qy[i]*this->y_val[j]); 99 96 iqr = cassign(0.0, qr); 100 97 ephase = cplx_exp(iqr); 101 98 102 99 //Let's multiply pixel(atomic) volume here 103 ephase = rcmult( vol_pix[j], ephase);100 ephase = rcmult(this->vol_pix[j], ephase); 104 101 //up_up 105 if ( inspin > 0.0 &&outspin > 0.0){102 if (this->inspin > 0.0 && this->outspin > 0.0){ 106 103 comp_sld = cassign(b_sld.uu, 0.0); 107 104 temp_fi = cplx_mult(comp_sld, ephase); … … 109 106 } 110 107 //down_down 111 if ( inspin < 1.0 &&outspin < 1.0){108 if (this->inspin < 1.0 && this->outspin < 1.0){ 112 109 comp_sld = cassign(b_sld.dd, 0.0); 113 110 temp_fi = cplx_mult(comp_sld, ephase); … … 115 112 } 116 113 //up_down 117 if ( inspin > 0.0 &&outspin < 1.0){114 if (this->inspin > 0.0 && this->outspin < 1.0){ 118 115 comp_sld = cassign(b_sld.re_ud, b_sld.im_ud); 119 116 temp_fi = cplx_mult(comp_sld, ephase); … … 121 118 } 122 119 //down_up 123 if ( inspin < 1.0 &&outspin > 0.0){120 if (this->inspin < 1.0 && this->outspin > 0.0){ 124 121 comp_sld = cassign(b_sld.re_du, b_sld.im_du); 125 122 temp_fi = cplx_mult(comp_sld, ephase); … … 129 126 130 127 if (i == 0){ 131 count += vol_pix[j];128 count += this->vol_pix[j]; 132 129 } 133 130 } … … 149 146 * Also assumes there is no polarization: No dependency on spin 150 147 */ 151 void GenI :: genicom(int npoints, double *q, double *I_out){152 //npoints is given negative for angular averaging 148 void genicom(GenI* this, int npoints, double *q, double *I_out){ 149 //npoints is given negative for angular averaging 153 150 // Assumes that q doesn't have qz component and sld_n is all real 154 151 //double Pi = 4.0*atan(1.0); 155 int is_sym = 0;152 int is_sym = this->n_pix < 0; 156 153 double qr = 0.0; 157 154 double sumj; 158 155 double sld_j = 0.0; 159 156 double count = 0.0; 160 if (n_pix < 0 ){ 161 is_sym = 1; 162 n_pix = n_pix * -1; 163 } 157 int n_pix = is_sym ? -this->n_pix : this->n_pix; 164 158 //Assume that pixel volumes are given in vol_pix in A^3 unit 165 159 // Loop over q-values and multiply apply matrix 166 160 for(int i=0; i<npoints; i++){ 167 sumj =0.0; 161 sumj =0.0; 168 162 for(int j=0; j<n_pix; j++){ 169 163 //Isotropic: Assumes all slds are real (no magnetic) … … 171 165 if (is_sym == 1){ 172 166 // approximation for a spherical symmetric particle 173 qr = sqrt( x_val[j]*x_val[j]+y_val[j]*y_val[j]+z_val[j]*z_val[j])*q[i];167 qr = sqrt(this->x_val[j]*this->x_val[j]+this->y_val[j]*this->y_val[j]+this->z_val[j]*this->z_val[j])*q[i]; 174 168 if (qr > 0.0){ 175 169 qr = sin(qr) / qr; 176 sumj += sldn_val[j] *vol_pix[j] * qr;170 sumj += this->sldn_val[j] * this->vol_pix[j] * qr; 177 171 } 178 172 else{ 179 sumj += sldn_val[j] *vol_pix[j];173 sumj += this->sldn_val[j] * this->vol_pix[j]; 180 174 } 181 175 } … … 184 178 //pragma omp parallel for 185 179 for(int k=0; k<n_pix; k++){ 186 sld_j = sldn_val[j] * sldn_val[k] * vol_pix[j] *vol_pix[k];187 qr = ( x_val[j]-x_val[k])*(x_val[j]-x_val[k])+188 ( y_val[j]-y_val[k])*(y_val[j]-y_val[k])+189 ( z_val[j]-z_val[k])*(z_val[j]-z_val[k]);180 sld_j = this->sldn_val[j] * this->sldn_val[k] * this->vol_pix[j] * this->vol_pix[k]; 181 qr = (this->x_val[j]-this->x_val[k])*(this->x_val[j]-this->x_val[k])+ 182 (this->y_val[j]-this->y_val[k])*(this->y_val[j]-this->y_val[k])+ 183 (this->z_val[j]-this->z_val[k])*(this->z_val[j]-this->z_val[k]); 190 184 qr = sqrt(qr) * q[i]; 191 185 if (qr > 0.0){ … … 198 192 } 199 193 if (i == 0){ 200 count += vol_pix[j];194 count += this->vol_pix[j]; 201 195 } 202 196 } -
src/sas/sascalc/calculator/c_extensions/sld2i_module.c
r1d014cb r0957bb3a 4 4 #include <Python.h> 5 5 #include <stdio.h> 6 #include <sld2i.h h>6 #include <sld2i.h> 7 7 8 8 #if PY_MAJOR_VERSION < 3 … … 35 35 void 36 36 del_sld2i(PyObject *obj){ 37 GenI* sld2i = static_cast<GenI *>(PyCapsule_GetPointer(obj, "GenI")); 38 delete sld2i; 39 return; 37 GenI* sld2i = (GenI *)(PyCapsule_GetPointer(obj, "GenI")); 38 PyMem_Free((void *)sld2i); 40 39 } 41 40 … … 43 42 * Create a GenI as a python object by supplying arrays 44 43 */ 45 PyObject * new_GenI(PyObject * , PyObject *args) {44 PyObject * new_GenI(PyObject *self, PyObject *args) { 46 45 PyObject *x_val_obj; 47 46 PyObject *y_val_obj; … … 79 78 OUTVECTOR(mz_val_obj, mz_val, n_x); 80 79 OUTVECTOR(vol_pix_obj, vol_pix, n_x); 81 GenI* sld2i = new GenI(n_pix,x_val,y_val,z_val,sldn_val,mx_val,my_val,mz_val,vol_pix,inspin,outspin,stheta); 80 GenI* sld2i = PyMem_Malloc(sizeof(GenI)); 81 if (sld2i != NULL) { 82 initGenI(sld2i, n_pix,x_val,y_val,z_val,sldn_val,mx_val,my_val,mz_val,vol_pix,inspin,outspin,stheta); 83 } 82 84 return PyCapsule_New(sld2i, "GenI", del_sld2i); 83 85 } … … 86 88 * GenI the given input (2D) according to a given object 87 89 */ 88 PyObject * genicom_inputXY(PyObject * , PyObject *args) {90 PyObject * genicom_inputXY(PyObject *self, PyObject *args) { 89 91 int npoints; 90 92 PyObject *qx_obj; … … 106 108 107 109 // Set the array pointers 108 void *temp = PyCapsule_GetPointer(gen_obj, "GenI"); 109 GenI* s = static_cast<GenI *>(temp); 110 111 s->genicomXY(npoints, qx, qy, I_out); 110 GenI* sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI"); 111 112 genicomXY(sld2i, npoints, qx, qy, I_out); 112 113 //return PyCObject_FromVoidPtr(s, del_genicom); 113 114 return Py_BuildValue("i",1); … … 117 118 * GenI the given 1D input according to a given object 118 119 */ 119 PyObject * genicom_input(PyObject * , PyObject *args) {120 PyObject * genicom_input(PyObject *self, PyObject *args) { 120 121 int npoints; 121 122 PyObject *q_obj; … … 134 135 135 136 // Set the array pointers 136 void *temp = PyCapsule_GetPointer(gen_obj, "GenI"); 137 GenI* s = static_cast<GenI *>(temp); 138 139 s->genicom(npoints, q, I_out); 137 GenI *sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI"); 138 139 genicom(sld2i, npoints, q, I_out); 140 140 //return PyCObject_FromVoidPtr(s, del_genicom); 141 141 return Py_BuildValue("i",1);
Note: See TracChangeset
for help on using the changeset viewer.
