! 7. Subroutines.					   File: subroutine7.f90
!------------------------------------------------------------------------------!

PROGRAM great_circle
	IMPLICIT NONE

	! Program to evaluate the distance along the earth's surface between 
	! two places whose latitudes and longitudes are (lat1,long1) and 
	! (lat2,long2). The earth's radius is declared in miles in the 
	! PARAMETER statement and the distance between these two places is 
	! output in miles, assuming the earth is perfectly spherical.

	! variable declarations:
	REAL :: lat1,long1,lat2,long2
	REAL :: surface_distance,direct_distance

	! inputs:
	CALL input_data(lat1,long1,lat2,long2)


	! evaluate distance along earth's surface between these two points
	CALL great_circle_distance(lat1,long1,lat2,long2,surface_distance, &
                                   direct_distance)

	! output this distance
	CALL output_data(lat1,long1,lat2,long2,surface_distance,direct_distance)

END PROGRAM great_circle


SUBROUTINE input_data(lat1,long1,lat2,long2)
	IMPLICIT NONE

	! prompts the user for input

	! dummy argument declarations:
	REAL, INTENT(OUT) :: lat1,long1,lat2,long2
	REAL, INTRINSIC :: SIN,COS,SQRT,ASIN

	PRINT *,  'Enter the lat. and long. (in degrees) of your starting point'
	READ *,  lat1,long1
	PRINT *,  'Enter the lat. and long. (in degrees) of your destination'
	READ *,  lat2,long2

END SUBROUTINE input_data


