%--------------------------                            
% program winkler.m
%--------------------------                            
 clear

%---------------- 
% INPUT DATA
%---------------- 
  load MESH   -ASCII
  load QUAD   -ASCII

% ---------------- 
% Define MESH Data
% ---------------- 
  NUMNP = MESH(1,1); 
  NUMEL = NUMNP-1;
  NUMEQ = 2*NUMNP;
  IB    = 4;
  for i=1:NUMNP;
     XORD(i)   = MESH(i+1,1);
     NPcode(i) = MESH(i+1,2);
  end

% ---------------- 
% Define QUAD Data
% ---------------- 
  NQPTS = QUAD(1,1);
  for i=1:NQPTS
     GPTS(i)=QUAD(i+1,1);
     GWTS(i)=QUAD(i+1,2);
  end

% --------------------------- 
% General initialization
% --------------------------- 
  for I=1:NUMEQ
     LHS(I)=0.0;
     RHS(I)=0.0;
     for J=1:4
        SK(I,J)=0.0; 
     end
  end
  for I=1:NUMNP
     NPBC(I)=0.0;
  end

  INITIAL

   NPBC

% -------------------------------- 
% FORMATION OF STIFFNESS MATRIX   
% and RIGHT HAND SIDE
% -------------------------------- 

  for I=1:NUMEL
% ----------------------------
% Initialize Element Variables
% ----------------------------
    for J=1:4
      for K=1:4
         S1(J,K)=0.0; 
         S2(J,K)=0.0; 
       end
       Qe(J)=0.0; 
    end

% ----------------------------
% Calculate element coordinate
% information.
% ----------------------------
     Xa = XORD(I);
     Xb = XORD(I+1);
     RL=Xb-Xa; 
     DxDu=RL/2.0; 
     DuDx=2.0/RL;   

% ------------------------- 
% Begin Gaussian Quadrature 
% ------------------------- 
  for J=1:NQPTS
     u  = GPTS(J);
     Wt = GWTS(J);
%    --------------------
%    Global coordinate of 
%    current Gauss point
%    --------------------
     Xg = (Xa+Xb)/2  + (RL/2.0)*u;

%    ----------------------------
%    INCLUDE COEF.m 
%    Defines: EIx, Kx, and Qx
%    ----------------------------
     COEF

%    ------------------------------
%    Calculate  shape functions wrt 
%    x at current Gauss point
%    ------------------------------
     No(1) = SF(0,1,u);
     No(2) = SF(0,2,u)*(DxDu);
     No(3) = SF(0,3,u);
     No(4) = SF(0,4,u)*(DxDu);

     Nxx(1) = SF(2,1,u)*(DuDx)^2;
     Nxx(2) = SF(2,2,u)*(DuDx)^2*(DxDu);
     Nxx(3) = SF(2,3,u)*(DuDx)^2;
     Nxx(4) = SF(2,4,u)*(DuDx)^2*(DxDu);


% ----------------------------
% Element stiffness matrices
% ----------------------------
    Wt = GWTS(J);
    for K=1:4
      for L=1:4
        S1(K,L)=S1(K,L) + ...
          Wt*Nxx(K)*EIx*Nxx(L)*DxDu;

        S2(K,L)=S2(K,L) + ...
          Wt*No(K)*Kx*No(L)*DxDu;
      end   
      Qe(K)=Qe(K) + Wt*No(K)*Qx*DxDu; 
   end   
  end   
% --------- Quadrature now complete 


%   --------------------------
%   Assemble element matrices 
%   into global matrix 
%   -------------------------- 
    K1=2*I-2;
    for K=1:4
      K1=K1+1; 
      L1=0; 
      for L=K:4
        L1=L1+1; 
        SK(K1,L1)=SK(K1,L1) + ...
                 S1(K,L)+S2(K,L); 
      end
      RHS(K1) =  RHS(K1) + Qe(K);
    end
   end
