% /*M-FILE Script SSTS__tether_MMET_III_replot MMM SSTSLAB */
% /*==================================================================================================
%  Simple Space Tether Simulation Laboratory Toolbox for Matlab 7.x
%
%  Copyright 2007 The SxLAB Family - Yi Chen - leo.chen.yi@gmail.com
% ====================================================================================================
%File description:
%       To replot the data in workspace, which from SGA__lib_tether_MMET_II.mdl
%       Orbits are cycles that tether goes around the earth,
%                                          theta = Orbits*(2*pi)
%       Fig-1: alpha ~ Orbits
%       Fig-2: psi   ~ Orbits
%       Fig-3: R     ~ Orbits
%
%        R(theta) = rp*(1+e)/(1+e*cos(theta))
%            where,
%                  rp - distance to periapsis, 6890km in default
%                   e - the magnitude of the eccentricity vector,which belong to [0,1]
%===================================================================================================
%  See Also:   SSTS__plot_position3d
%              SSTS__plot_position
%              SSTS__tether_dumbbell_replot
%              SSTS__tether_MMET_I_replot
%              SSTS__tether_MMET_III_replot
%              SSTS__tether_MMET_III_replot
%===================================================================================================
%
%===================================================================================================
%Revision -
%Date        Name    Description of Change email                 Location
%30-Jan-2007 Yi Chen Initial version       leo.chen.yi@gmail.com Glasgow
%HISTORY$
%==================================================================================================*/

% SSTS__tether_MMET_III_replot Begin

clear
home
close('all');
%set parameters
Fs = 1000;       % Sampling frequency
% t = (0:Fs)/Fs; % One second worth of samples
nfft=1024;
window = hamming(nfft);
noverlap=256;
dflag='none';

%  1 - on
% !1 - off
comet_is_on = 0;

%load data from mat file
load('SSTS__tether_MMET_III.mat');

%set variables according to the order in SSTS__tether_MMET_III.mdl
% mat file always have simulation time as the first line
Time              = SSTS__tether_MMET_III(1,:);

% Psi_upper         = SSTS__tether_MMET_III(2,:);
% Vel_Psi_upper     = SSTS__tether_MMET_III(3,:);
% Acc_Psi_upper     = SSTS__tether_MMET_III(4,:);
% 
% R                 = SSTS__tether_MMET_III(5,:);
% Vel_R             = SSTS__tether_MMET_III(6,:);
% Acc_R             = SSTS__tether_MMET_III(7,:);
% 
% Theta             = SSTS__tether_MMET_III(8,:);
% Vel_Theta         = SSTS__tether_MMET_III(9,:);
% Acc_Theta         = SSTS__tether_MMET_III(10,:);
% 
% dis_upper         = SSTS__tether_MMET_III(11,:);  
% vel_upper         = SSTS__tether_MMET_III(12,:);  
% acc_upper         = SSTS__tether_MMET_III(13,:);  
% 
% Orbits            = SSTS__tether_MMET_III(14,:);
% 
% Psi_lower         = SSTS__tether_MMET_III(15,:);
% Vel_Psi_lower     = SSTS__tether_MMET_III(16,:);
% Acc_Psi_lower     = SSTS__tether_MMET_III(17,:);
% 
% dis_lower         = SSTS__tether_MMET_III(18,:);  
% vel_lower         = SSTS__tether_MMET_III(19,:);  
% acc_lower         = SSTS__tether_MMET_III(20,:);  


%%%%%%%%%%%%%%%%%%%%%%%%%%%% (I) Upper Tether %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%2 - Psi_upper %%%%%%%%%%%%%%%%%%%%
%fig-2 plot psi
Psi_upper         = SSTS__tether_MMET_III(2,:);
Vel_Psi_upper     = SSTS__tether_MMET_III(3,:);
Acc_Psi_upper     = SSTS__tether_MMET_III(4,:);
h=figure(2000);

%Psi_upper
subplot(2,2,1);
grid on
hold on
plot(Time,Psi_upper);
xlabel('Time (Sec.)');
ylabel('\psi (rad)');
title('Upper Tehter \psi');

%Vel_Psi_upper
subplot(2,2,2);
grid on
hold on
plot(Time,Vel_Psi_upper);
xlabel('Time (Sec.)');
ylabel('Vel. of \psi(rad/s)');
title('Upper Tehter Vel. of \psi');

%Accl_Psi_upper
subplot(2,2,3);
grid on
hold on
plot(Time,Acc_Psi_upper);
xlabel('Time (Sec.)');
ylabel('Acc. of \psi(rad/s^2)');
title('Upper Tehter Acc. of \psi');

