PROGRAM Beam_Deflection
!-----------------------------------------------------------------------
! This program calculates deflections in a cantilever beam
! under a given load. Three different loading conditions are
! analyzed:
!   Case #1 : a single point load at the free end of the beam
!   Case #2 : a single point load at an interior point
!   Case #3 : load uniformly distributed along the beam
!
! Variables used are:
!   Elasticity : modulus of elasticity (Pa)
!   Inertia    : moment of inertia (m**4)
!   BeamLength : length of the beam (m)
!   NumPoints  : number of divisions at which deflections are found
!   Delta_x    : distance between points (m)
!   Option     : option selected by the user
!   Load       : load on the beam (N)
!   x          : distance from the free end of the beam (m)
!   a          : distances along the beam (m)
!   I          : control variable used in DO loop
!   Deflection : deflection at distance x (m)
!   Case_1     : internal function to compute deflection in Case #1
!   Case_2     :     "        "    "     "         "      " Case #2
!   Case_3     :     "        "    "     "         "      " Case #3
!
! Input:  Elasticity, Inertia, BeamLength, NumPoints, Option, Load
! Output: A table of deflections at points along the beam
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: NumPoints, Option, I
  REAL :: Elasticity, Inertia, BeamLength, Delta_x, Load, x, a, &
          Deflection
          

  ! Get beam information
  PRINT *, "Enter"
  WRITE (*, '(6X, A)', ADVANCE = "NO") "Modulus of elasticity (Pa): "
  READ *, Elasticity
  WRITE (*, '(6X, A)', ADVANCE = "NO") "Moment of inertia (m**4): "
  READ *, Inertia
  WRITE (*, '(6X, A)', ADVANCE = "NO") "Length of the beam (m): "
  READ *, BeamLength

  ! Repeat the following until user selects option 0 to stop
  DO
     PRINT *
     PRINT *, "The following options are available:"
     PRINT *, "   0 : Stop"
     PRINT *, "   1 : End Load"
     PRINT *, "   2 : Intermediate Load"
     PRINT *, "   3 : Uniform load"
     WRITE (*, '(1X, A)', ADVANCE = "NO") "Select one: "
     READ *, Option
 
     IF (Option < 0 .OR. Option > 3) THEN
        PRINT *, "Not a valid option"
 
     ELSE IF (Option == 0) THEN
        PRINT *, "Done processing for this beam"
        EXIT
 
     ELSE
        PRINT *,  "Enter total load for cases 1 and 2, unit load for case 3."
        WRITE (*, '(1X, A)', ADVANCE = "NO") "Load (N)? "
        READ *, Load
        IF (Option == 2) THEN
           WRITE (*, '(1X, A)', ADVANCE = "NO") "Enter the distance a (m):"
           READ *, a
        END IF
  
        WRITE (*, '(1X, A)', ADVANCE = "NO") "Enter number of points to use: "
        READ *, NumPoints
        Delta_x = BeamLength / REAL(NumPoints)
  
        PRINT *
        PRINT *, "Distance (m)    Deflection (m)"
        PRINT *, "=============================="
  
        x = 0.0
        DO I = 1, NumPoints
           SELECT CASE (Option)
              CASE (1)
                Deflection = &
                    Case_1(Elasticity, Inertia, BeamLength, Load, x)
              CASE (2)
                Deflection = &
                    Case_2(Elasticity, Inertia, BeamLength, Load, a, x)
              CASE (3)
                Deflection = &
                    Case_3(Elasticity, Inertia, BeamLength, Load, x)
           END SELECT
  
           PRINT '(1X, F7.2, E20.6)', x, Deflection
  
           x = x + Delta_x
  
        END DO
  
     END IF
 
  END DO

CONTAINS

  !- Case_1 -----------------------------------------------------------
  ! Computes deflections for a single point load at the free end
  ! of the beam.  Local variable used:
  !   Temp  : temporary variable used in calculating deflection
  !
  ! Accepts: Elasticity, Inertia, BeamLength, Load, x (see main program)
  ! Returns: Deflection
  !---------------------------------------------------------------------
  
  FUNCTION Case_1(Elasticity, Inertia, BeamLength, Load, x)
  
    REAL :: Case_1, Temp
    REAL, INTENT(IN) :: Elasticity, Inertia, BeamLength, Load, x
  
    Temp = x**3 - 3.0 * BeamLength**2 * x + 2.0 * BeamLength**3
    Case_1 = (-Load / (6.0 * Elasticity * Inertia)) * Temp
  
  END FUNCTION Case_1
  
  
  !- Case_2 ------------------------------------------------------------
  ! Computes deflections for a single point load at an interior*
  ! point of the beam.  Local variables used:
  !   b    : distance along the beam
  !   Temp : temporary variable used in calculating deflection
  !          (see main program)
  ! Accepts: Elasticity, Inertia, BeamLength, Load, a, x
  ! Returns: Deflection
  !---------------------------------------------------------------------
  
  FUNCTION Case_2(Elasticity, Inertia, BeamLength, Load, a, x)
  
    REAL :: Case_2, b, Temp
    REAL, INTENT(IN) :: Elasticity, Inertia, BeamLength, Load, a,  x
  
    b = BeamLength - a
    IF (x < b) THEN
       Temp = -a**3 + 3.0 * a**2 * BeamLength - 3.0 * a**2 * x
    ELSE
       Temp = (x - b)**3 - 3.0 * a**2 * (x - b) + 2.0 * a**3
    END IF
    Case_2 = (-Load / (6.0 * Elasticity * Inertia)) * Temp
  
  END FUNCTION Case_2
  
  
  !- Case_3 ------------------------------------------------------------
  ! Computes deflections for a load uniformly distributed along
  ! the beam.  Local variables used:
  !   TotalLoad  : total load
  !   Temp       : temporary variable used in calculating deflection
  !
  ! Accepts: Elasticity, Inertia, BeamLength, Load, x (see main program)
  ! Returns: Deflection
  !---------------------------------------------------------------------
  
  FUNCTION Case_3(Elasticity, Inertia, BeamLength, Load, x)
  
    REAL :: Case_3, Temp, TotalLoad
    REAL, INTENT(IN) :: Elasticity, Inertia, BeamLength, Load, x
  
    TotalLoad = Load * BeamLength
    Temp = x**4 - 4.0 * BeamLength**3 * x + 3.0 * BeamLength**4
    Case_3 = (-TotalLoad / (24.0 * Elasticity * Inertia * BeamLength))&
             * Temp
  
  END FUNCTION Case_3

END PROGRAM Beam_Deflection