%

%			3-D Orbit Visualization

%

%			   Michael Mulford

%			      ESM 4984

%			      7/24/98

%			

%

%  Description:

%

%	This program will prompt the user for the orbital elements used in 

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

%  the central planet, the orbiting body and its path, as well as several 

%  important vectors.  The length units are nondimensionalzed based on the

%  radius of the planet that the body is orbiting.  Several viewing options

%  are available in order to insure the optimal visualization for each

%  individual user.  



clear all



% Prompt the user to input the orbital elements

% Press enter to use the default values



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

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

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

  Or = input('Please enter the longitude of the ascending node (degrees) : ');

  wr = input('Please enter the argument of periapsis (degrees) : ');

  ir = input('Please enter the inclination (degrees) : ');



  if isempty(a) a = 2;  end

  if isempty(e) e = .3; end

  if isempty(nuar) nuar = 33;  end

  if isempty(Or) Or = 33;  end

  if isempty(wr) wr = 66;  end

  if isempty(ir) ir = 30;  end



% Prompt the user for a visualization preference



    qu = input('Would you like to see the XYZ axes drawn at (0,0,0)?  (y/n) : ','s');



% Change the angles to radians



  cc = 3.14159265358979/180;

  nua = nuar*cc;	O = Or*cc;	w = wr*cc;	i = ir*cc;



% Create an array of values for the orbit generation loop



  nu = 0:1:360;

  nur = nu*cc;



% Create a rotation matrix that will transform vectors 

% from the PQW to the IJK reference frame



  r11 = cos(O)*cos(w) - sin(O)*sin(w)*cos(i);

  r12 = -cos(O)*sin(w) - sin(O)*cos(w)*cos(i);

  r13 = sin(O)*sin(i);

  r21 = sin(O)*cos(w) + cos(O)*sin(w)*cos(i);

  r22 = -sin(O)*sin(w) + cos(O)*cos(w)*cos(i);

  r23 = -cos(O)*sin(i);

  r31 = sin(w)*sin(i);

  r32 = cos(w)*sin(i);

  r33 = cos(i);

  R = [r11 r12 r13;r21 r22 r23;r31 r32 33];



% Use a 'for' loop to generate the position of the body throughout its orbit

% and store the values in arrays



  for m=1:length(nu)

	rr(m) = a*(1-e*e)/(1+e*cos(nur(m)));

	rx(m) = rr(m)*cos(nur(m));

	ry(m) = rr(m)*sin(nur(m));

	v(m) = (2*(1/rr(m) - 1/(2*a)))^.5;

	ri(m) = r11*rx(m) + r12*ry(m);

	rj(m) = r21*rx(m) + r22*ry(m);

	rz(m) = r31*rx(m) + r32*ry(m);

  end



% Compute the epoch position



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

  rax = ra*cos(nua);

  ray = ra*sin(nua);

  va = (2*(1/ra - 1/(2*a)))^.5;

  ria = r11*rax + r12*ray;

  rja = r21*rax + r22*ray;

  rza = r31*rax + r32*ray;



% Draw the planet as a unit sphere  



  colormap(copper);

  sphere

  hold on



% Draw P,Q and h axes



  f = a/2+1.5; j = a/2+1;

  h = f*[r13,r23,r33];

  Q = f*[r12,r22,r32];

  E = f*[r11,r21,r31];

  h1= [0, h(1)]; h2 = [0, h(2)]; h3 = [0, h(3)];

  E1= [0, E(1)]; E2 = [0, E(2)]; E3 = [0, E(3)];

  Q1= [0, Q(1)]; Q2 = [0, Q(2)]; Q3 = [0, Q(3)];

  plot3(h1,h2,h3,'y')

  plot3(E1,E2,E3,'r')

  plot3(Q1,Q2,Q3,'g')



% Set the plotting limits



  set(gca,'xlim',[-j,j])

  set(gca,'ylim',[-j,j])

  set(gca,'zlim',[-j,j])



% Draw a partial 3-D box in the background



  gx1= [1.3*j,-j]; gx2=[-j,-j]; gx3=[-j,-j];

  gy1= [-j,-j]; gy2=[-j,j]; gy3=[-j,-j];

  gz1= [-j,-j]; gz2=[-j,-j]; gz3=[-j,j];

  gu1= [1.3*j,-j]; gu2=[-j,-j]; gu3=[j,j];

  gs1= [-j,-j]; gs2=[-j,j]; gs3=[j,j];

  gt1= [-j,-j]; gt2=[j,j]; gt3=[-j,j];

  plot3(gx1,gx2,gx3,'w:')

  plot3(gy1,gy2,gy3,'w:')

  plot3(gz1,gz2,gz3,'w:')

  plot3(gu1,gu2,gu3,'w:')

  plot3(gs1,gs2,gs3,'w:')

  plot3(gt1,gt2,gt3,'w:')



% Draw the orbit and the true epoch



  plot3(ri,rj,rz)

  plot3(ria,rja,rza,'wo')



% Set some viewing parameters



  axis('equal')

  set(gca,'AspectRatio',[1,.8])

  view(-30,30)

  set(gca,'xdir','reverse')

  set(gca,'ydir','reverse')



% Title the axes



  set(gca,'xlabel',text(0,0,'I'))

  set(gca,'ylabel',text(0,0,'J'))

  set(gca,'zlabel',text(0,0,'K'))

 

% Option to draw the IJK coordinate axes on the sphere



  if qu == 'y'

	gi1= [0,1.7]; gi2=[0,0]; gi3=[0,0];

	gj1= [0,0]; gj2=[0,1.7]; gj3=[0,0];

	gk1= [0,0]; gk2=[0,0]; gk3=[0,1.7];

	plot3(gi1,gi2,gi3,'w')

	plot3(gj1,gj2,gj3,'w')

	plot3(gk1,gk2,gk3,'w')

  end



% Optional command that can be used to animate the orbit path

%

% comet3(ri,rj,rz)