Geometry/basis_map_3d.f90 at master · lguez/Geometry · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
subroutine basis_map_3d ( u, v, a, ierror )

!*****************************************************************************80
!
!! BASIS_MAP_3D computes the matrix which maps one basis to another in 3D.
!
!  Discussion:
!
!    As long as the column vectors U1, U2 and U3 are linearly independent,
!    a matrix A will be computed that maps U1 to V1, U2 to V2, and
!    U3 to V3, where V1, V2 and V3 are the columns of V.
!
!    Depending on the values of the vectors, A may represent a
!    rotation, reflection, dilation, projection, or a combination of these
!    basic linear transformations.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    16 February 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, double precision U(3,3), the columns of U are the three
!    "domain" or "preimage" vectors, which should be linearly independent.
!
!    Input, double precision V(3,3), the columns of V are the three
!    "range" or "image" vectors.
!
!    Output, double precision A(3,3), a matrix with the property that
!    A * U1 = V1, A * U2 = V2 and A * U3 = V3.
!
!    Output, integer IERROR, error flag.
!    0, no error occurred.
!    nonzero, the matrix [ U1 | U2 | U3 ] is exactly singular.
!
  implicit none

  double precision a(3,3)
  double precision b(3,3)
  double precision c(3,3)
  double precision det
  integer ierror
  double precision u(3,3)
  double precision v(3,3)

  ierror = 0
!
!  Compute C = the inverse of [ U1 | U2 | U3 ].
!
  b(1:3,1:3) = u(1:3,1:3)

  call r8mat_inverse_3d ( b, c, det )

  if ( det == 0.0D+00 ) then
    ierror = 1
    return
  end if
!
!  A = [ V1 | V2 | V3 ] * inverse [ U1 | U2 | U3 ].
!
  a(1:3,1:3) = matmul ( v(1:3,1:3), c(1:3,1:3) )

  return
end