%------------------------------------------------------------
%     This is a MATLAB M-file for the solution of a beam on
%     an elastic foundation. 
%
%     The user must provide the following m-files, 
%        BeamData.m   - Provides beam data
%        loading.m    - Provides beam loading
%     Examples of these files is provided.
%
%     CE560 Advanced Mechanics of Solids
%     Colorado State University
%     Fort Collins, Colorado
%     Erik Thompson
%     November 9, 2003
%------------------------------------------------
       clear
%      --------------
%      Get Beam Data
%      --------------
       BeamData
       
       dx=rL/(INNI-1) ;
       rkEI=rk/EI ;

       x=-dx ;
       for i=1:INNI  ;
           x=x+dx ;
	   Y(i)=0.0;
           XP(i)=x ; 
	   BM(i)=0.0 ;
	   VS(i)=0.0 ;
       end ;

%      ------------------------------
%      Make stencil
%      ------------------------------
       p(1)=+1.0 ;
       p(2)=-4.0 ;
       p(3)=+6.0 + rkEI*(dx^4) ;
       p(4)=-4.0 ;
       p(5)=+1.0 ;


%      ------------------------------
%      Create Matrix
%      ------------------------------
       for i=1:INNI ;
	  for j=1:5 ;
	     c=i-3+j;
	     if c > 0
		if c<= INNI
	           A(i,c)=p(j) ;
		end
	     end
          end ;
          f(i)=0.0 ;
       end ;

%      ------------------------------
%      Determine RHS (loads)
%      nL = number of Load points
%      ------------------------------
       x=-dx ;
       c=(dx^4)/EI ;
       nL=INNI ;
       for i=1:INNI ;
          x=x+dx ;
 	  f(i)= -loading(x,rL,i,nL)*c ;
       end ;
       wL = f(1)/c;
       wR = f(INNI)/c;

%      ------------------------------
%      Boundary Conditions
%
%      Left hand end of beam.
%      ------------------------------
       if Lbc(1) ==  1 ;
%         -------------------
%         deflection is known
%         -------------------
	  A(1,1)=A(1,1)*1.0e+06 ;
	  f(1)=A(1,1)*uL ;
	  if Lbc(2) ==  1 ;
%            --------------
%            slope is known
%            --------------
	     A(2,2)=A(2,2)+1.0 ;
	     f(2)=f(2)+sL*(2.0*dx) ;
	  else ;
%            ---------------
%            moment is known
%            ---------------
             A(2,2)=A(2,2)-1.0 ;
	     f(2)=f(2) - 2.0*uL - mL*(dx^2)/EI ;
	  end ;
       else ;
%         -------------------
%         Shear is known
%         -------------------
	  if Lbc(2) ==  1 ;
%            --------------
%            slope is known
%            --------------
             A(1,2)=-8.0 ;
	     A(1,3)= 2.0 ;
	     A(2,2)= A(2,2)+1.0 ;
             f(1) = f(1) + vL*(2.0*dx^3/EI) - 2.0*sL*(2.0*dx)  ;
             f(2) = f(2) + sL*(2.0*dx)  ;
	  else ;
%            ---------------
%            moment is known
%            ---------------
             A(1,1)=A(1,1)-4.0 ;
	     A(1,2)=-4.0 ;
	     A(1,3)= 2.0 ;
	     A(2,2)=A(2,2)-1.0 ;
	     A(2,1)=-2.0 ;
	     f(1) = f(1) + vL*(2.0*dx^3/EI) + 2.0*mL*(dx^2/EI) ;
	     f(2) = f(2) - mL*(dx^2/EI) ;
	  end ;
       end ;

%      ------------------------------
%      Right hand end of beam.
%      ------------------------------
       Nf = INNI;
       if Rbc(1) ==  1 ;
%         -------------------
%         deflection is known
%         -------------------
	  A(Nf,Nf)=A(Nf,Nf)*1.0e+06 ;
	  f(Nf)=A(Nf,Nf)*uR ;
	  if Rbc(2) ==  1 ;
%            --------------
%            slope is known
%            --------------
	     A(Nf-1,Nf-1)=A(Nf-1,Nf-1)+1.0 ;
	     f(Nf-1)=f(Nf-1)-sR*(2.0*dx) ;
	  else ;
%            ---------------
%            moment is known
%            ---------------
             A(Nf-1,Nf-1)=A(Nf-1,Nf-1)-1.0 ;
	     f(Nf-1)=f(Nf-1) - 2.0*uR - mR*(dx^2)/EI ;
	  end ;
       else ;
%         -------------------
%         Shear is known
%         -------------------
	  if Rbc(2) ==  1 ;
%            --------------
%            slope is known
%            --------------
             A(Nf,Nf-1)=-8.0 ;
	     A(Nf,Nf-2)= 2.0 ;
	     A(Nf-1,Nf-1)= A(Nf-1,Nf-1)+1.0 ;
             f(Nf) = f(Nf) - vR*(2.0*dx^3/EI) + 2.0*sR*(2.0*dx)  ;
             f(Nf-1) = f(Nf-1) + sR*(2.0*dx)  ;
	  else ;
