! 14. The Do Loop.					File: doloop14.f90
!------------------------------------------------------------------------------!

MODULE initial
	IMPLICIT NONE
CONTAINS
	SUBROUTINE setup(lambdamin,lambdamax,delta,r0)
		IMPLICIT NONE

		! example: Rayleigh scattering						
		! This program evaluates the extent (the differential 
		! scattering cross section) to which non-polarized light 
		! is scattered through an angle theta by colloidal particles. 
		! It is measured by Rayleigh's ratio:
		!    R = 8*(alpha)**2 (pi/lambda)**4 (1 + (cos(THETA))**2)
		! in units of metres squared. alpha is the polarisability of 
		! the particle, lambda is the wavelength of the light. The 
		! inverse dependence on wavelength means that blue light is 
		! scattered away from the incident beam more than red. This 
		! is why the sky looks blue 
		! Run the program for wavelength range 400nm (deep violet) 
		! to 700nm (deep red)

		! declarations
		REAL, INTENT(OUT) :: lambdamin	! minimum wavelength selected
		REAL, INTENT(OUT) :: lambdamax	! maximum wavelength selected
		REAL, INTENT(OUT) :: delta	! wavelength increment
		REAL, INTENT(OUT) :: r0
		REAL :: theta		! scattering angle
		REAL :: alpha		! polarizability of particle
		REAL :: pi

		! inputs
		! read in values for scattering angle, initial wavelength, 
		! wavelength increment and number of points
		! all lengths expressed in nanometers (nm)

		PRINT *, 'This program evaluates the Rayleigh ratio R'
		PRINT *, 'for a given range of wavelengths and'
		PRINT *, 'scattering angle'
		PRINT *
		PRINT *, 'Specify scattering angle as a fraction of pi'
		READ *, theta
		PRINT *, 'Specify a minimum wavelength in nanometres'
		PRINT *, '(suggest 400nm, deep violet)'
		READ *, lambdamin
		PRINT *, 'Specify wavelength increment in nanometres'
		PRINT *, '(suggest 20nm)'
		READ *, delta
		PRINT *, 'Specify maximum wavelength in nanometres'
		PRINT *, '(suggest 700nm, deep red)'
		READ *, lambdamax

		PRINT *
		PRINT *, 'The results are written both to screen and saved'
		PRINT *, 'in the file: fort.31'

		pi = 4.*ATAN(1.)
		alpha = 6.e-3
		theta = theta*pi
		r0 = 8.*alpha*alpha*(pi**4)
		r0 = r0*(1. + COS(theta)**2)

		! echo input values to screen
		PRINT *
		PRINT *, 'Minimum wavelength is   ',lambdamin
		PRINT *, 'Maximum wavelength is   ',lambdamax
		PRINT *, 'Wavelength increment is ',delta
		PRINT *, 'Scattering angle is     ',theta,' radians'
		PRINT *

	END SUBROUTINE setup
END MODULE initial


MODULE lambda
	IMPLICIT NONE
CONTAINS
	SUBROUTINE ratio(lambdamin,lambdamax,delta,r0)
		IMPLICIT NONE

		INTEGER :: n, i
		REAL, INTENT(IN) :: lambdamin,lambdamax,delta,r0
		REAL :: r		! Rayleigh's ratio
		REAL :: lambda		! wavelength of light

		! number of values to be calculated
		n = INT( (lambdamax - lambdamin)/delta )


		!now we come to the do loop
		DO i = 0,n
		   lambda = lambdamin + delta*i
		   r = r0/(lambda**4)

		   ! output each value to the screen
		   WRITE(6,*) ' Lambda = ',lambda, ' ratio = ',r
		   ! also write each value to the file fort.10; plot these
		   WRITE (31,*) lambda, r
		END DO
	END SUBROUTINE ratio
END MODULE lambda


PROGRAM rayleigh
	USE initial
	USE lambda
	IMPLICIT NONE

	! note how concise this main program unit is.

	! declarations
	REAL :: lambdamin,lambdamax,delta,r0

	! provide input values and calculate lambda independent part of ratio
	CALL setup(lambdamin,lambdamax,delta,r0)

	! generate Rayleigh ratio for a range og wavelengths
	CALL ratio(lambdamin,lambdamax,delta,r0)

END PROGRAM rayleigh


!! 1.	The general form of the DO construct is
!!
!!		DO  count = initial,final,increment
!!			block of statements
!!		END DO
!!	where
!!
!!	count 	  = variable counting the passes through the loop
!!	initial   = starting value of count
!!	final     = final value of count
!!	increment = increment in loop variable, count
!!
!!	These MUST be integer variables and declared as such.
!!
!!	If 'increment' is omitted the increment is by default equal to 1.
!!	This is the case in the above example.

!! 2.	The number of trips through the loop is calculated to be
!!
!!		MAX((final-initial+increment)/increment,0)
!!
!!	This calculation is performed once before the looping begins to 
!!	calculate how many trips are to be made.

!! 3.	DO-loops may count downwards by selecting a negative increment. 
!!	In this case, final must be less than initial.

