home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / simtel / sigm / vols000 / vol024 / complex.lib < prev    next >
Text File  |  1984-04-29  |  4KB  |  188 lines
  1. (***************************************************************
  2. *
  3. *    HERE ARE SOME MORE ITEMS FOR YOUR LIBRARY, FOR THOSE
  4. *  WHO USE THESE OR ANY OTHER PORTIONS OF OUR LIBRARY ROUTINES
  5. *  IT WOULD BE GREATLY APPRECIATED IF YOU WOULD SEND US YOUR 
  6. *  UPDATES OR MOD'S. IN FACT, SINCE WE OFFERED OURS WHY NOT SEND
  7. *  US YOURS.--editor
  8. {++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
  9. {+        COMPLEX LIBRARY                    +}
  10. {++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
  11. {
  12. Routines in this library:
  13.  
  14. CREAD        -Enter a complex number
  15. CWRITE        -Write a complex number
  16. MAG        -Computes the modulus of a complex number
  17. ADD        -Adds two complex numbers
  18. SUB        -Subtracts two complex numbers
  19. MULT        -Multiplies a real with a complex
  20. PRODUCT        -Product of two complex numbers
  21. QUOTIENT    -Quotient of two complex numbers
  22. CCOS        -Cosine of a complex
  23. POLAR        -Writing a complex into polar form
  24. CLN        -Natural logarithm of a complex
  25. SIGN        -Changes the sign of a complex
  26. CHECK        -Checks to see if the function argument is outside range
  27.  
  28. }
  29. {++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
  30.  
  31. type
  32.   complex = record
  33.           re, im : real
  34.         end;
  35.   S$255 = string 255;
  36.  
  37. PROCEDURE HALT(message: S$255);EXTERNAL;
  38.  
  39.  
  40. {++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
  41.  
  42. procedure cread(var z: complex);
  43. begin
  44.   read(z.re, z.im);
  45. end;
  46.  
  47. procedure cwrite(var z: complex);
  48. begin
  49.   writeln('(', z.re, ',', z.im, ')');
  50. end;
  51.  
  52. function mag(var z: complex): real;
  53. { computes the modulus of a complex number }
  54. begin
  55.   mag := sqrt( sqr(z.re) + sqr(z.im) );
  56. end;
  57.  
  58. procedure add( u, v: complex;
  59.           var w: complex);
  60. begin
  61.   w.re := u.re + v.re;
  62.   w.im := u.im + v.im;
  63. end { add };
  64.  
  65. procedure sub(u, v: complex;
  66.          var w: complex);
  67. begin
  68.   w.re := u.re - v.re;
  69.   w.im := u.im - v.im;
  70. end { sub };
  71.  
  72. procedure mult(    a: real;
  73.            z: complex;
  74.            var w: complex);
  75. { Multiplies a real with a complex }
  76. begin
  77.   w.re := a * z.re;
  78.   w.im := a * z.im;
  79. end { mult };
  80.  
  81. procedure product( u, v: complex;
  82.           var w: complex);
  83. begin
  84.   w.re := (u.re * v.re) - (u.im * v.im);
  85.   w.im := (u.re * v.im) + (u.im * v.re);
  86. end { product };
  87.  
  88. procedure quotient( u, v: complex;
  89.            var w: complex);
  90. const
  91.   sqrtwo = 1.414213562373095; { square root of 2 }
  92. var
  93.   vr, vi,
  94.   a, b,
  95.   x1, x2,
  96.   y1, y2,
  97.   root   : real;
  98. begin
  99.   vr := abs(v.re);
  100.   vi := abs(v.im);
  101.   root := sqrtwo * sqrt(vr) * sqrt(vi);
  102.   a := vr + vi + root;
  103.   b := vr + vi - root;
  104.   if (a = 0.0) or (b = 0.0) then
  105.      HALT('W: dividing by 0 in procedure quotient');
  106.   x1 := u.re / a;
  107.   x2 := v.re / b;
  108.   y1 := u.im / a;
  109.   y2 := v.im / b;
  110.   w.re := x1 * x2 + y1 * y2;
  111.   w.im := x2 * y1 - x1 * y2;
  112. end { quotient };
  113.  
  114. procedure ccos(     z: complex;
  115.         var c: complex);
  116. { Cosine of a complex }
  117. var
  118.   ep, em, p, m: real;
  119. begin
  120.   ep := exp(z.im);
  121.   em := 1.0 / ep;
  122.   p := ep + em;
  123.   m := em - ep;
  124.   c.re := 0.5 * p * cos(z.re);
  125.   c.im := 0.5 * m * sin(z.re);
  126. end { ccos };
  127.  
  128. procedure polar(    u: complex;
  129.         var v: complex);
  130. { Writing a complex into polar form }
  131. const
  132.   halfpi = 1.570796326795; { pi / 2.0 }
  133. begin
  134.   if (u.re = 0.0) and (u.im = 0.0) then
  135.      HALT('W: conversion of 0 in procedure polar');
  136.  
  137.   if (u.re = 0.0) and (u.im <> 0) then
  138.     begin
  139.       v.re := mag(u);
  140.       v.im := halfpi; {pi / 2.0}
  141.     end
  142.   else
  143.     begin
  144.       v.re := mag(u);
  145.       v.im := arctan(u.im / u.re);
  146.     end;
  147. end { polar };
  148.  
  149. procedure cln(    z: complex;
  150.           var c: complex);
  151. { Natural logarithm of a complex }
  152. var
  153.   p: complex;
  154. begin
  155.   polar(z,p);
  156.   c.re := ln(p.re);
  157.   c.im := p.im;
  158. end { cln };
  159.  
  160. procedure sign(u: complex; var v: complex);
  161. { Changes the sign of a complex }
  162. begin
  163.   v.re := -u.re;
  164.   v.im := -u.im;
  165. end { sign };
  166.  
  167. procedure check(z: complex);
  168. { Checks to see if the function argument is outside range }
  169. var
  170.   a, b: real;
  171. begin
  172.   a := abs(z.re);
  173.   b := abs(z.im);
  174.   if   ((a < 1.0E-5) and (b < 1.0E-5))
  175.     or ((b <> 0.0) and (b < 1.0E-5)) then
  176.     begin
  177.       write('W: small argument which causes exponent error = ');
  178.       cwrite(z);
  179.       HALT(' ');
  180.     end;
  181.   if b > 50.0 then
  182.     begin
  183.       write('W: argument with imaginary part outside range = ');
  184.       cwrite(z);
  185.       HALT(' ');
  186.     end;
  187. end { check };
  188.