function nma_runProbeSimulation
%runs simulation of probe starting from some position and velosity vector for some delta time

%
% runs a simple simulation of a probe starting from some
% position and velosity vector for some delta time. 
%
% Author: Nasser Abbasi
% April 30, 2003.

% change history
% name   date   changes
% ----- -----   --------
% nma    050303  add ecc. vector to plot
%

TRUE=1;
FALSE=0;

EARTH_CANONICAL=1;
ENGLISH=2;
METRIC=3;

r0_vec = input('Enter initial position vector r0  [DU] >');
v0_vec = input('Enter initial velosity vector v0  [DU/TU]  >');
mu     = input('Enter mu (gravitational parameter) >');
%plotMode=input('plot mode? [1=show orbit trace, 2=show radial vector] >');
simulationTime = input('How long to run the simulation? [hrs]  >');
nSteps=input('number of steps? >');

% support other units later
%units=input('Which units? [1=earth canonical, 2=English, 3=Metric] >');

%
% covert all to canonical, then flip back at the end if needed
%
%TC=1;
%if(units ~=EARTH_CANONICAL)
%    switch(units)
%        case ENGLISH
%           simulationTime=


TC=1.239446309*10^-3;
simulationTime = simulationTime*60*60*TC;

% CIRCLE
% r0_vec=[0 1 .2];
% v0_vec=[.9 0 .123];

% hyparabola?
% r0_vec=[0 1.1 0 ];
% v0_vec=[1 0 1.4];

% good ELLIPS on the x-y plane only. nice sim do with no pause
%r0_vec=[1.2 0 0];
%v0_vec=[sqrt(2) 0 0]


% ELLIPS on the x-y plane only. nice sim do with no pause
%r0_vec=[1.2 0 0];
%v0_vec=[.1 1 0];

% ELLIPS on the x-y plane only. nice sim do with no pause
% this shows more clearly kepler second law. 
% use mu=5 to see it better
%r0_vec=[1.2 .9 0];
%v0_vec=[.1 1 0];


%50 hrs, 1000 steps  *******************
%r0_vec=[1.2 .9 .2];
%v0_vec=[.1 1 0];

% problem 1, Dr Orme given
%r0_vec=[0 1.1 0];
%v0_vec=[sqrt(2) 0 0 ];


%fprintf('r0 = %6.6f I + %6.6f J + %6.6f K\n',r0_vec(1),r0_vec(2),r0_vec(3));
%fprintf('v0 = %6.6f I + %6.6f J + %6.6f K\n',v0_vec(1),v0_vec(2),v0_vec(3));

fprintf('\nProbe is moving .....\n\n');


time_ = linspace(0,simulationTime,nSteps);
timeStepDuration = time_(2)-time_(1);

rx = zeros(length(time_),1);
ry = rx;
rz = rx;

vx = rx;
vy = rx;
vz = rx;

rx(1)=r0_vec(1);
ry(1)=r0_vec(2);
rz(1)=r0_vec(3);
vx(1)=v0_vec(1);
vy(1)=v0_vec(2);
vz(1)=v0_vec(3);

for(i=2:length(time_)-1)
    [r_vec,v_vec] = nma_moveProbe(r0_vec,v0_vec,timeStepDuration,mu);    
    
    rx(i)=r_vec(1);
    ry(i)=r_vec(2);
    rz(i)=r_vec(3);
    
    vx(i)=v_vec(1);
    vy(i)=v_vec(2);
    vz(i)=v_vec(3);
    
    r0_vec=r_vec;
    v0_vec=v_vec;    
end

figure;

% enable double buffering to reduce flickering.
set(gcf,'DoubleBuffer','on')
plot3(0,0,0,'*');
hold on;


if( abs(min(rx)-max(rx))<eps)
    startX=-Inf; endX=+Inf;
else
    startX=min(rx); endX=max(rx);
end

if( abs(min(ry)-max(ry))<eps)
    startY=-Inf; endY=+Inf;
else
    startY=min(ry); endY=max(ry);
end

if( abs(min(rz)-max(rz))<eps)
    startZ=-Inf; endZ=+Inf;
else
    startZ=min(rz); endZ=max(rz);
end

o=nma_getOrbitParams(r0_vec,v0_vec,mu);

periodInMin  = o.period*13.44686457;
periodHrs    = floor(periodInMin/60);
periodMins   = periodInMin-(periodHrs*60);

S2=sprintf('E=%6.2f, e=%2.6f, h=%6.2f, p=%6.2f, a=%6.3f,  Period=%6.3f DU [%4.0f hr:%2.0f min]', ...
    o.energy,norm(o.e_vec),norm(o.h_vec), o.p, o.a, o.period,periodHrs,periodMins);

    
h_radialLine = 0;
h_probe=0;
timeThisCycle=0;

plot3(0,0,0,'*');
hold on;
plot3( [0;o.rp_vec(1)] ,  [0;o.rp_vec(2)] ,  [0;o.rp_vec(3)],'k-');
text(o.rp_vec(1)/2 ,  o.rp_vec(2)/2 ,  o.rp_vec(3)/2,'rp');

plot3( [0;o.ra_vec(1)] ,  [0;o.ra_vec(2)] ,  [0;o.ra_vec(3)],'m-');
%plot3( o.ra_vec(1) ,  o.ra_vec(2) ,  o.ra_vec(3),'md');
text( o.ra_vec(1)/2 ,  o.ra_vec(2)/2 ,  o.ra_vec(3)/2,'ra');

view(45,45);
xlabel('x-axis');
ylabel('y-axis');
zlabel('z-axis');
keepPlottingTrace=TRUE;

for(i=1:length(rx)-1)
    
    if(h_radialLine ~= 0 )
        delete(h_radialLine);
        h_radialLine=0;
    end
    
    h_radialLine=plot3([0; rx(i)] , [0;ry(i)] , [0;rz(i)],'r'); % shows radial vector from origin.    
    
    if(keepPlottingTrace)
        plot3( rx(i) ,  ry(i) ,  rz(i),'+');
    end    
    
    if(h_probe ~=0)
        delete(h_probe);
    end    
    
    h_probe=plot3( rx(i) ,  ry(i) ,  rz(i),'o');
    
    currentTime  = i*timeStepDuration;
    currentTime  = (currentTime/TC);  %back to secs
    nhrs         = floor(currentTime/(60*60));
    nmin         = floor(  (currentTime- nhrs*60*60)/60 );
    nsec         = currentTime - ( (nhrs*60*60)+nmin*60 );
    
    S1=sprintf('t=%.0f:%2.0f:%2.0f ,r=%6.3f , r=[%6.3f %6.3f  %6.3f], v=%6.3f , v=[%6.3f  %6.3f  %6.3f]',...
        nhrs,nmin,nsec,norm([rx(i) ry(i) rz(i)]),rx(i),ry(i),rz(i),norm([vx(i) vy(i) vz(i)]),vx(i),vy(i),vz(i));
    
    title({S1 S2}); 
    axis([startX endX startY endY startZ endZ ]);
    grid on;
    
    drawnow;
    
    % to slow the simulation :)
    %    for(j=1:2*1000000)
    %        k=j+2;
    %    end
    
    timeThisCycle = timeThisCycle + timeStepDuration;
    if(abs(real(o.period)) >eps)
        if(timeThisCycle >= real(o.period) )
            timeThisCycle=0;
            keepPlottingTrace=FALSE;
        end
    end
end