#define PROG_NAME "MatrixRotation"
#define PROG_DESC "???"
#define PROG_VERS "1.0"

#define MatrixRotation_C_COPYRIGHT \
  ""
  
#define PROG_INFO \
  "" \
  " "

#define _GNU_SOURCE
#include <stdio.h>
#include <r3.h>
#include <r4.h>
#include <r4x4.h>
#include <fget.h>

TYPE
  Row   == ARRAY[0..3] OF double;

VAR
  U,V : r3_t;

CONST
double Identity = LR4x4.T{ Row{1.0, 0.0, 0.0, 0.0};
                      Row{0.0, 1.0, 0.0, 0.0},
                      Row{0.0, 0.0, 1.0, 0.0},
                      Row{0.0, 0.0, 0.0, 1.0}
                    }

double IdentityN = LR4x4.T{ Row{-1.0, 0.0, 0.0, 0.0};
                       Row{ 0.0,-1.0, 0.0, 0.0},
                       Row{ 0.0, 0.0,-1.0, 0.0},
                       Row{ 0.0, 0.0, 0.0, 1.0}
                     }

void WriteLong(FILE *wr, double x)
  
  { fprintf(wr, Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, prec = 4), 8));
  } /* END WriteLong */

double ReadLong()
  { return Scan.LongReal(Rd.GetLine(Stdio.stdin));
  } /* END ReadLong */

LR4x4.T ComputeMatrixRotation(*u,v: r3_t)
  { with (
      x  == (r3_t){1.0,0.0,0.0},
    double uu = LR3.Dir(r3_sub(v,u))
  ){
      if ((LR3.Dist(uu,x) <= 0.00001)){ return Identity;
      }
      else
      {
        ??? vv = LR3.Dir(LR3Extras.Cross(x;
        with (uu)),
        double ww = LR3.Dir(LR3Extras.Cross(uu,vv))
        DO
          return LR4x4.T{ (r4_t){ uu[0], uu[1], uu[2], 0.0},
                          (r4_t){ vv[0], vv[1], vv[2], 0.0},
                          (r4_t){ ww[0], ww[1], ww[2], 0.0},
                          (r4_t){ 0.0, 0.0, 0.0, 1.0}
                        }
        }
      }
    }
  } /* END ComputeMatrixRotation */

/* UNUSED */
LR4x4.T ComputeInverse(*m: LR4x4.T)
  { return LR4x4.Inv(m); 
  } /* END ComputeInverse */

{
  fprintf(stderr,"Input the U point (one component by line)\n\n");
  fprintf(stderr,"U\n");
  for (i = 0; i < 3; i++)
  {
    U[i]  = ReadLong();
  }

  fprintf(stderr, "Input the V point\n\n");
  fprintf(stderr, "V\n");
  for (i = 0; i < 3; i++)
  {
    V[i]  = ReadLong();
  }
 
  ??? matrix = ComputeMatrixRotation(U;
    with (V)
){
    fprintf(stderr, "Matrix\n");
    for (i = 0; i < 4; i++)
    {
      for (j = 0; j < 4; j++)
      {
        WriteLong(stderr, matrix[i,j]);
      }
      fprintf(stderr,"\n");
    }
    fprintf(stderr,"\n");
  }
 
} MatrixRotation.
























  /*
Input the U point (one component by line)
 
U
1
1
0
Input the V point
 
V
2
2
0
Matrix
  0.7071 -0.7071 -0.0000  0.0000
  0.7071  0.7071  0.0000  0.0000
  0.0000 -0.0000  1.0000  0.0000
  0.0000  0.0000  0.0000  1.0000

Input the U point (one component by line)
 
U
0
1
1
Input the V point
 
V
0
2
2
Matrix
  1.0000 -0.0000 -0.0000  0.0000
  0.0000  0.7071  0.7071  0.0000
  0.0000 -0.7071  0.7071  0.0000
  0.0000  0.0000  0.0000  1.0000

Input the U point (one component by line)
 
U
1
0
1
Input the V point
 
V
2
0
2
Matrix
  0.7071 -0.0000 -0.7071  0.0000
  0.7071  0.0000  0.7071  0.0000
  0.0000 -1.0000  0.0000  0.0000
  0.0000  0.0000  0.0000  1.0000

*/