c Cantilever beam bending analysis
c
c Usage:
c	beam < inputfile > outputfile
c
c Description:
c This program models a simple cantilever I-beam subjected to a
c static concentrated transverse load at its free end.  The z axis of
c the coordinate system lies along the axis of the beam.  The program
c computes the maximum bending stress and the maximum deflection.

	program main
c	============
c	local variables:
c	h	beam height, mm
c	w	flange width, mm
c	tf	flange thickness, mm
c	tw	web thickness, mm
c	l	beam length, mm
c	fx	load on free end in x direction, N
c	fy	load on free end in y direction, N
c	sfx	max bending stress due to fx, MPa
c	sfy	max bending stress due to fy, MPa
c	s	max stress in z direction stress, MPa
c	ix	2nd area moment of inertia about x axis, mm^4
c	iy	2nd area moment of inertia about y axis, mm^4
c	a	area, mm^2
c	hw	height of web, mm
c	dx	max deflection in x direction, mm
c	dy	max deflection in y direction, mm
c	d	max deflection, mm
c	e	Young's modulus (modulus of elasticity), MPa
c	z	location along z axis of beam, mm (0 at wall)
c	i	loop index
c	ns	number of stations along beam at which to report results

	implicit none
	double precision h,w,tf,tw,l,fx,fy,e
	double precision sfx,sfy,s,ix,iy,a,hw,dx,dy,d,z
	integer i,ns
	character*80 title

c	** set number of stations at which results should be reported
	ns = 10

c	** read input file
	call GetInp(title,w,h,tf,tw,l,fx,fy,e)

c	** compute geometric values
	hw = h - 2*tf
	a = hw*tw + 2*tf*w
	ix = tw*(hw**3)/12 + 2*(w*(tf**3)/12 + w*tf*((h-tf)/2)**2)
	iy = hw*(tw**3)/12 + 2*tf*(w**3)/12

c	** print intermediate results
	write(6,1)
	write(6,2) title
	write(6,3)
	write(6,4)
	write(6,5)
	write(6,6) l,h,w,tw,tf
	write(6,7)
	write(6,8)
	write(6,9)
	write(6,10) a,ix,iy
	write(6,11) fx,fy

c	** maximum stress
c
c	PJM 10-May-2000 ESIis05139  [for Linux]
c	Change '0.0' to '0.0d0' to force arg to be passed as double
c
	call Stress(0.0d0,w,h,l,fx,fy,ix,iy,e,sfx,sfy,s,dx,dy,d)
	write(6,12) s

c	** maximum deflection
	call Stress(l,w,h,l,fx,fy,ix,iy,e,sfx,sfy,s,dx,dy,d)
	write(6,13) dx,dy
	write(6,14) d

c	** intermediate stresses and deflections
	write(6,15)
	write(6,16)
	write(6,17)
	do 100 i = 0, ns
		if (i .eq. 0) then
			z = 0
		else if (i .eq. ns) then
			z = l
		else
			z = i*(l/ns)
		endif
		call Stress(z,w,h,l,fx,fy,ix,iy,e,sfx,sfy,s,dx,dy,d)
		write(6,18) z,s,dx,dy
100	continue

	stop

1	format('BEAM BENDING ANALYSIS')
2	format(/a)
3	format(/'GEOMETRY')
4	format(7x,'LENGTH',7x,'HEIGHT',7x,'WIDTH ',7x,'  WEB ',
     +		7x,'FLANGE')
5	format(7x,' (MM) ',7x,' (MM) ',7x,' (MM) ',7x,' (MM) ',
     +		7x,' (MM) ')
6	format(5(1x,f12.3))
7	format(/'SECTION PROPERTIES')
8	format('        AREA           IXX             IYY')
9	format('       (MM^2)         (MM^4)          (MM^4)')
10	format(f15.6,2(x,f15.6))
11	format(/'LOADS AT FREE END (N): ',2(x,f12.3))
12	format(/'MAXIMUM EFFECTIVE STRESS (MPA):',x,f15.6)
13	format('DEFLECTION AT FREE END (MM):   ',2(x,f15.6))
14	format('MAXIMUM DEFLECTION (MM):       ',x,f15.6)
15	format(/'STRESSES AND DEFLECTIONS')
16	format(/'        Z          STRESS        DX           DY')
17	format('       (MM)         (MPA)       (MM)         (MM)')
18	format(f12.3,x,f12.3,2(x,f12.6))

	end

	subroutine GetInp(title,w,h,tf,tw,l,fx,fy,e)
c	============================================
	implicit none
	character*(*) title
	double precision w,h,tf,tw,l,fx,fy,e

	read(5,1,err=999) title
	read(5,2,err=999) l,w,h,tf,tw
	read(5,3,err=999) fx,fy
	read(5,4,err=999) e

	return

999	write(6,'(''Input error'')')
	stop

1	format(a80)
2	format(5f15.8)
3	format(2f15.8)
4	format(f15.0)

	end

	subroutine Stress(z,w,h,l,fx,fy,ix,iy,e,sfx,sfy,s,dx,dy,d)
c	==========================================================
	implicit none
	double precision z,w,h,l,fx,fy,ix,iy,e
	double precision sfx,sfy,s,dx,dy,d

c	** compute deflection components at free end and stress components
c	** at wall using elementary method
	if (iy .gt. 0) then
		dx = fx*(z**2)*(3*l-z)/(6*e*iy)
		sfx = fx*(l-z)*(w/2)/iy
	else
		write(6,'(''Iyy .eq. 0; stress & deflection set to 0'')')
		dx = 0
		sfx = 0
	endif
	if (ix .gt. 0) then
		dy = fy*(z**2)*(3*l-z)/(6*e*ix)
		sfy = fy*(l-z)*(h/2)/ix
	else
		write(6,'(''Ixx .eq. 0; stress & deflection set to 0'')')
		dy = 0
		sfy = 0
	endif

c	** max deflection is Euclidean combination of deflections
	d = dsqrt(dx*dx + dy*dy)

c	** max bending stress is sum of max bending stresses
c	** note that all bending stresses are normal to cross section
	s = dabs(sfx) + dabs(sfy)

	return
	end
