Geometry/conv3d.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
73
74
75
76
77
78
79
80
subroutine conv3d ( axis, theta, n, cor3, cor2 )

!*****************************************************************************80
!
!! CONV3D converts 3D data to a 2D projection.
!
!  Discussion:
!
!    A "presentation angle" THETA is used to project the 3D point
!    (X3D, Y3D, Z3D) to the 2D projection (XVAL,YVAL).
!
!    If AXIS = 'X':
!
!      X2D = Y3D - sin ( THETA ) * X3D
!      Y2D = Z3D - sin ( THETA ) * X3D
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    24 February 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character AXIS, the coordinate axis to be projected.
!    AXIS should be 'X', 'Y', or 'Z'.
!
!    Input, double precision THETA, the presentation angle in degrees.
!
!    Input, integer N, the number of points.
!
!    Input, double precision COR3(3,N), the 3D points.
!
!    Output, double precision COR2(2,N), the 2D projections.
!
  implicit none

  integer n

  character axis
  double precision cor2(2,n)
  double precision cor3(3,n)
  double precision degrees_to_radians
  double precision stheta
  double precision theta

  stheta = sin ( degrees_to_radians ( theta ) )

  if ( axis == 'X' .or. axis == 'x' ) then

    cor2(1,1:n) = cor3(2,1:n) - stheta * cor3(1,1:n)
    cor2(2,1:n) = cor3(3,1:n) - stheta * cor3(1,1:n)

  else if ( axis == 'Y' .or. axis == 'y' ) then

    cor2(1,1:n) = cor3(1,1:n) - stheta * cor3(2,1:n)
    cor2(2,1:n) = cor3(3,1:n) - stheta * cor3(2,1:n)

  else if ( axis == 'Z' .or. axis == 'z' ) then

    cor2(1,1:n) = cor3(1,1:n) - stheta * cor3(3,1:n)
    cor2(2,1:n) = cor3(2,1:n) - stheta * cor3(3,1:n)

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CONV3D - Fatal error!'
    write ( *, '(a)' ) '  Illegal coordinate index = "' // axis // '".'
    stop 1

  end if

  return
end