function MRSMES1DIV(op)
close all
global n L q P E J obc wb1 wb2;
if nargin==0
    op = 1;
end
if op==1
    n = 33;
    L = 2;
    q = 1;
    P = 1;
    E = 1;
    J = 1;    
    obc = 1;
    wb1 = 1;
    wb2 = 1;
end
wyk = 0;

if op==2
    wb1 = 1;
end
if op==3
    wb1 = 2;
end
if op==4
    wb1 = 3;
end
if op==5
    wb1 = 4;
end
if op==6
    wb2 = 1;
end
if op==7
    wb2 = 2;
end
if op==8
    wb2 = 3;
end
if op==9
    wb2 = 4;
end
if op==10
    obc = 1;
end
if op==11
    obc = 2;
end
if op==12
    obc = 3;
end
if op==13
    z = inputdlg('PODAJ LICZBÊ WÊZ£ÓW n = ...','SIATKA',1,{num2str(n)});
    wyb = deal(z{:});
    n = eval(wyb);
end
if op==14
    z = inputdlg('PODAJ MODU£ YOUNG''A E = ...','MATERIA£',1,{num2str(E)});
    wyb = deal(z{:});
    E = eval(wyb);
end
if op==15
    z = inputdlg('PODAJ MODU£ BEZW£ADNOŒCI I = ...','GEOMETRIA',1,{num2str(J)});
    wyb = deal(z{:});
    J = eval(wyb);
end
if op==16
    wyk = 1;
end

h = L/(n-1);
X = 0:h:L;
EJ = E*J;

Mac = [0 0 1 0 0; 1/12 -2/3 0 2/3 -1/12; -1/12 4/3 -5/2 4/3 -1/12; -0.5 1 0 -1 0.5; 1 -4 6 -4 1];
for i=1:5
    Mac(i,:) = Mac(i,:)./h^(i-1);
end
A = zeros(n+4,n+4);
B = zeros(n+4,1);

if obc==3
    nr = 0;
    dmin = L;
    for i=1:n
        d = abs(X(i) - L/2);
        if d<dmin
            dmin = d;
            nr = i;
        end
    end
    nr
end

for i = 3:n+2
    A(i,(i-2):(i+2)) = EJ*Mac(5,:);
    B(i) = fun(X(i-2),L,obc);
end

wb = [1 0 1 0;1 1 0 0;0 1 0 1;0 0 1 1];

il = 0;
for i=1:4
    if wb(wb1,i) == 1
        il = il + 1;
        A(il,1:5) = Mac(i,:);
    end
end
il = n+2;
for i=1:4
    if wb(wb2,i) == 1
        il = il + 1;
        A(il,n:(n+4)) = Mac(i,:);
    end
end
if obc==3
    B(nr+2) = P/h;
end
wyz1 = det(A);
w = A^-1 * B;
hf = figure(1);
plot(X,w(3:n+2),'b-')
hold on
grid on
view(0,-90)

N = n-1;
Nss = 2*n;
bc = zeros(4,2);
il = 0;
for i=1:2
    if wb(wb1,i) == 1
        il = il + 1;
        bc(il,1) = i;
    end
end
for i=1:2
    if wb(wb2,i) == 1
        il = il + 1;
        bc(il,1) = Nss-2+i;
    end
end
bc(il+1:4,:) = [];
KG=zeros(Nss,Nss);
f=zeros(Nss,1);
for i=1:N
    Kel = EJ*[12/h^3 6/h^2 -12/h^3 6/h^2; 6/h^2 4/h -6/h^2 2/h; -12/h^3 -6/h^2 12/h^3 -6/h^2; 6/h^2 2/h -6/h^2 4/h]; 
    xs = 0.5*(X(i) + X(i+1)); 
    fel = fun(xs,L,obc)*[h/2; h^2/12; h/2; -h^2/12];    
    [KG,f] = assem([i 2*i-1 2*i 2*(i+1)-1 2*(i+1)],KG,Kel,f,fel);
end
if obc==3
    f(2*nr-1) = P;
