ssts__tether_mmet_i_replot.m

空间绳系卫星的仿真软件

% /*M-FILE Script SSTS__tether_MMET_I_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 SSTS__tether_MMET_I.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_III_replot
%===================================================================================================
%
%===================================================================================================
%Revision -
%Date        Name    Description of Change email                 Location
%30-Nov-2006 Yi Chen Initial version       [e-mail] Glasgow
%HISTORY$
%==================================================================================================*/

% SSTS__tether_MMET_I_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_I.mat');

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

Alpha_upper       = SSTS__tether_MMET_I(2,:);
Vel_Alpha_upper   = SSTS__tether_MMET_I(3,:);
Acc_Alpha_upper   = SSTS__tether_MMET_I(4,:);

Psi_upper         = SSTS__tether_MMET_I(5,:);
Vel_Psi_upper     = SSTS__tether_MMET_I(6,:);
Acc_Psi_upper    = SSTS__tether_MMET_I(7,:);

R                 = SSTS__tether_MMET_I(8,:);
Vel_R             = SSTS__tether_MMET_I(9,:);
Acc_R             = SSTS__tether_MMET_I(10,:);

Alpha_lower       = SSTS__tether_MMET_I(11,:);
Vel_Alpha_lower   = SSTS__tether_MMET_I(12,:);
Acc_Alpha_lower   = SSTS__tether_MMET_I(13,:);

Psi_lower         = SSTS__tether_MMET_I(14,:);
Vel_Psi_lower     = SSTS__tether_MMET_I(15,:);
Acc_Psi_lower    = SSTS__tether_MMET_I(16,:);

Orbits            = SSTS__tether_MMET_I(17,:);

Theta             = SSTS__tether_MMET_I(18,:);
Vel_Theta         = SSTS__tether_MMET_I(19,:);
Acc_Theta         = SSTS__tether_MMET_I(20,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%% (I) Upper Tether %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%  1 - Alpha upper %%%%%%%%%%%%%%%%%%%%
%fig-1 plot alpha
figure(1000)
hold on
grid on

%alpha
subplot(2,2,1);
hold on
grid on
plot( Time,Alpha_upper );
xlabel('Time (Sec.)');
ylabel('\alpha (rad)');
Title('Upper Tehter \alpha');

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

%Acc alpha
subplot(2,2,3);
hold on
grid on
plot(Time,Acc_Alpha_upper);
xlabel('Time (Sec.)');
ylabel('Acc. of \alpha (rad/s^2)');
Title('Upper Tehter Acc. of \alpha');

% alpha_Dalpha
subplot(2,2,4);
hold on
grid on
plot(Alpha_upper,Vel_Alpha_upper);
xlabel('\alpha(rad) ');
ylabel('d\alpha/dt(rad)');
Title('Phase Plane for Upper Tehter \alpha');
if( comet_is_on == 1 )
    comet(Alpha_upper,Vel_Alpha_upper);
end

%Alpha PSD
figure(1100)
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,Alpha_upper,'Fs',Fs);
Title('PSD of Upper Tehter \alpha');
% 
% subplot(2,1,2);
% % Hs=spectrum.periodogram;
% psd(Hs,Alpha_upper,'Fs',Fs);
% % axis([xmin xmax ymin ymax])
% axis([0  0.1 -100  0]);

%%%%%%%%%%%%%%%%%%%%2 - Psi_upper %%%%%%%%%%%%%%%%%%%%
%fig-2 plot psi
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


% Psi_upper PSD
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]);
%
%%%%%%%%%%%%%%%%%%%% 3 - R %%%%%%%%%%%%%%%%%%%%
%fig-3 plot R
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

% R PSD
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]);
%
%
%%%%%%%%%%%%%%%%%% 4 - Theta %%%%%%%%%%%%%%%%%%%%
%fig-1 plot Theta
figure(4000)

%Theta
subplot(2,2,1);
hold on
grid on
plot( Time,Theta );
xlabel('Time (Sec.)');
ylabel('\theta (rad)');
Title('Ang. of Tehter \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 Tehter \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 Tehter \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 Upper Tehter \theta');

if( comet_is_on == 1 )
    comet(Theta,Vel_Theta);
end

