% The script file pwoprt03 approximates damped free vibrations by the
% expression for linear not damped oscillations and damped vibrations.
% The iterations converge.

g=9.81; T1=0.1; f=1/T1; omega1=2*pi*f; amax=5; dt=T1/100; eta=0; 
a=amax; b=0; f=f; et=1/20/T1; % values used in final analysis

t=0:dt:4*T1-dt; t=t'; N=4*T1/dt;  % one point less 
E=exp(-eta*t); AE=amax*E;
H= exp(-i*omega1*t); RH=amax*real(H); 
Z=AE.*H; RZ=real(Z);
plot(t,RZ,'k',t,AE,'r',t,-AE,'r',t,RH,'b')
title('Displacements in not damped vibrations, envelops ')
pause
close
ae=0; be=0; fe=1.05*f; etae=0; r=0;

VALUES=[r, ae, be, fe]

for r=1:8

Y=RZ(1:N); T=0:dt:(N-1)*dt; T=T'; R=T;
E=ones(N,1); C=cos(2*pi*fe*T); S=sin(2*pi*fe*T); 
EC=E.*C; ES=E.*S; REC=R.*EC; RES=R.*ES; RREC=R.*REC; RRES=R.*RES;
E02=Y'*Y;
EJ=Y-ae*EC-be*ES; EJ2=EJ'*EJ; WE0=sqrt(EJ2/E02);

EA=2*EJ'*(-EC);
EAA=2*EC'*EC; EAB=2*ES'*EC; 
EAF=2*[-ae*(-2*pi*RES)-be*2*pi*REC]'*(-EC)+2*EJ'*(2*pi*RES);
EAD=2*[-ae*(-REC)-be*(-RES)]'*(-EC)+2*EJ'*(REC);

EB=2*EJ'*(-ES);
EBA=EAB; EBB=2*(ES'*ES);
EJF=[-ae*(-2*pi*RES)-be*2*pi*REC];
EBF=2*EJF'*(-ES)+2*EJ'*(-2*pi*REC);
EJD=[-ae*(-REC)-be*(-RES)];
EBD=2*EJD'*(-ES)+2*EJ'*(RES);

EF=2*EJ'*EJF;
EFA=EAF; EFB=EBF;
EFF=2*EJF'*EJF+2*EJ'*[-ae*(-4*pi*pi*RREC)-be*(-4*pi*pi*RRES)]; %?
EFD=2*EJD'*EJF+2*EJ'*(-ae*RRES+be*RREC)*2*pi; %?
   
ED=2*EJ'*EJD;
EDA=EAD; EDB=EBD; EDF=EFD;
EDD=2*EJD'*EJD+2*EJ'*[ae*RREC+be*RRES];

% W2=[EAA,EAB;EBA,EBB]; IW2=inv(W2); P2=-[EA;EB];
% X2=IW2*P2 

W3=[EAA,EAB,EAF;EBA,EBB,EBF;EFA,EFB,EFF]; IW3=inv(W3);
P3=-[EA;EB;EF];
X3=IW3*P3; X3=X3';
%RESULT=[r,X3]
ae=ae+X3(1); be=be+X3(2);fe=fe+X3(3);
EJ=Y-ae*EC-be*ES; EJ2=EJ'*EJ; WE0=sqrt(EJ2/E02);
ITERATIONS=[r,ae,be,fe,WE0]

end
%ae1=ae; be1=be; fe1=fe;

g=9.81; T1=0.1; f=1/T1; omega1=2*pi*f; amax=5; dt=T1/100; 
eta=1/20/T1;

t=0:dt:4*T1-dt; t=t'; N=4*T1/dt;  % one point less 
E=exp(-eta*t); AE=amax*E;
H= exp(-i*omega1*t); RH=amax*real(H); 
Z=AE.*H; RZ=real(Z);
plot(t,RZ,'k',t,AE,'r',t,-AE,'r',t,RH,'b')
title('Displacements in damped vibrations envelops and periodic component')
pause
close
yp=amax;
yk=amax*exp(-eta*4*T1)*cos(2*pi*f*4*T1);

ae=0; be=0; fe=f; etae=0.5; r=0;

VALUESB=[r, ae, be, fe,etae]

for r=1:8

Y=RZ(1:N); T=0:dt:(N-1)*dt; T=T'; R=T;
E=ones(N,1); C=cos(2*pi*fe*T); S=sin(2*pi*fe*T); 
EC=E.*C; ES=E.*S; REC=R.*EC; RES=R.*ES; RREC=R.*REC; RRES=R.*RES;
E02=Y'*Y;
EJ=Y-ae*EC-be*ES; EJ2=EJ'*EJ; WE0=sqrt(EJ2/E02);

EA=2*EJ'*(-EC);
EAA=2*EC'*EC; EAB=2*ES'*EC; 
EAF=2*[-ae*(-2*pi*RES)-be*2*pi*REC]'*(-EC)+2*EJ'*(2*pi*RES);
EAD=2*[-ae*(-REC)-be*(-RES)]'*(-EC)+2*EJ'*(REC);

EB=2*EJ'*(-ES);
EBA=EAB; EBB=2*(ES'*ES);
EJF=[-ae*(-2*pi*RES)-be*2*pi*REC];
EBF=2*EJF'*(-ES)+2*EJ'*(-2*pi*REC);
EJD=[-ae*(-REC)-be*(-RES)];
EBD=2*EJD'*(-ES)+2*EJ'*(RES);

EF=2*EJ'*EJF;
EFA=EAF; EFB=EBF;
EFF=2*EJF'*EJF+2*EJ'*[-ae*(-4*pi*pi*RREC)-be*(-4*pi*pi*RRES)]; %?
EFD=2*EJD'*EJF+2*EJ'*(-ae*RRES+be*RREC)*2*pi; %?
   
ED=2*EJ'*EJD;
EDA=EAD; EDB=EBD; EDF=EFD;
EDD=2*EJD'*EJD+2*EJ'*[ae*RREC+be*RRES];

% W2=[EAA,EAB;EBA,EBB]; IW2=inv(W2); P2=-[EA;EB];
% X2=IW2*P2 

W3=[EAA,EAB,EAF;EBA,EBB,EBF;EFA,EFB,EFF]; IW3=inv(W3);
P3=-[EA;EB;EF];
X3=IW3*P3; X3=X3';
%RESULT=[r,X3]
ae=ae+X3(1); be=be+X3(2);fe=fe+X3(3);
EJ=Y-ae*EC-be*ES; EJS=EJ'*EJ; WE0=sqrt(EJS/E02);
ITERATIONS=[r,ae,be,fe,WE0]
%r=r+1;

end
ae2=ae; be2=be; fe2=fe;

YC2=ae2*cos(2*pi*fe2*T)+be2*sin(2*pi*fe2*T);
plot(T,Y,'k',T,YC2,'r')
title('The vibrations and the approximation by three values a, b, f')
pause
close

dy=yp-yk; eta3=dy/4/T1/yp; AE=(ae+dy/2)*(1-eta3*T); 
YC3=AE.*cos(2*pi*fe2*T);
plot(T,Y,'k',T,YC3,'r')
title('The vibrations and the approximation by four values a,b,f,eta')
pause
close
EJ3=Y-YC3; EJS=EJ3'*EJ3; WE0=sqrt(EJS/E02);
VALUES=[a,b,f,et]
APPROXIMATION=[ae+dy/2,be,fe2,eta3,WE0]


















 
