ssts__tether_dumbbell_replot.m

空间绳系卫星的仿真软件

% /*M-FILE Script SSTS__tether_dumbbell_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_dumbbell.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
%
%===================================================================================================
%
%===================================================================================================
%Revision -
%Date        Name    Description of Change email                 Location
%27-Nov-2006 Yi Chen Initial version       [e-mail] Glasgow
%04-Dec-2006 Yi Chen add comet and PSD     [e-mail] Glasgow
%HISTORY$
%==================================================================================================*/

% SSTS__tether_dumbbell_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_dumbbell.mat');

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

Alpha       = SSTS_tether_dumbbell(2,:);
DAlpha_Dt   = SSTS_tether_dumbbell(3,:);
DDAlpha_DDt = SSTS_tether_dumbbell(4,:);

Psi         = SSTS_tether_dumbbell(5,:);
DPsi_Dt     = SSTS_tether_dumbbell(6,:);
DDPsi_DDt   = SSTS_tether_dumbbell(7,:);

R           = SSTS_tether_dumbbell(8,:);
DR_Dt       = SSTS_tether_dumbbell(9,:);
DDR_DDt     = SSTS_tether_dumbbell(10,:);

Orbits      = SSTS_tether_dumbbell(11,:);

Theta       = SSTS_tether_dumbbell(12,:);
DTheta_Dt   = SSTS_tether_dumbbell(13,:);
DDTheta_DDt = SSTS_tether_dumbbell(14,:);

%%%%%%%%%%%%%%%%%%  1 - Alpha %%%%%%%%%%%%%%%%%%%%
%fig-1 plot alpha
figure(1000)

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

%Dalpha_Dt
subplot(2,2,2);
hold on
grid on
plot(Time,DAlpha_Dt);
xlabel('Time (Sec.)');
ylabel('d\alpha/dt');

%DDalpha_DDt
subplot(2,2,3);
hold on
grid on
plot(Time,DDAlpha_DDt);
xlabel('Time (Sec.)');
ylabel('d(d\alpha/dt)/dt');

% alpha_Dalpha
subplot(2,2,4);
hold on
grid on
plot(Alpha,DAlpha_Dt);
xlabel('\alpha ');
ylabel('d\alpha/dt');

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

%%%%%%%%%%%%%%%%%%%%2 - Psi %%%%%%%%%%%%%%%%%%%%
%fig-2 plot psi
figure(2000)

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

%DPsi_Dt
subplot(2,2,2);
grid on
hold on
plot(Time,DPsi_Dt);
xlabel('Time (Sec.)');
ylabel('d\psi/dt');

%DDPsi_DDt
subplot(2,2,3);
grid on
hold on
plot(Time,DDPsi_DDt);
xlabel('Time (Sec.)');
ylabel('d(d\psi/dt)/dt');

% Psi_DPsi
subplot(2,2,4);
grid on
hold on
plot(Psi,DPsi_Dt);
xlabel('\Psi');
ylabel('d\Psi/dt ');
if( comet_is_on == 1 )
    comet(Psi,DPsi_Dt);
end


% Psi 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,'Fs',Fs);
Title('PSD of \psi');
% subplot(2,1,2);
% % Hs=spectrum.periodogram;
% psd(Hs,Psi,'Fs',Fs);
% % axis([xmin xmax ymin ymax])
% axis([0  0.1 -100  0]);
%
%%%%%%%%%%%%%%%%%%%% 3 - R %%%%%%%%%%%%%%%%%%%%
%fig-3 plot R
figure(3000)
%Psi
subplot(2,2,1);
grid on
hold on
plot(Time,R);
xlabel('Time (Sec.)');
ylabel('R (m)');

%DR_Dt
subplot(2,2,2);
grid on
hold on
plot(Time,DR_Dt);
xlabel('Time (Sec.)');
ylabel('dR/dt (m/s)');

%DDR_DDt
subplot(2,2,3);
grid on
hold on
plot(Time,DDPsi_DDt);
xlabel('Time (Sec.)');
ylabel('d(dR/dt)/dt (m/s^2)');

% R_DR
subplot(2,2,4);
grid on
hold on
plot(R,DR_Dt);
xlabel('R');
ylabel('dR/dt ');
if( comet_is_on == 1 )
    comet(R,DR_Dt);
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 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)');

%DTheta_Dt
subplot(2,2,2);
hold on
grid on
plot(Time,DTheta_Dt);
xlabel('Time (Sec.)');
ylabel('d\theta/dt');

%DDTheta_DDt
subplot(2,2,3);
hold on
grid on
plot(Time,DDTheta_DDt);
xlabel('Time (Sec.)');
ylabel('d(d\theta/dt)/dt');

% Theta_DTheta
subplot(2,2,4);
hold on
grid on
plot(Theta,DTheta_Dt);
xlabel('\theta ');
ylabel('d\theta/dt');

if( comet_is_on == 1 )
    comet(Theta,DTheta_Dt);
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 \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)');

%%%%%%%%%%%%%%%%%%%Other -  Alpha vs. Psi %%%%%%%%%%%%%%%%%%%%
% 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,DPsi_Dt,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]);


% alpha vs. psi
figure(6000)

subplot(3,1,1);
grid on
hold on
plot(Alpha,Psi);
xlabel('\alpha (rad)');
ylabel('\psi(rad)');
title('Angle');

subplot(3,1,2);
grid on
hold on
plot(DAlpha_Dt,DPsi_Dt);
xlabel('d\alpha/dt ');
ylabel('d\psi/dt');
title('Velocity');

subplot(3,1,3);
grid on
hold on
plot(DDAlpha_DDt,DDPsi_DDt);
xlabel('d(d\alpha/dt)/dt ');
ylabel('d(d\psi/dt)/dt');
title('Acceleration');

if( comet_is_on == 1 )
    %display comet plot
    subplot(3,1,1);
    comet(Alpha,Psi);

    subplot(3,1,2);
    comet(DAlpha_Dt,DPsi_Dt);

    subplot(3,1,3);
    comet(DDAlpha_DDt,DDPsi_DDt);
end



home
clear

% % SSTS__tether_dumbbell_replot End