!! 4.	Nested DO-loops. One DO-loop may contain another DO-loop. Such a 
!!	nested DO-loop must be completely contained within the outer DO-loop. 
!!	No overlapping is permitted. DO-loops may be repeatedly nested.
!!
!!	Example: Extract, study and run this program.
!!
!!	Problem: To find all three-digit numbers equal to the sum of the 
!!		 cubes of their digits
!!
!!	PROGRAM htu
!!		IMPLICIT NONE
!!
!!		!declarations
!!		INTEGER :: hundred, ten, unit
!!		LOGICAL :: truefalse
!!
!!		! inputs
!!		PRINT *, 'The following numbers are equal to'
!!		PRINT *, 'the sums of the cubes of their digits:'
!!		PRINT *
!!
!!		! run through the numbers from 100 to 999 
!!
!!		DO hundred =1,9
!!	   	   DO ten = 0,9
!!	      	      DO unit = 0,9
!!			 truefalse = (100*hundred + 10*ten + unit) &
!!			    == (hundred**3 + ten**3 + unit**3)
!!		 	 IF (truefalse) PRINT *, hundred, ten, unit
!!		      END DO
!!		   END DO
!!		END DO
!!
!!	END PROGRAM htu 
!!

!! 5.	The value of the loop index must NOT not be altered within the loop. 
!!	Its value, however, may be used within the loop.

!! 6.	In many instances, it will not be known in advance how many times a 
!!	loop needs to be executed. In this case, we can use a form of loop 
!!	which is not count-controlled:
!!
!!			DO
!!			   block of statements
!!			END DO
!!
!!	However, some criterion must be included to decide when to finish 
!!	cycling round the loop. Otherwise we have an INFINITE loop. Great 
!!	care must be taken to avoid such eventualities and the use of this 
!!	form of the DO construct is to be avoided if possible or used only 
!!	when it is certain that infinite loops will not arise.
!!
!!	Its use often arises in iterative schemes or in summing series where 
!!	we want to continue until a certain accuracy is achieved. Suppose 
!!	we are summing a series. We will continue to add terms to this series 
!!	until the magnitude of the latest term is so small that it does not 
!!	alter the current answer to the accuracy that is required. Then, all 
!!	we need do is to specify a small number, 'epsilon' say, and include 
!!	a transfer of control line in the loop as follows:
!!
!!		DO
!!		   ....
!!		   IF (term < epsilon) EXIT
!!		   ...
!!		END DO
!!
!!	If 'term < epsilon' then the EXIT statement will transfer control to 
!!	the statement immediately following the END DO statement of the loop 
!!	within which the  EXIT statement lies. Be careful if there are nested 
!!	loops.

!! 7.	STOP
!!
!!	Occasionally, a situation will arise where termination of a program is 
!!	required before the END statement has been reached. Perhaps an error 
!!	has been detected, or a wrong input value given, or .... The insertion 
!!	of STOP will terminate execution of the program immediately.

!! 8.	RETURN
!!
!!	If, within a procedure, we wish execution of the procedure to be 
!!	terminated and control returned to the program unit which called or 
!!	referenced it then we need only include the single statement: RETURN.

!! 9.	Exercises
!!
!!	(a) Identify the errors in the following program fragment:
!!
!!			INTEGER :: row, col, nrows, ncols
!!			DO row =1,nrows
!!			   row = row +1
!!			      DO col =1,ncols
!!			         PRINT *, row, col
!!			END DO
!!
!!	(b) WARNING: INFINITE LOOP!
!!	Do not run this program unless you know how to halt a running loop!
!!
!!	PROGRAM infinite
!!	   IMPLICIT NONE
!!
!!	   ! Before you run this program, make sure you know how to 
!!	   ! interrupt it: CTRL C. Beware of infinite loops!!!!!
!!	   DO
!!	      PRINT *, 'A confession:'
!!	      PRINT *, 'Physics is fun; it''s the physicist I can''t stand!'
!!	   END DO
!!	END PROGRAM infinite
!!
!!
!!	(c) The series representation of sin(x) is 
!!
!!		sin(x) = x - x**3/3! + x**5/5! - x**7/7! + ....
!!
!!	Write a function subprogram to evaluate this up to 'nmax' terms where 
!!	'nmax' is a given parameter. Then write a main program unit to obtain 
!!	sin(x) for a specified range of x. Also, compare your answer with the 
!!	intrinsic function SIN(x).
!!
!!	(d) A series for ATAN(X) is
!!
!!		arctan(x) = x - x**3/3 + x**5/5 - x**7/7 + ....
!!
!!	Write a program to calculate ATAN(X), but calling it by a name other
!!	than the Fortran intrinsic function. Hence check the value of pi using
!!
!!			pi = 4.0*ATAN(1.0)
!!
!!	(e) The sawtooth function 
!!
!!		y(x) = x/pi		0<=x<pi
!!		     = x/pi - 2.	pi<=x<=2*pi
!!
!!	can be approximated by the Fourier series
!!
!!		y(x) = sum of terms a_k from k=1 to k=N
!!
!!		a_k = (2/pi)*((-1)**(k-1)) sin(kx)/k
!!
!!	Write a program to sum the first N (to be read in) terms of this
!!	approximation, and have your program find the maximum deviation of
!!	the approximation from the true function. Read in M, the number of 
!!	values of x (between 0 to 2*pi) you wish y to be calculated for .
!!
!!	Plot both the exact function and the Fourier approximation using 
!!	GNUPLOT.

!! 10.	Read Chapter 6 of Ellis

!!	END OF FILE: doloop14.f90