SUBROUTINE great_circle_distance(lat1,long1,lat2,long2,surface_distance, &
				 direct_distance)
	IMPLICIT NONE

	! evaluates the distance along the earth's surface between two positions
	! with latitude and longitude (lat1,long1) and (lat2,long2)

	! dummy argument declarations:
	REAL,INTENT(IN)  :: lat1,long1,lat2,long2
	REAL,INTENT(OUT) :: surface_distance,direct_distance

	! local variables and constants:
	REAL,PARAMETER :: pi=3.14159265,radius=3957
	REAL,PARAMETER :: factor=pi/180.0       ! conversion: degrees to radians
	REAL :: theta1,fi1      ! lat. and long. in radians of starting position
	REAL :: theta2,fi2      ! lat. and long. in radians of destination
	REAL :: x1,y1,z1        ! rectangular coordinates of starting position
	REAL :: x2,y2,z2        ! rectangular coordinates of destination

	! convert angular coordinates to radians:
	theta1 = lat1*factor
	theta2 = lat2*factor
	fi1 = long1*factor
	fi2 = long2*factor

	! convert spherical coordinates to rectangular:
	x1 = radius*cos(theta1)*sin(fi1)
	x2 = radius*cos(theta2)*sin(fi2)
	y1 = radius*cos(theta1)*cos(fi1)
	y2 = radius*cos(theta2)*cos(fi2)
	z1 = radius*sin(theta1)
	z2 = radius*sin(theta2)

	! evaluate direct distance (i.e. through the earth) between 
	! (x1,y1,z1) and (x2,y2,z2):
	direct_distance = sqrt( (x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2 )

	! evaluate distance along earth's surface:
	surface_distance = 2 * radius*asin(0.5*direct_distance/radius)

END SUBROUTINE great_circle_distance


SUBROUTINE output_data(la1,lo1,la2,lo2,s_d,d_d)
	IMPLICIT NONE

	! outputs results

	! dummy argument declarations:
	REAL, INTENT(IN) :: la1,lo1,la2,lo2
	REAL, INTENT(IN) :: s_d,d_d

	PRINT *, 'distance along earth''s surface between (',la1,',',lo1,') &
        & and (',la2,',',lo2,') is ',s_d,' miles'
	PRINT *, 'The DIRECT distance is ',d_d,' miles'

END SUBROUTINE output_data


!! 1.	The above example consists of the main program and 3 subroutines.
!!	Observe how simple the main program now looks.
!!
!! 2.	The CALL statement.
!!
!!	In contrast to the FUNCTION subprogram, a subroutine is invoked
!!	by a CALL statement:
!!
!!		CALL  subroutine_name (argument_list)
!!
!!	The subroutine_name is the identifier of the subroutine. It is
!! 	NOT associated with any returned value and hence need not be 
!!	declared.
!!
!!	The argument_list is an optional list of ACTUAL arguments, of the
!!	form (expr1,expr2,,....). Values are sent to and returned from the 
!!	subroutine by means of this list. Note that the parentheses () are
!!	required with this list unless there are no arguments in which case
!!	they may be omitted. An argument may be an expression.
!!
!! 3.	The subroutine syntax.
!!
!!	This is like any other program except for the header line and the
!!	END statement which returns control to the calling program. Hence
!!
!!		SUBROUTINE  subroutine_name(arguments)
!!		...
!!		...
!!		END SUBROUTINE subroutine_name
!!
!!	When there are no arguments, the parentheses are omitted.
!!	The outputs of a subroutine are returned to the caller by means of
!!	the argument list. Indeed, an output argument may replace an input 
!!	argument.
!!
!!	A subroutine may invoke other subroutines or functions, but may not
!!	usually invoke itself (see text for a discussion of recursive 
!!	subroutines). Subroutine names follow the rules for variables (31 
!!	characters, first a letter, etc.).
!!
!!	Remember that variable names used in the subroutine are LOCAL to that
!!	subroutine and have no meaning outside it. The same name may frequently
!!	be used in different programs to stand for quite different variables.
!! 	So, be careful
!!
!!	Also, be aware that the arguments in the subroutine are DUMMY variables
!!	having no values associated with them. It is only when the subroutine 
!!	is called that specific ACTUAL values are passed to them either as a 
!!	value or as an actual variable.
!!
!!	The parameters passed to a subroutine need not have the same names as
!!	those in the subroutine heading line. For example, in the above program
!! 	we have:
!!
!!	CALL output_data(lat1,long1,lat2,long2,surface_distance,direct_distance)
!!
!!					and 	
!!
!!	SUBROUTINE output_data(la1,lo1,la2,lo2,s_d,d_d)
!!
!!	However the TYPE of each correpsonding element MUST match, i.e.
!!	`lat1' and `la1' must be of the same type (REAL in this case), and the
!!	same for `long1' and `lo1', `lat2' and `la2', etc.

!! 4.	Exercises.
!!
!!	1. Write a program which evaluates BOTH roots of a quadratic equation.
!!	Subroutines should be called to
!!		a) input the coefficients of the equation
!!		b) evaluate the roots of the equation
!!		c) output the value of the roots
!!		d) display an 'end of program' message on the screen
!!
!!	2. Write a subroutine, called SWAP, to interchange the values of two
!!	real variables. Then use SWAP to sort three numbers into ascending 
!!	order.
!!
!!	3.  Euler's method for the solution of 1st order odes:
!!
!!	We wish to solve
!!
!!			dy/dx =f(x,y)
!!
!!	Since the derivative is the limit, as h tends to zero, of
!!
!!			{y(x+h) - y(x)} / h
!!
!!	use this as an approximation to the derivative keeping h "sufficiently"
!!	small. 
!!
!!	Let  x_n = n*h  and  y_n = y(x_n) , then Euler's method is
!!
!!			y_(n+1) = y_n + h*f(x_n,y_n)
!!
!!	with a given initial condition  y_0 = a .
!!
!!	Write a subroutine to carry out this method, the function f(x,y) being
!!	supplied from a FUNCTION subprogram. Have the values for each (x_n,y_n)
!!	written to the file fort.100
!!
!!	Write a main program which invokes this subroutine to solve an ode 
!!	for which you know the answer. Compare the numerical and analytic
!!	solutions using gnuplot.

!! 5.	Read Chapter 4, $4.1-4.6 of Ellis.

!! End of file: subroutine7.f90