% Psi_DPsi
subplot(2,2,4);
grid on
hold on
plot(Psi_upper,Vel_Psi_upper);
xlabel('\psi (rad)');
ylabel('d\psi/dt (rad/s) ');
title('Phase Plane for Upper Tehter \psi');

if( comet_is_on == 1 )
    comet(Psi_upper,Vel_Psi_upper);
end
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_psi_upper', 'jpg');
close('all');

% Psi_upper PSD
h=figure(2100);

grid on
hold on

% PSD method - I
% Pxx_DAlpha_Dt = periodogram( DAlpha_Dt );
% % Create a PSD data object.
% Hpsd_DAlpha_Dt = dspdata.psd(Pxx_DAlpha_Dt,'Fs',Fs);
% % Plot the PSD data object.
% plot( Hpsd_DAlpha_Dt );

% PSD method - II
% subplot(2,1,1);
Hs=spectrum.periodogram;
psd(Hs,Psi_upper,'Fs',Fs);
title('PSD of Upper Tehter \psi');
% subplot(2,1,2);
% % Hs=spectrum.periodogram;
% psd(Hs,Psi_upper,'Fs',Fs);
% % axis([xmin xmax ymin ymax])
% axis([0  0.1 -100  0]);
%
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_PSD_psi_upper', 'jpg');
close('all');
%%clear

%%%%%%%%%%%%%%%%%%%% 3 - R %%%%%%%%%%%%%%%%%%%%

R                 = SSTS__tether_MMET_III(5,:);
Vel_R             = SSTS__tether_MMET_III(6,:);
Acc_R             = SSTS__tether_MMET_III(7,:);

%fig-3 plot R
h=figure(3000);
%R
subplot(2,2,1);
grid on
hold on
plot(Time,R);
xlabel('Time (Sec.)');
ylabel('R (m)');

%Vel_R
subplot(2,2,2);
grid on
hold on
plot(Time,Vel_R);
xlabel('Time (Sec.)');
ylabel('Vel. of R (m/s)');

%Acc_R
subplot(2,2,3);
grid on
hold on
plot(Time,Acc_R);
xlabel('Time (Sec.)');
ylabel('Acc. of R (m/s^2)');

% R_DR
subplot(2,2,4);
grid on
hold on
plot(R,Vel_R);
xlabel('R (m)');
ylabel('dR/dt (m/s) ');
title('Phase Plane for R');
if( comet_is_on == 1 )
    comet(R,Vel_R);
end
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_R', 'jpg');
close('all');

% R PSD
h=figure(3100);

grid on
hold on

% R method - I
% Pxx_DAlpha_Dt = periodogram( DAlpha_Dt );
% % Create a PSD data object.
% Hpsd_DAlpha_Dt = dspdata.psd(Pxx_DAlpha_Dt,'Fs',Fs);
% % Plot the PSD data object.
% plot( Hpsd_DAlpha_Dt );

% PSD method - II
% subplot(2,1,1);
Hs=spectrum.periodogram;
psd(Hs,R,'Fs',Fs);
title('PSD of Upper Tehter R');
% subplot(2,1,2);
% % Hs=spectrum.periodogram;
% psd(Hs,R,'Fs',Fs);
% % axis([xmin xmax ymin ymax])
% axis([0  0.1 -100  0]);
%
%
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_PSD_R', 'jpg');
close('all');
%clear


%%%%%%%%%%%%%%%%%% 4 - Theta %%%%%%%%%%%%%%%%%%%%
Theta             = SSTS__tether_MMET_III(8,:);
Vel_Theta         = SSTS__tether_MMET_III(9,:);
Acc_Theta         = SSTS__tether_MMET_III(10,:);

%fig-1 plot Theta
h=figure(4000);

%Theta
subplot(2,2,1);
hold on
grid on
plot( Time,Theta );
xlabel('Time (Sec.)');
ylabel('\theta (rad)');
title('Ang.of \theta');

%Vel_Theta
subplot(2,2,2);
hold on
grid on
plot(Time,Vel_Theta);
xlabel('Time (Sec.)');
ylabel('d\theta/dt(rad/s)');
title('Vel. of \theta');

%Acc_Theta
subplot(2,2,3);
hold on
grid on
plot(Time,Acc_Theta);
xlabel('Time (Sec.)');
ylabel('d(d\theta/dt)/dt(rad/s^2)');
title('Acc.of \theta');

