syms a1 a2 a3 a4 a5 x y h l q

psi = a1 * x^2 + a2 * x^2 * y + a3 * y^3 + a4 * x^2 * y^3 + a5 * y^5;


nab4 = diff(psi,'x',4) + diff(psi,'y',4) + 2*diff(diff(psi,'x',2),'y',2);
a5 = solve(nab4,'a5');
psi = subs(psi,'a5',a5);

sx = diff(psi,'y',2);
sy = diff(psi,'x',2);
txy = -diff(diff(psi,'x'),'y');


eq1 = subs(sy,'y',h/2) + q;
eq2 = subs(sy,'y',-h/2);
eq3 = subs(txy,'y',h/2);


a1 = solve(eq1+eq2, 'a1');
a2 = solve(eq3,'a2');
a4 = solve(subs(subs(eq2,'a1',a1),'a2',a2),'a4');
a2 = subs(a2,'a4',a4);


sx = subs(subs(subs(sx,'a1',a1),'a2',a2),'a4',a4);
sy = subs(subs(subs(sy,'a1',a1),'a2',a2),'a4',a4);
txy = subs(subs(subs(txy,'a1',a1),'a2',a2),'a4',a4);



int(subs(sx,'x',l/2),'y',-h/2,h/2)

a3 = solve(int(subs(sx,'x',l/2)*y,'y',-h/2,h/2),'a3');


sx = subs(sx,'a3',a3);
sy = subs(sy,'a3',a3);
txy = subs(txy,'a3',a3);





l = 3;
h = 1;
Nx = 100;
Ny = 20;
q = 1e6;

xa = linspace(-l/2, l/2, Nx);
ya = linspace(-h/2, h/2, Ny);

sx_ = zeros(Nx,Ny);
sy_ = zeros(Nx,Ny);
txy_ = zeros(Nx,Ny);

for j=1:Nx
    for k=1:Ny
        x = xa(j);
        y = ya(k);
        sx_(j,k) = eval(sx);
        sy_(j,k) = eval(sy);
        txy_(j,k) = eval(txy);
    end
end


figure(1)
pcolor(xa'*ones(1,Ny),ones(Nx,1)*ya,sx_)
axis equal
colorbar
title('\sigma_x')
xlabel('x')
ylabel('y')

figure(2)
pcolor(xa'*ones(1,Ny),ones(Nx,1)*ya,sy_)
axis equal
colorbar
title('\sigma_y')
xlabel('x')
ylabel('y')


figure(3)
pcolor(xa'*ones(1,Ny),ones(Nx,1)*ya,txy_)
axis equal
colorbar
title('\tau_{xy}')
xlabel('x')
ylabel('y')




figure(4)
plot(sx_(Nx/2,:),ya, -1/8*q*l^2 / (1/12*h^3)*ya, ya,'r')
xlabel('\sigma_x(l/2)')
ylabel('y')

figure(5)
plot(sy_(Nx/2,:),ya)
xlabel('\sigma_y(l/2)')
ylabel('y')

figure(6)
plot(txy_(Nx/4,:),ya, -q*l/4 * 1/2*(h^2/4 - ya.^2) / (1/12*h^3),ya,'r')
xlabel('\tau_{xy}(l/4)')
ylabel('y')



E = 210e9;
ni = 0.3;

ez_ = -ni/E * (sx_ + sy_);


eScale = 1e5;
figure(7)
surf(xa'*ones(1,Ny),ones(Nx,1)*ya,ez_*eScale)
title('\epsilon_z')
ylabel('y')
xlabel('x')
axis equal



s = 0.1;
eScale = 1e5;
figure(8)
surf(xa'*ones(1,Ny),ones(Nx,1)*ya, (ez_*eScale + 1)*s / 2)
hold on
surf(xa'*ones(1,Ny),ones(Nx,1)*ya, -(ez_*eScale + 1)*s / 2)
hold off
ylabel('y')
xlabel('x')
axis equal




sI = zeros(Nx,Ny);
sII = zeros(Nx,Ny);
nI = zeros(Nx,Ny,2);
nII = zeros(Nx,Ny,2);

for j=1:Nx
    for k=1:Ny
        T = [sx_(j,k) txy_(j,k)
             txy_(j,k) sy_(j,k)];
        [N,L] =  eig(T);
        sI(j,k) = L(2,2);
        sII(j,k) = L(1,1);
        nI(j,k,:) = N(:,2);
        nII(j,k,:) = N(:,1);
    end
end

        
figure(10)
axis equal
pcolor(xa'*ones(1,Ny),ones(Nx,1)*ya,sI)
colorbar
title('\sigma_I')
xlabel('x')
ylabel('y')

figure(11)
axis equal
pcolor(xa'*ones(1,Ny),ones(Nx,1)*ya,sII)
colorbar
title('\sigma_{II}')
xlabel('x')
ylabel('y')
        

X = xa'*ones(1,Ny);
Y = ones(Nx,1)*ya;

sIp = zeros(Nx,Ny);
sIm = zeros(Nx,Ny);
sIIp = zeros(Nx,Ny);
sIIm = zeros(Nx,Ny);

for j=1:Nx
    for k=1:Ny
        if sI(j,k) >=0
            sIp(j,k) = sI(j,k);
        else
            sIm(j,k) = sI(j,k);
        end            
        if sII(j,k) >=0
            sIIp(j,k) = sII(j,k);
        else
            sIIm(j,k) = sII(j,k);
        end            
    end
end

sscale = .5;
figure(12)
axis equal
mesh(X,Y,zeros(Nx,Ny))
hold on
quiver(X,Y, nI(:,:,1).*sIp,nI(:,:,2).*sIp,sscale,'color','b')
quiver(X,Y, -nI(:,:,1).*sIp,-nI(:,:,2).*sIp,sscale,'color','b')
quiver(X,Y, nI(:,:,1).*sIm,nI(:,:,2).*sIm,sscale,'color','r')
quiver(X,Y, -nI(:,:,1).*sIm,-nI(:,:,2).*sIm,sscale,'color','r')
hold off
view(2)

figure(13)
axis equal
mesh(X,Y,zeros(Nx,Ny))
hold on
quiver(X,Y, nII(:,:,1).*sIIp,nII(:,:,2).*sIIp,sscale,'color','b')
quiver(X,Y, -nII(:,:,1).*sIIp,-nII(:,:,2).*sIIp,sscale,'color','b')
quiver(X,Y, nII(:,:,1).*sIIm,nII(:,:,2).*sIIm,sscale,'color','r')
quiver(X,Y, -nII(:,:,1).*sIIm,-nII(:,:,2).*sIIm,sscale,'color','r')
hold off
view(2)



fscale = 0.5;
figure(14)
axis equal
mesh(X,Y,zeros(Nx,Ny))
hold on
quiver(xa, ones(1,Nx)*h/2, txy_(:,Ny)', sy_(:,Ny)', fscale,'color','b')
quiver(xa, ones(1,Nx)*-h/2, -txy_(:,1)', -sy_(:,1)', fscale,'color','b')
quiver(ones(1,Ny)*l/2, ya, sx_(Nx,:), txy_(Nx,:), fscale,'color','b')
quiver(ones(1,Ny)*-l/2, ya, -sx_(1,:), -txy_(1,:), fscale,'color','b')
hold off
view(2)
