home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / modelers / geomview / source.lha / Geomview / src / lib / gprim / discgrp / projective.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-02-16  |  4.3 KB  |  198 lines
  1. #include "options.h"
  2. #include "complex.h"
  3. #include "projective.h"
  4.  
  5. /*    sl2c_to_proj() converts sl2c_matrices to proj_matrices.  It is based on
  6. an explanation provided by Craig Hodgson of a program written by
  7. Diane Hoffoss which in turn was based on an algorithm explained to her
  8. by Bill Thurston.  */
  9.  
  10. void sl2c_to_proj(s, p)
  11. sl2c_matrix    s;
  12. proj_matrix    p;
  13. {
  14.     int            j;        /* which column of p    */
  15.     sl2c_matrix    ad_s,    /* s* = adjoint of s    */
  16.                 fs,        /* f(s) = s m s*        */
  17.                 temp;
  18.     static sl2c_matrix    m[4] = {{{{ 0.0, 0.0},{ 0.0, 1.0}},
  19.                                  {{ 0.0,-1.0},{ 0.0, 0.0}}},
  20.  
  21.                                 {{{ 0.0, 0.0},{ 1.0, 0.0}},
  22.                                  {{ 1.0, 0.0},{ 0.0, 0.0}}},
  23.  
  24.                                 {{{-1.0, 0.0},{ 0.0, 0.0}},
  25.                                  {{ 0.0, 0.0},{ 1.0, 0.0}}},
  26.  
  27.                                 {{{ 1.0, 0.0},{ 0.0, 0.0}},
  28.                                  {{ 0.0, 0.0},{ 1.0, 0.0}}}};
  29.  
  30.     for (j=0; j<4; j++) {
  31.         sl2c_adjoint(s, ad_s);
  32.         sl2c_mult(s, m[j], temp);
  33.         sl2c_mult(temp, ad_s, fs);
  34.         p[0][j] = fs[0][1].imag;
  35.         p[1][j] = fs[0][1].real;
  36.         p[2][j] = 0.5 * (fs[1][1].real - fs[0][0].real);
  37.         p[3][j] = 0.5 * (fs[1][1].real + fs[0][0].real);
  38.     }
  39.  
  40.     return;
  41. }
  42.  
  43.  
  44. void proj_to_sl2c(p, s)
  45. proj_matrix    p;
  46. sl2c_matrix    s;
  47. {
  48.     double    t2,
  49.             t3,
  50.             aa,
  51.             bb;
  52.  
  53.     /* Notation:  The entries of s are a, b, c, d (as is standard).    */
  54.     /* The complex conjugate of a is written as a' (read "a-bar").    */
  55.  
  56.     /* Outline of algorithm:  Write down the four matrices {sM0s*,    */
  57.     /* sM1s*, sM2s*, sM3s*} in terms of {a, b, c, d}, and find        */
  58.     /* expressions for {2a'a, 2a'b, 2a'c, 2a'd} as differences of    */
  59.     /* matrix entries.  Express these differences as functions of    */
  60.     /* the entries of p.  The matrix (2a'a, 2a'b; 2a'c, 2a'd) is a    */
  61.     /* multiple of the desired matrix (a, b; c, d), so when we        */
  62.     /* normalize the former we get the latter (we'd have to            */
  63.     /* normalize in any case, since this scheme can find            */
  64.     /* (a, b; c, d) only up to multiplication by a complex number    */
  65.     /* of modulus one).                                                */
  66.  
  67.     /* a may be zero, but a and b can't both be zero.  We'll use    */
  68.     /* the one with the bigger norm.                                */
  69.  
  70.     t2 = p[3][2] - p[2][2]; /* 00 entry of F(M2)    */
  71.     t3 = p[3][3] - p[2][3]; /* 00 entry of F(M3)    */
  72.     aa = t3 - t2;            /* aa = 2 * |a|^2        */
  73.     bb = t3 + t2;            /* bb = 2 * |b|^2        */
  74.     if (aa > bb) {
  75.         s[0][0].real =   aa;
  76.         s[0][0].imag =   0.0;
  77.         s[0][1].real =   p[3][1] - p[2][1];
  78.         s[0][1].imag =   p[3][0] - p[2][0];
  79.         s[1][0].real =   p[1][3] - p[1][2];
  80.         s[1][0].imag =   p[0][2] - p[0][3];
  81.         s[1][1].real =   p[0][0] + p[1][1];
  82.         s[1][1].imag =   p[1][0] - p[0][1];
  83.     }
  84.     else {
  85.         s[0][0].real =   p[3][1] - p[2][1];
  86.         s[0][0].imag =   p[2][0] - p[3][0];
  87.         s[0][1].real =   bb;
  88.         s[0][1].imag =   0.0;
  89.         s[1][0].real =   p[1][1] - p[0][0];
  90.         s[1][0].imag = - p[0][1] - p[1][0];
  91.         s[1][1].real =   p[1][3] + p[1][2];
  92.         s[1][1].imag = - p[0][2] - p[0][3];
  93.     }
  94.  
  95.     sl2c_normalize(s);
  96.  
  97.     return;
  98. }
  99.  
  100.  
  101. void proj_mult(a, b, product)
  102. proj_matrix    a,
  103.             b,
  104.             product;
  105. {
  106.     register int    i,
  107.                     j,
  108.                     k;
  109.     register double    sum;
  110.     proj_matrix        temp;
  111.  
  112.     for (i=0; i<4; i++)
  113.         for (j=0; j<4; j++) {
  114.             sum = 0.0;
  115.             for (k=0; k<4; k++)
  116.                 sum += a[i][k] * b[k][j];
  117.             temp[i][j] = sum;
  118.         }
  119.  
  120.     for (i=0; i<4; i++)
  121.         for (j=0; j<4; j++)
  122.             product[i][j] = temp[i][j];
  123.  
  124.     return;
  125. }
  126.  
  127.  
  128. void proj_copy(a, b)
  129. proj_matrix    a,
  130.             b;
  131. {
  132.     register int    i,
  133.                     j;
  134.  
  135.     for (i=0; i<4; i++)
  136.         for (j=0; j<4; j++)
  137.             a[i][j] = b[i][j];
  138.  
  139.     return;
  140. }
  141.  
  142.  
  143. /* proj_invert() assumes the matrix a is nonsingular, as will always    */
  144. /* be the case for matrices representing isometries of H^3.                */
  145.  
  146. void proj_invert(m, m_inv)
  147. proj_matrix    m,
  148.             m_inv;
  149. {
  150.     int            i, j, k;
  151.     double        scratch[4][8],
  152.                 *a[4],
  153.                 *temp;
  154.  
  155.     /* set up */
  156.     for (i=4; --i>=0; ) {
  157.         for (j=4; --j>=0; ) {
  158.             scratch[i][j]    = m[i][j];
  159.             scratch[i][j+4]    = (i == j) ? 1.0 : 0.0;
  160.         }
  161.         a[i] = scratch[i];
  162.     }
  163.  
  164.     /* do the forward part of Gaussian elimination */
  165.     for (j=0; j<4; j++) {
  166.         /* find the best pivot */
  167.         for (i=j+1; i<4; i++)
  168.             if (fabs(a[i][j]) > fabs(a[j][j])) {
  169.                 temp = a[i];
  170.                 a[i] = a[j];
  171.                 a[j] = temp;
  172.             }
  173.  
  174.         /* normalize row j */
  175.         for (i=j+1; i<8; i++)
  176.             a[j][i] /= a[j][j];
  177.  
  178.         /* clear the lower part of column j */
  179.         for (i=j+1; i<4; i++)
  180.             for (k=j+1; k<8; k++)
  181.                 a[i][k] -= a[i][j] * a[j][k];
  182.     }
  183.  
  184.     /* do the back substitution */
  185.     for (j=4; --j>=0; )
  186.         for (i=j; --i>=0; )
  187.             for (k=4; k<8; k++)
  188.                 a[i][k] -= a[i][j] * a[j][k];
  189.             
  190.     /* copy the answer into m_inv */
  191.     for (i=4; --i>=0; )
  192.         for (j=4; --j>=0; )
  193.             m_inv[i][j] = a[i][j+4];
  194.  
  195.     return;
  196. }
  197.  
  198.