Changeset 144e032a in sasview for src/sas/sascalc
- Timestamp:
- Nov 14, 2017 10:57:07 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:
- f0ce6e2
- Parents:
- a57ae07
- Location:
- src/sas/sascalc/calculator
- Files:
-
- 6 edited
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
ra1daf86 r144e032a 112 112 #endif // NEED_ERF 113 113 114 complex cassign(real, imag) 115 double real, imag; 116 { 117 complex x; 118 x.re = real; 119 x.im = imag; 120 return x; 121 } 122 123 124 complex cplx_add(x,y) 125 complex x,y; 126 { 127 complex z; 128 z.re = x.re + y.re; 129 z.im = x.im + y.im; 130 return z; 131 } 132 133 complex rcmult(x,y) 134 double x; 135 complex y; 136 { 137 complex z; 138 z.re = x*y.re; 139 z.im = x*y.im; 140 return z; 141 } 142 143 complex cplx_sub(x,y) 144 complex x,y; 145 { 146 complex z; 147 z.re = x.re - y.re; 148 z.im = x.im - y.im; 149 return z; 150 } 151 152 153 complex cplx_mult(x,y) 154 complex x,y; 155 { 156 complex z; 157 z.re = x.re*y.re - x.im*y.im; 158 z.im = x.re*y.im + x.im*y.re; 159 return z; 160 } 161 162 complex cplx_div(x,y) 163 complex x,y; 164 { 165 complex z; 166 z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im); 167 z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im); 168 return z; 169 } 170 171 complex cplx_exp(b) 172 complex b; 173 { 174 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 { 175 154 double br,bi; 176 155 br=b.re; 177 156 bi=b.im; 178 z.re = exp(br)*cos(bi); 179 z.im = exp(br)*sin(bi); 180 return z; 181 } 182 183 184 complex cplx_sqrt(z) //see Schaum`s Math Handbook p. 22, 6.6 and 6.10 185 complex z; 186 { 187 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 { 188 164 double zr,zi,x,y,r,w; 189 165 … … 193 169 if (zr==0.0 && zi==0.0) 194 170 { 195 c.re=0.0; 196 c.im=0.0; 197 return c; 198 } 199 else 200 { 171 c->re=0.0; 172 c->im=0.0; 173 } else { 201 174 x=fabs(zr); 202 175 y=fabs(zi); … … 205 178 r=y/x; 206 179 w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))); 207 } 208 else 209 { 180 } else { 210 181 r=x/y; 211 182 w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))); … … 213 184 if (zr >=0.0) 214 185 { 215 c.re=w; 216 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); 217 191 } 218 else 219 { 220 c.im=(zi >= 0) ? w : -w; 221 c.re=zi/(2.0*c.im); 222 } 223 return c; 224 } 225 } 226 227 complex cplx_cos(b) 228 complex b; 229 { 230 complex zero,two,z,i,bi,negbi; 231 zero = cassign(0.0,0.0); 232 two = cassign(2.0,0.0); 233 i = cassign(0.0,1.0); 234 bi = cplx_mult(b,i); 235 negbi = cplx_sub(zero,bi); 236 z = cplx_div(cplx_add(cplx_exp(bi),cplx_exp(negbi)),two); 237 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); 238 206 } 239 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
ra1daf86 r144e032a 54 54 polar_sld b_sld; 55 55 double qr = 0.0; 56 complex iqr = cassign(0.0, 0.0);57 complex ephase = cassign(0.0, 0.0);58 complex comp_sld = cassign(0.0, 0.0);59 60 complex sumj_uu;61 complex sumj_ud;62 complex sumj_du;63 complex sumj_dd;64 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; 65 65 66 66 double count = 0.0; 67 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); 68 72 69 73 //Assume that pixel volumes are given in vol_pix in A^3 unit … … 77 81 for(i=0; i<npoints; i++){ 78 82 //I_out[i] = 0.0; 79 sumj_uu = cassign(0.0, 0.0);80 sumj_ud = cassign(0.0, 0.0);81 sumj_du = cassign(0.0, 0.0);82 sumj_dd = cassign(0.0, 0.0);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); 83 87 //printf("i: %d\n", i); 84 88 //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); 85 89 86 90 for(j=0; j<this->n_pix; j++){ 87 //printf("j: %d\n", j);88 91 if (this->sldn_val[j]!=0.0 89 92 ||this->mx_val[j]!=0.0 … … 91 94 ||this->mz_val[j]!=0.0) 92 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], this->sldn_val[j],98 cassign(&temp_fi, 0.0, 0.0); 99 cal_msld(&b_sld, 0, qx[i], qy[i], this->sldn_val[j], 96 100 this->mx_val[j], this->my_val[j], this->mz_val[j], 97 101 this->inspin, this->outspin, this->stheta); 98 102 qr = (qx[i]*this->x_val[j] + qy[i]*this->y_val[j]); 99 iqr = cassign(0.0, qr);100 ephase = cplx_exp(iqr);103 cassign(&iqr, 0.0, qr); 104 cplx_exp(&ephase, iqr); 101 105 102 106 //Let's multiply pixel(atomic) volume here 103 ephase = rcmult(this->vol_pix[j], ephase);107 rcmult(&ephase, this->vol_pix[j], ephase); 104 108 //up_up 105 109 if (this->inspin > 0.0 && this->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);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 115 if (this->inspin < 1.0 && this->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);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 121 if (this->inspin > 0.0 && this->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);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 127 if (this->inspin < 1.0 && this->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 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){ -
src/sas/sascalc/calculator/sas_gen.py
ra1daf86 r144e032a 133 133 sldn = copy.deepcopy(self.data_sldn) 134 134 sldn -= self.params['solvent_SLD'] 135 model = mod.new_GenI((1 if self.is_avg else 0), 136 pos_x, pos_y, pos_z, 137 sldn, self.data_mx, self.data_my, 138 self.data_mz, self.data_vol, 139 self.params['Up_frac_in'], 140 self.params['Up_frac_out'], 141 self.params['Up_theta']) 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) 142 148 if len(qy): 143 149 qx, qy = _vec(qx), _vec(qy)
Note: See TracChangeset
for help on using the changeset viewer.