% Theta_DTheta
subplot(2,2,4);
hold on
grid on
plot(Theta,Vel_Theta);
xlabel('\theta(rad)');
ylabel('d\theta/dt(rad/s)');
title('Phase Plane for Tehter \theta');

if( comet_is_on == 1 )
    comet(Theta,Vel_Theta);
end
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_theta_upper', 'jpg');
close('all');

%Theta PSD
h=figure(4100);
grid on
hold on

% PSD method - I
% Pxx_DAlpha_Dt = periodogram( DAlpha_Dt );
% % Create a PSD data object.
% Hpsd_DAlpha_Dt = dspdata.psd(Pxx_DAlpha_Dt,'Fs',Fs);
% % Plot the PSD data object.
% plot( Hpsd_DAlpha_Dt );

% PSD method - II
% subplot(2,1,1);
Hs=spectrum.periodogram;
psd(Hs,Theta,'Fs',Fs);
title('PSD of Upper Tehter \theta');
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_PSD_theta_upper', 'jpg');
close('all');
%clear
% 
% subplot(2,1,2);
% % Hs=spectrum.periodogram;
% psd(Hs,Theta,'Fs',Fs);
% % axis([xmin xmax ymin ymax])
% axis([0  0.1 -100  0]);

%%%%%%%%%%%%%%%%%%%5 - Orbits %%%%%%%%%%%%%%%%%%%
Orbits            = SSTS__tether_MMET_III(14,:);

h=figure(5000);
hold on
grid on
plot( Time,Orbits );
xlabel('Time (Sec.)');
ylabel('Orbits (Cycles)');
title('Orbit ~ time');
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_Orbit', 'jpg');
close('all');
%clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%% (II) Lower Tether %%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%2 - Psi_lower %%%%%%%%%%%%%%%%%%%%
% 
Psi_lower         = SSTS__tether_MMET_III(15,:);
Vel_Psi_lower     = SSTS__tether_MMET_III(16,:);
Acc_Psi_lower     = SSTS__tether_MMET_III(17,:);

%fig-2 plot psi
h=figure(7000);

%Psi_upper
subplot(2,2,1);
grid on
hold on
plot(Time,Psi_lower);
xlabel('Time (Sec.)');
ylabel(' \psi (rad)');
title('Lower Outrigger Ang. of \psi');

%Vel_Psi_upper
subplot(2,2,2);
grid on
hold on
plot(Time,Vel_Psi_lower);
xlabel('Time (Sec.)');
ylabel('Vel. of \psi (rad/s)');
title('Lower Outrigger Vel. of \psi');

%Accl_Psi_upper
subplot(2,2,3);
grid on
hold on
plot(Time,Acc_Psi_lower);
xlabel('Time (Sec.)');
ylabel('Acc. of \psi(rad/s^2)');
title('Lower Outrigger Acc. of \psi');

% Psi_DPsi
subplot(2,2,4);
grid on
hold on
plot(Psi_lower,Vel_Psi_lower);
xlabel('\Psi (rad)');
ylabel('d\Psi/dt (rad/s) ');
title('Phase Plane for Lower Tehter \psi');
if( comet_is_on == 1 )
    comet(Psi_lower,Vel_Psi_lower);
end
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_psi_lower', 'jpg');
close('all');

% Psi_lower PSD
h=figure(7100);

grid on
hold on

% PSD method - I
% Pxx_DAlpha_Dt = periodogram( DAlpha_Dt );
% % Create a PSD data object.
% Hpsd_DAlpha_Dt = dspdata.psd(Pxx_DAlpha_Dt,'Fs',Fs);
% % Plot the PSD data object.
% plot( Hpsd_DAlpha_Dt );

% PSD method - II
% subplot(2,1,1);
Hs=spectrum.periodogram;
psd(Hs,Psi_lower,'Fs',Fs);
title('PSD of Lower Outrigger \psi');
% subplot(2,1,2);
% % Hs=spectrum.periodogram;
% psd(Hs,Psi_upper,'Fs',Fs);
% % axis([xmin xmax ymin ymax])
% axis([0  0.1 -100  0]);
%
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_PSD_psi_lower', 'jpg');
close('all');
%clear

