lib/complex.c

   1 /* complex - Complex number functions                                  */
   2 /* Copyright (c) 1990, by Luke Tierney                                 */
   3
   4 #include "xmath.h"
   5 #include "complex.h"
   6
   7 static double phase(c)
   8         Complex c;
   9 {
  10   double phi;
  11
  12   if (c.real == 0.0) {
  13     if (c.imag > 0.0) phi = PI / 2;
  14     else if (c.imag == 0.0) phi = 0.0;
  15     else phi = -PI / 2;
  16   }
  17   else {
  18     phi = atan(c.imag / c.real);
  19     if (c.real < 0.0) {
  20       if (c.imag > 0.0) phi += PI;
  21       else if (c.imag < 0.0) phi -= PI;
  22       else phi = PI;
  23     }
  24   }
  25   return(phi);
  26 }
  27
  28 double modulus(c)
  29         Complex c;
  30 {
  31   return(sqrt(c.real * c.real + c.imag * c.imag));
  32 }
  33
  34 Complex cart2complex(real, imag)
  35         double real, imag;
  36 {
  37   Complex val;
  38   val.real = real;
  39   val.imag = imag;
  40   return(val);
  41 }
  42
  43 static Complex polar2complex(mod, phi)
  44         double mod, phi;
  45 {
  46   Complex val;
  47   double cs, sn;
  48
  49   if (phi == 0) {
  50     cs = 1.0;
  51     sn = 0.0;
  52   }
  53   else if (phi == PI / 2) {
  54     cs = 0.0;
  55     sn = 1.0;
  56   }
  57   else if (phi == PI) {
  58     cs = -1.0;
  59     sn = 0.0;
  60   }
  61   else if (phi == -PI / 2) {
  62     cs = 0.0;
  63     sn = -1.0;
  64   }
  65   else {
  66     cs = cos(phi);
  67     sn = sin(phi);
  68   }
  69   val.real = mod * cs;
  70   val.imag = mod * sn;
  71   return(val);
  72 }
  73
  74 Complex cadd(c1, c2)
  75         Complex c1, c2;
  76 {
  77   return(cart2complex(c1.real + c2.real, c1.imag + c2.imag));
  78 }
  79
  80 Complex csub(c1, c2)
  81         Complex c1, c2;
  82 {
  83   return(cart2complex(c1.real - c2.real, c1.imag - c2.imag));
  84 }
  85
  86 Complex cmul(c1, c2)
  87         Complex c1, c2;
  88 {
  89   double m1, m2, p1, p2;
  90
  91   m1 = modulus(c1);
  92   p1 = phase(c1);
  93   m2 = modulus(c2);
  94   p2 = phase(c2);
  95   return(polar2complex(m1 * m2, p1 + p2));
  96 }
  97
  98 Complex cdiv(c1, c2)
  99         Complex c1, c2;
 100 {
 101   double m1, m2, p1, p2;
 102
 103   m1 = modulus(c1);
 104   p1 = phase(c1);
 105   m2 = modulus(c2);
 106   p2 = phase(c2);
 107   checkfzero(m2);
 108   return(polar2complex(m1 / m2, p1 - p2));
 109 }
 110
 111 static checkfzero(x)
 112      double x;
 113 {
 114   if (x == 0.0) xlfail("division by zero");
 115 }