home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume1 / 8707 / 75 / complex.c next >
Encoding:
C/C++ Source or Header  |  1990-07-13  |  5.6 KB  |  293 lines
  1. /*
  2.  *    Complex number library translated from Pascal source listing
  3.  *    in Dr. Dobbs October, 1984 "Simple Calculations with Complex Numbers"
  4.  *    by David Clark.
  5.  *
  6.  *    To insure that the same structure could be used as both an 'argument'
  7.  *    and 'result' immediate tranfer of structure members is done to a
  8.  *    local variable - in many of the functions this may not be necessary,
  9.  *    but it was done in all of the following functions for consistency
  10.  *
  11.  *    Author: R. Hellman
  12.  *    Date:    02/21/86
  13.  *
  14.  *    Released to the Public Domain
  15.  */
  16.  
  17. #include <malloc.h>
  18. #include <math.h>
  19.  
  20. #define TINY    2.2e-308    /* tiny value */
  21. #define LOGHUGE    709.778        /* natural log of huge value */
  22.  
  23. typedef    struct  {
  24.     double    re;
  25.     double    im;
  26. } complex;
  27.  
  28. /*
  29.  *    Initializes storage for each address of an array of complex pointers
  30.  *    *f[] with 'arraysz' number of elements
  31.  */
  32.  
  33. initcmplx(f, arraysz)
  34. complex    *f[];
  35. int    arraysz;
  36. {
  37.     int    i;
  38.  
  39.     for (i = 0; i < arraysz; i++)
  40.         if ((f[i] = (complex *)malloc(sizeof(complex))) == NULL)  {
  41.             puts("Error: complex number allocation\n");
  42.             exit(1);
  43.         }
  44. }
  45.  
  46. /* adds arg1 and arg2 and returns sum in result */
  47.  
  48. cadd(result, arg1, arg2)
  49. complex    *result, *arg1, *arg2;
  50. {
  51.     complex    targ1, targ2;
  52.  
  53.     targ1.re = arg1->re;
  54.     targ1.im = arg1->im;
  55.     targ2.re = arg2->re;
  56.     targ2.im = arg2->im;
  57.  
  58.     result->re = targ1.re + targ2.re;
  59.     result->im = targ1.im + targ2.im;
  60. }
  61.  
  62. /* subtracts arg1 and arg2 and returns sum in result */
  63.  
  64. csub(result, arg1, arg2)
  65. complex    *result, *arg1, *arg2;
  66. {
  67.     complex    targ1, targ2;
  68.  
  69.     targ1.re = arg1->re;
  70.     targ1.im = arg1->im;
  71.     targ2.re = arg2->re;
  72.     targ2.im = arg2->im;
  73.  
  74.     result->re = targ1.re - targ2.re;
  75.     result->im = targ1.im - targ2.im;
  76. }
  77.  
  78. /* multiplies arg1 and arg2 and returns product in result */
  79.  
  80. cmult(result, arg1, arg2)
  81. complex    *result, *arg1, *arg2;
  82. {
  83.     complex    targ1, targ2;
  84.  
  85.     targ1.re = arg1->re;
  86.     targ1.im = arg1->im;
  87.     targ2.re = arg2->re;
  88.     targ2.im = arg2->im;
  89.  
  90.     result->re = (targ1.re * targ2.re) - (targ1.im * targ2.im);
  91.     result->im = (targ1.im * targ2.re) + (targ1.re * targ2.im);
  92. }
  93.  
  94. /* divides arg1 and arg2 and returns quotient in result */
  95.  
  96. cdiv(result, arg1, arg2)
  97. complex    *result, *arg1, *arg2;
  98. {
  99.     double denom;
  100.     complex targ1,targ2;
  101.  
  102.     targ1.re = arg1->re;
  103.     targ1.im = arg1->im;
  104.     targ2.re = arg2->re;
  105.     targ2.im = arg2->im;
  106.     denom = (targ2.re * targ2.re) + (targ2.im * targ2.im);
  107.  
  108.     result->re = (targ1.re * targ2.re + targ1.im * targ2.im) / denom;
  109.     result->im = (targ1.im * targ2.re - targ1.re * targ2.im) / denom;
  110. }
  111.  
  112. /* converts a complex number arg in rectangular form to polar form */
  113.  
  114. polar(arg, modulus, amplitude)
  115. complex    *arg;
  116. double    *modulus, *amplitude;
  117. {
  118.     complex    targ;
  119.  
  120.     targ.re = arg->re;
  121.     targ.im = arg->im;
  122.  
  123.     if (fabs(targ.re) < TINY)
  124.         targ.re = 0.0;
  125.  
  126.     if (fabs(targ.im) < TINY)
  127.         targ.im = 0.0;
  128.  
  129.     *modulus = sqrt((targ.re * targ.re) + (targ.im * targ.im));
  130.  
  131.     if (targ.im == 0.0)
  132.         *amplitude = 0.0;
  133.     else if (targ.re == 0.0)  {
  134.         if (targ.im > 0.0)
  135.             *amplitude = M_PI_2;
  136.         else
  137.             *amplitude = -M_PI_2;
  138.     }
  139.     else if ((log(fabs(targ.im)) - log(fabs(targ.re))) > LOGHUGE) {
  140.         if (targ.re > 0.0)  {
  141.             if (targ.im > 0.0)
  142.                 *amplitude = M_PI_2;
  143.             else
  144.                 *amplitude = -M_PI_2;
  145.         }
  146.         else if (targ.im > 0.0)
  147.             *amplitude = -M_PI_2;
  148.         else
  149.             *amplitude = M_PI_2;
  150.         }
  151.     }
  152.     else
  153.         *amplitude = atan(targ.im / targ.re);
  154. }
  155.  
  156. /* raises arg to the positive integral power and returns answer in result */
  157.  
  158. ctopower(result, arg, power)
  159. complex    *result, *arg;
  160. int    power;
  161. {
  162.     int    i;
  163.     double    modulus, amplitude, newmod, powamp;
  164.     complex targ;
  165.  
  166.     targ.re = arg->re;
  167.     targ.im = arg->im;
  168.  
  169.     if (power == 0)  {
  170.         result->re = 1.0;
  171.         result->im = 0.0;
  172.     }
  173.     else {
  174.         polar(&targ, &modulus, &litude);
  175.         newmod = 1.0;
  176.  
  177.         if (power > 0)
  178.             for (i = 0; i < power; i++)
  179.                 newmod = newmod * modulus;
  180.         else
  181.             for (i = 0; i < abs(power); i++)
  182.                 newmod = newmod / modulus;
  183.  
  184.         powamp = power * amplitude;
  185.  
  186.         result->re = newmod * cos(powamp);
  187.         result->im = newmod * sin(powamp);
  188.     }
  189. }
  190.  
  191. /* raises e to the arg and returns the answer in result */
  192.  
  193. cexp(result, arg)
  194. complex    *result, *arg;
  195. {
  196.     double    expo;
  197.     complex    targ;
  198.  
  199.     targ.re = arg->re;
  200.     targ.im = arg->im;
  201.     expo = exp(targ.re);
  202.  
  203.     result->re = expo * cos(targ.im);
  204.     result->im = expo * sin(targ.im);
  205. }
  206.  
  207. /* takes the natural logarithm of arg and returns the answer in result */
  208.  
  209. cln(result, arg)
  210. complex    *result, *arg;
  211. {
  212.     double    modulus, amplitude;
  213.     complex    targ;
  214.     
  215.     targ.re = arg->re;
  216.     targ.im = arg->im;
  217.  
  218.     polar(&targ, &modulus, &litude);
  219.  
  220.     result->re = log(modulus);
  221.     result->im = amplitude;
  222. }
  223.  
  224. /* raise complex number arg1 to complex power arg2 and return answer in result */
  225.  
  226. ctoc(result, arg1, arg2)
  227. complex *result, *arg1, *arg2;
  228. {
  229.     complex    logpart, expo, targ1, targ2;
  230.  
  231.     targ1.re = arg1->re;
  232.     targ1.im = arg1->im;
  233.     targ2.re = arg2->re;
  234.     targ2.im = arg2->im;
  235.  
  236.     cln(&logpart, &targ1);
  237.     cmult(&expo, &targ2, &logpart);
  238.     cexp(result, &expo);
  239. }
  240.  
  241. /* takes the sine of arg and returns it in result */
  242.  
  243. csin(result, arg)
  244. complex *result, *arg;
  245. {
  246.     complex    targ, exp1, exp2, part1,
  247.         part2, sum, divisor;
  248.  
  249.     targ.re = arg->re;
  250.     targ.im = arg->im;
  251.     exp1.re = -targ.im;
  252.     exp1.im = targ.re;
  253.  
  254.     cexp(&part1, &exp1);
  255.  
  256.     exp2.re = targ.im;
  257.     exp2.im = -targ.re;
  258.  
  259.     cexp(&part2, &exp2);
  260.     csub(&sum, &part1, &part2);
  261.  
  262.     divisor.re = 0.0;
  263.     divisor.im = 2.0;
  264.  
  265.     cdiv(result, &sum, &divisor);
  266. }
  267.  
  268. /* takes the cosine of arg and returns it in result */
  269.  
  270. ccos(result, arg)
  271. complex    *result, *arg;
  272. {
  273.     complex    targ, exp1, exp2, part1,
  274.         part2, sum, divisor;
  275.  
  276.     targ.re = arg->re;
  277.     targ.im = arg->im;
  278.     exp1.re = -targ.im;
  279.     exp1.im = targ.re;
  280.  
  281.     cexp(&part1, &exp1);
  282.  
  283.     exp2.re = targ.im;
  284.     exp2.im = -targ.re;
  285.  
  286.     cexp(&part2, &exp2);
  287.     cadd(&sum, &part1, &part2);
  288.  
  289.     divisor.re = 2.0;
  290.     divisor.im = 0.0;
  291.  
  292.     cdiv(result, &sum, &divisor);
  293.