Changes in / [f0ce6e2:aab23e3] in sasview
- Location:
- src/sas/sascalc/calculator
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
src/sas/sascalc/calculator/c_extensions/libfunc.c
r144e032a rb6c8abe 95 95 // spintheta: angle (anti-clock-wise) between neutron spin(up) and x axis 96 96 // Note: all angles are in degrees. 97 void cal_msld(polar_sld *p_sld,int isangle, double qx, double qy, double bn,97 polar_sld cal_msld(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 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; 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; 132 133 133 134 //No mag means no further calculation 134 if (isangle>0) 135 if (isangle>0){ 135 136 if (m_max < 1.0e-32){ 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; 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; 213 213 } 214 214 -
src/sas/sascalc/calculator/c_extensions/libfunc.h
r144e032a rb6c8abe 18 18 //double gamln(double x); 19 19 20 void cal_msld(polar_sld*,int isangle, double qx, double qy, double bn, double m01, double mtheta1,20 polar_sld cal_msld(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
r144e032a ra1daf86 112 112 #endif // NEED_ERF 113 113 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 { 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; 154 175 double br,bi; 155 176 br=b.re; 156 177 bi=b.im; 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 { 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; 164 188 double zr,zi,x,y,r,w; 165 189 … … 169 193 if (zr==0.0 && zi==0.0) 170 194 { 171 c->re=0.0; 172 c->im=0.0; 173 } else { 195 c.re=0.0; 196 c.im=0.0; 197 return c; 198 } 199 else 200 { 174 201 x=fabs(zr); 175 202 y=fabs(zi); … … 178 205 r=y/x; 179 206 w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))); 180 } else { 207 } 208 else 209 { 181 210 r=x/y; 182 211 w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))); … … 184 213 if (zr >=0.0) 185 214 { 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); 215 c.re=w; 216 c.im=zi/(2.0*w); 191 217 } 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); 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; 206 238 } 207 239 -
src/sas/sascalc/calculator/c_extensions/librefl.h
r144e032a r9e531f2 5 5 double re; 6 6 double im; 7 } Cplx;7 } complex; 8 8 9 9 typedef struct { 10 Cplx a;11 Cplx b;12 Cplx c;13 Cplx d;10 complex a; 11 complex b; 12 complex c; 13 complex d; 14 14 } matrix; 15 15 16 void cassign(Cplx*,double real, double imag);16 complex cassign(double real, double imag); 17 17 18 void cplx_add(Cplx*, Cplx x,Cplx y);18 complex cplx_add(complex x,complex y); 19 19 20 void rcmult(Cplx*, double x,Cplx y);20 complex rcmult(double x,complex y); 21 21 22 void cplx_sub(Cplx*, Cplx x,Cplx y);22 complex cplx_sub(complex x,complex y); 23 23 24 void cplx_mult(Cplx*, Cplx x,Cplx y);24 complex cplx_mult(complex x,complex y); 25 25 26 void cplx_div(Cplx*, Cplx x,Cplx y);26 complex cplx_div(complex x,complex y); 27 27 28 void cplx_exp(Cplx*, Cplx b);28 complex cplx_exp(complex b); 29 29 30 void cplx_sqrt(Cplx*, Cplx z);30 complex cplx_sqrt(complex z); 31 31 32 void cplx_cos(Cplx*, Cplx b);32 complex cplx_cos(complex 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
r144e032a ra1daf86 54 54 polar_sld b_sld; 55 55 double qr = 0.0; 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;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; 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);72 68 73 69 //Assume that pixel volumes are given in vol_pix in A^3 unit … … 81 77 for(i=0; i<npoints; i++){ 82 78 //I_out[i] = 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);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); 87 83 //printf("i: %d\n", i); 88 84 //q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); // + qz[i]*qz[i]); 89 85 90 86 for(j=0; j<this->n_pix; j++){ 87 //printf("j: %d\n", j); 91 88 if (this->sldn_val[j]!=0.0 92 89 ||this->mx_val[j]!=0.0 … … 94 91 ||this->mz_val[j]!=0.0) 95 92 { 96 // printf("i,j: %d,%d\n", i,j);97 93 //anisotropic 98 cassign(&temp_fi,0.0, 0.0);99 cal_msld(&b_sld,0, qx[i], qy[i], this->sldn_val[j],94 temp_fi = cassign(0.0, 0.0); 95 b_sld = cal_msld(0, qx[i], qy[i], this->sldn_val[j], 100 96 this->mx_val[j], this->my_val[j], this->mz_val[j], 101 97 this->inspin, this->outspin, this->stheta); 102 98 qr = (qx[i]*this->x_val[j] + qy[i]*this->y_val[j]); 103 cassign(&iqr,0.0, qr);104 cplx_exp(&ephase,iqr);99 iqr = cassign(0.0, qr); 100 ephase = cplx_exp(iqr); 105 101 106 102 //Let's multiply pixel(atomic) volume here 107 rcmult(&ephase,this->vol_pix[j], ephase);103 ephase = rcmult(this->vol_pix[j], ephase); 108 104 //up_up 109 105 if (this->inspin > 0.0 && this->outspin > 0.0){ 110 c assign(&comp_sld,b_sld.uu, 0.0);111 cplx_mult(&temp_fi,comp_sld, ephase);112 cplx_add(&sumj_uu,sumj_uu, temp_fi);106 comp_sld = cassign(b_sld.uu, 0.0); 107 temp_fi = cplx_mult(comp_sld, ephase); 108 sumj_uu = cplx_add(sumj_uu, temp_fi); 113 109 } 114 110 //down_down 115 111 if (this->inspin < 1.0 && this->outspin < 1.0){ 116 c assign(&comp_sld,b_sld.dd, 0.0);117 cplx_mult(&temp_fi,comp_sld, ephase);118 cplx_add(&sumj_dd,sumj_dd, temp_fi);112 comp_sld = cassign(b_sld.dd, 0.0); 113 temp_fi = cplx_mult(comp_sld, ephase); 114 sumj_dd = cplx_add(sumj_dd, temp_fi); 119 115 } 120 116 //up_down 121 117 if (this->inspin > 0.0 && this->outspin < 1.0){ 122 c assign(&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);118 comp_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); 125 121 } 126 122 //down_up 127 123 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 } 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 132 129 133 130 if (i == 0){ -
src/sas/sascalc/calculator/sas_gen.py
r144e032a ra1daf86 133 133 sldn = copy.deepcopy(self.data_sldn) 134 134 sldn -= self.params['solvent_SLD'] 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) 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']) 148 142 if len(qy): 149 143 qx, qy = _vec(qx), _vec(qy)
Note: See TracChangeset
for help on using the changeset viewer.