end
[u,r]= solveq(KG,f,bc);
plot(X,u(1:2:Nss),'r-');
legend('MRS','MES');
title(['W^M^R^S_m_a_x = ' num2str(max(w(3:n+2)),6) ', W^M^E^S_m_a_x = ' num2str(max(u(1:2:Nss)),6)])
p = line([0 L],[0 0]);
set(p,'Color','k');
set(p,'LineWidth',[10]);
if wyk==1
    hf2 = figure(2);    
    set(hf2,'Name','WYKRESY SI£ PRZEKROJOWYCH');
    set(hf2,'NumberTitle','off');    
    S = zeros(3,4);
    Pos = zeros(3,4);
    for k=1:3
        sp = subplot(3,2,2*k-1);
        hold on
        grid on
        s = zeros(n,1);
        for i=3:n+2
            s(i-2) = Mac(k+1,:) * w((i-2:i+2))*EJ;
        end
        plot(X,s,'b-');
        for i=1:n
            plot([X(i) X(i)],[0 s(i)],'b-');
        end
        nap = 'K¥T UGIÊCIA';
        if k==2
            nap = 'MOMENT ZGINAJ¥CY';
        end
        if k==3
            nap = 'SI£A POPRZECZNA';
        end
        title(['MRS: ' nap]);
        S(k,:) = axis;
        subplot(3,2,2*k);
        hold on
        grid on
        title(['MES: ' nap]);
     %   Pos(k,:) = get(sp,'OuterPosition');
    end
    no = 2;
    s = zeros(3*N,no); 
    for i=1:N
        uel = [0 u(2*i-1) u(2*i) 0 u(2*(i+1)-1) u(2*(i+1))];
        z = beam2s_kat([X(i) X(i+1)],[0 0],[E 1 J],uel,[0 fun(0.5*(X(i) + X(i+1)),L,obc)],no);
        s(3*i-2,:) = -z(:,1)';
        s(3*i,:) = z(:,2)';
        s(3*i-1,:) = -z(:,3)';
        for k=1:3
            subplot(3,2,2*k);
            eldia2_test([X(i) X(i+1)],[0 0],s(3*i-(3-k),:)',[4 4],1);
        end
    end                    
    for k=1:3
        sp = subplot(3,2,2*k);
      %  set(sp,'OuterPosition',Pos(k,:));      
        axis(S(k,:)); 
    end
end
eps = 10^-16;
if rcond(A)<eps
    msgbox(str2mat('Z£E PODPARCIE BELKI!'),'B£¥D','error','modal');
end
set(hf,'MenuBar','none');
Pos = get(hf,'Position');
set(hf,'Name','ANALIZA BELEK ZGINANYCH MRS i MES');
set(hf,'NumberTitle','off');
set(hf,'Position',[50 50 1.5*Pos(3) 1.5*Pos(4)]);
et = str2mat('PODPARCIE &LEWE','>SWOBODNE','>WSPORNIK','>PRZESUWNE','>BRAK','PODPARCIE &PRAWE','>SWOBODNE','>WSPORNIK','>PRZESUWNE','>BRAK',...
    '&OBCI¥¯ENIE','>CI¥G£E NA CA£EJ BELCE','>CI¥G£E NA ŒRODKU BELKI [0.25*L,0.75*L]','>SI£A SKUPIONA',...
    '&SIATKA','&MATERIA£','&GEOMETRIA','&WYKRESY','&KONIEC');
wyw = str2mat('',strcat('MRSMES1DIV(2)'),strcat('MRSMES1DIV(3)'),strcat('MRSMES1DIV(4)'),strcat('MRSMES1DIV(5)'),...
    '',strcat('MRSMES1DIV(6)'),strcat('MRSMES1DIV(7)'),strcat('MRSMES1DIV(8)'),strcat('MRSMES1DIV(9)'),...
    '',strcat('MRSMES1DIV(10)'),strcat('MRSMES1DIV(11)'),strcat('MRSMES1DIV(12)'),strcat('MRSMES1DIV(13)'),strcat('MRSMES1DIV(14)'),strcat('MRSMES1DIV(15)'),...
    strcat('MRSMES1DIV(16)'),strcat('close all'));
m = makemenu(hf,et,wyw);

function [z] = fun(x,L,obc)
q = 1;
P = 1;
z = 0;
if obc==1
    z = q;
end
if obc==2
    z = 0;
    if x>L/4 & x<3*L/4
        z = q;
    end
    if x==L/4 | x==3*L/4
        z = q/2;
    end
end

function [sfac]=eldia2_test(ex,ey,es,plotpar,sfac,eci)
%eldia2(ex,ey,es,plotpar,sfac)
%eldia2(ex,ey,es,plotpar,sfac,eci)
%[sfac]=eldia2(ex,ey,es)
%[sfac]=eldia2(ex,ey,es,plotpar)
%
%-------------------------------------------------------------
% PURPOSE: 
%  Draw the section force diagrams of a two dimensional beam element.
%  	
% INPUT:
%	ex = [ x1 x2 ]
%	ey = [ y1 y2 ]	element node coordinates.
%
%	es = [  X1;
%		X2;
%		... ] 	vector containing the section force
%			in Nbr evaluation points along the beam.
%	
%    plotpar=[ linecolor, elementcolor] 
%
%             linecolor=1 -> black    elementcolor=1 -> black
%                       2 -> blue                    2 -> blue
%                       3 -> magenta                 3 -> magenta
%                       4 -> red                     4 -> red
%
%	sfac = [scalar]	scale factor.
%
%	eci = [	x1;
%		x2;
%		... ]	local x-coordinates of the evaluation points (Nbr).
%       If not given, the evaluation points are assumed to be uniformly
%       distributed.
%	
%-------------------------------------------------------------

% LAST MODIFIED: O Dahlblom  2004-10-01
% Copyright (c)  Division of Structural Mechanics and
%                Division of Solid Mechanics.
%                Lund Institute of Technology
%-------------------------------------------------------------

 if ~((nargin==3)|(nargin==4)|(nargin==5)|(nargin==6))
    disp('??? Wrong number of input arguments!')
    return
 end
 a=size(ex); b=size(ey);
 
 if ~((a-b)==[0 0])
    disp('??? Check size of coordinate input arguments!')
    return
 end

 c=size(es);
 Nbr=c(1);
 b=[ex(2)-ex(1);ey(2)-ey(1)];
 Length=sqrt(b'*b);
 n=b/Length;
 if nargin==3
    plotpar=[2 1]
 end
 if((nargin==3)|(nargin==4))
    sfac=(0.2*Length)/max(abs(es));
 end
 if((nargin==3)|(nargin==4)|(nargin==5))
    eci=[0:Length/(Nbr-1):Length]'; 
 end   
 if (size(plotpar)==[1 2])
    if plotpar(1)==1 ; s1='-k';
    elseif plotpar(1)==2 ; s1='-b';
    elseif plotpar(1)==3 ; s1='-m';
    elseif plotpar(1)==4 ; s1='-r';
    else disp('??? Error in variable plotpar(1)!');
       return;
    end
    if plotpar(2)==1 ; s2='-k';
    elseif plotpar(2)==2 ; s2='-b';
    elseif plotpar(2)==3 ; s2='-m';
    elseif plotpar(2)==4 ; s2='-r';
    else disp('??? Error in variable plotpar(2)!');
       return;
    end
 else
    disp('??? Check size of "plotpar" input argument!')
    return
 end
 a=size(eci); 
 if ~(c(1)==a(1)|(c(2)==1))
    disp('??? Check size of "es" or "eci" input arguments!')
    return 
 end 
       
 es=es*sfac;
% From local x-coordinates to global coordinates of the beam element.
 A=zeros(Nbr,2);
 A(1,1)=ex(1);
 A(1,2)=ey(1);
 for i=2:Nbr
	A(i,1)=A(1,1)+eci(i)*n(1);
	A(i,2)=A(1,2)+eci(i)*n(2);
 end

 B=A;
 for i=1:Nbr
	A(i,1)=A(i,1)+es(i)*n(2);
	A(i,2)=A(i,2)-es(i)*n(1);
 end

% ************* plot commands *******************
 hold on
 %axis equal
 % Plots diagram.
 plot(A(:,1),A(:,2),s1);
 plot([ex(1) A(1,1)],[ey(1) A(1,2)],s1);
 plot([ex(2) A(Nbr,1)],[ey(2) A(Nbr,2)],s1);
% Plots stripes in diagram.
 for i=1:Nbr
	plot([B(i,1) A(i,1)],[B(i,2) A(i,2)],s1);
 end
% Plots element.	
 a=plot(ex',ey',s2);
 set(a,'LineWidth',[2]);
 hold off;

 function [es,edi,eci]=beam2s_kat(ex,ey,ep,ed,eq,n)
% es=beam2s(ex,ey,ep,ed)
% es=beam2s(ex,ey,ep,ed,eq)
% [es,edi,eci]=beam2s(ex,ey,ep,ed,eq,n)
%---------------------------------------------------------------------
%    PURPOSE
%      Compute section forces in two dimensional beam element (beam2e). 
% 
%    INPUT:  ex = [x1 x2]
%            ey = [y1 y2]     element node coordinates
%
%            ep = [E A I]     element properties,
%                              E:  Young's modulus
%                              A:  cross section area
%                              I:  moment of inertia
%
%            ed = [u1 ... u6] element displacements
%
%            eq = [qx qy]     distributed loads, local directions 
%
%            n : number of evaluation points ( default=2 )
%          
%    OUTPUT: es = [ N1 V1 M1 ;  section forces, local directions, in 
%                   N2 V2 M2 ;  n points along the beam, dim(es)= n x 3
%                   .........]  
%           
%            edi = [ u1 v1 ;    element displacements, local directions,
%                    u2 v2 ;    in n points along the beam, dim(es)= n x 2
%                   .......]    
%
%            eci = [ x1  ;      local x-coordinates of the evaluation 
%                    x2 ;       points, (x1=0 and xn=L)
%                    ...]
%-------------------------------------------------------------------------
 
% LAST MODIFIED: K Persson    1995-08-23
% Copyright (c)  Division of Structural Mechanics and
%                Department of Solid Mechanics.
%                Lund Institute of Technology
%-------------------------------------------------------------
  EA=ep(1)*ep(2); EI=ep(1)*ep(3);
  b=[ ex(2)-ex(1); ey(2)-ey(1) ];
  L=sqrt(b'*b);

  if length(ed(:,1)) > 1 
     disp('Only one row is allowed in the ed matrix !!!')
     return
   end
  
  qx=0; qy=0;  if nargin>4;  qx=eq(1); qy=eq(2); end 
    
  ne=2;        if nargin>5;  ne=n; end;
     
  C=[0   0   0    1   0   0;
     0   0   0    0   0   1;
     0   0   0    0   1   0;
     L   0   0    1   0   0;
     0   L^3  L^2 0   L   1;
     0 3*L^2 2*L  0   1   0];

  n=b/L;

  G=[n(1) n(2)  0    0    0   0;
    -n(2) n(1)  0    0    0   0;
      0    0    1    0    0   0;
      0    0    0   n(1) n(2) 0;
      0    0    0  -n(2) n(1) 0;
      0    0    0    0    0   1];
 
  M=inv(C)*(G*ed'-[0 0 0 -qx*L^2/(2*EA) qy*L^4/(24*EI) qy*L^3/(6*EI)]' );
  
  A=[M(1) M(4)]';  B=[M(2) M(3) M(5) M(6)]';
  
  x=[0:L/(ne-1):L]';   zero=zeros(size(x));    one=ones(size(x));
  
  u=[x one]*A-(x.^2)*qx/(2*EA);
  du=[one zero]*A-x*qx/EA;
  v=[x.^3 x.^2 x one]*B+(x.^4)*qy/(24*EI); 
  dv=[3*x.^2 2*x one zero]*B+(x.^3)*qy/(6*EI);
  d2v=[6*x 2*one zero zero]*B+(x.^2)*qy/(2*EI);
  d3v=[6*one zero zero zero]*B+x*qy/EI;
  
  N=EI*dv; M=EI*d2v; V=-EI*d3v; 
  es=[N V M];
  edi=[u v];
  eci=x;
%--------------------------end--------------------------------
 

