! 13. The Case Construct.				File: case13.f90
!------------------------------------------------------------------------------!

MODULE out
	IMPLICIT NONE
CONTAINS
	SUBROUTINE output(x_1,x_2,x_3)
		IMPLICIT NONE

		! this subroutine prints out the values of x_1,x_2,x_3 

		! dummy argument declarations
		REAL, INTENT(IN) :: x_1
		COMPLEX, INTENT(IN) :: x_2,x_3

		! outputs
		PRINT *,'The first root is ',x_1
		PRINT *
		PRINT *,'The second root is ',x_2
		PRINT *
		PRINT *,'The third root is ',x_3

	END SUBROUTINE output
END MODULE out


PROGRAM cubic
	USE out
	IMPLICIT NONE

	! this program obtains the roots of a cubic equation where possible

	! declarations
	INTEGER :: which
	REAL :: a0,a1,a2,p,q,r,s,s1,s2,t,z1
	REAL, PARAMETER :: epsilon=1.0e-6
	COMPLEX :: z2,z3

	! input
	PRINT *, 'This program calculates, when possible, the roots of '
	PRINT *, 'the cubic equation:'
	PRINT *
	PRINT *, '   z**3 + a2*z**2 + a1*z + a0 = 0  '
	PRINT *
	PRINT *, 'Type in values for a0, a1, a2 '
	PRINT *
	READ *, a0,a1,a2

	! intermediate variable definitions
	q = a1/3 -a2*a2/9
	r = (a1*a2-3*a0)/6 - a2**3/27
	s = q**3+r*r
	which = s/epsilon	! an integer
  
	! Find roots 
	SELECT CASE (which)
	CASE (1:)
		! s>0, one real root and pair of complex conjugate roots
		t = SQRT(s)

		! take the cube root of r+t; but watch out for negative r+t!!!
		s1 = SIGN( (ABS(r+t))**(1.0/3.0) ,r+t)
		! take the cube root of r-t; but watch out for negative r-t!!!
		s2 = SIGN( (ABS(r-t))**(1.0/3.0) ,r-t)

		! the roots are
		z1 = s1+s2 -a2/3
		z2 = CMPLX(-(s1+s2)/2-a2/3, SQRT(3.0)*(s1-s2)/2)
		z3 = CMPLX(-(s1+s2)/2-a2/3, -SQRT(3.0)*(s1-s2)/2)

		! output roots
		CALL output(z1,z2,z3)

	CASE (0)
		! s=0, all roots real and at least two are zero
		s1 = SIGN( (ABS(r))**(1.0/3.0) ,r)
		s2 = s1

		! the roots are
		z1 = 2*s1-a2/3
		z2 = -s1-a2/3
		z3 = z2

		! output roots
		CALL output(z1,z2,z3)

	CASE (:-1)
		! s<0, all roots real but cannot be found by this method
		PRINT *, 'The equation has three real roots '
		PRINT *, 'They cannot be found by this method '

	END SELECT

END PROGRAM cubic

!! 1.	The multiple blockIF construct works through the decision criteria in 
!!	the order in which they appear. Some occasions require decisions to be 
!!	made in order, and particularly when there is overlap between criteria.
!!	On the other hand, the decision criteria are often mutually exclusive 
!!	with no overlap. To cater for this case Fortran90 provides an 
!!	alternative form of selection called the CASE construct. Its syntax is
!!
!!		SELECT CASE (case expression)
!!		CASE (case selector)
!!			block of statements
!!		CASE (case selector)
!!			block of statements
!!		.......
!!		......
!!		END SELECT
!!
!!	a) Remember: the decision criteria must not overlap.
!!	b) In a blockIF, the expression in the IF() must be a logical 
!!	expression. In the CASE construct, the case expression in the CASE() 
!!	may be a logical expression, an integer expression, or a character 
!!	expression (but NOT a real expression). Note that they may be 
!!	expressions and not just single integer constants or variables.
!!	c) When the case expression is evaluated, the block of statements 
!!	following the CASE with this value will executed.
!!	d) A CASE DEFAULT statement may be included to cover cases where no 
!!	other match has been found.
!!	e) The case selector may take the forms:
!!
!!			single_value
!!			low_value:
!!			:high_value
!!			low_value:high_value
!!
!!	or a list of any combination of these, where the colon indicates a 
!!	range of values.

!! 2.	Exercise
!!
!!	Use the case construct to write a program which either adds, subtracts, 
!!	multiplies or divides two integers. The program should ask you to choose
!!	which option you wish. Each of the operations should be performed in 
!!	a separate function subprogram. The functions should be encapsulated in 
!!	a module.

!!	End of file: case13.f90	
