! 18. Matrices and Vectors.				File: matrix18.f90
!------------------------------------------------------------------------------!

MODULE rotation
	IMPLICIT NONE
CONTAINS
	FUNCTION rotate(angle)
		IMPLICIT NONE

		! this creates the 2x2 rotation matrix

		! declarations
		REAL, DIMENSION(2,2) :: rotate
		REAL, INTENT(IN) :: angle
		! local variables
		REAL :: c,s

		c = COS(angle)
		s = SIN(angle)

		rotate = RESHAPE( (/c,-s,s,c/),(/2,2/) )

	END FUNCTION rotate


	SUBROUTINE input(theta,vec_r)
		IMPLICIT NONE

		! declarations
		REAL, INTENT(OUT) :: theta
		REAL, DIMENSION(2), INTENT(OUT) :: vec_r

		! inputs
		PRINT *
		PRINT *,'    ROTATION OF (X,Y) AXES   '
		PRINT *
		PRINT *,'This program calculates the new (x,y) coordinates '
		PRINT *,'of a point in the plane as a result of an axes '
		PRINT *,'rotation through an angle theta.'
		PRINT *
		PRINT *,'Specify the angle of rotation in degrees'
		PRINT *,'between 0 and 360'
		PRINT *
		READ *, theta
		PRINT *
		PRINT *,'Specify the (x,y) coordinates in original frame'
		PRINT *
		READ *, vec_r
		PRINT *

	END SUBROUTINE input


	SUBROUTINE output(theta,vec_r,vec_rprime)
		IMPLICIT NONE

		! declarations
		REAL, INTENT(IN) :: theta
		REAL, DIMENSION(2), INTENT(IN) :: vec_r,vec_rprime

		! output
		PRINT *
		PRINT *,'The rotation angle is (in degrees) ', theta
		PRINT *
		PRINT *,'The original coordinates are ',vec_r
		PRINT *
		PRINT *,'The transformed coordinates are ',vec_rprime
		PRINT *

	END SUBROUTINE output


	SUBROUTINE check(angle0,vec_1,vec_2)
		IMPLICIT NONE

		! declarations
		REAL, PARAMETER :: epsilon=1.0e-4
		REAL :: distance_1,distance_2,dotprod,angle,pi
		REAL, DIMENSION(2), INTENT(IN) :: vec_1,vec_2
		REAL, INTENT(IN) :: angle0

		pi = 4.0*ATAN(1.0)

		! evaluate the distances from the common origin
		! check that the angle between vec_1 and vec_2 = theta 
		! (test for equality in their cosines)
		distance_1 = SQRT( DOT_PRODUCT(vec_1,vec_1) )
		distance_2 = SQRT( DOT_PRODUCT(vec_2,vec_2) )
		dotprod = DOT_PRODUCT(vec_1,vec_2)
		! divide by magnitudes to get COS(theta)
		! remember definition of scalar product
		dotprod = dotprod/(distance_1*distance_2)
		! ensure that this is less than 1 in magnitude since COS() 
		! must be less than 1
		IF (ABS(dotprod) > 1.0+epsilon) THEN
			PRINT *,'Error! COS(theta) > 1'
			STOP
		ELSE IF (ABS(dotprod-1.0) < epsilon) THEN
			dotprod = 1.0
		ELSE IF (ABS(dotprod+1.0) < epsilon) THEN
			dotprod = -1.0 
		END IF

		IF (ABS(distance_1 - distance_2) < epsilon) THEN
			PRINT*
			PRINT*,'The distance from the origin is conserved'
			PRINT *
		ELSE
			PRINT*
			PRINT*,'The distance from the origin is NOT conserved'
			PRINT *
		END IF

		IF (ABS(dotprod-COS(angle0)) < epsilon) THEN		
			PRINT*
			PRINT*,'The angle of rotation checks.'
			PRINT *
		ELSE
			PRINT*
			PRINT*,'The angle of rotation does NOT check'
			PRINT *
		END IF

	END SUBROUTINE check

END MODULE rotation


PROGRAM axis
	USE rotation
	IMPLICIT NONE

	! declarations
	REAL :: theta,thetarad,pi
	REAL, DIMENSION(2) :: vec_r,vec_rprime

	pi = 4.0*ATAN(1.0)

	! input angle and coords
	CALL input(theta,vec_r)

	! convert theta to radians
	thetarad = pi*theta/180.0

	! evaluate new coords
	vec_rprime = MATMUL(rotate(thetarad),vec_r)

	! output new coords
	CALL output(theta,vec_r,vec_rprime)

	! check: distance from the origin conserved in this transformation
	! and angle between old and new vectors is the angle of rotation
	CALL check(thetarad,vec_r,vec_rprime)

END PROGRAM axis


!! 1.	Fortran90 will handle arrays with rank up to 7 (i.e. 7-dimensional 
!!	arrays). On most occasions in physics and mathematics, we deal with 
!!	vectors (rank-1 arrays) and matrices (rank-2 arrays). 
!!
!!	Arrays can be arguments in subroutine calls or function references.
!!
!!	Fortran90 possesses a number of intrinsic functions which operate on 
!!	arrays of any dimension. Three of these are restricted to vector and 
!!	matrix operations. These are:
!!
!!	DOT_PRODUCT(vector1,vector2)
!!		This gives the scalar (dot) product of two vectors (rank-1 
!!		arrays). The answer is a scalar. It is equal to the sum of 
!!		the products of vector(i)*vector(i). It applies to vectors 
!!		of any extent (not just 3D vectors).
!!
!!	MATMUL(matrix1,matrix2)
!!		This gives the product of the matrix1*matrix2. Matrix1 must 
!!		be of rank 2. Matrix2 can be of rank 1 or 2. 
!!		The second dimension of matrix1 must equal the first dimension 
!!		of matrix2 (usual condition to be able to multiply 2 matrices).
!!
!!	TRANSPOSE(matrix)
!!		Transposes the matrix (rank-2 array).

!! 2.	Other useful functions.
!!
!!	RESHAPE
!!		We have already met this function.
!!
!!	SIZE(array,dim)
!!		This integer function returns the EXTENT (number of elements) 
!!		of 'array' along the specified dimension 'dim'. If 'dim' is 
!!		absent then the TOTAL number of the elements in the array is 
!!		returned.
!!
!!	LBOUND(array,dim)
!!		This returns the lower index bound of the dimension 'dim'. 
!!		If 'dim' is absent, the result is a rank-1 array containing 
!!		all the lower index bounds.
!!
!!	UBOUND(array,dim)
!!		This returns the upper index bound of the dimension 'dim'. 
!!		If 'dim' is absent, the result is a rank-1 array containing 
!!		all the upper index bounds.
!!
!!	For other intrinsics (e.g. MAXVAL, MINVAL, PRODUCT, SUM) see Appendix 
!!	of textbook.

!! 3.	Array-valued functions.
!!
!!	Note that the result of the reference to the external function 'rotate'
!!	is an array. Three points must be remembered when dealing with array 
!!	valued functions:
!!
!!		(a) the function must have an explicit interface. This is 
!!		ensured when placed within a module.
!!		(b) the function TYPE must be declared WITHIN the function 
!!		subprogram (and not as previously for scalar valued functions).
!!		(c) the result of the reference is an explicit-shape array.

!!	END OF FILE: matrix18.f90
