! 6. External Functions. 				   File: function6.f90
!------------------------------------------------------------------------------!

PROGRAM	interference
	IMPLICIT NONE

	! this program illustrates the use of external functions

	! example: this program calculates the intensity of an 
	! interference pattern for multiple (n) narrow slits

	! variable declarations:		
	REAL :: d		! distance between neighbouring slits
	INTEGER :: n_slits	! number of narrow slits
	REAL :: lambda		! wavelength of light
	REAL :: d_over_l	! ratio of d to lambda
	REAL :: phi		! relative phase of interfering waves
	REAL :: theta		! angle between direction of field point 
				! and normal to the plane of the slits
	REAL :: intensity	! intensity for each theta
	INTEGER :: j
	REAL, PARAMETER :: pi=3.1415926

	! intrinsic functions invoked:
	REAL, INTRINSIC :: REAL, SIN

	! external functions invoked
	REAL, EXTERNAL :: sinsin

	! inputs
	PRINT *, 'This program calculates the interference pattern of'
	PRINT *, 'n narrow slits separated by a distance x times the ' 
	PRINT *, 'wavelength of the light as a function of angle, theta.'
	PRINT *

	PRINT *, 'Specify d_over_l, the separation between neighbouring slits'
	PRINT *, 'in multiples of the light wavelength'
	PRINT *
	READ *, d_over_l
	PRINT *
	PRINT *, 'How many slits?'
	READ *, n_slits

	! compute theta, intensity for a range of theta values.
	! this DO loop repeats this part of the program, incrementing
	! j by one each time, starting at j=0 and ending when j=1000

	DO j=0,1000

		theta = REAL(j)*pi/4000.0
		phi = 2.*pi*d_over_l*SIN(theta)
		intensity = sinsin(phi,n_slits)
 
		! output to fort.18
		! produces a file consisting of two columns of figures 
		! which can be submitted to Gnuplot	
		WRITE(18,*) theta, intensity

	END DO

END PROGRAM interference

REAL FUNCTION sinsin(x,n)
	IMPLICIT NONE

	! evaluates (sin(nx/2)/sin(x/2))**2

	! dummy argument declarations:
	INTEGER, INTENT(IN) :: n
	REAL, INTENT(IN) :: x

	! this blockIF construct distinguishes between the case when 
	! x=0 or x nonzero

	IF ( ABS(x) < 1.E-5 ) THEN
		sinsin = n*n
	ELSE
		sinsin = ( SIN(0.5*REAL(n)*x) / SIN(0.5*x) )**2
	END IF

END FUNCTION sinsin

!! 1.	Now run this and use Gnuplot (or other graphics package) to view it. 
!!	Try n_slits=6, d_over_l=10

!! 2.	FUNCTIONs only ever return a single value to the calling program
!!	unit. If more than one quantity is to be calculated then a 
!!	SUBROUTINE should be used (see next session).

!! 3.	The structure of a function subprogram is
!!
!!		type  FUNCTION  function_name(arguments)
!!		...
!!		function_name = ...
!!		...
!!		END FUNCTION function_name
!!
!!	where
!!
!!	a) type = INTEGER, REAL, COMPLEX, CHARACTER, LOGICAL
!!	b) the value computed is associated with the function-name
!!	c) a function may have zero or more arguments
!!	d) somewhere inside the subprogram there must exist a statement 
!!	of the form "function_name = some value"
!!	e) a function may invoke another function or subroutine but may not 
!!	usually invoke itself (see text for a discussion of recursive functions)
!!
!!	Remember, however, that the variables used within subprograms are 
!!	local to that program only. Therefore they MUST be declared (REAL, 
!!	INTEGER, etc.) as in the calling program. When control is returned 
!!	to the calling program, these local variables disappear.
!!
!!	The arguments are dummy arguments. They are assigned values only when
!!	the subprogram is called (and hence evaluated for specific given 
!!	argument values).

!! 4.	Function subprograms are of two types. There are those which are 
!!	provided and are contained in a Fortran library of functions. Such
!!	functions are called INTRINSIC. There are those which you create
!!	yourself or acquire which are called EXTERNAL. Indeed you may build
!!	your own library of functions.
!!
!!	It is good practice to declare all the functions used in a program
!!	to be INTRINSIC or EXTERNAL. This avoids certain kinds of programming
!!	error. For example, some compilers provide additional functions which
!!	you may be unaware of and hence you may accidently use a name already
!!	in the library.

!! 5.	Exercises
!!
!!	1. Write a program which evaluates the ratio of the two roots of a 
!!	quadratic equation. The user should be requested for the coefficients
!!	of the equation and these should then be passed to an EXTERNAL FUNCTION
!!	which returns the ratio of the two roots. This ratio should then be 
!!	output from the main program unit.
!!
!!	2. Write an external function 'gaussian(x)' to evaluate the Gaussian 
!!	probability distribution
!!			f(x) = EXP(-x**2/2)/SQRT(2.0*pi)
!!
!!	3. Write a main program unit to invoke the function 'gaussian' to 
!!	calculate the Maxwellian velocity distribution (in 1D) 
!!
!!			f(v) = EXP(-v**2/(2*vt**2))/(SQRT(2*pi)*vt)
!!
!!	for a 100 values of v ranging from -3*vt to +3*vt. The constant vt is 
!!	the thermal velocity given by vt = SQRT(k*T/m) where k is Boltzmann's 
!!	constant, T is the temperature in K, and m is the particle mass. 
!!	Execute the program for electrons at a temperature of 10,000K.

!! End of file: function6.f90