%            ---------------
%            moment is known
%            ---------------
             A(Nf,Nf)=A(Nf,Nf)-4.0 ;
	     A(Nf,Nf-1)=-4.0 ;
	     A(Nf,Nf-2)= 2.0 ;
	     A(Nf-1,Nf-1)=A(Nf-1,Nf-1)-1.0 ;
	     A(Nf-1,Nf)=-2.0 ;
	     f(Nf) = f(Nf) + vR*(2.0*dx^3/EI) + 2.0*mR*(dx^2/EI) ;
	     f(Nf-1) = f(Nf-1) - mR*(dx^2/EI) ;
	  end ;
       end ;
	      
%      ----------------
%      Solve equations 
%      ----------------
       Y = A\f';

%      --------------------------------
%      Calculate values at ghost points
%      Left side
%      --------------------------------
       if Lbc(2) == 1
          Ym1 = Y(2) - 2*dx*sL;
       else
	  Ym1   = (dx^2)*(mL/EI) + 2*Y(1) - Y(2);
       end
       if Lbc(1) == 1
          Ym2 = (dx^4)*(wL/EI) + 4*Ym1 - 6*Y(1) + 4*Y(2) -1*Y(3);
       else
          Ym2 = -2*(dx^3)*(vL/EI) + 2*Ym1 - 2*Y(2) + Y(3);
       end
%      --------------------------------
%      Calculate values at ghost points
%      Right side
%      --------------------------------
       if Rbc(2) == 1
          Yp1 = Y(INNI-1) + 2*dx*sR;
       else
	  Yp1   = (dx^2)*(mR/EI) + 2*Y(INNI) - Y(INNI-1);
       end
       if Rbc(1) == 1
          Yp2 = (dx^4)*(wR/EI) + 4*Yp1 - 6*Y(INNI) + 4*Y(INNI-1) -1*Y(INNI-2);
       else
          Yp2 =  2*(dx^3)*(vR/EI) + 2*Yp1 - 2*Y(INNI-1) + Y(INNI-2);
       end

%      --------------
%      Bending Moment
%      --------------
       ibgn=2 ;
       iend=INNI-1 ;
       fac=EI/(dx^2) ;
       for i=2:iend ;
	  BM(i) = (Y(i-1)-2.0*Y(i)+Y(i+1))*fac ;
       end ;
       BM(1) = (Ym1-2*Y(1)+Y(2))*fac;
       BM(INNI) = (Yp1-2*Y(INNI)+Y(INNI-1))*fac;

%      ----------- 
%      Shear 
%      -----------
       ibgn=3 ;
       iend=INNI-2 ;
       fac=EI/(2.0*dx^3) ;
       for i=3:iend ;
	  VS(i) = (-Y(i-2)+2.0*Y(i-1)-2.0*Y(i+1)+Y(i+2))*fac ;
       end ;
       VS(1) = (-Ym2+2.0*Ym1-2.0*Y(2)+Y(3))*fac ;
       VS(2) = (-Ym1+2.0*Y(1)-2.0*Y(3)+Y(4))*fac ;
       VS(INNI)   = (-Y(INNI-2)+2.0*Y(INNI-1)-2.0*Yp1+Yp2)*fac ;
       VS(INNI-1) = (-Y(INNI-3)+2.0*Y(INNI-2)-2.0*Y(INNI)+Yp1)*fac ;

%      ---------------------
%      Prepare for plotting
%      ---------------------
       EId = ['EI = ' num2str(EI)] ;
       RLd = ['Length = ' num2str(rL)];
       RKd = ['Foundation stiffness = ' num2str(rk)] ;
       Yxd = ['Max. Deflection = ' num2str(max(Y))] ;
       Ynd = ['Min. Deflection = ' num2str(min(Y))] ;
       Vxd = ['Max. Shear  = ' num2str(max(VS))] ;
       Vnd = ['Min. Shear  = ' num2str(min(VS))] ;
       Mxd = ['Max. Moment  = ' num2str(max(BM))] ;
       Mnd = ['Min. Moment  = ' num2str(min(BM))] ;

%      ------------------
%      Plot results
%      ------------------
       subplot(2,2,1) ;
       plot(XP(1:INNI),Y(1:INNI)) ;
       grid
       ylabel( ' Deflection ')
       xlabel( ' Distance along beam ')
       subplot(2,2,2) ;
       plot(XP(1:INNI),VS(1:INNI)) ;
       grid
       ylabel( ' Beam Shear ')
       xlabel( ' Distance along beam ')
       subplot(2,2,3) ;
       plot(XP(1:INNI),BM(1:INNI));
       grid
       ylabel( ' Bending Moment ');
       xlabel( ' Distance along beam ');
       subplot(2,2,4);
       brdrX = [0 100 100   0   0];
       brdrY = [0   0 100 100   0];
       plot(brdrX,brdrY)
       text(10,90, RLd ) ;
       text(10,80, EId ) ;
       text(10,70, RKd ) ;
       text(10,60, Yxd ) ;
       text(10,50, Ynd ) ;
       text(10,40, Vxd ) ;
       text(10,30, Vnd ) ;
       text(10,20, Mxd ) ;
       text(10,10, Mnd ) ;
       axis off ;
