! 3. Arithmetic Operations.				File: arithmetic3.f90
!------------------------------------------------------------------------------!

PROGRAM arithmetic
	IMPLICIT NONE

	! This session discusses how the computer does basic arithmetic 
	! with integer and real numbers.
	
	! Example: Physical constants in SI Units	

	! Variables:
	!	c		: speed of light
	!	e		: elementary charge (Coulomb)
	!	k		: Boltzmann constant (JK-1)
	!	m		: electron mass	
	!	h		: Planck constant (Js)
	!	eps0		: Permittivity of free space (Fm-1)
	!	Rydberg		: Rydberg constant (m-1)
	!	Bohradius	: Bohr radius
	!	el_radius	: classical electron radius
	!	lambda_ev	: wavelength associated with 1eV
	!	freq_ev		: frequency associated with 1eV 
	!	energy_ev	: energy associated with 1eV 
	!	temp_ev		: temperature associated with 1eV (K)

	! Variable declarations
	REAL, PARAMETER :: pi=3.1415926,c=2.9979e8,e=1.6022e-19
	REAL, PARAMETER :: k=1.3807E-23,m=9.1095e-31,h=6.6262e-34
	REAL, PARAMETER :: eps0=8.8542e-12
	REAL :: el_radius,lambda_ev,freq_ev,energy_ev,temp_ev
!	REAL :: Rydberg,Bohradius
	REAL :: Rydberg_2,Bohradius_2


       	! Evaluate derived constants

!	Rydberg=m*e**4/(8.0*(eps0**2)*c*(h**3))	!as per definition
	Rydberg_2=(m/h)*(e/eps0**2)*(e/h)**2*(e/c)/8.0 !alternative grouping
!	Bohradius=eps0*h**2/(pi*m*e**2) !as per definition
	Bohradius_2=(h/m)*(h/e)*(eps0/e)/pi !alternative grouping of terms
	el_radius=(e/(4.0*pi*eps0*(m*c**2)))*e
	lambda_ev=h*c/e
	freq_ev=e/h
	energy_ev=e
	temp_ev=e/k

	! Output
	PRINT *, 'Speed of light             =', c,  '  m per s'
	PRINT *, 'Boltzmann constant         =', k, '  J per K'
	PRINT *, 'Elementary charge          =', e,  '  C'
	PRINT *, 'Electron mass              =', m, '  kg'
	PRINT *, 'Planck constant            =', h, '  J per s'
	PRINT *, 'Permittivity of free space =', eps0, '  F per m'
	PRINT *
!	PRINT *, 'Rydberg constant          =',Rydberg, '  per m'
!	PRINT *, 'Bohr radius               =',Bohradius, '  m'
	PRINT *, 'Rydberg constant          =',Rydberg_2, '  per m'
	PRINT *, 'Bohr radius               =',Bohradius_2, '  m'
	PRINT *, 'Classical electron radius =',el_radius, '  m'
	PRINT *
	PRINT *, 'Wavelength associated with 1eV  =',lambda_ev, '  m'
	PRINT *, 'Frequency associated with 1eV   =',freq_ev, '  per s'
	PRINT *, 'Energy associated with 1eV      =',energy_ev, '  J'
	PRINT *, 'Temperature associated with 1eV =',temp_ev, '  K'

END PROGRAM arithmetic

