mag.cpp

   1 /* Magnitude of complex z
   2
   3         z               mag(z)
   4         -               ------
   5
   6         a               a
   7
   8         -a              a
   9
  10         (-1)^a          1
  11
  12         exp(a + i b)    exp(a)
  13
  14         a b             mag(a) mag(b)
  15
  16         a + i b         sqrt(a^2 + b^2)
  17
  18 Notes
  19
  20         1. Handles mixed polar and rectangular forms, e.g. 1 + exp(i pi/3)
  21
  22         2. jean-francois.debroux reports that when z=(a+i*b)/(c+i*d) then
  23
  24                 mag(numerator(z)) / mag(denominator(z))
  25
  26            must be used to get the correct answer. Now the operation is
  27            automatic.
  28 */
  29
  30 #include "stdafx.h"
  31 #include "defs.h"
  32
  33 void
  34 eval_mag(void)
  35 {
  36         push(cadr(p1));
  37         eval();
  38         mag();
  39 }
  40
  41 void
  42 mag(void)
  43 {
  44         save();
  45         p1 = pop();
  46         push(p1);
  47         numerator();
  48         yymag();
  49         push(p1);
  50         denominator();
  51         yymag();
  52         divide();
  53         restore();
  54 }
  55
  56 void
  57 yymag(void)
  58 {
  59         save();
  60         p1 = pop();
  61         if (isnegativenumber(p1)) {
  62                 push(p1);
  63                 negate();
  64         } else if (car(p1) == symbol(POWER) && equaln(cadr(p1), -1))
  65                 // -1 to a power
  66                 push_integer(1);
  67         else if (car(p1) == symbol(POWER) && cadr(p1) == symbol(E)) {
  68                 // exponential
  69                 push(caddr(p1));
  70                 real();
  71                 exponential();
  72         } else if (car(p1) == symbol(MULTIPLY)) {
  73                 // product
  74                 push_integer(1);
  75                 p1 = cdr(p1);
  76                 while (iscons(p1)) {
  77                         push(car(p1));
  78                         mag();
  79                         multiply();
  80                         p1 = cdr(p1);
  81                 }
  82         } else if (car(p1) == symbol(ADD)) {
  83                 // sum
  84                 push(p1);
  85                 rect(); // convert polar terms, if any
  86                 p1 = pop();
  87                 push(p1);
  88                 real();
  89                 push_integer(2);
  90                 power();
  91                 push(p1);
  92                 imag();
  93                 push_integer(2);
  94                 power();
  95                 add();
  96                 push_rational(1, 2);
  97                 power();
  98                 simplify_trig();
  99         } else
 100                 // default (all real)
 101                 push(p1);
 102         restore();
 103 }
 104
 105 #if SELFTEST
 106
 107 static char *s[] = {
 108
 109         "mag(a+i*b)",
 110         "(a^2+b^2)^(1/2)",
 111
 112         "mag(exp(a+i*b))",
 113         "exp(a)",
 114
 115         "mag(1)",
 116         "1",
 117
 118         "mag(-1)",
 119         "1",
 120
 121         "mag(1+exp(i*pi/3))",
 122         "3^(1/2)",
 123
 124         "mag((a+i*b)/(c+i*d))",
 125         "(a^2+b^2)^(1/2)/((c^2+d^2)^(1/2))",
 126
 127         "mag(exp(i theta))",
 128         "1",
 129
 130         "mag(exp(-i theta))",
 131         "1",
 132
 133         "mag((-1)^theta)",
 134         "1",
 135
 136         "mag((-1)^(-theta))",
 137         "1",
 138
 139         "mag(3*(-1)^theta)",
 140         "3",
 141
 142         "mag(3*(-1)^(-theta))",
 143         "3",
 144
 145         "mag(-3*(-1)^theta)",
 146         "3",
 147
 148         "mag(-3*(-1)^(-theta))",
 149         "3",
 150 };
 151
 152 void
 153 test_mag(void)
 154 {
 155         test(__FILE__, s, sizeof s / sizeof (char *));
 156 }
 157
 158 #endif