lib/complex.c

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