ssts__tether_mmet_ii_replot.m

空间绳系卫星的仿真软件

% /*M-FILE Script SSTS__tether_MMET_II_replot MMM SSTSLAB */
% /*==================================================================================================
%  Simple Space Tether Simulation Laboratory Toolbox for Matlab 7.x
%
%  Copyright 2007 The SxLAB Family - Yi Chen - [e-mail]
% ====================================================================================================
%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_II_replot
%              SSTS__tether_MMET_II_replot
%===================================================================================================
%
%===================================================================================================
%Revision -
%Date        Name    Description of Change email                 Location
%30-Nov-2006 Yi Chen Initial version       [e-mail] Glasgow
%HISTORY$
%==================================================================================================*/

% SSTS__tether_MMET_II_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_II.mat');

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

Psi_upper         = SSTS__tether_MMET_II(2,:);
Vel_Psi_upper     = SSTS__tether_MMET_II(3,:);
Acc_Psi_upper     = SSTS__tether_MMET_II(4,:);

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

Theta             = SSTS__tether_MMET_II(8,:);
Vel_Theta         = SSTS__tether_MMET_II(9,:);
Acc_Theta         = SSTS__tether_MMET_II(10,:);

Orbits            = SSTS__tether_MMET_II(11,:);

Psi_lower         = SSTS__tether_MMET_II(12,:);
Vel_Psi_lower     = SSTS__tether_MMET_II(13,:);
Acc_Psi_lower     = SSTS__tether_MMET_II(14,:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%% (I) Upper Tether %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%2 - Psi_upper %%%%%%%%%%%%%%%%%%%%
%fig-2 plot psi
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-II-JPG\SSTS__tether_MMET_II_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-II-JPG\SSTS__tether_MMET_II_replot_PSD_Upper_Tehter', 'jpg');
close('all');

%%%%%%%%%%%%%%%%%%%% 3 - R %%%%%%%%%%%%%%%%%%%%
%fig-3 plot R
h=figure(3000);
%Psi_upper
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-II-JPG\SSTS__tether_MMET_II_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-II-JPG\SSTS__tether_MMET_II_replot_PSD_R', 'jpg');
close('all');
%
%%%%%%%%%%%%%%%%%% 4 - Theta %%%%%%%%%%%%%%%%%%%%
%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-II-JPG\SSTS__tether_MMET_II_replot_theta', '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-II-JPG\SSTS__tether_MMET_II_replot_PSD_theta', 'jpg');
close('all');
% 
% 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 %%%%%%%%%%%%%%%%%%%
h=figure(5000);
hold on
grid on
plot( Time,Orbits );
xlabel('Time (Sec.)');
ylabel('Orbits (Cycles)');
title('Orbit ~ time');

saveas(h,'.\MMET-II-JPG\SSTS__tether_MMET_II_replot_orbit', 'jpg');
close('all');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% (II) Lower Tether %%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%Psi_lower
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_lower
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_lower
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-II-JPG\SSTS__tether_MMET_II_replot_psi_lower', 'jpg');
close('all');

% Psi_upper 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');

saveas(h,'.\MMET-II-JPG\SSTS__tether_MMET_II_replot_PSD_psi_lower', 'jpg');
close('all');
% subplot(2,1,2);
% % Hs=spectrum.periodogram;
% psd(Hs,Psi_upper,'Fs',Fs);
% % axis([xmin xmax ymin ymax])
% axis([0  0.1 -100  0]);
%

%%%%%%%%%%%%%%%%%%%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-II-JPG\SSTS__tether_MMET_II_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-II-JPG\SSTS__tether_MMET_II_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,Acc_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,Acc_Psi_upper);
end

saveas(h,'.\MMET-II-JPG\SSTS__tether_MMET_II_replot_acc_psi_upper_lower', 'jpg');
close('all');

home
clear
close('all');
disp('SSTS__tether_MMET_II_replot() done...')

% SSTS__tether_MMET_II_replot End