home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / FRACTAL.ZIP / COMPLEX.PAS next >
Encoding:
Pascal/Delphi Source File  |  1985-03-11  |  7.5 KB  |  292 lines
  1. { Complex Function Library for use with TURBO PASCAL include statement
  2.   based on D.D. Clark article in Dr. Dobb's Journal, Oct. 1984, with
  3.   modifications to make Polar actually work properly, etc.
  4.            J.M.Weiss, 1572 Peacock Ave., Sunnyvale, CA 94087
  5.                             3/6/84 }
  6.  
  7. TYPE Complex = RECORD
  8.                   re : Real;
  9.                   im : Real
  10.                END; {of RECORD}
  11.  
  12. {Following is for Turbo-87 Version with 8087 support }
  13. CONST Precision: Real=1.5e-16;
  14.       max_real:Real=1.67e308;
  15.       big_sqrt:Real=1.29e154;
  16.       closest:Real=4.19e-307;
  17.  
  18. {Following is for Turbo Version without 8087 support
  19. CONST Precision: Real=1.e-11;
  20.       max_real:Real=1.e38;
  21.       big_sqrt:Real=1.e19;
  22.       closest:Real=1.e-39; }
  23.  
  24. PROCEDURE Cadd( Arg1, Arg2: Complex; VAR Result: Complex);
  25. {adds Arg1 and Arg2 and returns sum in Result}
  26.  
  27. BEGIN {Cadd}
  28.    Result.re := Arg1.re + Arg2.re;
  29.    Result.im := Arg1.im + Arg2.im;
  30. END;{ of Cadd}
  31.  
  32. PROCEDURE Csub( Arg1, Arg2: Complex; VAR Result: Complex);
  33. {subtracts Arg2 and Arg1 and returns difference in Result}
  34.  
  35. BEGIN {Csub}
  36.    Result.re := Arg1.re - Arg2.re;
  37.    Result.im := Arg1.im - Arg2.im;
  38. END;{ of Csub}
  39.  
  40. PROCEDURE Cmult( Arg1, Arg2: Complex; VAR Result: Complex );
  41. {multiplies Arg1 by Arg2 and returns quotient in Result}
  42.  
  43. BEGIN {Cmult}
  44.    Result.re := Arg1.re*Arg2.re - Arg1.im*Arg2.im;
  45.    Result.im := Arg1.im*Arg2.re + Arg1.re*Arg2.im
  46. END; {of Cmult}
  47.  
  48. PROCEDURE Cdiv( Arg1, Arg2: Complex; VAR Result: Complex );
  49. {divides Arg1 by Arg2 and returns quotient in Result}
  50.  
  51. VAR Denom: Real;
  52.  
  53. BEGIN {Cdiv}
  54.    Denom := Sqr(Arg2.re) + Sqr(Arg2.im);
  55.    Result.re := (Arg1.re*Arg2.re + Arg1.im*Arg2.im)/Denom;
  56.    Result.im := (Arg1.im*Arg2.re - Arg1.re*Arg2.im)/Denom
  57. END; {of Cdiv}
  58.  
  59.  
  60.  
  61.  
  62.  
  63. PROCEDURE Cexp( Arg: Complex; VAR Result: Complex );
  64. {raises e to the arg and returns as result}
  65.  
  66. VAR   Expo: Real;
  67.  
  68. BEGIN {Cexp }
  69.    Expo := Exp(Arg.re);
  70.    Result.re := Expo*Cos(Arg.im);
  71.    Result.im := Expo*Sin(Arg.im)
  72. END; {of Cexp}
  73.  
  74.  
  75.  
  76.  
  77.  
  78. PROCEDURE Csin( Arg: Complex; VAR Result:  Complex );
  79. {takes the sine of arg and returns as result}
  80. {Warning 15 digits precision ==> sin(pi)=1.e-16}
  81.  
  82. VAR One,Exp1,Part1,Part2,Sum,Divisor:    Complex;
  83.  
  84. BEGIN {Csin}
  85.    One.re := 1.0;
  86.    One.im := 0.0;
  87.    Exp1.re := - Arg.im;        { z*i }
  88.    Exp1.im := Arg.re;
  89.    Cexp(Exp1, Part1);           { exp(zi) }
  90.    Cdiv(One,Part1,Part2);
  91.    Csub(Part1, Part2, Sum);
  92.    Divisor.re := 0.0;
  93.    Divisor.im := 2.0;
  94.    Cdiv(Sum, Divisor, Result);       { exp(zi) - exp(-zi))/(2i) }
  95.    IF Abs(Result.re) < Precision THEN Result.re := 0.0;
  96.    IF Abs(Result.im) < Precision THEN Result.im :=0.0;
  97. END; { of Csin }
  98.  
  99. PROCEDURE Ccos( Arg: Complex; VAR Result: Complex );
  100. {takes the cosine of arg and returns as result}
  101. {Precision warning: 15 digits:  cos(pi/2)=6.e-17}
  102.  
  103. VAR One,Exp1,Part1,Part2,Sum,Divisor:    Complex;
  104.  
  105. BEGIN {Ccos}
  106.    One.re := 1.0;
  107.    One.im := 0.0;
  108.    Exp1.re := - Arg.im;        { z*i }
  109.    Exp1.im := Arg.re;
  110.    Cexp(Exp1, Part1);           { exp(zi) }
  111.    Cdiv(One,Part1,Part2);
  112.    Cadd(Part1, Part2, Sum);
  113.    Divisor.re := 2.0;
  114.    Divisor.im := 0.0;
  115.    Cdiv(Sum, Divisor, Result);       { exp(zi) + exp(-zi))/2 }
  116.    IF Abs(Result.re) < Precision THEN Result.re := 0.0;
  117.    IF Abs(Result.im) < Precision THEN Result.im :=0.0;
  118. END; {of Ccos }
  119.  
  120. PROCEDURE Ctan (Arg: Complex; VAR Result: Complex);
  121.  
  122. VAR Arg1,Arg2: Complex;
  123.  
  124. BEGIN {Ctan}
  125.    Csin(Arg,Arg1);
  126.    Ccos(Arg,Arg2);
  127.    Cdiv( Arg1,Arg2, Result)
  128. END; {of Ctan}
  129.  
  130. PROCEDURE Cmplx (Arg1, Arg2: real; VAR Result: Complex);
  131.  
  132. BEGIN {Cmplx}
  133.    Result.re := Arg1;
  134.    Result.im := Arg2
  135. END; {of Cmplx}
  136.  
  137. PROCEDURE Conj(VAR Arg: Complex);
  138. {take complex conjugate of complex variable Arg}
  139.  
  140. BEGIN {Conj}
  141.    Arg.im := -Arg.im
  142. END; {of Conj}
  143.  
  144. PROCEDURE Rect(Modulus,Phase: Real; VAR Result: Complex);
  145. {convert modulus and phase to rectangular coordinate}
  146.  
  147. BEGIN {Rect}
  148.    Result.re := Modulus*Cos(Phase);
  149.    Result.im := Modulus*Sin(Phase)
  150. END; {of Rect}
  151.  
  152.  
  153.  
  154. PROCEDURE Csqr(Arg: Complex; VAR Result:Complex);
  155. {take complex square of Arg and return Result}
  156.  
  157. BEGIN {Csqr}
  158.    Cmult(Arg,Arg,Result)
  159. END; {of Csqr}
  160.  
  161. PROCEDURE Cscale(Scale: Real; Arg: Complex; VAR Result:Complex);
  162.  
  163. BEGIN {Cscale}
  164.    Result.re := Scale*Arg.re;
  165.    Result.im := Scale*Arg.im
  166. END; {of Cscale}
  167.  
  168. FUNCTION SafeDivide(x,y:Real) : Real;
  169. { divide x by y and intercept incipient underflow or overflow}
  170.  
  171. LABEL 1;
  172.  
  173. BEGIN  {SafeDivide}
  174.    if Abs(y) < 1.0 then        {check overflow}
  175.       if Abs(x) > Abs(y)*max_real then
  176.          SafeDivide := max_real
  177.       else
  178.          goto 1
  179.    else if abs(x)  < 1.0 then  {check underflow}
  180.       if Abs(y) > Abs(x)*max_real then
  181.          SafeDivide := 0.0
  182.       else
  183.          goto 1
  184.    else
  185. 1:    SafeDivide := x/y         {normal calculation}
  186. END;  {of SafeDivide}
  187.  
  188. FUNCTION ATan2(x,y:Real) : Real;
  189. {take arctangent of x/y and account for right angles and proper quadrant}
  190. VAR
  191.     a: Real;
  192. BEGIN {ATan2}
  193.    if y = 0.0 then     {check for exceptional cases with y=0.0}
  194.       if(x = 0.0) then  {if both arguments are zero, return zero}
  195.          a := 0.0
  196.       else
  197.          a := 0.5*pi   {otherwise its a right angle}
  198.    else begin
  199.       a := SafeDivide(x,y);  { try the divide, preventing unseemly errors}
  200.       if a>= 1.e308 then
  201.          a:= 0.5*pi
  202.       else begin
  203.          if a<> 0.0 then
  204.          begin
  205.             a:=Abs(a);
  206.             if a>big_sqrt then
  207.                a:=0.5*pi
  208.             else if a< closest then
  209.                a:=0.0
  210.             else
  211.                a:=ArcTan(a);
  212.          end;
  213.          if y < 0.0 then   {choose upper or lower half of plane}
  214.             a:= pi - a;
  215.       end;
  216.    end;
  217.    if x < 0.0 then         {choose left or right half of plane}
  218.       a := -a;
  219.    ATan2 := a
  220. end;  {of ATan2}
  221.  
  222. PROCEDURE Polar(Arg: Complex; VAR Modulus, Phase: Real);
  223. {converts a complex number arg in rectangular form to polar form}
  224.  
  225. BEGIN {Polar}
  226.    if Abs(Arg.re) < closest then
  227.       Arg.re := 0.0;
  228.    if Abs(Arg.im) < closest then
  229.       Arg.im := 0.0;
  230.    Modulus := Sqrt(Sqr(Arg.re) + Sqr(Arg.im));
  231.    Phase := ATan2(Arg.im,Arg.re)
  232. END; { of Polar}
  233.  
  234. PROCEDURE Ctopower( Arg: Complex; Power: Integer; VAR Result: Complex );
  235. {raises Arg to integral Power and returns the answer in Result}
  236.  
  237. VAR  i: Integer;
  238.      Modulus, Phase, NewMod, PowAmp: Real;
  239.  
  240. BEGIN {Ctopower}
  241.    IF Power = 0 THEN BEGIN
  242.       Result.re := 1.0;
  243.       Result.im := 0.0
  244.    END
  245.    ELSE BEGIN
  246.       Polar(Arg, Modulus, Phase);
  247.       NewMod := 1;
  248.       IF Power >0 THEN
  249.          FOR i := 1 TO Power DO NewMod := NewMod*Modulus
  250.       ELSE
  251.          FOR i := 1 TO Abs(Power) DO NewMod := NewMod/Modulus;
  252.       PowAmp := Power*Phase;
  253.       Result.re := NewMod*Cos(PowAmp);
  254.       Result.im := NewMod*Sin(PowAmp)
  255.    END
  256. END; { of Ctopower}
  257.  
  258. PROCEDURE Cln( Arg: Complex; VAR Result: Complex );
  259. {takes natural log of arg  as result}
  260.  
  261. VAR Modulus,Phase: Real;
  262.  
  263. BEGIN {Cln}
  264.    Polar(Arg, Modulus, Phase);
  265.    Result.re := Ln(Modulus);
  266.    Result.im := Phase
  267. END; { of Cln }
  268.  
  269. PROCEDURE CtoC( Arg1,Arg2: Complex; VAR Result: Complex );
  270. {takes complex number arg1 to complex power arg2}
  271.  
  272. VAR   LogPart, Expo: Complex;
  273.  
  274. BEGIN { CtoC }
  275.    Cln(Arg1, LogPart);
  276.    Cmult(Arg2,LogPart, Expo);
  277.    Cexp(Expo, Result )
  278. END; {of CtoC}
  279.  
  280. PROCEDURE Csqrt(Arg: Complex; VAR Result:Complex);
  281. {take complex square root of Arg and return Result}
  282.  
  283. VAR
  284.    Argmag,Argphs: Real;
  285.  
  286. BEGIN {Csqrt}
  287.    Polar(Arg,Argmag,Argphs);
  288.    Argmag:=Sqrt(Argmag);
  289.    Argphs:=0.5*Argphs;
  290.    Rect(Argmag,Argphs,Result)
  291. END; {of Csqrt}
  292.