home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / modelers / geomview / source.lha / Geomview / src / lib / gprim / discgrp / complex.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-06-29  |  3.7 KB  |  214 lines
  1. #include "options.h"
  2. #include "complex.h"
  3.  
  4. complex one        = {1.0, 0.0};
  5. complex zero    = {0.0, 0.0};
  6.  
  7. complex cplx_plus(z0, z1)
  8. complex z0, z1;
  9. {
  10.     complex sum;
  11.     
  12.     sum.real = z0.real + z1.real;
  13.     sum.imag = z0.imag + z1.imag;
  14.     return(sum);
  15. }
  16.  
  17.  
  18. complex cplx_minus(z0, z1)
  19. complex z0, z1;
  20. {
  21.     complex diff;
  22.     
  23.     diff.real = z0.real - z1.real;
  24.     diff.imag = z0.imag - z1.imag;
  25.     return(diff);
  26. }
  27.  
  28.  
  29. complex cplx_div(z0, z1)
  30. complex z0, z1;
  31. {
  32.     double    mod_sq;
  33.     complex    quotient;
  34.  
  35.     mod_sq            =  z1.real * z1.real  +  z1.imag * z1.imag;
  36.     quotient.real    = (z0.real * z1.real  +  z0.imag * z1.imag)/mod_sq;
  37.     quotient.imag    = (z0.imag * z1.real  -  z0.real * z1.imag)/mod_sq;
  38.     return(quotient);
  39. }
  40.  
  41.  
  42. complex cplx_mult(z0, z1)
  43. complex z0, z1;
  44. {
  45.     complex product;
  46.  
  47.     product.real    = z0.real * z1.real  -  z0.imag * z1.imag;
  48.     product.imag    = z0.real * z1.imag  +  z0.imag * z1.real;
  49.     return(product);
  50. }
  51.  
  52.  
  53. double modulus(z0)
  54. complex z0;
  55. {
  56.     return( sqrt(z0.real * z0.real  +  z0.imag * z0.imag) );
  57. }
  58.  
  59.  
  60. complex cplx_sqrt(z)
  61. complex z;
  62. {
  63.     double    mod,
  64.             arg;
  65.     complex    result;
  66.  
  67.     mod = sqrt(modulus(z));
  68.     if (mod == 0.0)
  69.         return(zero);
  70.     arg = 0.5 * atan2(z.imag, z.real);
  71.     result.real = mod * cos(arg);
  72.     result.imag = mod * sin(arg);
  73.     return(result);
  74. }
  75.  
  76.  
  77.  
  78. void sl2c_mult(a, b, product)
  79. sl2c_matrix    a,
  80.             b,
  81.             product;
  82. {
  83.     sl2c_matrix    temp;
  84.  
  85.     temp[0][0] = cplx_plus(cplx_mult(a[0][0], b[0][0]),  cplx_mult(a[0][1], b[1][0]));
  86.     temp[0][1] = cplx_plus(cplx_mult(a[0][0], b[0][1]),  cplx_mult(a[0][1], b[1][1]));
  87.     temp[1][0] = cplx_plus(cplx_mult(a[1][0], b[0][0]),  cplx_mult(a[1][1], b[1][0]));
  88.     temp[1][1] = cplx_plus(cplx_mult(a[1][0], b[0][1]),  cplx_mult(a[1][1], b[1][1]));
  89.     product[0][0] = temp[0][0];
  90.     product[0][1] = temp[0][1];
  91.     product[1][0] = temp[1][0];
  92.     product[1][1] = temp[1][1];
  93.     return;
  94. }
  95.  
  96.  
  97. void sl2c_copy(a, b)
  98. sl2c_matrix    a,
  99.             b;
  100. {
  101.     a[0][0] = b[0][0];
  102.     a[0][1] = b[0][1];
  103.     a[1][0] = b[1][0];
  104.     a[1][1] = b[1][1];
  105.     return;
  106. }
  107.  
  108.  
  109. /* normalizes a matrix to have determinant one */
  110. void sl2c_normalize(a)
  111. sl2c_matrix a;
  112. {
  113.     complex det,
  114.             factor;
  115.     double    arg,
  116.             mod;
  117.  
  118.     /* compute determinant */
  119.     det = cplx_minus(cplx_mult(a[0][0], a[1][1]),  cplx_mult(a[0][1], a[1][0]));
  120.     if (det.real == 0.0 && det.imag == 0.0) {
  121.         printf("singular sl2c_matrix\n");
  122.         exit(0);
  123.     }
  124.  
  125.     /* convert to polar coordinates */
  126.     arg = atan2(det.imag, det.real);
  127.     mod = modulus(det);
  128.  
  129.     /* take square root */
  130.     arg *= 0.5;
  131.     mod = sqrt(mod);
  132.  
  133.     /* take reciprocal */
  134.     arg = -arg;
  135.     mod = 1.0/mod;
  136.  
  137.     /* return to rectangular coordinates */
  138.     factor.real = mod * cos(arg);
  139.     factor.imag = mod * sin(arg);
  140.  
  141.     /* normalize matrix */
  142.     a[0][0] = cplx_mult(a[0][0], factor);
  143.     a[0][1] = cplx_mult(a[0][1], factor);
  144.     a[1][0] = cplx_mult(a[1][0], factor);
  145.     a[1][1] = cplx_mult(a[1][1], factor);
  146.  
  147.     return;
  148. }
  149.  
  150.  
  151. /* inverts a matrix;  assumes determinant is already one */
  152. void sl2c_invert(a, a_inv)
  153. sl2c_matrix a,
  154.             a_inv;
  155. {
  156.     complex    temp;
  157.  
  158.     temp = a[0][0];
  159.     a_inv[0][0] = a[1][1];
  160.     a_inv[1][1] = temp;
  161.  
  162.     a_inv[0][1].real = -a[0][1].real;
  163.     a_inv[0][1].imag = -a[0][1].imag;
  164.  
  165.     a_inv[1][0].real = -a[1][0].real;
  166.     a_inv[1][0].imag = -a[1][0].imag;
  167.  
  168.     return;
  169. }
  170.  
  171.  
  172. /* computes the square of the norm of a matrix */
  173. /* relies on the assumption that the sl2c matrix is stored in memory as 8 consecutive doubles */
  174. /* IS THIS RELIABLE? */
  175. double sl2c_norm_squared(a)
  176. sl2c_matrix a;
  177. {
  178.     register int    i;
  179.     register double    *p;
  180.     register double    sum;
  181.  
  182.     p = (double *) a;
  183.     sum = 0.0;
  184.     for (i=8; --i>=0; ) {
  185.         sum += *p * *p;
  186.         p++;
  187.     }
  188.     return(sum);
  189. }
  190.  
  191.  
  192. void sl2c_adjoint(a, ad_a)
  193. sl2c_matrix    a,
  194.             ad_a;
  195. {
  196.     complex    temp;
  197.  
  198.     /* transpose */
  199.     temp = a[0][1];
  200.     ad_a[0][0] = a[0][0];
  201.     ad_a[0][1] = a[1][0];
  202.     ad_a[1][0] = temp;
  203.     ad_a[1][1] = a[1][1];
  204.  
  205.     /* conjugate */
  206.     ad_a[0][0].imag = -ad_a[0][0].imag;
  207.     ad_a[0][1].imag = -ad_a[0][1].imag;
  208.     ad_a[1][0].imag = -ad_a[1][0].imag;
  209.     ad_a[1][1].imag = -ad_a[1][1].imag;
  210.  
  211.     return;
  212. }
  213.  
  214.