% demo_elp.m     December 8, 2012

% graphics display of the orbital elements
% of the Moon using function elp2000.m

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

% read leap seconds data file

readleap;

% unit conversion constants

dtr = pi / 180;

rtd = 180 / pi;

% begin simulation

clc; home;
   
fprintf('\nprogram demo_elp\n');
   
fprintf('\n< orbital elements of the moon >\n');

% request initial calendar date

fprintf('\nplease input a UTC calendar date\n');

[month, day, year] = getdate;

itry = 5;

for i = 1:1:itry
    
   fprintf('\nplease input the simulation period (days)\n');

   ndays = input('? ');
   
   if (ndays > 0)
       
      break;
      
   end
   
end

for i = 1:1:itry
    
   fprintf('\nplease input the graphics step size (days)\n');

   dtstep = input('? ');
   
   if (dtstep > 0)
       
      break;
      
   end 
   
end

for i = 1:1:itry
    
   fprintf('\nplease select the orbital element to plot\n');
   
   fprintf('\n  <1> semimajor axis');

   fprintf('\n  <2> eccentricity');

   fprintf('\n  <3> orbital inclination');

   fprintf('\n  <4> argument of perigee');

   fprintf('\n  <5> right ascension of the ascending node');

   fprintf('\n  <6> true anomaly');
   
   fprintf('\n  <7> apogee radius');
   
   fprintf('\n  <8> perigee radius');
   
   fprintf('\n  <9> geocentric distance\n\n');

   oetype = input('? ');
   
   if (oetype >= 1 && oetype <= 9)
       
      break;
      
   end
   
end

% compute initial utc julian date

jdutc = julian(month, day, year);

% compute tdt julian date

jdatei = utc2tdt(jdutc);

% create data

t = -dtstep;

npts = 0;

while (1)
    
   t = t + dtstep;
   
   tjd = jdatei + t;
   
   npts = npts + 1;
      
   % compute orbital elements
   
   oev = elp2000(tjd);
   
   % determine serial date number
   
   [month, day, year] = gdate(tjd);

   % serial date number

   x(npts) = datenum(year, month, day);

   switch oetype
       
   case 1
       
      % semimajor axis
      
      y(npts) = oev(1);
      
   case 2
       
      % eccentricity
      
      y(npts) = oev(2);
      
   case 3
       
      % inclination
      
      y(npts) = oev(3) * rtd;
      
   case 4
       
      % argument of perigee
      
      y(npts) = oev(4) * rtd;
      
   case 5
       
      % raan
      
      y(npts) = oev(5) * rtd;
      
   case 6
       
      % true anomaly
      
      y(npts) = oev(6) * rtd;
      
   case 7
       
      % apogee radius
      
      y(npts) = oev(1) * (1 + oev(2));
      
   case 8
       
      % perigee radius
      
      y(npts) = oev(1) * (1 - oev(2));
      
   case 9
       
      % geocentric distance
      
      y(npts) = oev(1) * (1 - oev(2)) / (1 + oev(2) * cos(oev(6))); 
      
   end
   
   % check for end of simulation
   
   if (t >= ndays)
       
      break;
      
   end 
   
end   

% create and label plot
   
plot(x, y, 'LineWidth', 2);

% label x axis with calendar dates

datetick('x', 2);

switch oetype
    
case 1
    
   title ('Semimajor Axis of the Moon', 'FontSize', 16);
   
   ylabel('Semimajor Axis (kilometers)', 'FontSize', 12);
   
case 2
    
   title ('Orbital Eccentricity of the Moon', 'FontSize', 16);
   
   ylabel('Orbital Eccentricity', 'FontSize', 12);
   
case 3
    
   title ('Orbital Inclination of the Moon', 'FontSize', 16);
   
   ylabel('Orbital Inclination (degrees)', 'FontSize', 12);
   
case 4
    
   title ('Argument of Perigee of the Moon', 'FontSize', 16);
   
   ylabel('Argument of Perigee (degrees)', 'FontSize', 12);
   
case 5
    
   title ('RAAN of the Moon', 'FontSize', 16);
   
   ylabel('RAAN (degrees)', 'FontSize', 12);
   
case 6
    
   title ('True Anomaly of the Moon', 'FontSize', 16); 
   
   ylabel('True Anomaly (degrees)', 'FontSize', 12');
   
case 7
    
   title ('Apogee Radius of the Moon', 'FontSize', 16);
   
   ylabel('Apogee Radius (kilometers)', 'FontSize', 12');
   
case 8
    
   title ('Perigee Radius of the Moon', 'FontSize', 16);
   
   ylabel('Perigee Radius (kilometers)', 'FontSize', 12');
   
case 9
    
   title ('Geocentric Distance of the Moon', 'FontSize', 16);
   
   ylabel('Geocentric Distance (kilometers)', 'FontSize', 12'); 
   
end   
   
xlabel('Calendar Date', 'FontSize', 12);
 
% turn on grid and zoom feature

grid;
   
zoom on;

% create eps graphics file with tiff preview

print -depsc -tiff -r300 demo_elp.eps