%

%			2-D Orbit Visualization

%

%			   Michael Mulford

%			      ESM 4984

%			      7/24/98

%			

%

%  Description:

%

%	This program will prompt the user for some parameters describing 

%  path of a celestial body, and will generate an image of: the central planet,

%  the orbiting body and its path, as well as separate plot of the velocity

%  of the orbiting body.  The length units are nondimensionalzed based on the

%  radius of the planet that the body is orbiting. 





  clear all



%  Prompt the user for a, e and nu

%  Press enter to choose the default values



  a = input('Please enter a semimajor axis length (DU) :  ');

  e = input('Please enter the eccentricity (0-0.99) :  ');

  nue = input('Please enter the true anomaly (degrees) :  ');

  nu = 0:1:360;



  if isempty(a) a = 2;  end

  if isempty(e) e = .3; end

  if isempty(nue) nue = 33;  end



% Change the angles to radians



  cc = 3.14159265358979/180;

  nue = nue*cc;	nur = nu*cc;



%  Generate some values to draw the circular planet



  n = linspace(0,2*pi,40);

  Y = sin(n);

  X = cos(n);



%  Calculate the epoch position



  ra = a*(1-e*e)/(1+e*cos(nue));

  rax = ra*cos(nue);

  ray = ra*sin(nue);



%  Use a 'for loop' to generate the path of the orbit in the 

%  ecliptic plane, also calculate the velocity at each point



  for i=1:length(nu)

	r(i) = a*(1-e*e)/(1+e*cos(nur(i)));

	rx(i) = r(i)*cos(nur(i));

	ry(i) = r(i)*sin(nur(i));

	v(i) = (2*(1/r(i) - 1/(2*a)))^.5;

  end



%  Plot the orbital path, the epoch position and the planet



  plot(rx,ry)

  hold on

  plot(rax,ray,'wo')

  hold on

  plot(X,Y)

  fill(X,Y,'b')

  axis('equal')

  

% Draw P and Q axes



  f = a/2 + 1;

  Px = [0,f];Py = [0,0];

  Qx = [0,0];Qy = [0,f];

  plot(Px,Py,'r')

  plot(Qx,Qy,'g')



  pause



% Ask the user if they want to see a graph of the velocity



  qu = input('Would you like to plot the velocity versus position?  (y/n) : ','s');



% Plot the velocity versus position



  if qu == 'y'

        figure

        axis('normal')

	plot(nu,v)

        set(gca,'xlim',[0,360])

	%title('Nondimensionalized Velocity of the Celestial Body versus Position (Degrees)')

  end



%  Option to see the orbit animated

%  comet(rx,ry)