home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / pascal / library / dos / complex / complex.inc
Encoding:
Text File  |  1987-11-08  |  4.9 KB  |  167 lines
  1. {***********************************************************************}
  2. {                                                                       }
  3. {  TURBO Pascal library of Complex Number routines adapted from :       }
  4. {  October 1984 Dr. Dobbs Journal by John Lucas 30 Sept. 1984           }
  5. {                                                                       }
  6. {-----------------------------------------------------------------------}
  7. {                                                                       }
  8. {  Global Declarations needed by this Library :                         }
  9. {                                                                       }
  10. {  type                                                                 }
  11. {     complex = record                                                  }
  12. {        re : real;                                                     }
  13. {        im : real;                                                     }
  14. {     end;                                                              }
  15. {                                                                       }
  16. {***********************************************************************}
  17.  
  18. procedure cadd (var result : complex; arg1,arg2 : complex);
  19. begin
  20.    result.re := arg1.re + arg2.re;
  21.    result.im := arg1.im + arg2.im
  22. end;
  23.  
  24. procedure csub (var result : complex; arg1,arg2 : complex);
  25. begin
  26.    result.re := arg1.re - arg2.re;
  27.    result.im := arg1.im - arg2.im
  28. end;
  29.  
  30. procedure cmult (var result : complex; arg1,arg2 : complex);
  31. begin
  32.    result.re := arg1.re * arg2.re - arg1.im * arg2.im;
  33.    result.im := arg1.im * arg2.re - arg1.re * arg2.im
  34. end;
  35.  
  36. procedure cdiv (var result : complex; arg1,arg2 : complex);
  37. var
  38.    denom : real;
  39. begin
  40.    denom := sqr(arg2.re) + sqr(arg2.im);
  41.    result.re := (arg1.re * arg2.re + arg1.im * arg2.im)/denom;
  42.    result.im := (arg1.im * arg2.re - arg2.re * arg2.im)/denom
  43. end;
  44.  
  45. procedure polar (arg : complex; var modulus,amplitude : real);
  46. const                                                      
  47.    lnmaxreal = 87.49823353;
  48.    piover2   = 1.570796327;
  49.    closest   = 1E-19;
  50. begin
  51.    with arg do
  52.    begin
  53.       if abs(re) < closest then
  54.          re := 0.0;
  55.       if abs(im) < closest then
  56.          im := 0.0;
  57.       modulus := sqrt(sqr(re) + sqr(im));
  58.       if im = 0.0 then
  59.          amplitude := 0.0
  60.       else
  61.          if re = 0.0 then
  62.             if im = 0.0 then
  63.                amplitude := piover2
  64.             else
  65.                amplitude := -piover2
  66.          else
  67.             if (ln(abs(im)) - ln(abs(re)) > lnmaxreal) then
  68.                if re > 0.0 then
  69.                   if im > 0.0 then
  70.                      amplitude := piover2
  71.                   else
  72.                      amplitude := -piover2
  73.                else
  74.                   if im > 0.0 then
  75.                      amplitude := -piover2
  76.                   else
  77.                      amplitude := piover2
  78.             else
  79.                amplitude := arctan(im/re)
  80.    end;
  81. end;
  82.  
  83. procedure ctopower (var result : complex; arg : complex; power : integer);
  84. var
  85.    i : integer;
  86.    modulus,amplitude,newmod,powamp : real;
  87. begin
  88.    if power = 0 then
  89.    begin
  90.       result.re := 1.0;
  91.       result.im := 0.0
  92.    end
  93.    else
  94.    begin
  95.       polar(arg,modulus,amplitude);
  96.       newmod := 1.0;
  97.       if power > 0 then
  98.          for i := 1 to power do
  99.             newmod := newmod * modulus
  100.       else
  101.          for i := 1 to abs(power) do
  102.             newmod := newmod/modulus;
  103.       powamp := power * amplitude;
  104.       result.re := newmod * cos(powamp);
  105.       result.im := newmod * sin(powamp)
  106.    end;
  107. end;
  108.  
  109. procedure cexp (var result : complex; arg : complex);
  110. var
  111.    expo : real;
  112. begin
  113.    expo := exp(arg.re);
  114.    result.re := expo * cos(arg.im);
  115.    result.im := expo * sin(arg.im)
  116. end;
  117.  
  118. procedure cln (var result : complex; arg : complex);
  119. var
  120.    modulus,amplitude : real;
  121. begin
  122.    polar(arg,modulus,amplitude);
  123.    result.re := ln(modulus);
  124.    result.im := amplitude
  125. end;
  126.  
  127. procedure ctoc (var result : complex; arg1,arg2 : complex);
  128. var
  129.    logpart,expo : complex;
  130. begin
  131.    cln(logpart,arg1);
  132.    cmult(expo,arg2,logpart);
  133.    cexp(result,expo)
  134. end;
  135.  
  136. procedure csin (var result : complex; arg : complex);
  137. var
  138.    exp1,exp2,part1,part2,sum,divisor : complex;
  139. begin
  140.    exp1.re := -arg.im;
  141.    exp1.im := arg.re;
  142.    cexp(part1,exp1);
  143.    exp2.re := arg.im;
  144.    exp2.im := -arg.re;
  145.    cexp(part2,exp2);
  146.    csub(sum,part1,part2);
  147.    divisor.re := 0.0;
  148.    divisor.im := 2.0;
  149.    cdiv(result,sum,divisor)
  150. end;
  151.  
  152. procedure ccos (var result : complex; arg : complex);
  153. var
  154.    exp1,exp2,part1,part2,sum,divisor : complex;
  155. begin
  156.    exp1.re := -arg.im;
  157.    exp1.im := arg.re;
  158.    cexp(part1,exp1);
  159.    exp2.re := arg.im;
  160.    exp2.im := -arg.re;
  161.    cexp(part2,exp2);
  162.    cadd(sum,part1,part2);
  163.    divisor.re := 2.0;
  164.    divisor.im := 0.0;
  165.    cdiv(result,sum,divisor)
  166. end;
  167.