%
% Solve plates by Moris Levy's single trigonometric series method
% Решава плочи чрез метода на Морис Леви на единичните тригонометрични редове
%
clear all
format compact
syms x y
% input data - Входни данни за плочата
% [m]  [m]  [m]      [Pa]    [-]    
a=6; b=2.5; t=0.1; E=10.e+9; nu=0.10; 
% Load function - Функция на натоварването и граници
% [N/m^2]                          Concentrated force - Съсредоточена сила
p=50.0e+3; xl=0; xu=a;  
% Degree of descretization - Степен на дискретизация на редовете и по осите
% m      x     y
nm=11; nx=24; ny=24; 
%==================================================================
mv=[1:2:nm]; 
Ds=E*t^3/(12*(1-nu^2));

for m=mv, 
        al(m)=pi*m/a;
        q(m)=2/a/Ds*p*single(int(sin(al(m)*x),x,0,a));
        Yb(m)=4*p*a^4/Ds/(m*pi)^5;
        th(m)=al(m)*b;
        A(m)=-Yb(m)*(sinh(th(m))+th(m)*cosh(th(m)))/(cosh(th(m))*sinh(th(m))+th(m));
        D(m)=Yb(m)*al(m)*sinh(th(m))/(cosh(th(m))*sinh(th(m))+th(m));
end

xi=[0:1./nx:1]; eta=[-1:1./2/ny:1];

for i=1:length(xi),
    xx=xi(i)*a;
    for j=1:length(eta),
        yy=eta(j)*b;
        w(i,j)=0;Mx(i,j)=0;My(i,j)=0;Mxy(i,j)=0;
        for m=mv,
            Ym=A(m)*cosh(al(m)*yy)+D(m)*yy*sinh(al(m)*yy)+Yb(m);
            Yim=al(m)*A(m)*sinh(al(m)*yy)+D(m)*(sinh(al(m)*yy)+al(m)*yy*cosh(al(m)*yy));
            Yiim=al(m)^2*A(m)*cosh(al(m)*yy)+D(m)*al(m)*(2*cosh(al(m)*yy)+al(m)*yy*sinh(al(m)*yy));
            w(i,j)=w(i,j)+Ym*sin(al(m)*xx);
            Mx(i,j)=Mx(i,j)-Ds*(nu*Yiim-al(m)^2*Ym)*sin(al(m)*xx);
            My(i,j)=My(i,j)-Ds*(Yiim-nu*al(m)^2*Ym)*sin(al(m)*xx);
            Mxy(i,j)=Mxy(i,j)-Ds*(1-nu)*al(m)*Yim*cos(al(m)*xx);
        end
    end
end

figure
surfc(eta*b,xi*a,w)
xlabel('y (m)'),ylabel('x (m)'),zlabel('w (m)')
figure
surfc(eta*b,xi*a,Mx)
xlabel('y (m)'),ylabel('x (m)'),zlabel('M_x (N)')
figure
surfc(eta*b,xi*a,My)
xlabel('y (m)'),ylabel('x (m)'),zlabel('M_y (N)')
figure
surfc(eta*b,xi*a,Mxy)
xlabel('y (m)'),ylabel('x (m)'),zlabel('M_{xy} (N)')

