%----------------------------------------------------------------
%       Phase Plane Trajectory(位相平面軌跡)
%  
%                   
%     q'' +  2 te w q' + w^2 q=0
%                  
%     x1:=q   x2:=q'
%     then
%     x'1=x2
%     x'2=-2 te w x2 - w^2 x1
% 
%       w=1  te=0.2 
%       x1(0)=-5  x2(0)=5
%
%       Notice!
%       call trajectory1_eq0.m
%
%            2006.3.31                          produced by kizawa
%--------------------------------------------------------------------------

clc
clear all
format long e

global ni;

%----setup---%
sampling=0.01;%sampling time of simulation 
%kaisu=1000;
%timespan=sampling*kaisu;% simulation time
%--- initial state ----
x1=-5;%[rad]
x2=+5;%[rad/s]        

x0=[x1  x2 ]; %initial condition


%---mechanical parameter---
wn=1;        
te =0.2; 


%-----------------  Non linear system------------------------------------------------
%---definition---   
%  y(1,1)   y(2,1)        
%           dy(1,1)        
%-------------------------------------------------------------------------------------

%  Y(1,1)     Y(2,1)               
%    x1          x2
%--------------------------------------------------------------------------

disp('It is under calculation now!!')
ni=1;
tend=10000;  %6000   応答時間に合わせて変える
flag=0;
while(ni <= tend) 
    
             [Tn,Yn]=ode45(@trajectory2_eq0,[0 sampling],x0,[],wn,te); %変えること
                    
                                                                                                    
           if ni==1
              YY(1,:)=Yn(1,1:2); 
              TT(1,:)=Tn(1,1); 
             % VIN(1,:)=Vin;
             
           else    
             YY(ni,:)=Yn(size(Yn,1),1:2);
             TT(ni,1)=sampling*ni; 
             %VIN(ni,1)=Vin;
         
           end
           %---- Up data------
           size(Yn,1);
           x0=Yn(size(Yn,1),1:2);
           ni=ni+1;
           
           
end 


figure(1)
plot(YY(:,1) ,YY(:,2));
ax1=findobj(gca,'Type','Axes','Visible','on');
set(ax1,'fontsize',12);  %目盛りを大きく
axis([-10  +10 -8  +8])
xlabel('y_1','fontsize',12);
ylabel('y_2','fontsize',12)
title('Phase Plane Trajectory','fontsize',12)
text(-5,6,'( y_1(0),y_2(0) )=(-5,5)' );
hold on
plot([-10 10],[0 0],'k')
plot([0 0],[-8 8],'k')

hold off








