home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / cebit_91 / biomorph / complex.c < prev    next >
Encoding:
Text File  |  1994-03-01  |  14.8 KB  |  518 lines
  1. /**********************************************************/
  2. /*                                                        */
  3. /*                       COMPLEX.C                        */
  4. /*    INCLUDE-Datei zum Rechnen mit komplexen Zahlen.     */
  5. /*                                                        */
  6. /*                  Sprache: Turbo C 2.0                  */
  7. /*                                                        */
  8. /*          (C) 1991 Ralf Baethke-Franke  & toolbox       */
  9. /*                                                        */
  10. /**********************************************************/
  11.  
  12.  
  13. /**********************************************************/
  14. /*                                                        */
  15. /* Die Funktionen und ihre Berechnung:                    */
  16. /* -----------------------------------                    */
  17. /*                                                        */
  18. /* Es gilt grundsätzlich: z = x + iy                      */
  19. /*                        mit x = z.re und y = z.im       */
  20. /* Absolutbetrag:         r = abs_c(z)                    */
  21. /*                        r = |z| = √(x² + y²)            */
  22. /* Argument:              w =       atan(y/x)             */
  23. /*                                                        */
  24. /* 1. Grundrechenarten                                    */
  25. /* 1.1. Addition        z = add_c(z1,z2)                  */
  26. /*      z1 + z2 = (x1 + x2) + i(y1 + y2)                  */
  27. /* 1.2. Subtraktion     z = sub_c(z1,z2)                  */
  28. /*      z1 - z2 = (x1 - x2) - i(y1 - y2)                  */
  29. /* 1.3. Multiplikation  z = mul_c(z1,z2)                  */
  30. /*      z1 * z2 = x1x2 - y1y2 + i(x1y2 + y1x2)            */
  31. /* 1.4. Division        z = div_c(z1,z2)                  */
  32. /*        z1    x1x2 + y1y2     x1y2 - y1x2               */
  33. /*        --  = ----------- + i -----------               */
  34. /*        z2        |z|²            |z|²                  */
  35. /*                                                        */
  36. /* 2. Potenzen und Wurzeln                                */
  37. /* 2.1. Potenzen und Exponentialfunktion                  */
  38. /* 2.1.1. Quadrat         z = sqr_c(z1)                   */
  39. /*        z² = (x1² - y1²) + 2*i*x1y1                     */
  40. /* 2.1.2. reelle Potenz   z = pot_c(z1,n)                 */
  41. /*        z1ⁿ = r1ⁿ (cos(nw1) + isin(nw1))                */
  42. /* 2.1.3. komplexe Potenz z = pot_cc(z1,z2)               */
  43. /*        z1^z2 = exp(z2*ln(z1))                          */
  44. /* 2.1.4. Exponentialfkt. z = exp_c(z1)                   */
  45. /*        exp(z1) = exp(x1)*(cos(y1) + isin(y1))          */
  46. /* 2.2. Wurzeln und natürl. Logarithmus                   */
  47. /* 2.2.1. Quadratwurzel   z = sqrt_c(z1)                  */
  48. /*        √z1 =  ±√ ((r1 + x1)/2) ± i * √((r1 - x1)/2)    */
  49. /*  Es wird nur die erste der 4 mögl. Wurzeln berechnet   */
  50. /* 2.2.2. reelle Wurzel   z = rad_c(z1,n)                 */
  51. /*     ⁿ√z1 = z1^(1/n)                                    */
  52. /* 2.2.3. komplexe Wurzel z = rad_cc(z1,z2)               */
  53. /*        (z2)√z1 = z1^(1/z2)                             */
  54. /* 2.2.4. Logarithmus     z = ln_c(z1)                    */
  55. /*        ln(z1) = ln(r1) + i*w1                          */
  56. /* Es wird der Hauptzweig des nat. Logarithmus berechnet  */
  57. /*                                                        */
  58. /* 3. Trigonometrische Funktionen                         */
  59. /* 3.1. Sinus           z = sin_c(z1)                     */
  60. /*      sin(z1) = sin(x1)cosh(y1) + i * cos(x1)sinh(y1)   */
  61. /* 3.2. Cosinus         z = cos_c(z1)                     */
  62. /*      cos(z1) = cos(x1)cosh(y1) - i * sin(x1)sinh(y1)   */
  63. /* 3.3. Tangens         z = tan_c(z1)                     */
  64. /*      tan(z1) = sin(z1)/cos(z1)                         */
  65. /* 3.4. Cotangens       z = cot_c(z1)                     */
  66. /*      cot(z1) = cos(z1)/sin(z1)                         */
  67. /*                                                        */
  68. /* 4. Trigonometrische Umkehrfunktionen                   */
  69. /* 4.1. Arcussinus      z = asin_c(z1)                    */
  70. /*      asin(z1) = -i*ln(i*z1 + √(1 - z²))                */
  71. /* 4.2. Arcuscosinus    z = acos_c(z1)                    */
  72. /*      acos(z1) = -i*ln(z1 + √(z1² -1))                  */
  73. /* 4.3. Arcustangens    z = atan_c(z1)                    */
  74. /*                  1       1 + i*z1²                     */
  75. /*      atan(z1) = --- * ln(---------)                    */
  76. /*                 2*i      1 - i*z1²                     */
  77. /* 4.4. Arcuscotangens  z = acot_c(z1)                    */
  78. /*                 -1       i*z1² + 1                     */
  79. /*      acot(z1) = --- * ln(---------)                    */
  80. /*                 2*i      i*z1² - 1                     */
  81. /*                                                        */
  82. /* 5. Hyperbelfunktionen                                  */
  83. /* 5.1. Hyperbelsinus     z = sinh_c(z1)                  */
  84. /*      sinh(z1) = sinh(x1)cos(y1) + i*cosh(x1)sin(y1)    */
  85. /* 5.2. Hyperbelcosinus   z = cosh_c(z1)                  */
  86. /*      cosh(z1) = cosh(x1)cos(y1) + i*sinh(x1)sin(y1)    */
  87. /* 5.3. Hyperbeltangens   z = tanh_c(z1)                  */
  88. /*      tanh(z1) = sinh(z1) / cosh(z1)                    */
  89. /* 5.4. Hyperbelcotangens z = coth_c(z1)                  */
  90. /*      coth(z1) = cosh(z1) / sinh(z1)                    */
  91. /*                                                        */
  92. /* 6. Umkehrungen der Hyperbelfunktionen                  */
  93. /* 6.1. Arcus-Hyperbelsinus     z = asinh_c(z1)           */
  94. /*      asinh(z1) = ln(z1 + √(z1² + 1))                   */
  95. /* 6.2. Arcus-Hyperbelcosinus   z = acosh_c(z1)           */
  96. /*      acosh(z1) = ln(z1 + √(z1² - 1))                   */
  97. /* 6.3. Arcus-Hyperbeltangens   z = atanh_c(z1)           */
  98. /*                  1      1 + z                          */
  99. /*      atanh(z1) = - * ln(-----)                         */
  100. /*                  2      1 - z                          */
  101. /* 6.4. Arcus-Hyperbelcotangens z = atanh_c(z1)           */
  102. /*                  1      z + 1                          */
  103. /*      atanh(z1) = - * ln(-----)                         */
  104. /*                  2      z - 1                          */
  105. /*                                                        */
  106. /**********************************************************/
  107.  
  108. #include <math.h>
  109.  
  110. #define PI M_PI
  111. #define PI_HALBE PI/2
  112. #define ZWEI_PI 2 * PI
  113. #define INF 9e307
  114.  
  115. typedef struct
  116. {
  117.    double re, im;
  118. } compl;
  119.  
  120. /******************* Funktionsprototypen *****************/
  121.  
  122. double arg_opt ( double arg );
  123. double sqr ( double x );
  124. int sgn ( double i );
  125. int sign ( double i );
  126. double abs_c ( compl z );
  127. double arg_c ( compl z );
  128. compl konj_kompl ( compl z );
  129. compl add_c ( compl z1, compl z2 );
  130. compl sub_c ( compl z1, compl z2 );
  131. compl mul_c ( compl z1, compl z2 );
  132. compl div_c ( compl z1, compl z2 );
  133. compl sqr_c ( compl z1 );
  134. compl sqrt_c ( compl z1 );
  135. compl pot_c ( compl z1, double p);
  136. compl rad_c ( compl z1, double r );
  137. compl pot_cc ( compl z1, compl z2 );
  138. compl rad_cc ( compl z1, compl z2 );
  139. compl exp_c ( compl z1 );
  140. compl ln_c ( compl z1 );
  141. compl sin_c ( compl z1 );
  142. compl cos_c(compl z1);
  143. compl tan_c ( compl z1 );
  144. compl cot_c ( compl z1 );
  145. compl asin_c ( compl z1 );
  146. compl acos_c ( compl z1 );
  147. compl atan_c ( compl z1 );
  148. compl acot_c ( compl z1 );
  149. compl sinh_c ( compl z1 );
  150. compl cosh_c ( compl z1 );
  151. compl tanh_c ( compl z1 );
  152. compl coth_c ( compl z1 );
  153. compl asinh_c ( compl z1 );
  154. compl acosh_c ( compl z1 );
  155. compl atanh_c ( compl z1 );
  156. compl acoth_c ( compl z1 );
  157.  
  158. /********************* Hilfsfunktionen ********************/
  159.  
  160. double arg_opt ( double arg )
  161. /* dient zur Vermeidung von zu großen / kleinen Argu-     */
  162. /* menten bei Winkelfunktionen                            */
  163. {
  164.    if ( arg > ZWEI_PI )
  165.                    arg -= floor ( arg / ZWEI_PI ) * ZWEI_PI;
  166.    if ( arg < - ZWEI_PI )
  167.                    arg += floor ( arg / ZWEI_PI ) * ZWEI_PI;
  168.    return ( arg );
  169. }
  170.  
  171. double sqr ( double x )     /* Quadrat einer reellen Zahl */
  172. {
  173.    return( x * x );
  174. }
  175.  
  176. int sgn ( double i )     /* Vorzeichen einer reellen Zahl */
  177. {                        /* sgn(0) = +1                   */
  178.    int sgn;
  179.    sgn = ( i < 0 ) ? -1 : 1;
  180.    return ( sgn );
  181. }
  182.  
  183. int sign ( double i )    /* Vorzeichen einer reellen Zahl */
  184. {                        /* sign(0) = 0                   */
  185.    int sign;
  186.    sign = ( !i ) ? 0 : sgn ( i );
  187.    return ( sign );
  188. }
  189.  
  190. double abs_c ( compl z )   /* Betrag einer komplexen Zahl */
  191. {
  192.    return ( hypot ( z.re, z.im ) );
  193. }
  194.  
  195. double arg_c ( compl z )    /* Argument einer kompl. Zahl */
  196. {
  197.    double arg;
  198.    if ( !z.re ) arg = sign ( z.im ) * PI_HALBE;
  199.    else arg = ( z.re > 0 ) ? atan ( z.im / z.re ) :
  200.                    atan ( z.im / z.re ) + sgn ( z.im ) * PI;
  201.    return ( arg );
  202. }
  203.  
  204. compl konj_kompl ( compl z )  /* konjugiert komplexe Zahl */
  205. {                             /* zu z -> z'               */
  206.    z.im = -z.im;
  207.    return ( z );
  208. }
  209.  
  210. /************** Grundrechenarten: +, -, *, / **************/
  211.  
  212. compl add_c ( compl z1, compl z2 )
  213. {
  214.    compl z;
  215.    z.re = z1.re + z2.re;
  216.    z.im = z1.im + z2.im;
  217.    return ( z );
  218. }
  219.  
  220. compl sub_c ( compl z1, compl z2 )
  221. {
  222.    compl z;
  223.    z.re = z1.re - z2.re;
  224.    z.im = z1.im - z2.im;
  225.    return ( z );
  226. }
  227.  
  228. compl mul_c ( compl z1, compl z2 )
  229. {
  230.    compl z;
  231.    z.re = z1.re * z2.re - z1.im * z2.im;
  232.    z.im = z1.re * z2.im + z1.im * z2.re;
  233.    return ( z );
  234. }
  235.  
  236. compl div_c ( compl z1, compl z2 )
  237. {
  238.    compl z;
  239.    if ( abs_c ( z2 ) == 0 )
  240.    {
  241.       z.re =  INF * sign ( z2.re );
  242.       z.im = -INF * sign ( z2.im );
  243.    }
  244.    else
  245.    {
  246.       z.re = ( z1.re * z2.re + z1.im * z2.im ) /
  247.                                        sqr ( abs_c ( z2 ) );
  248.       z.im = ( z1.im * z2.re - z1.re * z2.im ) /
  249.                                        sqr ( abs_c ( z2 ) );
  250.    }
  251.    return ( z );
  252. }
  253.  
  254. /***************** Potenzen und Wurzeln *******************/
  255.  
  256. compl sqr_c ( compl z1 )
  257. {
  258.    compl z;
  259.    z.re = sqr ( z1.re ) - sqr ( z1.im );
  260.    z.im = 2 * z1.re * z1.im;
  261.    return ( z );
  262. }
  263.  
  264. compl sqrt_c ( compl z1 )
  265. {
  266.    compl z;
  267.    z.re = sqrt ( ( abs_c ( z1 ) + z1.re ) / 2 );
  268.    z.im = sqrt ( ( abs_c ( z1 ) - z1.re ) / 2 );
  269.    return(z);
  270. }
  271.  
  272. compl pot_c ( compl z1, double p)
  273. {
  274.    compl z;
  275.    double betr, arg;
  276.    if ( p == 1 ) z = z1;
  277.    else if ( p == -1 )
  278.    {
  279.       z.re = 1;
  280.       z.im = 0;
  281.       z = div_c ( z, z1 );
  282.    }
  283.    else
  284.    {
  285.       betr = pow ( abs_c ( z1 ), p );
  286.       arg  = arg_opt ( p * arg_c ( z1 ) );
  287.       z.re = betr * cos ( arg );
  288.       z.im = betr * sin ( arg );
  289.    }
  290.    return ( z );
  291. }
  292.  
  293. compl rad_c ( compl z1, double r )
  294. {
  295.    compl z;
  296.    z = pot_c ( z1, 1 / r );
  297.    return ( z );
  298. }
  299.  
  300. compl pot_cc ( compl z1, compl z2 )
  301. {
  302.    compl z;
  303.    z = exp_c ( mul_c ( z2, ln_c ( z1 ) ) );
  304.    return ( z );
  305. }
  306.  
  307. compl rad_cc ( compl z1, compl z2 )
  308. {
  309.    compl z;
  310.    z.re = 1; z.im = 0;
  311.    z = pot_cc ( z1, div_c ( z, z2 ) );
  312.    return ( z );
  313. }
  314.  
  315. compl exp_c ( compl z1 )
  316. {
  317.    compl z;
  318.    z1.im = arg_opt ( z1.im );
  319.    z1.re = exp ( z1.re );
  320.    z.re  = z1.re * cos ( z1.im );
  321.    z.im  = z1.re * sin ( z1.im );
  322.    return ( z );
  323. }
  324.  
  325. compl ln_c ( compl z1 )
  326. {
  327.    compl z;
  328.    z.re = log ( abs_c ( z1 ) );
  329.    z.im = arg_c ( z1 );
  330.    return(z);
  331. }
  332.  
  333. /************** Trigonometrische Funktionen ***************/
  334.  
  335. compl sin_c ( compl z1 )
  336. {
  337.    compl z;
  338.    z1.re = arg_opt ( z1.re );
  339.    z.re = sin ( z1.re ) * cosh ( z1.im );
  340.    z.im = cos ( z1.re ) * sinh ( z1.im );
  341.    return(z);
  342. }
  343.  
  344. compl cos_c(compl z1)
  345. {
  346.    compl z;
  347.    z1.re = arg_opt ( z1.re );
  348.    z.re =  cos ( z1.re ) * cosh ( z1.im );
  349.    z.im = -sin ( z1.re ) * sinh ( z1.im );
  350.    return(z);
  351. }
  352.  
  353. compl tan_c ( compl z1 )
  354. {
  355.    compl z;
  356.    z = div_c ( sin_c ( z1 ), cos_c ( z1 ) );
  357.    return ( z );
  358. }
  359.  
  360. compl cot_c ( compl z1 )
  361. {
  362.    compl z;
  363.    z = div_c ( cos_c ( z1 ), sin_c ( z1 ) );
  364.    return ( z );
  365. }
  366.  
  367. /****** Umkehrungen der trigonometrischen Funktionen ******/
  368.  
  369. compl asin_c ( compl z1 )
  370. {
  371.    compl z, c1, c2, c3;
  372.    c1.re = 1; c1.im = 0;
  373.    c2.re = 0; c2.im = 1;
  374.    c3.re = 0; c3.im = -1;
  375.    z = sqr_c ( z1 );
  376.    z = sub_c ( c1, z);
  377.    z = sqrt_c ( z );
  378.    z = add_c ( mul_c ( c2, z1 ), z );
  379.    z = ln_c ( z );
  380.    z = mul_c ( c3, z );
  381.    return ( z );
  382. }
  383.  
  384. compl acos_c ( compl z1 )
  385. {
  386.    compl z, c1, c2;
  387.    c1.re = 1; c1.im = 0;
  388.    c2.re = 0; c2.im = -1;
  389.    z = sqr_c ( z1 );
  390.    z = sub_c ( z, c1);
  391.    z = sqrt_c ( z );
  392.    z = add_c ( z1, z );
  393.    z = ln_c ( z );
  394.    z = mul_c ( c2, z );
  395.    return ( z );
  396. }
  397.  
  398. compl atan_c ( compl z1 )
  399. {
  400.    compl z, c1, c2, c3, a1, a2;
  401.    c1.re = 0; c1.im = 1;
  402.    c2.re = 1; c2.im = 0;
  403.    c3.re = 0; c3.im = 2;
  404.    z = mul_c ( c1, z1 );
  405.    a1 = add_c ( c2, z );
  406.    a2 = sub_c ( c2, z );
  407.    z = div_c ( a1, a2 );
  408.    z = ln_c ( z );
  409.    a1 = div_c ( c2, c3 );
  410.    z = mul_c ( a1, z );
  411.    return ( z );
  412. }
  413.  
  414. compl acot_c ( compl z1 )
  415. {
  416.    compl z, c1, c2, c3, a1, a2;
  417.    c1.re = 0; c1.im = 1;
  418.    c2.re = 1; c2.im = 0;
  419.    c3.re = 0; c3.im = -2;
  420.    z = mul_c ( c1, z1 );
  421.    a1 = add_c ( z, c2 );
  422.    a2 = sub_c ( z, c2 );
  423.    z = div_c ( a1, a2 );
  424.    z = ln_c ( z );
  425.    a1 = div_c ( c2, c3 );
  426.    z = mul_c ( a1, z );
  427.    return ( z );
  428. }
  429.  
  430. /******************* Hyperbel-Funktionen ******************/
  431.  
  432. compl sinh_c ( compl z1 )
  433. {
  434.    compl z;
  435.    z1.im = arg_opt ( z1.im );
  436.    z.re = sinh ( z1.re ) * cos ( z1.im );
  437.    z.im = cosh ( z1.re ) * sin ( z1.im );
  438.    return ( z );
  439. }
  440.  
  441. compl cosh_c ( compl z1 )
  442.  
  443. {
  444.    compl z;
  445.    z1.im = arg_opt ( z1.im );
  446.    z.re = cosh ( z1.re ) * cos ( z1.im );
  447.    z.im = sinh ( z1.re ) * sin ( z1.im );
  448.    return ( z );
  449. }
  450.  
  451. compl tanh_c ( compl z1 )
  452. {
  453.    compl z;
  454.    z = div_c ( sinh_c ( z1 ), cosh_c ( z1 ) );
  455.    return ( z );
  456. }
  457.  
  458. compl coth_c ( compl z1 )
  459. {
  460.    compl z;
  461.    z = div_c ( cosh_c ( z1 ), sinh_c ( z1 ) );
  462.    return ( z );
  463. }
  464.  
  465. /*********** Umkehrungen der Hyperbel-Funktionen **********/
  466.  
  467. compl asinh_c ( compl z1 )
  468. {
  469.    compl z, c;
  470.    c.re = 1; c.im = 0;
  471.    z = sqr_c ( z1 );
  472.    z = add_c ( z, c );
  473.    z = sqrt_c ( z );
  474.    z = add_c ( z1, z );
  475.    z = ln_c ( z );
  476.    return ( z );
  477. }
  478.  
  479. compl acosh_c ( compl z1 )
  480. {
  481.    compl z, c;
  482.    c.re = -1; c.im = 0;
  483.    z = sqr_c ( z1 );
  484.    z = add_c ( z, c );
  485.    z = sqrt_c ( z );
  486.    z = add_c ( z1, z );
  487.    z = ln_c ( z );
  488.    return ( z );
  489. }
  490.  
  491. compl atanh_c ( compl z1 )
  492. {
  493.    compl z, c1, c2, a1, a2;
  494.    c1.re = 1; c1.im = 0;
  495.    c2.re = 2; c2.im = 0;
  496.    a1 = add_c ( c1, z1 );
  497.    a2 = sub_c ( c1, z1 );
  498.    z = div_c ( a1, a2 );
  499.    z = ln_c ( z );
  500.    z = div_c ( z, c2 );
  501.    return ( z );
  502. }
  503.  
  504. compl acoth_c ( compl z1 )
  505. {
  506.    compl z, c1, c2, a1, a2;
  507.    c1.re = 1; c1.im = 0;
  508.    c2.re = 2; c2.im = 0;
  509.    a1 = add_c ( z1, c1 );
  510.    a2 = sub_c ( z1, c1 );
  511.    z = div_c ( a1, a2 );
  512.    z = ln_c ( z );
  513.    z = div_c ( z, c2 );
  514.    return ( z );
  515. }
  516.  
  517. /******************* Ende von COMPLEX.C ******************/
  518.