% Priklad na nalezeni optimalni pozice podpor na spojitem nosniku o dvou
% polich a previslem konci
function [X,Y,Z] = spojak()
clear, clc, close all;

%%% Pro symbolicke reseni
% syms x1 x2 L EI ;
% syms f fi2 fi3 ;

% meze
lb = [0.2,0.75] ;
ub = [0.55,1.0];

x1 = 0.39520915418298 ;
x2 = 0.817732328824383 ;

F = moments([x1,x2]) ;

%% Optimum
% L = 12, f = 0.3*0.5*25 ;
% x =   0.3952    0.8175
% fval =    8.9923
[x,fval] = fmincon(@fmoments,[0.15,0.55],[],[],[],[],lb,ub)

%% Vykresleni plochy a optim
%
[X,Y] = meshgrid(lb(1):0.005:ub(1), lb(2):0.005:ub(2)); 

[m n]=size(X) ;

X = reshape(X,1,[]) ;
Y = reshape(Y,1,[]) ;
Z = zeros(1,length(X)) ;

for i=1:length(X)
    M = moments([X(i),Y(i)]) ;
    Z(i) = max(abs(M));
    Mx(i) = max(M);
    Mn(i) = min(M);
end

XX = reshape(X,m,n) ;
YY = reshape(Y,m,n) ;
ZZ = reshape(Z,m,n) ;

% vsechna reseni metodou Nelder-Mead z CSV souboru
% importovane z Excelu, ulozene do matice Matlabu
load('nm.mat') ;

subplot(2,2,1)
% surf(XX,YY,ZZ)
[~,h] = contour(XX,YY,reshape(abs(Mx+Mn),m,n),9) ; hold on ;
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
plot([0.5481], [0.8321],'-ko', 'MarkerEdgeColor','k',...
                'MarkerFaceColor','r',...
                'MarkerSize',6)
title('abs(M_{max}+M_{min}), optimum [0.5481,0.8321]')
subplot(2,2,2)
[~,h] = contour(XX,YY,ZZ,20) ; hold on ;
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
plot([0.3952], [0.8175],'-ko', 'MarkerEdgeColor','k',...
                'MarkerFaceColor','r',...
                'MarkerSize',6)
title('max(abs(M_{max};M_{min})), optimum [0.3952,0.8175]')
%%%%%%%%%%%%%%%
plot(nm(:,1), nm(:,2),'ko', 'MarkerEdgeColor','k',...
                'MarkerFaceColor','r',...
                'MarkerSize',6)

subplot(2,2,3)
[~,h] = contour(XX,YY,reshape(Mx,m,n),20) ;
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
title('M_{max}')
subplot(2,2,4)
[~,h] = contour(XX,YY,reshape(Mn,m,n),20) ;
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
title('M_{min}')
hold on ;

end

%% Definice cilove funkce
%
function M=fmoments(X)

M = max(abs(moments(X))) ; 

end

%% Vypocet prubehu momentu deformacni metodou
%
function M=moments(X)

x1 = X(1) ;
x2 = X(2) ;

% delka
L = 12 ;
% vlastni tiha
f = 0.3*0.5*25 ;
% tuhost
EI = 20 ;

L1 = x1*L ;
L2 = x2*L - L1 ;
L3 = L - L1 - L2 ;

%%% Sestaveni rovnic v symbolickem zapisu
% R1 =
%  
% (f*(L*x1 - L*x2)^2)/12 - (4*EI*fi2)/(L*x1 - L*x2) - (2*EI*fi3)/(L*x1 - L*x2) - (L^2*f*x1^2)/8 + (3*EI*fi2)/(L*x1)
%  
%  
% R2 =
%  
% (f*(L - L*x2)^2)/2 - (f*(L*x1 - L*x2)^2)/12 - (2*EI*fi2)/(L*x1 - L*x2) - (4*EI*fi3)/(L*x1 - L*x2)

% Prava strana
F = [ - (f*(L*x1 - L*x2)^2)/12 + (L^2*f*x1^2)/8; -(f*(L - L*x2)^2)/2 + (f*(L*x1 - L*x2)^2)/12 ] ;
% matice tuhosti
K = [ - (4*EI)/(L*x1 - L*x2) + (3*EI)/(L*x1) , - (2*EI)/(L*x1 - L*x2) ;
    - (2*EI)/(L*x1 - L*x2), - (4*EI)/(L*x1 - L*x2)] ;

% deformace - natoceni
fi = K\F ;

fi2 = fi(1) ;
fi3 = fi(2) ;

% koncove momenty
M21 = -f*L1*L1/8 + 3*EI*fi2/L1 ;
M23 = f*L2*L2/12 + 4*EI*fi2/L2 + 2*EI*fi3/L2 ;
M32 = -f*L2*L2/12 + 2*EI*fi2/L2 + 4*EI*fi3/L2 ;
M34 = f*L3*L3/2 ;

%%% kontrola rovnovahy
% R1 = M21 + M23 
% R2 = M32 + M34

% extremy momentu
Z1 = (f*L1*L1/2 + M21)/L1 ;
X1 = Z1/f ;
Mmax1 = Z1*X1/2 ;

Z2 = (f*L2*L2/2 + M23 + M32)/L2 ;

if Z2>0
    X2 = Z2/f ;
    Mmax2 = -M23 + Z2*X2/2 ;
else
    Mmax2 = 0 ;
end

M = [Mmax1 Mmax2 M21 M32 ] ; 

% M = max(abs([Mmax1 Mmax2 M21 M32 ])) ; 
% M = abs(max([Mmax1 Mmax2])+ min([M21 M32 ])) ; 

end