Changeset f54e82cf in sasview for src/sas/sascalc/calculator/c_extensions/librefl.c
- Timestamp:
- Jul 31, 2018 2:36:52 AM (6 years ago)
- Branches:
- ESS_GUI, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc
- Children:
- d9b7197
- Parents:
- 3067196
- git-author:
- Torin Cooper-Bennun <torin.cooper-bennun@…> (07/30/18 08:35:02)
- git-committer:
- Torin Cooper-Bennun <torin.cooper-bennun@…> (07/31/18 02:36:52)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/sas/sascalc/calculator/c_extensions/librefl.c
r4c29e4d rf54e82cf 103 103 #endif // NEED_ERF 104 104 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; 105 void cassign(Cplx *x, double real, double imag) 106 { 107 x->re = real; 108 x->im = imag; 109 } 110 111 112 void cplx_add(Cplx *z, Cplx x, Cplx y) 113 { 114 z->re = x.re + y.re; 115 z->im = x.im + y.im; 116 } 117 118 void rcmult(Cplx *z, double x, Cplx y) 119 { 120 z->re = x*y.re; 121 z->im = x*y.im; 122 } 123 124 void cplx_sub(Cplx *z, Cplx x, Cplx y) 125 { 126 z->re = x.re - y.re; 127 z->im = x.im - y.im; 128 } 129 130 131 void cplx_mult(Cplx *z, Cplx x, Cplx y) 132 { 133 z->re = x.re*y.re - x.im*y.im; 134 z->im = x.re*y.im + x.im*y.re; 135 } 136 137 void cplx_div(Cplx *z, Cplx x, Cplx y) 138 { 139 z->re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im); 140 z->im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im); 141 } 142 143 void cplx_exp(Cplx *z, Cplx b) 144 { 166 145 double br,bi; 167 146 br=b.re; 168 147 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; 148 z->re = exp(br)*cos(bi); 149 z->im = exp(br)*sin(bi); 150 } 151 152 153 void cplx_sqrt(Cplx *c, Cplx z) //see Schaum`s Math Handbook p. 22, 6.6 and 6.10 154 { 179 155 double zr,zi,x,y,r,w; 180 156 … … 184 160 if (zr==0.0 && zi==0.0) 185 161 { 186 c.re=0.0; 187 c.im=0.0; 188 return c; 189 } 190 else 191 { 162 c->re=0.0; 163 c->im=0.0; 164 } else { 192 165 x=fabs(zr); 193 166 y=fabs(zi); … … 196 169 r=y/x; 197 170 w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r))); 198 } 199 else 200 { 171 } else { 201 172 r=x/y; 202 173 w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r))); … … 204 175 if (zr >=0.0) 205 176 { 206 c.re=w; 207 c.im=zi/(2.0*w); 177 c->re=w; 178 c->im=zi/(2.0*w); 179 } else { 180 c->im=(zi >= 0) ? w : -w; 181 c->re=zi/(2.0*c->im); 208 182 } 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; 183 } 184 } 185 186 void cplx_cos(Cplx *z, Cplx b) 187 { 188 // cos(b) = (e^bi + e^-bi)/2 189 // = (e^b.im e^-i bi.re) + e^-b.im e^i b.re)/2 190 // = (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 191 // = 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 192 // = (e^b.im + 1/e^b.im)/2 cos(b.re) + (-e^b.im + 1/e^b.im)/2 sin(b.re) i 193 // = cosh(b.im) cos(b.re) - sinh(b.im) sin(b.re) i 194 double exp_b_im = exp(b.im); 195 z->re = 0.5*(+exp_b_im + 1.0/exp_b_im) * cos(b.re); 196 z->im = -0.5*(exp_b_im - 1.0/exp_b_im) * sin(b.re); 229 197 } 230 198
Note: See TracChangeset
for help on using the changeset viewer.