



%% not interesting- dkip this part




clear all
close all
clc
 
%constant
lib2Kg=0.329;
lib2N=lib2Kg*9.807;
feet2m=0.3048;

ei=2.5*10^6*lib2N*feet2m^2;
gj=0.4*10^6*lib2N*feet2m^2;
b=(40+80)*2*feet2m;
clalpha=2*pi;
c=8*feet2m;

ac=0.25*c; %[m]
ea=0.26*c; %[m]
n=100; %Number of elements
l=b/2/n; %[m]
le=b/2/n;
rho=1.225; %[kg/m^3]
%ang=20*pi/180;
econ=ea-ac;
m=6*lib2Kg/feet2m;%[kg];
EI=ei;
s=c*l;
Length=b/2; 
mw=1.5;%[kg]
d=0.01*c;%[m] !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!to gain analitical solution- d=0
 
 
K=3*EI/Length;
M=m*Length;

p0=K*2;






%% start watching from here!!!

w_gag=5;
w_struct_1=sqrt(K/M); T_struct_1=2*pi/w_struct_1;
w_struct_3=sqrt(K/M)*5; T_struct_3=2*pi/w_struct_3;
w_gust=w_struct_1/2; T_gust=2*pi/w_gust;
v=12;


dt=min(T_struct_3/10,T_gust/10);
T=max(T_gust*5,T_struct_1*5);


size_sig=2^nextpow2(T/dt);

for i=1:size_sig
    t(i)=dt*i-dt;
    if t(i)<T_gust
        gust(i)=0.5*w_gag*(1-cos(w_gust*t(i)));
        w_dot(i)=0.5*w_gag*w_gust*sin(w_gust*t(i));
    else
        gust(i)=0;
        w_dot(i)=0;
    end
end

s_kussner=v*t/c*2;
kussner=(s_kussner.^2+s_kussner)./(s_kussner.^2+2.82*s_kussner+0.8);
cnv=conv(kussner,w_dot)*dt;
ghust_load=cnv(1:length(t));
F=ghust_load;



signal=F;

fs=1/dt; % fs is the sampling frequency
nfft=length(signal);
y=fft(signal,nfft)/nfft;
y=conj(y);
%y=signal*dftmtx(nfft)/nfft;




y(1)=y(1)/2;
yy=2*sqrt(y.*conj(y)); % the 2 factor here is in order to get the amplitude peak-to-peak
Amp=yy(1:nfft/2); % here I get rid of all frequencies above half the sampling frequency
Amp=2*abs(y(1:nfft/2));
phase=atan(real(y(1:nfft/2))./imag(y(1:nfft/2)));
f=(0:(nfft/2)-1)/nfft*fs;
plot(f,Amp);


size=100;


f=f(1:size);
mekad=y(1:size)*2;
figure(2)
plot(f,real(mekad(1:size)),f,imag(mekad(1:size)),f,abs(mekad(1:size)))
legend ('imaginarry part' ,'real part','amplitude','location','best')
grid
w=2*pi*f;


figure(3) %plot signal


for i=1:length(phase)
    if imag(y(i))==0
        if real(y(i))<pi
        phase(i)=pi/2;
        else
        phase(i)=-pi/2*0; 
        end
    end
end



for i=1:length(t)
    
        FF(i)=real(mekad)*cos(w'*t(i))+imag(mekad)*sin(w'*t(i));
        %FF(i)=Amp*cos(w'*t(i)+phase'/2*2);
        
end

plot(t,gust,t,F,t,FF)
legend('orugunal gust gust','gust after kussner', 'gust after kusner filtered and returned' ,'location','best')
grid
%s=(0:0.01:10);
%kusner=(s.^2+s)./(s.^2+2.82*s+0.8);
%figure (4)
%plot(s,kusner)


