%
% Solve plates by double trigonometric series
% Решава плочи чрез двойни тригонометрични редове
%
clear all
format compact
syms x y
% input data - Входни данни за плочата
% [m]  [m]  [m]      [Pa]    [-]    
a=6; b=5;t=0.1; E=10.e+9; nu=0.10; 
% Load function - Функция на натоварването и граници
% [N/m^2]                          Concentrated force - Съсредоточена сила
p=10.0e+3; xl=2.2; xu=4.2; yl=1.1; yu=4.0; c=a/2; d=b/2; 
% Degree of descretization - Степен на дискретизация на редовете и по осите
% m     n      x     y
nm=11; nn=11; nx=24; ny=24; 
%==================================================================
mv=[1:1:nm]; nv=[1:1:nn];
D=E*t^3/(12*(1-nu^2));
for m=mv, for n=nv,
        q(m,n)=4/(a*b)*single(int(sin(m*pi*x/a)*int(p*sin(n*pi*y/b),y,yl,yu),x,xl,xu));
        %q(m,n)=4*p/(a*b)*sin(m*pi*c/a)*sin(n*pi*d/b); % Съсредоточена сила
        A(m,n)=q(m,n)/(D*pi^4)/(m^2/a^2+n^2/b^2)^2;
    end, end

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

for i=1:length(xi),
    for j=1:length(eta),
        w(i,j)=0;Mx(i,j)=0;My(i,j)=0;Mxy(i,j)=0;
        for m=mv,
            for n=nv,
                w(i,j)=w(i,j)+A(m,n)*sin(m*pi*xi(i))*sin(n*pi*eta(j));
                Mx(i,j)=Mx(i,j)+A(m,n)*D*pi^2*(m^2/a^2+nu*n^2/b^2)*sin(m*pi*xi(i))*sin(n*pi*eta(j));
                My(i,j)=My(i,j)+A(m,n)*D*pi^2*(nu*m^2/a^2+n^2/b^2)*sin(m*pi*xi(i))*sin(n*pi*eta(j));
                Mxy(i,j)=Mxy(i,j)-A(m,n)*D*pi^2*(1-nu)*m*n*cos(m*pi*xi(i))*cos(n*pi*eta(j));
            end
        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)')