!! 1.	We confine attention here to the two numeric data types REAL and 
!!	INTEGER. The arithmetic operators are:
!!
!!		+	-	*	/	**	=
!!
!!	These are, respectively, addition, subtraction, multiplication, 
!!	division, exponentiation, and the assignment operator.
!!
!!	In fact, there is a set of such operators for each of the numeric data 
!!	types. Although the symbols used are the same in each case, the actual 
!!	operations carried out by the computer differ with the type of data 
!!	being operated on. 
!!
!!	Arithmetic expressions can be constructed from the data types and the 
!!	operators +,-,*,/,** much as in ordinary algebra. However, Fortran is 
!!	only a 'pseudo-mathematical' language. You should appreciate in 
!!	expressions just how the computer goes about reducing it from some
!!	complex form involving many operators and possibly different data 
!!	types to a final variable of a particular data type.
!!
!!	Points to note are:
!!
!!	# all operators must be shown explicitly
!!	# operators can never appear next to one another
!!	# parenthesis should be used in expressions to enforce clarity or 
!!	avoid ambiguity [e.g. a*b+c is different from a*(b+c)  ]
!!	# when the sequence of operations is not clarified by parenthesis, the
!!	'strength' of operations is
!!			- all exponentiations performed first
!!			- next all multiplications and divisions
!!			- finally, all additions and subtractions
!!
!!	# within a sequence of multiplications and/or divisions, or additions 
!!	and/or subtractions, again not fully specified by parenthesis, the 
!!	operations are carried out LEFT to RIGHT.
!!		e.g. a/b*c  means  (a/b)*c  and not  a/(b*c)
!!		e.g. i-j+k  means  (i-j)+k  and not  i-(j+k)
!!	# successive exponentiations, not specified by parenthesis, are 
!!	performed RIGHT to LEFT.
!!		e.g. a**b**c  means  a**(b**c)  and not  (a**b)**c
!!	Bear these points in mind. BUT, remember that you want your program 
!!	to be easily and unambiguously readable by yourself AND others.
!!
!!	A note on exponentiations: Fortran distinguishes between integer powers
!!	and real powers. Integer powers are regarded as repeated multiplication
!!	(e.g. x**3 is evaluated as x*x*x). Real powers, such as x**0.3, are 
!!	evaluated by the algorithm
!!				x**0.3 = antilog(0.3*log(x))
!!	Thus, x**2 and x**2.0 are not the same. They are evaluated in different
!!	ways, the latter being more time consuming and less accurate.

!! 2.	The operation type
!!
!!	It is important to know what data type is produced by each operation.
!!	The type of value produced by an operation is determined by the type 
!!	of its operands.
!!
!!	Operations take two forms:
!!
!!	operand1  	operator	operand2	(binary operator)
!!			operator	operand		(unary operator)
!!
!!	examples of the former:		x+y, 2.0/3.0, x**3
!!	examples of the latter:		-count, +3.0
!!	These are evaluated like the binary operations (0-count) and (0.0+3.0)
!!	respectively. Therefore, all arithmetic expressions may be thought of
!!	as binary operations.	


!! 3.	Integer arithmetic.
!!
!!	Integer addition/subtraction/multiplication are performed exactly 
!!	provided at no stage in a calculation an integer becomes larger than 
!!	the computer can handle. 
!!
!!	Integer division is performed as if the numbers were real but with 
!!	any fractional part being discarded: the result must also be an 
!!	integer. This leads to some surprising results if one is unaware 
!!	as to how the computer treats integer division. For example, you 
!!	might expect the expression 2/3 to have the value 0.666667. However, 
!!	2/3 has the value 0, rather unexpectedly. The type of value produced 
!!	by an operation is determined not only by the nature of the operation 
!!	but also by the type of its operands. It is important to know what 
!!	data type is produced by each operation. 
!!
!!	The program recognises 2 and 3 as integers and therefore the / in
!!	2/3 as meaning integer division. Therefore 2/3=0 
!!	If the expression had been 2.0/3.0 then the program recognises 2.0 
!!	and 3.0 as real numbers and / as real number division. The answer 
!!	then would be 2.0/3.0=0.666667
!!
!!	Examples: 	12/8=1		12.0/8.0=1.50000
!!			257/258=0	257.0/258.0=0.996124
!!			5/3*6 = 6	5.0/3.0*6.0 = 10.0
!!	but note:	5*6/3 = 6*5/3 = 5*(6/3) = (6/3)*5 = 10  as it should.
!!	Parentheses avoid ambiguity!
!!
!!	Exercise:  Write a short program to demonstrate the truth of the 
!!	above. Read in integers int_1,int_2 and print out the integer 
!!	int_1/int_2. Integers may be converted to real form by use of the 
!!	'intrinsic function' REAL (of which more later) used as follows:
!!
!!			r_1=REAL(int_1)
!!			r_2=REAL(int_2)
!!
!!	This converts integer int_1 to real r_1  (e.g. 4  to  4.0)
!!	Print out r_1,r_2, and r_1/r_2 for comparison.