%  ---- Global Matrices are assembled


% -------------------------
% BOUNDARY CONDITIONS
% -------------------------
  B = 1.0E+06;
  for I=1:NUMNP
   if NPBC(I) == 1 | NPBC(I) == 3
     I1 = 2*I-1; 
     SK(I1,1)=SK(I1,1)*B;
     RHS(I1)=LHS(I1)*SK(I1,1);
   end
   if NPBC(I) == 2 | NPBC(I) == 3
     I2=2*I;
     SK(I2,1)=SK(I2,1)*B;
     RHS(I2)=LHS(I2)*SK(I2,1);
   end
  end


% -------------------------
% CALL EQUATION SOLVER
% -------------------------
  LHS = sGAUSS(SK,RHS,NUMEQ,IB);

% -------------------------------------- 
% Nodal values for w and dw/dx % are now 
% in LHS.  Use these values to calculate 
% shear and moment at center of each element.
%
% First calculate shape functions
% and their derivatives.  Ic is
% counter for number elements used.
% ------------------------------------ 
  Ic = 0;
  for I=1:NUMEL
    Xa = XORD(I);
    Xb = XORD(I+1);
    RL=Xb-Xa; 
    Xg=(Xa+Xb)/2;
    DuDx=2.0/RL;   
    DxDu=RL/2.0;   

    No(1) = SF(0,1,0);
    No(2) = SF(0,2,0)*(DxDu);
    No(3) = SF(0,3,0);
    No(4) = SF(0,4,0)*(DxDu);

    Nx(1) = SF(1,1,0)*DuDx;;
    Nx(2) = SF(1,2,0)*(DxDu)*DuDx;
    Nx(3) = SF(1,3,0)*DuDx;
    Nx(4) = SF(1,4,0)*(DxDu)*DuDx;

    Nxx(1)   = SF(2,1,0)*DuDx^2;
    Nxx(2)   = SF(2,2,0)*DxDu*DuDx^2;
    Nxx(3)   = SF(2,3,0)*DuDx^2;
    Nxx(4)   = SF(2,4,0)*DxDu*DuDx^2;

    Nxxx(1)   = SF(3,1,0)*DuDx^3;
    Nxxx(2)   = SF(3,2,0)*DxDu*DuDx^3;
    Nxxx(3)   = SF(3,3,0)*DuDx^3;
    Nxxx(4)   = SF(3,4,0)*DxDu*DuDx^3;

% --------------------------- 
%  Calculate deflection, slope
%  shear and moment at center
%  of element.
% --------------------------- 
    COEF

    W=0;
    S=0;
    V=0;
    M=0;
    for J=1:4
       W = W + No(J)*LHS(2*I-2+J);
       S = S + Nx(J)*LHS(2*I-2+J);
       M = M + EIx*Nxx(J)*LHS(2*I-2+J);
       V = V + EIx*Nxxx(J)*LHS(2*I-2+J)...
                +dEIx*Nxx(J)*LHS(2*I-2+J);
    end
   
% --------------------------- 
%  Save center values except 
%  for a pin-element.
% --------------------------- 
    if EIx  > 0
      Ic = Ic + 1;
      ANS(Ic,1)  =  Xg;
      ANS(Ic,2)  =  W;
      ANS(Ic,3)  =  S;
      ANS(Ic,4)  =  M;
      ANS(Ic,5)  =  V;
    end

  end
% --------------------------- 
% Solution is now complete.
% Save values.
% --------------------------- 

  save ANS ANS -ASCII

 NOTE = ...
 ['Solution is saved in file:   '; ...
  '        ANS                  '; ...
  'where the columns  are:      '; ...
  '                             '; ...
  '[ x | w | dw/dx | M | V ]    '; ...
  '                             '; ...
  'To plot, e.g.  w vs x,  type '; ...
  'plot(ANS(:,1),ANS(:,2))      ']


