-- beam2.adb  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.

with Ada.Text_IO; use Ada.Text_IO;
with Real_Arrays; use Real_Arrays;
with Simeq;

procedure Beam2 is
  E : constant := 10.0; -- arbitrary constants,
  I : constant :=  9.0; -- values and units can be found
  W : constant :=  7.0;
  L : constant :=  8.0;
  x, y, h, slope : long_float;
  n : Integer;


  procedure Difeq is
    -- difference equation solution
    n : constant := 7; -- number of equations, zeros on front and back
    A : Real_Matrix(1..n,1..n);
    U : Real_Vector(1..n);    -- A U := F  F is p(x)
    F : Real_Vector(1..n);   -- possibly corrected for boundaries
    C : Real_Vector(1..5) := (1.0, -4.0, 6.0, -4.0, 1.0);
    x, y, h : Long_Float;
    ii, j : Integer;
  begin

    --  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(-25y(x) +48y(x+h) -36y(x+2h) +16y(x+3h) -3y(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

    h := L/Long_float(n+1); -- must count h"s for ends
    Put_Line("Fourth order difference equation solution");
    Put_Line("boundary dy/dx:=0 fourth order");
    Put_Line("n="&Integer'Image(n)&", h="&Long_Float'Image(h));
    -- C(1), C(2), C(3), C(4), C(5)
    for i in 1..n loop
      for j in 1..n loop
        A(i,j) := 0.0;
      end loop;
    end loop;
    for ii in  1..n loop
      F(ii) := -W*h*h*h*h/(L*E*I);
    end loop;
    -- 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;   -- Dirchilet
    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;    -- Dirchilet
    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

    for i in 3..n-2 loop   -- all inner rows
      for j in i-2..i-2+4 loop
        A(i,j) := C(j-i+3);
      end loop;
    end loop;

    -- debug print of simultaneous equations to be solved
    Put_Line("solve A y := F for y ");
    Put_Line("                           A                            F ");
    for i in 1..n loop
      for j in 1..n loop
        Put_Line(Long_Float'Image(A(i,j)));
      end loop;
      Put_Line(Long_Float'Image(F(i)));
    end loop;


    U := simeq(A, F);
    Put_Line("fourth order difference equations are exact within");
    Put_Line("numerical accuracy when solution is fourth order.");
    Put_Line("x:=0.0  y:=0.0  boundary");
    for j in 1..n loop
      x := Long_float(j)*H;
      y := U(j);
      Put_Line("x:="&Long_Float'Image(x)&" y:= "&Long_Float'Image(y));
    end loop;
    Put_Line("x:="&Long_Float'image(Long_Float(n+1)*h)&"  y:=0.0  boundary");
  end Difeq;

begin

  Put_Line("beam2.adb running");
  n := 9;
  h := L/Long_float(n-1);
  Put_Line("E="&Long_Float'Image(E)&" I="&Long_Float'Image(I)&
           " W="&Long_Float'Image(W));
  Put_Line("L="&Long_Float'Image(L)&" h:="&Long_Float'Image(h)&
           " n:="&integer'Image(n));
  Put_Line("Analytic Solution.");
  for j in 0..n-1 loop
    x := Long_Float(j)*h;
    y := (-W* x*x *(x-L)*(x-L))/(24.0*L*E*I);
    Put_Line("beam x:="&Long_Float'Image(X)&"  y:="&Long_Float'Image(Y));
    slope := ((-W*x*x*x)/(6.0*L) + (W*x*x)/4.0 - (W*L*x)/12.0)/(E*I);
    Put_Line("                       slope:="&Long_Float'Image(Slope));
  end loop;

  Put_Line("Numerical Solution.");
  Difeq; -- above
end Beam2;