!! 4.	Mixed-mode expressions
!!
!!	The computer requires the operands in binary arithmetic operations 
!!	+-*/ to be of the same type. If operands of a given operator are of 
!!	different types, Fortran 'promotes' one to be of the same type as the 
!!	other. Fortran has a 'hierarchy' of types:
!!
!!	# COMPLEX is higher than REAL 
!!	# REAL is higher than INTEGER
!!
!!	Involved arithmetic expressions can be broken down into a sequence of 
!!	sub-expressions which are binary operations. If such sub-expressions 
!!	are of mixed type, the rule is that the operand lower in the hierarchy
!!	gets automatically promoted to the higher of the two types before the 
!!	operation is carried out. For example, in 5/4.0, first the integer 5 is
!!	promoted to the real 5.0 then the division is performed as real division
!!	to obtain 1.250000 .
!!
!!	Examples:
!!
!!	expression		form			value/type
!!	__________		____			__________
!!
!!	3*4.0			integer*real		12.0 (real)
!!	2.0e2 + 1		real+integer		2.01e2 (real)
!!	3/2.0e2			integer/real		1.5e-2 (real)
!!	(0.0,1.0) - 2.0		complex-real		(-2.0,1.0) (complex)
!!	2 + (2.0,1.0)		integer+complex		(4.0,1.0) (complex)
!!
!!	BE VERY CAREFUL IN EVALUATING MIXED-MODE EXPRESSIONS!!!!!
!!
!!	E.g  	4 + 2/4.0 - 3.0*(4 - 6/5) = -4.50000
!!		4 + 2/4 - 3*(4.0 -6.0/5) = -4.40000
!!	Make sure you understand the differences in the above two expressions.

!! 5.	The assignment operator =
!!
!!	The syntax for this is
!!
!!			variable_name = expression
!!
!!	It is not the same as the algebraic equality sign. First the 
!!	expression is evaluated and then its value then placed in the 
!!	storage cell called 'variable_name'. 
!!
!!	To emphasize this, consider the Fortran statement
!!
!!			count = count + 1
!!
!!	where count is an integer. This is clearly invalid in algebra, 
!!	but quite in order in Fortran. It states: take the present value 
!!	stored in count then add 1 to it, finally placing this new value 
!!	in count and replacing the previous value.
!!
!!	The expression on the right hand side of the assignment will always 
!!	be evaluated as given by the above hierarchy of types.
!!	If this value is of different type from the type which 'variable_name' 
!!	is declared to be, then it is automatically converted before it is 
!!	placed in the storage cell.
!!
!!	Remember that only a single variable name may appear on the left of the
!!	assignment operator.

