Changeset 144e032a in sasview for src/sas/sascalc/calculator/c_extensions/librefl.c

Timestamp:

Nov 14, 2017 10:57:07 AM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

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

Message:

tinycc doesn't return structures, so must pass return structure as pointer

File:

• src/sas/sascalc/calculator/c_extensions/librefl.c (modified) (4 diffs)

src/sas/sascalc/calculator/c_extensions/librefl.c

@@ -112 +112 @@
 #endif // NEED_ERF
 
-complex cassign(real, imag)
-        double real, imag;
-{
-        complex x;
-        x.re = real;
-        x.im = imag;
-        return x;
-}
-
-
-complex cplx_add(x,y)
-        complex x,y;
-{
-        complex z;
-        z.re = x.re + y.re;
-        z.im = x.im + y.im;
-        return z;
-}
-
-complex rcmult(x,y)
-        double x;
-    complex y;
-{
-        complex z;
-        z.re = x*y.re;
-        z.im = x*y.im;
-        return z;
-}
-
-complex cplx_sub(x,y)
-        complex x,y;
-{
-        complex z;
-        z.re = x.re - y.re;
-        z.im = x.im - y.im;
-        return z;
-}
-
-
-complex cplx_mult(x,y)
-        complex x,y;
-{
-        complex z;
-        z.re = x.re*y.re - x.im*y.im;
-        z.im = x.re*y.im + x.im*y.re;
-        return z;
-}
-
-complex cplx_div(x,y)
-        complex x,y;
-{
-        complex z;
-        z.re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);
-        z.im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);
-        return z;
-}
-
-complex cplx_exp(b)
-        complex b;
-{
-        complex z;
+void cassign(Cplx *x, double real, double imag)
+{
+        x->re = real;
+        x->im = imag;
+}
+
+
+void cplx_add(Cplx *z, Cplx x, Cplx y)
+{
+        z->re = x.re + y.re;
+        z->im = x.im + y.im;
+}
+
+void rcmult(Cplx *z, double x, Cplx y)
+{
+        z->re = x*y.re;
+        z->im = x*y.im;
+}
+
+void cplx_sub(Cplx *z, Cplx x, Cplx y)
+{
+        z->re = x.re - y.re;
+        z->im = x.im - y.im;
+}
+
+
+void cplx_mult(Cplx *z, Cplx x, Cplx y)
+{
+        z->re = x.re*y.re - x.im*y.im;
+        z->im = x.re*y.im + x.im*y.re;
+}
+
+void cplx_div(Cplx *z, Cplx x, Cplx y)
+{
+        z->re = (x.re*y.re + x.im*y.im)/(y.re*y.re + y.im*y.im);
+        z->im = (x.im*y.re - x.re*y.im)/(y.re*y.re + y.im*y.im);
+}
+
+void cplx_exp(Cplx *z, Cplx b)
+{
         double br,bi;
         br=b.re;
         bi=b.im;
-        z.re = exp(br)*cos(bi);
-        z.im = exp(br)*sin(bi);
-        return z;
-}
-
-
-complex cplx_sqrt(z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10
-        complex z;
-{
-        complex c;
+        z->re = exp(br)*cos(bi);
+        z->im = exp(br)*sin(bi);
+}
+
+
+void cplx_sqrt(Cplx *c, Cplx z)    //see Schaum`s Math Handbook p. 22, 6.6 and 6.10
+{
         double zr,zi,x,y,r,w;
 
@@ -193 +169 @@
         if (zr==0.0 && zi==0.0)
         {
-    c.re=0.0;
-        c.im=0.0;
-    return c;
-        }
-        else
-        {
+                c->re=0.0;
+                c->im=0.0;
+        } else {
                 x=fabs(zr);
                 y=fabs(zi);
@@ -205 +178 @@
                         r=y/x;
                         w=sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
-                }
-                else
-                {
+                } else {
                         r=x/y;
                         w=sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
@@ -213 +184 @@
                 if (zr >=0.0)
                 {
-                        c.re=w;
-                        c.im=zi/(2.0*w);
+                        c->re=w;
+                        c->im=zi/(2.0*w);
+                } else {
+                        c->im=(zi >= 0) ? w : -w;
+                        c->re=zi/(2.0*c->im);
                 }
-                else
-                {
-                        c.im=(zi >= 0) ? w : -w;
-                        c.re=zi/(2.0*c.im);
-                }
-                return c;
-        }
-}
-
-complex cplx_cos(b)
-        complex b;
-{
-        complex zero,two,z,i,bi,negbi;
-        zero = cassign(0.0,0.0);
-        two = cassign(2.0,0.0);
-        i = cassign(0.0,1.0);
-        bi = cplx_mult(b,i);
-        negbi = cplx_sub(zero,bi);
-        z = cplx_div(cplx_add(cplx_exp(bi),cplx_exp(negbi)),two);
-        return z;
+        }
+}
+
+void cplx_cos(Cplx *z, Cplx b)
+{
+        // cos(b) = (e^bi + e^-bi)/2
+        //        = (e^b.im e^-i bi.re) + e^-b.im e^i b.re)/2
+        //        = (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
+        //        = 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
+        //        = (e^b.im + 1/e^b.im)/2 cos(b.re) + (-e^b.im + 1/e^b.im)/2 sin(b.re) i
+        //        = cosh(b.im) cos(b.re) - sinh(b.im) sin(b.re) i
+        double exp_b_im = exp(b.im);
+        z->re = 0.5*(+exp_b_im + 1.0/exp_b_im) * cos(b.re);
+        z->im = -0.5*(exp_b_im - 1.0/exp_b_im) * sin(b.re);
 }