% FEM, galerkinov metod
% jna: Ay'' + By' + Cy + D = 0, y(0) = 0, y'(1) = 1
% sa simbolickim promenljivama
% trapezno krilo

clear all, 
clc

% ulazni podaci
ag = 4.774;
Cmac = -0.04;
b2 = 5;
e = 0.25;
V = 50;
rho = 1.225;
G = 2.769e10;
t = 0.001;
alpha = 4*pi/180;

% jedan konacni element
Ne = 1;
a = 0;
b = b2;
x = linspace(a,b,10*Ne+1);

% pomocna promenljiva za resavanje integrala
syms xs

% karakteristike krila
l = 0.8-0.06*xs;
J = 0.0257*t*(0.8-0.06*xs)^3;
dJdy = -0.00463*t*(0.8-0.06*xs)^2;

% koeficijenti dif. j-ne
A = G*J;
B = G*dJdy;
C = rho*V^2/2*ag*l^2*e;
D = rho*V^2/2*l^2*(ag*alpha*e + Cmac);

M = 2;
% matrica koeficijenata
K = zeros(M, M);
% matrica slobodnih clanova
f = zeros(M, 1);

% proracun matrica koeficijenata i slobodnih clanova
for i = 1:M
    for j = 1:M
        K(i,j) = int((-A*j*(j+1)/b2^2*(1-xs/b2)^(j-1) + B*(j+1)/b2*(1-xs/b2)^j +...
            C*(1-(1-xs/b2)^(j+1))) * (1-(1-xs/b2)^(i+1)), a, b);
    end
    f(i,1) = int(-D * (1-(1-xs/b2)^(i+1)), a, b);
end

% koeficijenti
A = K\f;
% resenje
y = 0;
for i = 1:M
    y = y + A(i)*(1-(1-x/b2).^(i+1));
end

% prikaz rezultata
figure,
set(gca,'FontSize',20)
plot(x,y*180/pi,'-o','LineWidth',3);
xlabel('y');
ylabel('\theta [\circ]');
title(['M = ', int2str(M), ', V = ', int2str(V), 'm/s'])
axis([0 b2 -1000 1000]);

y(10*Ne+1)*180/pi