!! 7.	Finite number representations
!!
!!	Real arithmetic:
!!
!!	Physicists are well-used to the idea of limited accuracy in data, with
!!	results which are correct to so many significant figures. Computers 
!!	also suffer from limited accuracy. Any data type has only a finite size
!!	storage cell allocated to it. This means that the ordinary rules of 
!!	mathematics are not obeyed when applied to arithmetic operations with 
!!	real numbers of finite length. For example, consider 
!!
!!		0.400000 + 12345678. - 12345677.
!!
!!	According to the rules this is evaluated left to right. The first sum
!!	will give the result 12345678. , with the 0.4 being entirely lost. The
!!	subtraction of the final number then yields 1.0000000 which is 
!!	incorrect. We have incurred a loss of significance by adding a small 
!!	number to a large one.
!!
!!	Question: What is the result if we enforce, using parenthesis, a 
!!	different order:
!!
!!		0.400000 + (12345678. - 12345677.)
!!
!!	All computers have maximum (positive and negative) real numbers, 
!!	which on our system is typically 1.E38 . Thus, if any calculation 
!!	goes out of range we obtain a 'floating-point overflow' - usually 
!!	indicated on most systems by a 'run-time' error'. This can often be 
!!	avoided by suitable groupings of terms to keep quantities within range.
!!
!!	Exercise: Consider the above program 'arithmetic' . 
!!	Comment out the lines involving  Rydberg_2  and  Bohradius_2 and 
!!	activate those involving Rydberg and Bohradius.  
!!	Compile and execute the program. Explain the difference in output.
!!
!!	To find out the largest and smallest numbers which your machine will
!!	allow, use the intrinsic functions HUGE(x) and TINY(x). Write a 
!!	program with the following lines:
!!
!!		REAL :: u,v
!!		INTEGER :: i
!!
!!		u = HUGE(1.0)
!!		v = TINY(1.0)
!!		i = HUGE(1)
!!
!!	Now have these variables printed out and you will see the largest and 
!!	smallest real and integer numbers represented by your machine. If the 
!!	large values are exceeded on either the positive or negative sides 
!!	then an 'overflow' will occur, your program will cease execution and 
!!	an appropriate message will appear on your terminal. If any real number
!!	is less in magnitude than v=TINY(1.0) an 'underflow' will occur; the 
!!	program will continue to execute using the value zero.
!!	Beware of division by zero: a common source of error! 
!!
!!	[There is another 'small' number which is also important: the 'machine
!!	epsilon'. This is the smallest real positive number which when added 
!!	to 1.0 gives an answer different from 1.0]

!! 8.	The PARAMETER attribute
!!
!!	Rather than repeat a frequently used constant (such as pi=3.1415926)
!!	at various locations in a program with all the attendant likelihood 
!!	of error, it is preferable to define a NAMED CONSTANT by use of the 
!!	PARAMETER ATTRIBUTE. Observe how this is done for the basic physical 
!!	constants in the program above.
!!
!!	Note also that if a certain combination of constants occurs repeatedly 
!!	then it is sensible to declare a named constant to represent it.
!!	For example, physicists often need to use hbar=h/2pi where h is 
!!	Planck's constant. Then we could include the line
!!
!!		REAL, PARAMETER :: pi=3.1415926, h=6.6262e-34,hbar=h/(2.0*pi)
!!
!!	and use hbar rather than having to calculate h/2pi each time it arises.

!! 9.	Exercises 3.
!!
!!	a) Write correct Fortran statements, assuming all quantities real, for:
!!		i) (x**2-y**2)/(x(y-1))
!!		ii) x - x**3/3! + x**5/5!
!!		iii) [3x(z+1) y**2 + yz/2]/(1+x+x**3)
!!		iv) C = 4 Pi Epsilon ab /(b-a)	[capacitance of conc. spheres]
!!
!!	b) What value is stored by the following assignment statements,
!!	assuming X is real and I integer:
!!		i)	X = 12 - (8+7/3)/3
!!		ii)	I = 15/6*6/15
!!		iii)	I = 3.6*(14.0**(3/7))
!!		iv)	I = 1./6. + 1./6. + 1./6. + 1./6. + 1./6. + 1./6.
!!		vi)	Write a program to check your answers. 		
!!
!!	c) Identify the errors in the following assignment statements:
!!		i)	3x = beta - gamma
!!		ii)	3.0 = 9.0/3.0
!!		iii)	3.0 = weight
!!		iv)	zinc = density*-2.0
!!		v)	x=y=z
!!
!!	d) Write a program to evaluate and print out:
!!		i)	y=1.0E+30,  	x=(1.0/y)**10
!!		ii)	y=1.0E10,	x=y+1./y
!!		iii)	y=3.2E+02,	x=(1/2)*y
!!	Explain these answers.

!! 11.	Now read Chapter 3 (omitting $3.7) of the textbook.

!!	End of file: arithmetic3.f90