%%%%%%%%%%%%%%%%%%%III -  Alpha vs. Psi_upper %%%%%%%%%%%%%%%%%%%%
% alpha Transfer function
% figure(3000)
% grid on
% hold on
% [ H_vel_alpha , ft ]=tfe(Orbits,DAlpha_Dt,nfft,Fs,window,noverlap,dflag);
% [ H_vel_psi , ft ]  =tfe(Orbits,Vel_Psi_upper,nfft,Fs,window,noverlap,dflag);
% 
% loglog(ft,abs(H_vel_alpha),':r',ft,abs(H_vel_psi),'b');
% xlabel('Frequency(Hz.)');
% ylabel('| (d\alpha/d\theta)/(d\psi/d\theta) | ( dB )');
% legend('(d\alpha/d\theta)','d\psi/d\theta');
% title('| (d\alpha/d\theta)/(d\psi/d\theta) | Transfer Function Estimate')
% axis([0  50 0  40]);


% %psi
% h=figure(9300);
% % subplot(3,1,1);
% grid on
% hold on
% plot(Psi_upper,Psi_lower);
% xlabel('Upper Ang. of \psi (rad)');
% ylabel('Lower Ang. of \psi (rad)');
% title('Upper vs. Lower Ang. of \psi');
% saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_ang_psi_upper_lower', 'jpg');
% close('all');
% 
% h=figure(9400);
% % subplot(3,1,2);
% grid on
% hold on
% plot(Vel_Psi_upper,Vel_Psi_lower);
% xlabel('Upper Vel. of \psi (rad/s)');
% ylabel('Lower Vel. of \psi (rad/s)');
% title('Upper vs. Lower Vel. of \psi '); 
% saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_vel_psi_upper_lower', 'jpg');
% close('all');
% 
% h=figure(9500);
% % subplot(3,1,3);
% grid on
% hold on
% plot(Acc_Psi_upper,Acc_Psi_lower);
% xlabel('Upper Acc. of \psi');
% ylabel('Lower Acc. of \psi');
% title('Upper vs. Lower Acc. of \psi ');
% 
% if( comet_is_on == 1 )
%     %display comet plot
%     %upper vs. lower
%     subplot(3,2,1);
%     comet(Psi_upper,Psi_lower);
% 
%     subplot(3,2,2);
%     comet(Vel_Psi_upper,Vel_Psi_lower);
% 
%     subplot(3,2,3);
%     comet(Acc_Psi_upper,Accl_Psi_lower);
%     
%     %upper
% %     subplot(3,2,1);
% %     comet(Alpha_upper,Psi_upper);
% % 
% %     subplot(3,2,2);
% %     comet(Vel_Alpha_upper,Vel_Psi_upper);
% % 
% %     subplot(3,2,3);
% %     comet(Acc_Alpha_upper,Accl_Psi_upper);
% end
% saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_acc_psi_upper_lower', 'jpg');
% close('all');
% %clear

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% IV axial motion %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time domain
dis_upper         = SSTS__tether_MMET_III(11,:);  
vel_upper         = SSTS__tether_MMET_III(12,:);  
acc_upper         = SSTS__tether_MMET_III(13,:);  

%displacement 
h=figure(100);
grid on
hold on
plot(Time,dis_upper,'-.',Time,dis_lower,'--');
xlabel('Time (Sec.)');
ylabel('Displacement (m)');
title('Axial Displacement ');
legend('upper','lower');
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_axial_dis', 'jpg');
close('all');

%Velocity
h=figure(200);
grid on
hold on
plot(Time,vel_upper,'-.',Time,vel_lower,'--');
xlabel('Time (Sec.)');
ylabel('Velocity (m/s)');
title('Axial Velocity ');
legend('upper','lower');
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_axial_vel', 'jpg');
close('all');

%Acceleration
h=figure(300);
grid on
hold on
plot(Time,acc_upper,'-.',Time,acc_lower,'--');
xlabel('Time (Sec.)');
ylabel('Acceleration (m/s^2)');
title('Axial Acceleration ');
legend('upper','lower');
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_axial_acc', 'jpg');
 close('all');
 
 % Frequency domain
 %upper
h=figure(400);
[H_acc_upper,ft]=tfe(Acc_Psi_upper,acc_upper,nfft,Fs,window,noverlap,dflag);
%lower
[H_acc_lower,ft]=tfe(Acc_Psi_lower,acc_lower,nfft,Fs,window,noverlap,dflag);


loglog(ft,abs(H_acc_upper),':r' ,ft,abs(H_acc_lower),'b');
xlabel('frequency(Hz.)');
ylabel('|axial acceleration / acceleration of \psi|(dB.)');
title('Transformation function of  ''axial acceleration / acceleration of \psi''');
legend('upper','lower');
saveas(h,'.\MMET-III-JPG\SSTS__tether_MMET_III_replot_axial_acc_feq', 'jpg');
close('all');


home
clear


% SSTS__tether_MMET_III_replot End