eigenmath/mag.cpp

159 lines
2.0 KiB
C++
Raw Permalink Normal View History

2005-10-12 01:25:25 +02:00
/* Magnitude of complex z
z mag(z)
- ------
a a
-a a
(-1)^a 1
exp(a + i b) exp(a)
a b mag(a) mag(b)
a + i b sqrt(a^2 + b^2)
2006-08-09 19:57:45 +02:00
Notes
2005-10-12 01:25:25 +02:00
2006-08-09 19:57:45 +02:00
1. Handles mixed polar and rectangular forms, e.g. 1 + exp(i pi/3)
2. jean-francois.debroux reports that when z=(a+i*b)/(c+i*d) then
mag(numerator(z)) / mag(denominator(z))
must be used to get the correct answer. Now the operation is
automatic.
2005-10-12 01:25:25 +02:00
*/
#include "stdafx.h"
#include "defs.h"
void
eval_mag(void)
{
push(cadr(p1));
eval();
mag();
}
void
mag(void)
2006-08-09 19:57:45 +02:00
{
save();
p1 = pop();
push(p1);
numerator();
yymag();
push(p1);
denominator();
yymag();
divide();
restore();
}
void
yymag(void)
2005-10-12 01:25:25 +02:00
{
save();
p1 = pop();
if (isnegativenumber(p1)) {
push(p1);
negate();
} else if (car(p1) == symbol(POWER) && equaln(cadr(p1), -1))
// -1 to a power
push_integer(1);
else if (car(p1) == symbol(POWER) && cadr(p1) == symbol(E)) {
// exponential
push(caddr(p1));
real();
exponential();
} else if (car(p1) == symbol(MULTIPLY)) {
// product
push_integer(1);
p1 = cdr(p1);
while (iscons(p1)) {
push(car(p1));
mag();
multiply();
p1 = cdr(p1);
}
} else if (car(p1) == symbol(ADD)) {
// sum
push(p1);
rect(); // convert polar terms, if any
p1 = pop();
push(p1);
real();
push_integer(2);
power();
push(p1);
imag();
push_integer(2);
power();
add();
push_rational(1, 2);
power();
simplify_trig();
} else
// default (all real)
push(p1);
restore();
}
2007-05-08 16:57:30 +02:00
#if SELFTEST
2005-10-09 23:24:06 +02:00
static char *s[] = {
2005-10-08 04:20:34 +02:00
2005-10-09 23:24:06 +02:00
"mag(a+i*b)",
"(a^2+b^2)^(1/2)",
2005-10-08 04:20:34 +02:00
2005-10-09 23:24:06 +02:00
"mag(exp(a+i*b))",
"exp(a)",
2005-10-08 04:20:34 +02:00
2005-10-09 23:24:06 +02:00
"mag(1)",
"1",
2005-10-08 04:20:34 +02:00
2005-10-09 23:24:06 +02:00
"mag(-1)",
"1",
2005-10-08 04:20:34 +02:00
2005-10-09 23:24:06 +02:00
"mag(1+exp(i*pi/3))",
"3^(1/2)",
2006-08-09 19:57:45 +02:00
"mag((a+i*b)/(c+i*d))",
"(a^2+b^2)^(1/2)/((c^2+d^2)^(1/2))",
"mag(exp(i theta))",
"1",
"mag(exp(-i theta))",
"1",
"mag((-1)^theta)",
"1",
"mag((-1)^(-theta))",
"1",
"mag(3*(-1)^theta)",
"3",
"mag(3*(-1)^(-theta))",
"3",
"mag(-3*(-1)^theta)",
"3",
"mag(-3*(-1)^(-theta))",
"3",
2005-10-09 23:24:06 +02:00
};
2005-10-08 04:20:34 +02:00
void
2005-10-09 23:24:06 +02:00
test_mag(void)
2005-10-08 04:20:34 +02:00
{
2005-10-09 23:24:06 +02:00
test(__FILE__, s, sizeof s / sizeof (char *));
2005-10-08 04:20:34 +02:00
}
2007-05-08 16:57:30 +02:00
#endif