- Timestamp:
- Dec 6, 2017 11:04:17 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:
- a64f8f4, b6b81a3, 0a88623, 66285f5
- Parents:
- f53d684 (diff), bc4f5bb (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent. - git-author:
- Paul Kienzle <pkienzle@…> (12/06/17 11:04:17)
- git-committer:
- GitHub <noreply@…> (12/06/17 11:04:17)
- Location:
- src/sas
- Files:
-
- 3 added
- 12 deleted
- 27 edited
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
src/sas/sascalc/calculator/c_extensions/libfunc.c
rb6c8abe r144e032a 95 95 // spintheta: angle (anti-clock-wise) between neutron spin(up) and x axis 96 96 // Note: all angles are in degrees. 97 polar_sld cal_msld(int isangle, double qx, double qy, double bn,97 void cal_msld(polar_sld *p_sld, int isangle, double qx, double qy, double bn, 98 98 double m01, double mtheta1, double mphi1, 99 99 double spinfraci, double spinfracf, double spintheta) … … 124 124 double my = 0.0; 125 125 double mz = 0.0; 126 polar_sld p_sld; 127 p_sld.uu = sld; 128 p_sld.dd = sld; 129 p_sld.re_ud = 0.0; 130 p_sld.im_ud = 0.0; 131 p_sld.re_du = 0.0; 132 p_sld.im_du = 0.0; 126 double uu = sld; 127 double dd = sld; 128 double re_ud = 0.0; 129 double im_ud = 0.0; 130 double re_du = 0.0; 131 double im_du = 0.0; 133 132 134 133 //No mag means no further calculation 135 if (isangle>0) {134 if (isangle>0) { 136 135 if (m_max < 1.0e-32){ 137 p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu; 138 p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd; 139 return p_sld; 140 } 141 } 142 else{ 143 if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){ 144 p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu; 145 p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd; 146 return p_sld; 147 } 148 } 149 150 //These are needed because of the precision of inputs 151 if (in_spin < 0.0) in_spin = 0.0; 152 if (in_spin > 1.0) in_spin = 1.0; 153 if (out_spin < 0.0) out_spin = 0.0; 154 if (out_spin > 1.0) out_spin = 1.0; 155 156 if (q_x == 0.0) q_angle = pi / 2.0; 157 else q_angle = atan(q_y/q_x); 158 if (q_y < 0.0 && q_x < 0.0) q_angle -= pi; 159 else if (q_y > 0.0 && q_x < 0.0) q_angle += pi; 160 161 q_angle = pi/2.0 - q_angle; 162 if (q_angle > pi) q_angle -= 2.0 * pi; 163 else if (q_angle < -pi) q_angle += 2.0 * pi; 164 165 if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){ 166 m_perp = 0.0; 167 } 168 else { 169 m_perp = m_max; 170 } 171 if (is_angle > 0){ 172 m_phi *= pi/180.0; 173 m_theta *= pi/180.0; 174 mx = m_perp * cos(m_theta) * cos(m_phi); 175 my = m_perp * sin(m_theta); 176 mz = -(m_perp * cos(m_theta) * sin(m_phi)); 177 } 178 else{ 179 mx = m_perp; 180 my = m_phi; 181 mz = m_theta; 182 } 183 //ToDo: simplify these steps 184 // m_perp1 -m_perp2 185 m_perp_x = (mx) * cos(q_angle); 186 m_perp_x -= (my) * sin(q_angle); 187 m_perp_y = m_perp_x; 188 m_perp_x *= cos(-q_angle); 189 m_perp_y *= sin(-q_angle); 190 m_perp_z = mz; 191 192 m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta)); 193 m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta)); 194 m_sigma_z = (m_perp_z); 195 196 //Find b 197 p_sld.uu -= m_sigma_x; 198 p_sld.dd += m_sigma_x; 199 p_sld.re_ud = m_sigma_y; 200 p_sld.re_du = m_sigma_y; 201 p_sld.im_ud = m_sigma_z; 202 p_sld.im_du = -m_sigma_z; 203 204 p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu; 205 p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd; 206 207 p_sld.re_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.re_ud; 208 p_sld.im_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.im_ud; 209 p_sld.re_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.re_du; 210 p_sld.im_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.im_du; 211 212 return p_sld; 136 uu = sqrt(sqrt(in_spin * out_spin)) * uu; 137 dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * dd; 138 } 139 } 140 else if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){ 141 uu = sqrt(sqrt(in_spin * out_spin)) * uu; 142 dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * dd; 143 } else { 144 145 //These are needed because of the precision of inputs 146 if (in_spin < 0.0) in_spin = 0.0; 147 if (in_spin > 1.0) in_spin = 1.0; 148 if (out_spin < 0.0) out_spin = 0.0; 149 if (out_spin > 1.0) out_spin = 1.0; 150 151 if (q_x == 0.0) q_angle = pi / 2.0; 152 else q_angle = atan(q_y/q_x); 153 if (q_y < 0.0 && q_x < 0.0) q_angle -= pi; 154 else if (q_y > 0.0 && q_x < 0.0) q_angle += pi; 155 156 q_angle = pi/2.0 - q_angle; 157 if (q_angle > pi) q_angle -= 2.0 * pi; 158 else if (q_angle < -pi) q_angle += 2.0 * pi; 159 160 if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){ 161 m_perp = 0.0; 162 } 163 else { 164 m_perp = m_max; 165 } 166 if (is_angle > 0){ 167 m_phi *= pi/180.0; 168 m_theta *= pi/180.0; 169 mx = m_perp * cos(m_theta) * cos(m_phi); 170 my = m_perp * sin(m_theta); 171 mz = -(m_perp * cos(m_theta) * sin(m_phi)); 172 } 173 else{ 174 mx = m_perp; 175 my = m_phi; 176 mz = m_theta; 177 } 178 //ToDo: simplify these steps 179 // m_perp1 -m_perp2 180 m_perp_x = (mx) * cos(q_angle); 181 m_perp_x -= (my) * sin(q_angle); 182 m_perp_y = m_perp_x; 183 m_perp_x *= cos(-q_angle); 184 m_perp_y *= sin(-q_angle); 185 m_perp_z = mz; 186 187 m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta)); 188 m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta)); 189 m_sigma_z = (m_perp_z); 190 191 //Find b 192 uu -= m_sigma_x; 193 dd += m_sigma_x; 194 re_ud = m_sigma_y; 195 re_du = m_sigma_y; 196 im_ud = m_sigma_z; 197 im_du = -m_sigma_z; 198 199 uu = sqrt(sqrt(in_spin * out_spin)) * uu; 200 dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * dd; 201 202 re_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * re_ud; 203 im_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * im_ud; 204 re_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * re_du; 205 im_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * im_du; 206 } 207 p_sld->uu = uu; 208 p_sld->dd = dd; 209 p_sld->re_ud = re_ud; 210 p_sld->im_ud = im_ud; 211 p_sld->re_du = re_du; 212 p_sld->im_du = im_du; 213 213 } 214 214 -
src/sas/sascalc/calculator/c_extensions/libfunc.h
rb6c8abe r144e032a 18 18 //double gamln(double x); 19 19 20 polar_sld cal_msld(int isangle, double qx, double qy, double bn, double m01, double mtheta1,20 void cal_msld(polar_sld*, int isangle, double qx, double qy, double bn, double m01, double mtheta1, 21 21 double mphi1, double spinfraci, double spinfracf, double spintheta); 22 22 -
src/sas/sascalc/calculator/c_extensions/librefl.c
r4c29e4d r144e032a 7 7 #include <stdio.h> 8 8 #include <stdlib.h> 9 #if defined (_MSC_VER)9 #if defined _MSC_VER || defined __TINYCC__ 10 10 #define NEED_ERF 11 11 #endif … … 21 21 22 22 23 #ifdef _WIN32 23 #ifdef __TINYCC__ 24 # ifdef isnan 25 # undef isnan 26 # endif 27 # ifdef isfinite 28 # undef isfinite 29 # endif 30 # define isnan(x) (x != x) 31 # define isfinite(x) (x != INFINITY && x != -INFINITY) 32 #elif defined _WIN32 24 33 # include <float.h> 25 34 # if !defined __MINGW32__ || defined __NO_ISOCEXT … … 30 39 # define isinf(x) (!_finite(x) && !_isnan(x)) 31 40 # endif 32 # ifndef finite33 # define finite(x) _finite(x)41 # ifndef isfinite 42 # define isfinite(x) _finite(x) 34 43 # endif 35 44 # endif … … 84 93 double erf(double x) 85 94 { 86 if (! finite(x)) {95 if (!isfinite(x)) { 87 96 if (isnan(x)) return x; /* erf(NaN) = NaN */ 88 97 return (x>0 ? 1.0 : -1.0); /* erf(+-inf) = +-1.0 */ … … 94 103 double erfc(double x) 95 104 { 96 if (! finite(x)) {105 if (!isfinite(x)) { 97 106 if (isnan(x)) return x; /* erfc(NaN) = NaN */ 98 107 return (x>0 ? 0.0 : 2.0); /* erfc(+-inf) = 0.0, 2.0 */ … … 103 112 #endif // NEED_ERF 104 113 105 complex cassign(real, imag) 106 double real, imag; 107 { 108 complex x; 109 x.re = real; 110 x.im = imag; 111 return x; 112 } 113 114 115 complex cplx_add(x,y) 116 complex x,y; 117 { 118 complex z; 119 z.re = x.re + y.re; 120 z.im = x.im + y.im; 121 return z; 122 } 123 124 complex rcmult(x,y) 125 double x; 126 complex y; 127 { 128 complex z; 129 z.re = x*y.re; 130 z.im = x*y.im; 131 return z; 132 } 133 134 complex cplx_sub(x,y) 135 complex x,y; 136 { 137 complex z; 138 z.re = x.re - y.re; 139 z.im = x.im - y.im; 140 return z; 141 } 142 143 144 complex cplx_mult(x,y) 145 complex x,y; 146 { 147 complex z; 148 z.re = x.re*y.re - x.im*y.im; 149 z.im = x.re*y.im + x.im*y.re; 150 return z; 151 } 152 153 complex cplx_div(x,y) 154 complex x,y; 155 { 156 complex z; 157 z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im); 158 z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im); 159 return z; 160 } 161 162 complex cplx_exp(b) 163 complex b; 164 { 165 complex z; 114 void cassign(Cplx *x, double real, double imag) 115 { 116 x->re = real; 117 x->im = imag; 118 } 119 120 121 void cplx_add(Cplx *z, Cplx x, Cplx y) 122 { 123 z->re = x.re + y.re; 124 z->im = x.im + y.im; 125 } 126 127 void rcmult(Cplx *z, double x, Cplx y) 128 { 129 z->re = x*y.re; 130 z->im = x*y.im; 131 } 132 133 void cplx_sub(Cplx *z, Cplx x, Cplx y) 134 { 135 z->re = x.re - y.re; 136 z->im = x.im - y.im; 137 } 138 139 140 void cplx_mult(Cplx *z, Cplx x, Cplx y) 141 { 142 z->re = x.re*y.re - x.im*y.im; 143 z->im = x.re*y.im + x.im*y.re; 144 } 145 146 void cplx_div(Cplx *z, Cplx x, Cplx y) 147 { 148 z->re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im); 149 z->im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im); 150 } 151 152 void cplx_exp(Cplx *z, Cplx b) 153 { 166 154 double br,bi; 167 155 br=b.re; 168 156 bi=b.im; 169 z.re = exp(br)*cos(bi); 170 z.im = exp(br)*sin(bi); 171 return z; 172 } 173 174 175 complex cplx_sqrt(z) //see Schaum`s Math Handbook p. 22, 6.6 and 6.10 176 complex z; 177 { 178 complex c; 157 z->re = exp(br)*cos(bi); 158 z->im = exp(br)*sin(bi); 159 } 160 161 162 void cplx_sqrt(Cplx *c, Cplx z) //see Schaum`s Math Handbook p. 22, 6.6 and 6.10 163 { 179 164 double zr,zi,x,y,r,w; 180 165 … … 184 169 if (zr==0.0 && zi==0.0) 185 170 { 186 c.re=0.0; 187 c.im=0.0; 188 return c; 189 } 190 else 191 { 171 c->re=0.0; 172 c->im=0.0; 173 } else { 192 174 x=fabs(zr); 193 175 y=fabs(zi); … … 196 178 r=y/x; 197 179 w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))); 198 } 199 else 200 { 180 } else { 201 181 r=x/y; 202 182 w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))); … … 204 184 if (zr >=0.0) 205 185 { 206 c.re=w; 207 c.im=zi/(2.0*w); 186 c->re=w; 187 c->im=zi/(2.0*w); 188 } else { 189 c->im=(zi >= 0) ? w : -w; 190 c->re=zi/(2.0*c->im); 208 191 } 209 else 210 { 211 c.im=(zi >= 0) ? w : -w; 212 c.re=zi/(2.0*c.im); 213 } 214 return c; 215 } 216 } 217 218 complex cplx_cos(b) 219 complex b; 220 { 221 complex zero,two,z,i,bi,negbi; 222 zero = cassign(0.0,0.0); 223 two = cassign(2.0,0.0); 224 i = cassign(0.0,1.0); 225 bi = cplx_mult(b,i); 226 negbi = cplx_sub(zero,bi); 227 z = cplx_div(cplx_add(cplx_exp(bi),cplx_exp(negbi)),two); 228 return z; 192 } 193 } 194 195 void cplx_cos(Cplx *z, Cplx b) 196 { 197 // cos(b) = (e^bi + e^-bi)/2 198 // = (e^b.im e^-i bi.re) + e^-b.im e^i b.re)/2 199 // = (e^b.im cos(-b.re) + e^b.im sin(-b.re) i)/2 + (e^-b.im cos(b.re) + e^-b.im sin(b.re) i)/2 200 // = e^b.im cos(b.re)/2 - e^b.im sin(b.re)/2 i + 1/e^b.im cos(b.re)/2 + 1/e^b.im sin(b.re)/2 i 201 // = (e^b.im + 1/e^b.im)/2 cos(b.re) + (-e^b.im + 1/e^b.im)/2 sin(b.re) i 202 // = cosh(b.im) cos(b.re) - sinh(b.im) sin(b.re) i 203 double exp_b_im = exp(b.im); 204 z->re = 0.5*(+exp_b_im + 1.0/exp_b_im) * cos(b.re); 205 z->im = -0.5*(exp_b_im - 1.0/exp_b_im) * sin(b.re); 229 206 } 230 207 -
src/sas/sascalc/calculator/c_extensions/librefl.h
r9e531f2 r144e032a 5 5 double re; 6 6 double im; 7 } complex;7 } Cplx; 8 8 9 9 typedef struct { 10 complex a;11 complex b;12 complex c;13 complex d;10 Cplx a; 11 Cplx b; 12 Cplx c; 13 Cplx d; 14 14 } matrix; 15 15 16 complex cassign(double real, double imag);16 void cassign(Cplx*, double real, double imag); 17 17 18 complex cplx_add(complex x,complex y);18 void cplx_add(Cplx*, Cplx x,Cplx y); 19 19 20 complex rcmult(double x,complex y);20 void rcmult(Cplx*, double x,Cplx y); 21 21 22 complex cplx_sub(complex x,complex y);22 void cplx_sub(Cplx*, Cplx x,Cplx y); 23 23 24 complex cplx_mult(complex x,complex y);24 void cplx_mult(Cplx*, Cplx x,Cplx y); 25 25 26 complex cplx_div(complex x,complex y);26 void cplx_div(Cplx*, Cplx x,Cplx y); 27 27 28 complex cplx_exp(complex b);28 void cplx_exp(Cplx*, Cplx b); 29 29 30 complex cplx_sqrt(complex z);30 void cplx_sqrt(Cplx*, Cplx z); 31 31 32 complex cplx_cos(complex b);32 void cplx_cos(Cplx*, Cplx b); 33 33 34 34 double intersldfunc(int fun_type, double n_sub, double i, double nu, double sld_l, double sld_r); -
src/sas/sascalc/calculator/c_extensions/sld2i.c
rb523c0e r144e032a 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 is_avg, 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; 31 this->is_avg = is_avg; 37 32 this->n_pix = npix; 38 33 this->x_val = x; … … 47 42 this->outspin = out_spin; 48 43 this->stheta = s_theta; 49 } ;44 } 50 45 51 46 /** 52 47 * Compute 2D anisotropic 53 48 */ 54 void GenI :: genicomXY(int npoints, double *qx, double *qy, double *I_out){55 //npoints is given negative for angular averaging 49 void genicomXY(GenI* this, int npoints, double *qx, double *qy, double *I_out){ 50 //npoints is given negative for angular averaging 56 51 // Assumes that q doesn't have qz component and sld_n is all real 57 52 //double q = 0.0; … … 59 54 polar_sld b_sld; 60 55 double qr = 0.0; 61 complex iqr = cassign(0.0, 0.0);62 complex ephase = cassign(0.0, 0.0);63 complex comp_sld = cassign(0.0, 0.0);64 65 complex sumj_uu;66 complex sumj_ud;67 complex sumj_du;68 complex sumj_dd;69 complex temp_fi;56 Cplx iqr; 57 Cplx ephase; 58 Cplx comp_sld; 59 60 Cplx sumj_uu; 61 Cplx sumj_ud; 62 Cplx sumj_du; 63 Cplx sumj_dd; 64 Cplx temp_fi; 70 65 71 66 double count = 0.0; 72 //check if this computation is for averaging 67 int i, j; 68 69 cassign(&iqr, 0.0, 0.0); 70 cassign(&ephase, 0.0, 0.0); 71 cassign(&comp_sld, 0.0, 0.0); 73 72 74 73 //Assume that pixel volumes are given in vol_pix in A^3 unit … … 76 75 //int y_size = 0; //in Ang 77 76 //int z_size = 0; //in Ang 78 77 79 78 // Loop over q-values and multiply apply matrix 80 81 for(int i=0; i<npoints; i++){ 79 80 //printf("npoints: %d, npix: %d\n", npoints, this->n_pix); 81 for(i=0; i<npoints; i++){ 82 82 //I_out[i] = 0.0; 83 sumj_uu = cassign(0.0, 0.0);84 sumj_ud = cassign(0.0, 0.0);85 sumj_du = cassign(0.0, 0.0);86 sumj_dd = cassign(0.0, 0.0);87 //printf ("%d", i);83 cassign(&sumj_uu, 0.0, 0.0); 84 cassign(&sumj_ud, 0.0, 0.0); 85 cassign(&sumj_du, 0.0, 0.0); 86 cassign(&sumj_dd, 0.0, 0.0); 87 //printf("i: %d\n", i); 88 88 //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); 89 89 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 { 90 for(j=0; j<this->n_pix; j++){ 91 if (this->sldn_val[j]!=0.0 92 ||this->mx_val[j]!=0.0 93 ||this->my_val[j]!=0.0 94 ||this->mz_val[j]!=0.0) 95 { 96 // printf("i,j: %d,%d\n", i,j); 93 97 //anisotropic 94 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]);99 iqr = cassign(0.0, qr);100 ephase = cplx_exp(iqr);101 98 cassign(&temp_fi, 0.0, 0.0); 99 cal_msld(&b_sld, 0, qx[i], qy[i], this->sldn_val[j], 100 this->mx_val[j], this->my_val[j], this->mz_val[j], 101 this->inspin, this->outspin, this->stheta); 102 qr = (qx[i]*this->x_val[j] + qy[i]*this->y_val[j]); 103 cassign(&iqr, 0.0, qr); 104 cplx_exp(&ephase, iqr); 105 102 106 //Let's multiply pixel(atomic) volume here 103 ephase = rcmult(vol_pix[j], ephase);107 rcmult(&ephase, this->vol_pix[j], ephase); 104 108 //up_up 105 if ( inspin > 0.0 &&outspin > 0.0){106 c omp_sld = cassign(b_sld.uu, 0.0);107 temp_fi = cplx_mult(comp_sld, ephase);108 sumj_uu = cplx_add(sumj_uu, temp_fi);109 if (this->inspin > 0.0 && this->outspin > 0.0){ 110 cassign(&comp_sld, b_sld.uu, 0.0); 111 cplx_mult(&temp_fi, comp_sld, ephase); 112 cplx_add(&sumj_uu, sumj_uu, temp_fi); 109 113 } 110 114 //down_down 111 if ( inspin < 1.0 &&outspin < 1.0){112 c omp_sld = cassign(b_sld.dd, 0.0);113 temp_fi = cplx_mult(comp_sld, ephase);114 sumj_dd = cplx_add(sumj_dd, temp_fi);115 if (this->inspin < 1.0 && this->outspin < 1.0){ 116 cassign(&comp_sld, b_sld.dd, 0.0); 117 cplx_mult(&temp_fi, comp_sld, ephase); 118 cplx_add(&sumj_dd, sumj_dd, temp_fi); 115 119 } 116 120 //up_down 117 if ( inspin > 0.0 &&outspin < 1.0){118 c omp_sld = cassign(b_sld.re_ud, b_sld.im_ud);119 temp_fi = cplx_mult(comp_sld, ephase);120 sumj_ud = cplx_add(sumj_ud, temp_fi);121 if (this->inspin > 0.0 && this->outspin < 1.0){ 122 cassign(&comp_sld, b_sld.re_ud, b_sld.im_ud); 123 cplx_mult(&temp_fi, comp_sld, ephase); 124 cplx_add(&sumj_ud, sumj_ud, temp_fi); 121 125 } 122 126 //down_up 123 if (inspin < 1.0 && outspin > 0.0){ 124 comp_sld = cassign(b_sld.re_du, b_sld.im_du); 125 temp_fi = cplx_mult(comp_sld, ephase); 126 sumj_du = cplx_add(sumj_du, temp_fi); 127 } 128 127 if (this->inspin < 1.0 && this->outspin > 0.0){ 128 cassign(&comp_sld, b_sld.re_du, b_sld.im_du); 129 cplx_mult(&temp_fi, comp_sld, ephase); 130 cplx_add(&sumj_du, sumj_du, temp_fi); 131 } 129 132 130 133 if (i == 0){ 131 count += vol_pix[j];134 count += this->vol_pix[j]; 132 135 } 133 136 } … … 142 145 I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix 143 146 } 144 //printf ("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0],mz_val[0]);147 //printf("count = %d %g %g %g %g\n", count, this->sldn_val[0],this->mx_val[0], this->my_val[0], this->mz_val[0]); 145 148 } 146 149 /** … … 149 152 * Also assumes there is no polarization: No dependency on spin 150 153 */ 151 void GenI :: genicom(int npoints, double *q, double *I_out){152 //npoints is given negative for angular averaging 154 void genicom(GenI* this, int npoints, double *q, double *I_out){ 155 //npoints is given negative for angular averaging 153 156 // Assumes that q doesn't have qz component and sld_n is all real 154 157 //double Pi = 4.0*atan(1.0); 155 int is_sym = 0;156 158 double qr = 0.0; 157 159 double sumj; 158 160 double sld_j = 0.0; 159 161 double count = 0.0; 160 if (n_pix < 0 ){ 161 is_sym = 1; 162 n_pix = n_pix * -1; 163 } 162 int i, j, k; 163 164 164 //Assume that pixel volumes are given in vol_pix in A^3 unit 165 165 // Loop over q-values and multiply apply matrix 166 for(i nt i=0; i<npoints; i++){167 sumj =0.0; 168 for( int j=0; j<n_pix; j++){166 for(i=0; i<npoints; i++){ 167 sumj =0.0; 168 for(j=0; j<this->n_pix; j++){ 169 169 //Isotropic: Assumes all slds are real (no magnetic) 170 170 //Also assumes there is no polarization: No dependency on spin 171 if ( is_sym== 1){171 if (this->is_avg == 1){ 172 172 // 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];173 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 174 if (qr > 0.0){ 175 175 qr = sin(qr) / qr; 176 sumj += sldn_val[j] *vol_pix[j] * qr;176 sumj += this->sldn_val[j] * this->vol_pix[j] * qr; 177 177 } 178 178 else{ 179 sumj += sldn_val[j] *vol_pix[j];179 sumj += this->sldn_val[j] * this->vol_pix[j]; 180 180 } 181 181 } … … 183 183 //full calculation 184 184 //pragma omp parallel for 185 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]);185 for(k=0; k<this->n_pix; k++){ 186 sld_j = this->sldn_val[j] * this->sldn_val[k] * this->vol_pix[j] * this->vol_pix[k]; 187 qr = (this->x_val[j]-this->x_val[k])*(this->x_val[j]-this->x_val[k])+ 188 (this->y_val[j]-this->y_val[k])*(this->y_val[j]-this->y_val[k])+ 189 (this->z_val[j]-this->z_val[k])*(this->z_val[j]-this->z_val[k]); 190 190 qr = sqrt(qr) * q[i]; 191 191 if (qr > 0.0){ … … 198 198 } 199 199 if (i == 0){ 200 count += vol_pix[j];200 count += this->vol_pix[j]; 201 201 } 202 202 } 203 203 I_out[i] = sumj; 204 if ( is_sym == 1){204 if (this->is_avg == 1) { 205 205 I_out[i] *= sumj; 206 206 } 207 207 I_out[i] *= (1.0E+8 / count); //in cm (unit) / number; //to be multiplied by vol_pix 208 208 } 209 //printf 209 //printf("count = %d %g %g %g %g\n", count, sldn_val[0],mx_val[0], my_val[0], mz_val[0]); 210 210 } -
src/sas/sascalc/calculator/c_extensions/sld2i_module.c
r1d014cb ra1daf86 4 4 #include <Python.h> 5 5 #include <stdio.h> 6 #include <sld2i.hh>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 #if PY_MAJOR_VERSION < 3 38 GenI* sld2i = (GenI *)obj; 39 #else 40 GenI* sld2i = (GenI *)(PyCapsule_GetPointer(obj, "GenI")); 41 #endif 42 PyMem_Free((void *)sld2i); 40 43 } 41 44 … … 43 46 * Create a GenI as a python object by supplying arrays 44 47 */ 45 PyObject * new_GenI(PyObject * , PyObject *args) {48 PyObject * new_GenI(PyObject *self, PyObject *args) { 46 49 PyObject *x_val_obj; 47 50 PyObject *y_val_obj; … … 52 55 PyObject *mz_val_obj; 53 56 PyObject *vol_pix_obj; 54 Py_ssize_t n_x; 55 //PyObject rlimit_obj; 56 //PyObject npoints_obj; 57 //PyObject nrbins_obj; 58 //PyObject nphibins_obj; 59 int n_pix; 57 Py_ssize_t n_x, n_y, n_z, n_sld, n_mx, n_my, n_mz, n_vol_pix; 58 int is_avg; 60 59 double* x_val; 61 60 double* y_val; … … 69 68 double outspin; 70 69 double stheta; 71 72 if (!PyArg_ParseTuple(args, "iOOOOOOOOddd", &n_pix, &x_val_obj, &y_val_obj, &z_val_obj, &sldn_val_obj, &mx_val_obj, &my_val_obj, &mz_val_obj, &vol_pix_obj, &inspin, &outspin, &stheta)) return NULL; 73 OUTVECTOR(x_val_obj, x_val, n_x); 74 OUTVECTOR(y_val_obj, y_val, n_x); 75 OUTVECTOR(z_val_obj, z_val, n_x); 76 OUTVECTOR(sldn_val_obj, sldn_val, n_x); 77 OUTVECTOR(mx_val_obj, mx_val, n_x); 78 OUTVECTOR(my_val_obj, my_val, n_x); 79 OUTVECTOR(mz_val_obj, mz_val, n_x); 80 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); 82 return PyCapsule_New(sld2i, "GenI", del_sld2i); 70 PyObject *obj; 71 GenI* sld2i; 72 73 //printf("new GenI\n"); 74 if (!PyArg_ParseTuple(args, "iOOOOOOOOddd", &is_avg, &x_val_obj, &y_val_obj, &z_val_obj, &sldn_val_obj, &mx_val_obj, &my_val_obj, &mz_val_obj, &vol_pix_obj, &inspin, &outspin, &stheta)) return NULL; 75 INVECTOR(x_val_obj, x_val, n_x); 76 INVECTOR(y_val_obj, y_val, n_y); 77 INVECTOR(z_val_obj, z_val, n_z); 78 INVECTOR(sldn_val_obj, sldn_val, n_sld); 79 INVECTOR(mx_val_obj, mx_val, n_mx); 80 INVECTOR(my_val_obj, my_val, n_my); 81 INVECTOR(mz_val_obj, mz_val, n_mz); 82 INVECTOR(vol_pix_obj, vol_pix, n_vol_pix); 83 sld2i = PyMem_Malloc(sizeof(GenI)); 84 //printf("sldi:%p\n", sld2i); 85 if (sld2i != NULL) { 86 initGenI(sld2i,is_avg,(int)n_x,x_val,y_val,z_val,sldn_val,mx_val,my_val,mz_val,vol_pix,inspin,outspin,stheta); 87 } 88 obj = PyCapsule_New(sld2i, "GenI", del_sld2i); 89 //printf("constructed %p\n", obj); 90 return obj; 83 91 } 84 92 … … 86 94 * GenI the given input (2D) according to a given object 87 95 */ 88 PyObject * genicom_inputXY(PyObject * , PyObject *args) {89 int npoints;96 PyObject * genicom_inputXY(PyObject *self, PyObject *args) { 97 PyObject *gen_obj; 90 98 PyObject *qx_obj; 99 PyObject *qy_obj; 100 PyObject *I_out_obj; 101 Py_ssize_t n_qx, n_qy, n_out; 91 102 double *qx; 92 PyObject *qy_obj;93 103 double *qy; 94 PyObject *I_out_obj;95 Py_ssize_t n_out;96 104 double *I_out; 97 PyObject *gen_obj; 98 99 if (!PyArg_ParseTuple(args, "OiOOO", &gen_obj, &npoints, &qx_obj, &qy_obj, &I_out_obj)) return NULL; 100 OUTVECTOR(qx_obj, qx, n_out); 101 OUTVECTOR(qy_obj, qy, n_out); 105 GenI* sld2i; 106 107 //printf("in genicom_inputXY\n"); 108 if (!PyArg_ParseTuple(args, "OOOO", &gen_obj, &qx_obj, &qy_obj, &I_out_obj)) return NULL; 109 sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI"); 110 INVECTOR(qx_obj, qx, n_qx); 111 INVECTOR(qy_obj, qy, n_qy); 102 112 OUTVECTOR(I_out_obj, I_out, n_out); 113 //printf("qx, qy, I_out: %d %d %d, %d %d %d\n", qx, qy, I_out, n_qx, n_qy, n_out); 103 114 104 115 // Sanity check 105 //if(n_in!=n_out) return Py_BuildValue("i",-1); 106 107 // 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); 116 //if(n_q!=n_out) return Py_BuildValue("i",-1); 117 118 genicomXY(sld2i, (int)n_qx, qx, qy, I_out); 119 //printf("done calc\n"); 112 120 //return PyCObject_FromVoidPtr(s, del_genicom); 113 121 return Py_BuildValue("i",1); … … 117 125 * GenI the given 1D input according to a given object 118 126 */ 119 PyObject * genicom_input(PyObject * , PyObject *args) {120 int npoints;127 PyObject * genicom_input(PyObject *self, PyObject *args) { 128 PyObject *gen_obj; 121 129 PyObject *q_obj; 130 PyObject *I_out_obj; 131 Py_ssize_t n_q, n_out; 122 132 double *q; 123 PyObject *I_out_obj;124 Py_ssize_t n_out;125 133 double *I_out; 126 PyObject *gen_obj; 127 128 if (!PyArg_ParseTuple(args, "OiOO", &gen_obj, &npoints, &q_obj, &I_out_obj)) return NULL; 129 OUTVECTOR(q_obj, q, n_out); 134 GenI *sld2i; 135 136 if (!PyArg_ParseTuple(args, "OOO", &gen_obj, &q_obj, &I_out_obj)) return NULL; 137 sld2i = (GenI *)PyCapsule_GetPointer(gen_obj, "GenI"); 138 INVECTOR(q_obj, q, n_q); 130 139 OUTVECTOR(I_out_obj, I_out, n_out); 131 140 132 141 // Sanity check 133 //if(n_in!=n_out) return Py_BuildValue("i",-1); 134 135 // 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); 140 //return PyCObject_FromVoidPtr(s, del_genicom); 142 //if (n_q!=n_out) return Py_BuildValue("i",-1); 143 144 genicom(sld2i, (int)n_q, q, I_out); 141 145 return Py_BuildValue("i",1); 142 146 } -
src/sas/sascalc/calculator/sas_gen.py
rb58265c3 r144e032a 118 118 self.is_avg = is_avg 119 119 120 def _gen(self, x, y, i):120 def _gen(self, qx, qy): 121 121 """ 122 122 Evaluate the function … … 129 129 pos_y = self.data_y 130 130 pos_z = self.data_z 131 len_x = len(pos_x)132 131 if self.is_avg is None: 133 len_x *= -1134 132 pos_x, pos_y, pos_z = transform_center(pos_x, pos_y, pos_z) 135 len_q = len(x)136 133 sldn = copy.deepcopy(self.data_sldn) 137 134 sldn -= self.params['solvent_SLD'] 138 model = mod.new_GenI(len_x, pos_x, pos_y, pos_z, 139 sldn, self.data_mx, self.data_my, 140 self.data_mz, self.data_vol, 141 self.params['Up_frac_in'], 142 self.params['Up_frac_out'], 143 self.params['Up_theta']) 144 if y == []: 145 mod.genicom(model, len_q, x, i) 146 else: 147 mod.genicomXY(model, len_q, x, y, i) 135 # **** WARNING **** new_GenI holds pointers to numpy vectors 136 # be sure that they are contiguous double precision arrays and make 137 # sure the GC doesn't eat them before genicom is called. 138 # TODO: rewrite so that the parameters are passed directly to genicom 139 args = ( 140 (1 if self.is_avg else 0), 141 pos_x, pos_y, pos_z, 142 sldn, self.data_mx, self.data_my, 143 self.data_mz, self.data_vol, 144 self.params['Up_frac_in'], 145 self.params['Up_frac_out'], 146 self.params['Up_theta']) 147 model = mod.new_GenI(*args) 148 if len(qy): 149 qx, qy = _vec(qx), _vec(qy) 150 I_out = np.empty_like(qx) 151 #print("npoints", qx.shape, "npixels", pos_x.shape) 152 mod.genicomXY(model, qx, qy, I_out) 153 #print("I_out after", I_out) 154 else: 155 qx = _vec(qx) 156 I_out = np.empty_like(qx) 157 mod.genicom(model, qx, I_out) 148 158 vol_correction = self.data_total_volume / self.params['total_volume'] 149 return self.params['scale'] * vol_correction * i + \ 150 self.params['background'] 159 result = (self.params['scale'] * vol_correction * I_out 160 + self.params['background']) 161 return result 151 162 152 163 def set_sld_data(self, sld_data=None): … … 156 167 self.sld_data = sld_data 157 168 self.data_pos_unit = sld_data.pos_unit 158 self.data_x = sld_data.pos_x159 self.data_y = sld_data.pos_y160 self.data_z = sld_data.pos_z161 self.data_sldn = sld_data.sld_n162 self.data_mx = sld_data.sld_mx163 self.data_my = sld_data.sld_my164 self.data_mz = sld_data.sld_mz165 self.data_vol = sld_data.vol_pix169 self.data_x = _vec(sld_data.pos_x) 170 self.data_y = _vec(sld_data.pos_y) 171 self.data_z = _vec(sld_data.pos_z) 172 self.data_sldn = _vec(sld_data.sld_n) 173 self.data_mx = _vec(sld_data.sld_mx) 174 self.data_my = _vec(sld_data.sld_my) 175 self.data_mz = _vec(sld_data.sld_mz) 176 self.data_vol = _vec(sld_data.vol_pix) 166 177 self.data_total_volume = sum(sld_data.vol_pix) 167 178 self.params['total_volume'] = sum(sld_data.vol_pix) … … 180 191 :return: (I value) 181 192 """ 182 if x.__class__.__name__ == 'list':193 if isinstance(x, list): 183 194 if len(x[1]) > 0: 184 195 msg = "Not a 1D." 185 196 raise ValueError(msg) 186 i_out = np.zeros_like(x[0])187 197 # 1D I is found at y =0 in the 2D pattern 188 out = self._gen(x[0], [] , i_out)198 out = self._gen(x[0], []) 189 199 return out 190 200 else: … … 199 209 :Use this runXY() for the computation 200 210 """ 201 if x.__class__.__name__ == 'list': 202 i_out = np.zeros_like(x[0]) 203 out = self._gen(x[0], x[1], i_out) 204 return out 211 if isinstance(x, list): 212 return self._gen(x[0], x[1]) 205 213 else: 206 214 msg = "Q must be given as list of qx's and qy's" … … 214 222 where qx,qy are 1D ndarrays (for 2D). 215 223 """ 216 if qdist.__class__.__name__ == 'list': 217 if len(qdist[1]) < 1: 218 out = self.run(qdist) 219 else: 220 out = self.runXY(qdist) 221 return out 224 if isinstance(qdist, list): 225 return self.run(qdist) if len(qdist[1]) < 1 else self.runXY(qdist) 222 226 else: 223 227 mesg = "evalDistribution is expecting an ndarray of " 224 228 mesg += "a list [qx,qy] where qx,qy are arrays." 225 229 raise RuntimeError(mesg) 230 231 def _vec(v): 232 return np.ascontiguousarray(v, 'd') 226 233 227 234 class OMF2SLD(object): … … 1041 1048 self.line_z = line_z 1042 1049 1050 def _get_data_path(*path_parts): 1051 from os.path import realpath, join as joinpath, dirname, abspath 1052 # in sas/sascalc/calculator; want sas/sasview/test 1053 return joinpath(dirname(realpath(__file__)), 1054 '..', '..', 'sasview', 'test', *path_parts) 1055 1043 1056 def test_load(): 1044 1057 """ … … 1046 1059 """ 1047 1060 from mpl_toolkits.mplot3d import Axes3D 1048 current_dir = os.path.abspath(os.path.curdir) 1049 print(current_dir) 1050 for i in range(6): 1051 current_dir, _ = os.path.split(current_dir) 1052 tfile = os.path.join(current_dir, "test", "CoreXY_ShellZ.txt") 1053 ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf") 1054 if os.path.isfile(tfile): 1055 tfpath = tfile 1056 ofpath = ofile 1057 break 1061 tfpath = _get_data_path("1d_data", "CoreXY_ShellZ.txt") 1062 ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf") 1063 if not os.path.isfile(tfpath) or not os.path.isfile(ofpath): 1064 raise ValueError("file(s) not found: %r, %r"%(tfpath, ofpath)) 1058 1065 reader = SLDReader() 1059 1066 oreader = OMFReader() 1060 output = decode(reader.read(tfpath))1061 ooutput = decode(oreader.read(ofpath))1067 output = reader.read(tfpath) 1068 ooutput = oreader.read(ofpath) 1062 1069 foutput = OMF2SLD() 1063 1070 foutput.set_data(ooutput) … … 1088 1095 plt.show() 1089 1096 1097 def test_save(): 1098 ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf") 1099 if not os.path.isfile(ofpath): 1100 raise ValueError("file(s) not found: %r"%(ofpath,)) 1101 oreader = OMFReader() 1102 omfdata = oreader.read(ofpath) 1103 omf2sld = OMF2SLD() 1104 omf2sld.set_data(omfdata) 1105 writer = SLDReader() 1106 writer.write("out.txt", omf2sld.output) 1107 1090 1108 def test(): 1091 1109 """ 1092 1110 Test code 1093 1111 """ 1094 current_dir = os.path.abspath(os.path.curdir) 1095 for i in range(3): 1096 current_dir, _ = os.path.split(current_dir) 1097 ofile = os.path.join(current_dir, "test", "A_Raw_Example-1.omf") 1098 if os.path.isfile(ofile): 1099 ofpath = ofile 1100 break 1112 ofpath = _get_data_path("coordinate_data", "A_Raw_Example-1.omf") 1113 if not os.path.isfile(ofpath): 1114 raise ValueError("file(s) not found: %r"%(ofpath,)) 1101 1115 oreader = OMFReader() 1102 ooutput = decode(oreader.read(ofpath)) 1103 foutput = OMF2SLD() 1104 foutput.set_data(ooutput) 1105 writer = SLDReader() 1106 writer.write(os.path.join(os.path.dirname(ofpath), "out.txt"), 1107 foutput.output) 1116 omfdata = oreader.read(ofpath) 1117 omf2sld = OMF2SLD() 1118 omf2sld.set_data(omfdata) 1108 1119 model = GenSAS() 1109 model.set_sld_data(foutput.output) 1110 x = np.arange(1000)/10000. + 1e-5 1111 y = np.arange(1000)/10000. + 1e-5 1112 i = np.zeros(1000) 1113 model.runXY([x, y, i]) 1120 model.set_sld_data(omf2sld.output) 1121 x = np.linspace(0, 0.1, 11)[1:] 1122 return model.runXY([x, x]) 1114 1123 1115 1124 if __name__ == "__main__": 1125 #test_load() 1126 #test_save() 1127 #print(test()) 1116 1128 test() 1117 test_load() -
src/sas/sascalc/dataloader/file_reader_base_class.py
r20fa5fe r3053a4a 27 27 28 28 class FileReader(object): 29 # List of Data1D and Data2D objects to be sent back to data_loader30 output = []31 # Current plottable_(1D/2D) object being loaded in32 current_dataset = None33 # Current DataInfo object being loaded in34 current_datainfo = None35 29 # String to describe the type of data this reader can load 36 30 type_name = "ASCII" … … 43 37 # Able to import the unit converter 44 38 has_converter = True 45 # Open file handle46 f_open = None47 39 # Default value of zero 48 40 _ZERO = 1e-16 49 41 42 def __init__(self): 43 # List of Data1D and Data2D objects to be sent back to data_loader 44 self.output = [] 45 # Current plottable_(1D/2D) object being loaded in 46 self.current_dataset = None 47 # Current DataInfo object being loaded in 48 self.current_datainfo = None 49 # File path sent to reader 50 self.filepath = None 51 # Open file handle 52 self.f_open = None 53 50 54 def read(self, filepath): 51 55 """ … … 54 58 :param filepath: The full or relative path to a file to be loaded 55 59 """ 60 self.filepath = filepath 56 61 if os.path.isfile(filepath): 57 62 basename, extension = os.path.splitext(os.path.basename(filepath)) … … 96 101 self.current_datainfo = None 97 102 self.current_dataset = None 103 self.filepath = None 98 104 self.output = [] 99 105 -
src/sas/sascalc/dataloader/readers/cansas_reader.py
r61f329f0 r2469df7 82 82 83 83 def read(self, xml_file, schema_path="", invalid=True): 84 if schema_path != "" or invalid != True:84 if schema_path != "" or not invalid: 85 85 # read has been called from self.get_file_contents because xml file doens't conform to schema 86 86 _, self.extension = os.path.splitext(os.path.basename(xml_file)) … … 942 942 pos, "z", datainfo.sample.position.z, 943 943 {"unit": datainfo.sample.position_unit}) 944 if written == True:944 if written: 945 945 self.append(pos, sample) 946 946 … … 955 955 ori, "yaw", datainfo.sample.orientation.z, 956 956 {"unit": datainfo.sample.orientation_unit}) 957 if written == True:957 if written: 958 958 self.append(ori, sample) 959 959 … … 1002 1002 size, "z", datainfo.source.beam_size.z, 1003 1003 {"unit": datainfo.source.beam_size_unit}) 1004 if written == True:1004 if written: 1005 1005 self.append(size, source) 1006 1006 … … 1058 1058 size, "z", aperture.size.z, 1059 1059 {"unit": aperture.size_unit}) 1060 if written == True:1060 if written: 1061 1061 self.append(size, apert) 1062 1062 … … 1081 1081 written = written | self.write_node(det, "SDD", item.distance, 1082 1082 {"unit": item.distance_unit}) 1083 if written == True:1083 if written: 1084 1084 self.append(det, instr) 1085 1085 … … 1091 1091 written = written | self.write_node(off, "z", item.offset.z, 1092 1092 {"unit": item.offset_unit}) 1093 if written == True:1093 if written: 1094 1094 self.append(off, det) 1095 1095 … … 1103 1103 item.orientation.z, 1104 1104 {"unit": item.orientation_unit}) 1105 if written == True:1105 if written: 1106 1106 self.append(ori, det) 1107 1107 … … 1115 1115 item.beam_center.z, 1116 1116 {"unit": item.beam_center_unit}) 1117 if written == True:1117 if written: 1118 1118 self.append(center, det) 1119 1119 … … 1125 1125 written = written | self.write_node(pix, "z", item.pixel_size.z, 1126 1126 {"unit": item.pixel_size_unit}) 1127 if written == True:1127 if written: 1128 1128 self.append(pix, det) 1129 1129 self.write_node(det, "slit_length", item.slit_length, -
src/sas/sascalc/dataloader/readers/danse_reader.py
raf3e9f5 r2469df7 157 157 # Store all data 158 158 # Store wavelength 159 if has_converter == Trueand self.current_datainfo.source.wavelength_unit != 'A':159 if has_converter and self.current_datainfo.source.wavelength_unit != 'A': 160 160 conv = Converter('A') 161 161 wavelength = conv(wavelength, … … 164 164 165 165 # Store distance 166 if has_converter == Trueand detector.distance_unit != 'm':166 if has_converter and detector.distance_unit != 'm': 167 167 conv = Converter('m') 168 168 distance = conv(distance, units=detector.distance_unit) … … 170 170 171 171 # Store pixel size 172 if has_converter == Trueand detector.pixel_size_unit != 'mm':172 if has_converter and detector.pixel_size_unit != 'mm': 173 173 conv = Converter('mm') 174 174 pixel = conv(pixel, units=detector.pixel_size_unit) -
src/sas/sascalc/dataloader/readers/sesans_reader.py
r849094a r3053a4a 12 12 from ..file_reader_base_class import FileReader 13 13 from ..data_info import plottable_1D, DataInfo 14 from ..loader_exceptions import FileContentsException , DataReaderException14 from ..loader_exceptions import FileContentsException 15 15 16 16 # Check whether we have a converter available … … 18 18 try: 19 19 from sas.sascalc.data_util.nxsunit import Converter 20 except :20 except ImportError: 21 21 has_converter = False 22 22 _ZERO = 1e-16 … … 46 46 line = self.nextline() 47 47 params = {} 48 while not line.startswith("BEGIN_DATA"):48 while line and not line.startswith("BEGIN_DATA"): 49 49 terms = line.split() 50 50 if len(terms) >= 2: … … 63 63 raise FileContentsException("Wavelength has no units") 64 64 if params["SpinEchoLength_unit"] != params["Wavelength_unit"]: 65 raise FileContentsException("The spin echo data has rudely used " 66 "different units for the spin echo length " 67 "and the wavelength. While sasview could " 68 "handle this instance, it is a violation " 69 "of the file format and will not be " 70 "handled by other software.") 65 raise FileContentsException( 66 "The spin echo data has rudely used " 67 "different units for the spin echo length " 68 "and the wavelength. While sasview could " 69 "handle this instance, it is a violation " 70 "of the file format and will not be " 71 "handled by other software.") 71 72 72 73 headers = self.nextline().split() … … 86 87 87 88 if not data.size: 88 raise FileContentsException("{} is empty".format( path))89 raise FileContentsException("{} is empty".format(self.filepath)) 89 90 x = data[:, headers.index("SpinEchoLength")] 90 91 if "SpinEchoLength_error" in headers: -
src/sas/sascalc/fit/MultiplicationModel.py
r574adc7 r0a9cbc3 68 68 try: 69 69 multiplicity = p_model.multiplicity 70 except :70 except AttributeError: 71 71 multiplicity = 1 72 72 ## functional multiplicity of the model … … 76 76 self.non_fittable = p_model.non_fittable 77 77 self.multiplicity_info = [] 78 self.fun_list = {}78 self.fun_list = [] 79 79 if self.non_fittable > 1: 80 80 try: … … 82 82 self.fun_list = p_model.fun_list 83 83 self.is_multiplicity_model = True 84 except :84 except AttributeError: 85 85 pass 86 86 else: -
src/sas/sascalc/fit/qsmearing.py
r50fcb09 r2469df7 90 90 #print "data1D.dx[0]",data1D.dx[0],data1D.dxl[0] 91 91 # If we found resolution smearing data, return a QSmearer 92 if _found_resolution == True:92 if _found_resolution: 93 93 return pinhole_smear(data, model) 94 94 … … 113 113 break 114 114 # If we found slit smearing data, return a slit smearer 115 if _found_slit == True:115 if _found_slit: 116 116 return slit_smear(data, model) 117 117 return None -
src/sas/sascalc/pr/invertor.py
rd04ac05 r2469df7 222 222 elif name == 'est_bck': 223 223 value = self.get_est_bck() 224 if value == 1: 225 return True 226 else: 227 return False 224 return value == 1 228 225 elif name in self.__dict__: 229 226 return self.__dict__[name] … … 460 457 461 458 # If we need to fit the background, add a term 462 if self.est_bck == True:459 if self.est_bck: 463 460 nfunc_0 = nfunc 464 461 nfunc += 1 … … 506 503 507 504 # Keep a copy of the last output 508 if self.est_bck == False:505 if not self.est_bck: 509 506 self.out = c 510 507 self.cov = err … … 658 655 file.write("#slit_width=%g\n" % self.slit_width) 659 656 file.write("#background=%g\n" % self.background) 660 if self.est_bck == True:657 if self.est_bck: 661 658 file.write("#has_bck=1\n") 662 659 else: … … 738 735 elif line.startswith('#has_bck='): 739 736 toks = line.split('=') 740 if int(toks[1]) == 1: 741 self.est_bck = True 742 else: 743 self.est_bck = False 737 self.est_bck = int(toks[1]) == 1 744 738 745 739 # Now read in the parameters -
src/sas/sascalc/pr/num_term.py
ra1b8fee r2469df7 55 55 medi = 0 56 56 for i in range(dv): 57 if odd == True:57 if odd: 58 58 medi = osc[int(med)] 59 59 else: … … 98 98 new_osc3.append(self.osc_list[i]) 99 99 100 if flag9 == True:100 if flag9: 101 101 self.dataset = new_osc1 102 elif flag8 == True:102 elif flag8: 103 103 self.dataset = new_osc2 104 104 else: … … 141 141 div = len(nts) 142 142 tem = float(div) / 2.0 143 odd = self.is_odd(div) 144 if odd == True: 143 if self.is_odd(div): 145 144 nt = nts[int(tem)] 146 145 else: … … 148 147 return nt, self.alpha_list[nt - 10], self.mess_list[nt - 10] 149 148 except: 150 #TODO: check the logic above and make sure it doesn't 149 #TODO: check the logic above and make sure it doesn't 151 150 # rely on the try-except. 152 151 return self.nterm_min, self.invertor.alpha, '' -
src/sas/sasgui/guiframe/local_perspectives/plotting/Plotter1D.py
r7432acb r2469df7 827 827 on_Modify Plot Property_close 828 828 """ 829 if self.appD.okay_clicked == True:829 if self.appD.okay_clicked: 830 830 info = self.appD.get_current_values() 831 831 self.appearance_selected_plot.custom_color = \ -
src/sas/sasgui/perspectives/calculator/media/sas_calculator_help.rst
r5ed76f8 r1b67f3e 26 26 intensity from the particle is 27 27 28 .. image:: gen_i.png 28 .. math:: 29 30 I(\vec Q) = \frac{1}{V}\left| 31 \sum_j^N v_j \beta_j \exp(i\vec Q \cdot \vec r_j)\right|^2 29 32 30 33 Equation 1. … … 46 49 atomic structure (such as taken from a PDB file) to get the right normalization. 47 50 48 *NOTE! $\beta_j$displayed in the GUI may be incorrect but this will not51 *NOTE!* $\beta_j$ *displayed in the GUI may be incorrect but this will not 49 52 affect the scattering computation if the correction of the total volume V is made.* 50 53 … … 56 59 ^^^^^^^^^^^^^^^^^^^ 57 60 58 For magnetic scattering, only the magnetization component, $ M_\perp$,59 perpendicular to the scattering vector $ Q$ contributes to the magnetic61 For magnetic scattering, only the magnetization component, $\mathbf{M}_\perp$, 62 perpendicular to the scattering vector $\vec Q$ contributes to the magnetic 60 63 scattering length. 61 64 … … 64 67 The magnetic scattering length density is then 65 68 66 .. image:: dm_eq.png 69 .. math:: 70 71 \beta_M = \frac{\gamma r_0}{2 \mu_B}\sigma \cdot \mathbf{M}_\perp 72 = D_M\sigma \cdot \mathbf{M}_\perp 67 73 68 74 where the gyromagnetic ratio is $\gamma = -1.913$, $\mu_B$ is the Bohr … … 81 87 .. image:: gen_mag_pic.png 82 88 83 Now let us assume that the angles of the *Q* vector and the spin-axis (x')84 to the x-axis are $\phi$ and $\theta_\text{up}$ respectively (see above). Then,89 Now let us assume that the angles of the $\vec Q$ vector and the spin-axis ($x'$) 90 to the $x$-axis are $\phi$ and $\theta_\text{up}$ respectively (see above). Then, 85 91 depending upon the polarization (spin) state of neutrons, the scattering 86 92 length densities, including the nuclear scattering length density ($\beta_N$) … … 89 95 * for non-spin-flips 90 96 91 .. image:: sld1.png 97 .. math:: 98 \beta_{\pm\pm} = \beta_N \mp D_M M_{\perp x'} 92 99 93 100 * for spin-flips 94 101 95 .. image:: sld2.png 102 .. math:: 103 \beta_{\pm\mp} = - D_M(M_{\perp y'} \pm i M_{\perp z'}) 96 104 97 105 where 98 106 99 .. image:: mxp.png107 .. math:: 100 108 101 .. image:: myp.png 109 M_{\perp x'} &= M_{0q_x}\cos\theta_\text{up} + M_{0q_y}\sin\theta_\text{up} \\ 110 M_{\perp y'} &= M_{0q_y}\cos\theta_\text{up} - M_{0q_x}\sin\theta_\text{up} \\ 111 M_{\perp z'} &= M_{0z} \\ 112 M_{0q_x} &= (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\ 113 M_{0q_y} &= (M_{0y}\sin\phi - M_{0y}\cos\phi)\sin\phi 102 114 103 .. image:: mzp.png 104 105 .. image:: mqx.png 106 107 .. image:: mqy.png 108 109 Here the $M0_x$, $M0_y$ and $M0_z$ are the $x$, $y$ and $z$ 110 components of the magnetisation vector in the laboratory $xyz$ frame. 115 Here the $M_{0x}$, $M_{0y}$ and $M_{0z}$ are 116 the $x$, $y$ and $z$ components of the magnetisation vector in the 117 laboratory $x$-$y$-$z$ frame. 111 118 112 119 .. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ … … 148 155 uses the Debye equation below providing a 1D output 149 156 150 .. image:: gen_debye_eq.png 157 .. math:: 158 159 I(|\vec Q|) = \frac{1}{V}\sum_j^N v_j\beta_j \sum_k^N v_k \beta_k 160 \frac{\sin(|\vec Q||\vec r_j - \vec r_k|)}{|\vec Q||\vec r_j - \vec r_k|} 151 161 152 162 where $v_j \beta_j \equiv b_j$ is the scattering -
src/sas/sasgui/perspectives/calculator/model_editor.py
r1c206d9 r2469df7 332 332 msg = ("%s is not a valid Python name. Only alphanumeric \n" \ 333 333 "and underscore allowed" % self.name) 334 334 335 335 #Now check if the name already exists 336 336 if not self.overwrite_name and self.good_name: … … 344 344 msg = "Name exists already." 345 345 346 if self.good_name == False:346 if not self.good_name: 347 347 self.name_tcl.SetBackgroundColour('pink') 348 348 info = 'Error' -
src/sas/sasgui/perspectives/fitting/basepage.py
r3bd677b reee94bf 1643 1643 if item_page[2].__class__.__name__ == "ComboBox": 1644 1644 if item_page_info[2] in self.model.fun_list: 1645 fun_val = self.model.fun_list [item_page_info[2]]1645 fun_val = self.model.fun_list.index(item_page_info[2]) 1646 1646 self.model.setParam(item_page_info[1], fun_val) 1647 1647 if item_page[3] is not None: … … 1687 1687 selection = value 1688 1688 if value in self.model.fun_list: 1689 selection = self.model.fun_list [value]1689 selection = self.model.fun_list.index(value) 1690 1690 item_page[2].SetValue(selection) 1691 1691 self.model.setParam(param_name, selection) … … 3365 3365 if item[2].__class__.__name__ == "ComboBox": 3366 3366 if content[name][1] in self.model.fun_list: 3367 fun_val = self.model.fun_list [content[name][1]]3367 fun_val = self.model.fun_list.index(content[name][1]) 3368 3368 self.model.setParam(name, fun_val) 3369 3369 try: … … 3415 3415 if item[2].__class__.__name__ == "ComboBox": 3416 3416 if value[0] in self.model.fun_list: 3417 fun_val = self.model.fun_list [value[0]]3417 fun_val = self.model.fun_list.index(value[0]) 3418 3418 self.model.setParam(name, fun_val) 3419 3419 # save state -
src/sas/sasgui/perspectives/fitting/fitpage.py
r3bd677b rbfeb823 1726 1726 and not self.temp_multi_functional: 1727 1727 return None 1728 # Get the func name list 1729 list = self.model.fun_list 1730 if len(list) == 0: 1731 return None 1732 # build function (combo)box 1733 ind = 0 1734 while(ind < len(list)): 1735 for key, val in list.items(): 1736 if val == ind: 1737 fun_box.Append(key, val) 1738 break 1739 ind += 1 1728 for index, choice in enumerate(self.model.fun_list): 1729 fun_box.Append(choice, index) 1740 1730 1741 1731 def _on_select_accuracy(self, event): … … 1765 1755 value = fun_box.GetValue() 1766 1756 if value in self.model.fun_list: 1767 fun_val = self.model.fun_list [value]1757 fun_val = self.model.fun_list.index(value) 1768 1758 1769 1759 self.model.setParam(name, fun_val) -
src/sas/sasgui/perspectives/fitting/media/fitting.rst
r3bd677b rc926a97 18 18 19 19 Polarisation/Magnetic Scattering <magnetism/magnetism> 20 21 Oriented Particles <orientation/orientation> 20 22 21 23 Information on the SasView Optimisers <optimizer> -
src/sas/sasgui/perspectives/invariant/invariant_panel.py
r7432acb r2469df7 707 707 708 708 # reset power_out to default to get ready for another '_on_text' 709 if self.is_power_out == True:709 if self.is_power_out: 710 710 self.state.container = copy.deepcopy(self.inv_container) 711 711 self.state.timestamp = self._get_time_stamp() -
src/sas/sasgui/perspectives/invariant/invariant_state.py
r1fa4f736 r2469df7 655 655 : return: None 656 656 """ 657 if self.cansas == True:657 if self.cansas: 658 658 return self._read_cansas(path) 659 659 else: … … 763 763 """ 764 764 # Sanity check 765 if self.cansas == True:765 if self.cansas: 766 766 doc = self.write_toXML(datainfo, invstate) 767 767 # Write the XML document -
src/sas/sasgui/perspectives/pr/inversion_state.py
r1fa4f736 r2469df7 389 389 390 390 """ 391 if self.cansas == True:391 if self.cansas: 392 392 return self._read_cansas(path) 393 393 else: … … 505 505 """ 506 506 # Sanity check 507 if self.cansas == True:507 if self.cansas: 508 508 doc = self.write_toXML(datainfo, prstate) 509 509 # Write the XML document -
src/sas/sasgui/perspectives/pr/pr.py
rcb62bd5 r2469df7 407 407 y[i] = value 408 408 409 if self._normalize_output == True:409 if self._normalize_output: 410 410 y = y / total 411 411 dy = dy / total 412 elif self._scale_output_unity == True:412 elif self._scale_output_unity: 413 413 y = y / pmax 414 414 dy = dy / pmax … … 544 544 lines = buff.split('\n') 545 545 for line in lines: 546 if data_started == True:546 if data_started: 547 547 try: 548 548 toks = line.split() -
src/sas/sasgui/plottools/PlotPanel.py
ra1b8fee r2469df7 476 476 return 477 477 self.mousemotion = True 478 if self.leftdown == True and self.mousemotion == True:478 if self.leftdown and self.mousemotion: 479 479 ax = event.inaxes 480 480 if ax is not None: # the dragging is perform inside the figure -
src/sas/sasgui/plottools/fitDialog.py
r7432acb r2469df7 44 44 fitting and derives and displays specialized output parameters based 45 45 on the scale choice of the plot calling it. 46 46 47 47 :note1: The fitting is currently a bit convoluted as besides using 48 48 plottools.transform.py to handle all the conversions, it uses … … 55 55 This would considerably simplify the code and remove the need I think 56 56 for LineModel.py and possibly fittins.py altogether. -PDB 7/10/16 57 57 58 58 :note2: The linearized fits do not take resolution into account. This 59 59 means that for poor resolution such as slit smearing the answers will … … 142 142 """ 143 143 144 # set up sizers first. 144 # set up sizers first. 145 145 # vbox is the panel sizer and is a vertical sizer 146 146 # The first element of the panel is sizer which is a gridbagsizer … … 151 151 sizer = wx.GridBagSizer(5, 5) 152 152 sizer_button = wx.BoxSizer(wx.HORIZONTAL) 153 153 154 154 #size of string boxes in pixels 155 155 _BOX_WIDTH = 100 … … 395 395 sizer_button.Add(self.btClose, 0, 396 396 wx.LEFT | wx.RIGHT | wx.ADJUST_MINSIZE, 10) 397 397 398 398 vbox.Add(sizer) 399 self.static_line_1 = wx.StaticLine(self, -1) 399 self.static_line_1 = wx.StaticLine(self, -1) 400 400 vbox.Add(self.static_line_1, 0, wx.EXPAND, 0) 401 401 vbox.Add(sizer_button, 0, wx.EXPAND | wx.BOTTOM | wx.TOP, 10) … … 439 439 # makes transformation for y as a line to fit 440 440 if self.x != []: 441 if self.checkFitValues(self.xminFit) == True:441 if self.checkFitValues(self.xminFit): 442 442 # Check if the field of Fit Dialog contain values 443 443 # and use the x max and min of the user -
src/sas/sasgui/plottools/plottable_interactor.py
ra1b8fee r2469df7 166 166 from within the boundaries of an artist. 167 167 """ 168 if self._context_menu == True:168 if self._context_menu: 169 169 self._context_menu = False 170 170 evt.artist = self.marker … … 216 216 """ 217 217 if not evt.artist.__class__.__name__ == "AxesSubplot": 218 if self._context_menu == False:218 if not self._context_menu: 219 219 self.base.plottable_selected(None) 220 220 try: -
src/sas/sasgui/plottools/plottables.py
r2d9526d r2469df7 227 227 max_value = None 228 228 for p in self.plottables: 229 if p.hidden == True:229 if p.hidden: 230 230 continue 231 231 if p.x is not None: … … 1062 1062 Renders the plottable on the graph 1063 1063 """ 1064 if self.interactive == True:1064 if self.interactive: 1065 1065 kw['symbol'] = self.symbol 1066 1066 kw['id'] = self.id
Note: See TracChangeset
for help on using the changeset viewer.