%Theta PSD
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');
% 
% 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 %%%%%%%%%%%%%%%%%%%
figure(5000)
hold on
grid on
plot( Time,Orbits );
xlabel('Time (Sec.)');
ylabel('Orbits (Cycles)');
title('Orbit ~ time');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% (II) Lower Tether %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%  1 - Alpha lower %%%%%%%%%%%%%%%%%%%%
%fig-1 plot alpha
figure(6000)

%alpha
subplot(2,2,1);
hold on
grid on
plot( Time,Alpha_lower );
xlabel('Time (Sec.)');
ylabel('\alpha (rad)');
Title('Lower Outrigger Ang. of \alpha');

%Vel alpha
subplot(2,2,2);
hold on
grid on
plot(Time,Vel_Alpha_lower);
xlabel('Time (Sec.)');
ylabel('Vel. of \alpha(rad/s)');
Title('Lower Outrigger Vel. of \alpha');
%Acc alpha
subplot(2,2,3);
hold on
grid on
plot(Time,Acc_Alpha_lower);
xlabel('Time (Sec.)');
ylabel('Acc. of \alpha (rad/s^2)');
Title('Lower Outrigger Acc. of \alpha');

% alpha_Dalpha
subplot(2,2,4);
hold on
grid on
plot(Alpha_lower,Vel_Alpha_lower);
xlabel('\alpha(rad) ');
ylabel('d\alpha/dt(rad)');
Title('Phase Plane for lower Tehter \alpha');
if( comet_is_on == 1 )
    comet(Alpha_lower,Vel_Alpha_lower);
end

%Alpha PSD
figure(6100)

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,Alpha_lower,'Fs',Fs);
Title('Lower Outrigger \alpha');
% 
% subplot(2,1,2);
% % Hs=spectrum.periodogram;
% psd(Hs,Alpha_upper,'Fs',Fs);
% % axis([xmin xmax ymin ymax])
% axis([0  0.1 -100  0]);

%%%%%%%%%%%%%%%%%%%%2 - Psi_lower%%%%%%%%%%%%%%%%%%%%
%fig-2 plot psi
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


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

%%%%%%%%%%%%%%%%%%%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]);


% Upper alpha vs. psi
figure(8000)
% subplot(3,1,1);
grid on
hold on
plot(Alpha_upper,Psi_upper);
xlabel('Upper Ang. of \alpha (rad)');
ylabel('Upper Ang. of \psi   (rad)');
Title('Upper Tether Ang. of \alpha ~ \psi ');

figure(8100)
% subplot(3,1,2);
grid on
hold on
plot(Vel_Alpha_upper,Vel_Psi_upper);
xlabel('Upper Vel .of \alpha (rad/s) ');
ylabel('Upper Vel .of \psi   (rad/s)');
Title('Upper Tether Vel.of \alpha ~ \psi');

figure(8200)
% subplot(3,1,3);
grid on
hold on
plot(Acc_Alpha_upper,Acc_Psi_upper);
xlabel('Upper Acc.of \alpha (rad/s^2) ');
ylabel('Upper Acc.of \psi   (rad/s^2)');
Title('Upper Tether Acc.of \alpha ~ \psi ');

if( comet_is_on == 1 )
    %display comet plot
    subplot(3,1,1);
    comet(Alpha_upper,Psi_upper);

    subplot(3,1,2);
    comet(Vel_Alpha_upper,Vel_Psi_upper);

    subplot(3,1,3);
    comet(Acc_Alpha_upper,Accl_Psi_upper);
end

% Lower vs. upper alpha vs. psi
%alpha
figure(9000)
% subplot(3,1,1);
grid on
hold on
plot(Alpha_upper,Alpha_lower);
xlabel('Upper Ang. of \alpha');
ylabel('Lower Ang. of \alpha');
Title('Upper vs. Lower Ang. of \alpha');

figure(9100)
% subplot(3,1,2);
grid on
hold on
plot(Vel_Alpha_upper,Vel_Alpha_lower);
xlabel('Upper Vel.of \alpha (rad/s)');
ylabel('Lower Vel.of \alpha (rad/s)');
Title('Upper vs. Lower Vel.of \alpha'); 

figure(9200)
% subplot(3,1,3);
grid on
hold on
plot(Acc_Alpha_upper,Acc_Alpha_lower);
xlabel('Upper Acc. of \alpha (rad/s^2)');
ylabel('Lower Acc. of \alpha (rad/s^2)');
Title('Upper vs. Lower Acc. of \alpha' );

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

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

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

home
clear


% SSTS__tether_MMET_I_replot End