! beam2.f90  ordinary differential equation for a loaded beam 
!            handle Neumann boundary  dy/dx=0 at x=0 and x=L  
!
!  given a beam of length L, from 0 < x < L
!  with Young's Modulus of E
!  with moment of inertia I
!  with p(x) being the load density  e.g. force per unit length
!  with both ends fixed, meaning:
!       y=0 at x=0, y=0 at x=L, dy/dx=0 at x=0, dy/dx=0 at x=L
!  then the differential equation that defines the y position of
!  the beam at every x coordinate is
!
!      d^4 y
!   EI ----- = p(x)   with the convention for downward is negative
!      dx^4
!
!  for uniformly distributed force (weight) p(x) becomes  - W/L
!  This simple case can be integrated and solved analytically:
!
!      d^4 y 
!   EI ----- = -W/L
!      dx^4
!
!     d^3 y
!  EI ----- = -W/L x + A   ( A is constant of integration, value found later)
!     dx^3
!
!     d^2 y
!  EI ----- = -W/L x^2/2 + A x + B 
!     dx^2
!
!     dy
!  EI -- = -W/L x^3/6 + A x^2/2 + B x + C  we see C=0 because dy/dx=0 at x=0
!     dx
!
!  EI y = -W/L x^4/24 + A x^3/6 + B x^2/2 + D  we see D=0 because y=0 at x=0
!
!  substituting y=0 at x=L in the equation above gives
!
!     0 = -W/L L^4/24 + A L^3/6 + B L^2/2
!
!  substituting dy/dx=0 at x=L in the equation above the above gives
!
!     0 = -W/L L^3/6 + A L^2/2 + B L
!
!  solving two equations in two unknowns  A = W/2  B = - WL/12
!  then substituting for A and B in  EI y = ... and combining terms
!
!         W
!  y = --------  x^2 (x-L)^2
!      24 L EI
!
!  The known solution for a specific case is valuable to check your
!  programming of a numerical solution before computing a more
!  general case of p(x) where no closed form solution may exists.

module stuff
  double precision, parameter :: E = 10.0 ! arbitrary constants,
  double precision, parameter :: I =  9.0 ! values and units can be found
  double precision, parameter :: W =  7.0
  double precision, parameter :: L =  8.0
end module stuff

program beam2
  use stuff
  double precision :: x, y, h, slope
  integer :: j, n

  interface
    subroutine difeq()
    end subroutine difeq
  end interface

  print *, 'beam2.f90 running'
  n = 9
  h = L/(n-1)
  print *, 'E=', E, ' I=', I, ' W=', W
  print *, 'L=', L, ' h=', h, ' n=', n
  print *, 'Analytic Solution.'
  do j=0,n-1
    x = j*h
    y = (-W* x*x *(x-L)*(x-L))/(24.0*L*E*I)
    print *, 'beam x=', x, '  y=', y
    slope = ((-W*x*x*x)/(6.0*L) + (W*x*x)/4.0 - (W*L*x)/12.0)/(E*I)
    print *, '                       slope=', slope 
  end do

  print *, 'Numerical Solution.'
  call difeq() ! below 
end program beam2 


subroutine difeq()
  use stuff
  ! difference equation solution  
  double precision, dimension(100,100) :: A
  double precision, dimension(100) :: U  ! A U = F  F is p(x) 
  double precision, dimension(100) :: F  ! possibly corrected for boundaries 
  double precision, dimension(5) :: C
  data c /1.0d0, -4.0d0, 6.0d0, -4.0d0, 1.0d0/
  double precision :: x, y, h
  integer :: ii, j, n

  interface
    ! simeq.interface f90  solve simultaneous equations AX=Y 
    subroutine simeq(n, sz, A, Y, X)
      integer, intent(in) :: n  ! number of equations
      integer, intent(in) :: sz ! dimension of arrays
      double precision, dimension(sz,sz), intent(in) :: A
      double precision, dimension(sz), intent(in) :: Y
      double precision, dimension(sz), intent(out) :: X
    end subroutine simeq
  end interface

  !  d^4 y / dx^4 as a difference equation =
  !  1/h^4 ( 1 y(x-2h) - 4 y(x-h) + 6 y(x) - 4 y(x+h) + 1 y(x+2h) )
  !  1/h^4 moved to other side of equations (inner)
  !
  !  dy/dx at x = 1/12h(-25 y(x) +48 y(x+h) -36 y(x+2h) +16 y(x+3h) -3 y(x+4h))
  !  multiply both sides by 12h, but constant side is zero for this problem
  !  y(0) = 0 for this problem thus no constant contribution
  !  flip for x=L  -3, +16, -36, +48   y(L)=0 thus no -25
   
  n = 7 ! number of equations, zeros on front and back 
  h = L/(n+1) ! must count h's for ends 
  print *, 'Fourth order difference equation solution '
  print *, 'boundary dy/dx=0 fourth order'
  print '("n=",I3,", h=",E10.5,", C= ",5F6.2)', &
           n, h, C(1), C(2), C(3), C(4), C(5)
  do ii=1,n
    do j=1,n
      A(ii,j) = 0.0
    end do
  end do
  do ii=1,n
    F(ii) = -W*h*h*h*h/(L*E*I)
  end do
  ! handle first two and last two rows with boundary conditions 
  ii=1;   j=1;   A(ii,j) =  48.0
          j=2;   A(ii,j) = -36.0
          j=3;   A(ii,j) =  16.0
          j=4;   A(ii,j) =  -3.0
               F(ii)     =  0.0*12.0*h ! Neumann 
  ii=2;   j=1;   A(ii,j) = -4.0
          j=2;   A(ii,j) =  6.0
          j=3;   A(ii,j) = -4.0
          j=4;   A(ii,j) =  1.0    ! Dirichlet 
  ii=n-1; j=n-3; A(ii,j) =  1.0
          j=n-2; A(ii,j) = -4.0
          j=n-1; A(ii,j) =  6.0
          j=n;   A(ii,j) = -4.0    ! Dirichlet 
  ii=n;   j=n-3; A(ii,j) =  -3.0
          j=n-2; A(ii,j) =  16.0
          j=n-1; A(ii,j) = -36.0
          j=n;   A(ii,j) =  48.0
	       F(ii)     =  0.0*12.0*h ! Neumann 

  do ii=3,n-2   ! all inner rows 
    do j=ii-2, ii-2+4
      A(ii,j) = C(j-ii+3)
    end do
  end do

  ! debug print of simultaneous equations to be solved 
    print *, 'solve A y = F for y '
    print *, '                           A                            F '
    do ii=1,n
      do j=1,n
        print '(F7.3 )', A(ii,j)
      end do
      print *, F(ii)
    end do


  call simeq(n, 100, A, F, U)
  print *, 'fourth order difference equations are exact within'
  print *, 'numerical accuracy when solution is fourth order.'
  print *, 'x=0.0  y=0.0  boundary'
  do j=1,n
    x = (j)*h
    y = U(j)
    print *, 'x=', x,' y= ', y
  end do
  print *, 'x=', (n+1)*h, '  y=0.0  boundary'

end subroutine difeq