clear all
close all

%Tower Dimensions and Mechanical Properties
h=100; %height
dz=0.1;
z=0:dz:h; %height vector
E = 200 * 10^9; %Young's Modulus of Steel [Pa = N/m^2]
I = 0.5*( (pi/64) * ((4.5^4-4.4962^4) + (3.8^4-3.785^4)) ); %Average Moment of Inertia of Tapered Tower (diameter range: [4.5 3.8][m^4]
max_stress = 400 * 10^6; %Pa
c = 4.5/2; %maximum height of cross section from centroidal axis
max_moment = (max_stress * I) / c;

%Constants
cT = 0.6; %drag coefficient for tower
rho = 1.225; %density of air in kg/m^3
cD = 0.04; %drag coefficient for airfoil shape (SMALL) units: Fd / A
A = 350; %reference area of airfoil in m^2
f = 0;

%Wind speed vector
u = linspace(0,50,1001);

%Main Loop
for i=1:length(u)
    %Thrust Force due to blades
    fT(i) = 0.5*rho*cD*A.*(u(i))^2;
    Ft = [zeros(1,1000) fT(i)];
    
    %Distributed Wind Force
    fW(i) = 0.5*rho*cT*4.15*(u(i))^2;
    Fw = repmat(fW(i),1,1001);
    
    %Total Loading on Tower
    F_total = Ft + Fw;
    
    %Reactions
    Ay(i) = fT(i) + (fW(i) * h);
    MA(i) = fT(i)*h + 0.5*fW(i)*h^2;
    
    %Initial Values
    y4 = F_total/(E*I); %Load / EI
    y3 = Ay(i)/(E*I); %Inital shear at bottom / EI
    y2 = MA(i)/(E*I); %Initial bending moment at bottom / EI
    y1 = 0; %No rotation at fixed end of cantilever beam
    y = 0; %No deflection at fixed end of cantilever beam
    
    %Discretization Loops
    for k=2:length(F_total)
        y3(k) = y3(k-1) - y4(k-1)*dz;
        y2(k) = y2(k-1) - y3(k-1)*dz;
        y1(k) = y1(k-1) + y2(k-1)*dz;
        y(k) =   y(k-1) + y1(k-1)*dz;
    end
    
    %Plot
    if y2*E*I < max_moment
        %if mod(i+9,10)==0
            f = f + 1;
            plot(y,z,'-b')
            grid on
            axis([0 2.5 0 110])
            xlabel('Deflection (meters)')
            ylabel('Height (meters)')
            title(['wind speed = ', num2str(u(i)), ' m/s'])
            %saveas(gcf, ['deflection2_',sprintf('%03d',f),'.png'])
            pause(0.001)
        %end
    else
        plot(y2*E*I,z,'-r')
        grid on
        xlabel('Bending moment (N-m)')
        ylabel('Height (meterts)')
        title(['Bending Moment Diagram at Fracture; Wind speed = ', num2str(u(i)), ' m/s'])
        break
    end
end