function Control1_3_JMS(version,units)
% Solution to control1, 2023/11/03
% J.M.Soler

% Set data
if version=='A'
  T = 1.4;   % orbital period in years (yr)
  dp = 0.6;  % periastron in astronomical units (au)
  da = 1.8;  % afoastron in au
elseif version=='B'
  T = 1.9;
  dp = 0.3;
  da = 2.0;
elseif version=='C'
  T = 0.9;
  dp = 0.8;
  da = 1.5;
elseif version=='D'
  T = 2.6;
  dp = 0.2;
  da = 2.8;
else
  error('Unknown version')
end
nt = 1000; % number of orbital points

% Set units
yr = 3.15576e7;       % one year in s
au = 1.495978707e11;  % one astronomical unit in m
if nargin<2
  units = 'au';
end
if strcmp(units,'au')
  unitLength = 1;
  unitTime = 1;
elseif strcmp(units,'MKS')
  unitLength = 1/au;
  unitTime = 1/yr;
else
  error('Unknown units')
end
unitEnergy = (unitLength/unitTime)^2;         % unit of energy/mass

% Find orbital constants
[a,b,c,e,h] = ParametrosElipse_JMS(dp,da);
GM = (2*pi)^2*a^3/T^2; % gravitational constant times atracting mass (au^3/yr^2)

% Find orbital trajectory
theta = linspace(0,4*pi,nt);
[AE,t,r,R] = OrbitaEliptica_JMS(h,e,T,theta);

% Plot orbital trajectory
r1 = dp*(1+e/2);
r2 = dp*(1+e);
r3 = dp/(1-e);
ir = find( dp<=r & r<=r1 ); % red points
iy = find( r1<=r & r<r2 );  % yellow points
ig = find( r2<=r & r<r3 );  % green points
ib = find( r3<=r & r<=da ); % blue points

figure(11)
plot( R(1,ir)/unitLength,R(2,ir)/unitLength,'r.', ...
      R(1,iy)/unitLength,R(2,iy)/unitLength,'y.', ...
      R(1,ig)/unitLength,R(2,ig)/unitLength,'g.', ...
      R(1,ib)/unitLength,R(2,ib)/unitLength,'b.', 'MarkerSize',5 )
axis equal
if strcmp(units,'au')
  xlabel('x (a.u.)')
  ylabel('y (a.u.)')
else
  xlabel('x (m)')
  ylabel('y (m)')
end
set(gca,'FontSize',15)
grid on

% Find and plot orbit's area
figure(12)
A0 = pi*a*b;
A = 0.5*integral(theta,r.^2);
plot(t/T,A/A0,'LineWidth',1.5)
xlabel('t/T')
ylabel('A/A_0')
set(gca,'FontSize',15)
grid on

% Find and plot orbital velocity
DthetaDt = derivative(t,theta);
DrDt = derivative(t,r);
v = sqrt( DrDt.^2 + (r.*DthetaDt).^2 );

figure(13)
plot(t/unitTime,DthetaDt*unitTime,'LineWidth',1.5)
if strcmp(units,'au')
  xlabel('time (yr)')
  ylabel('d\theta/dt (rad/yr)')
else
  xlabel('time (s)')
  ylabel('d\theta/dt (rad/s)')
end
set(gca,'FontSize',15)
grid on

figure(14)
plot( t/unitTime,DrDt/(unitLength/unitTime), ...
      t/unitTime,v/(unitLength/unitTime), 'LineWidth',1.5 )
legend('radial velocity','total velocity')
if strcmp(units,'au')
  xlabel('time (yr)')
  ylabel('dr/dt (au/yr)')
else
  xlabel('time (s)')
  ylabel('dr/dt (m/s)')
end
set(gca,'FontSize',15)
grid on

% Find and plot orbital energies
ek = 0.5*v.^2;  % kinetic energy per unit mass of orbiting body
ep = -GM./r;
et = ek + ep;

figure(15)
plot( t/unitTime,ek/unitEnergy, ...
      t/unitTime,ep/unitEnergy, ...
      t/unitTime,et/unitEnergy, 'LineWidth',1.5 )
legend('E_k/M','E_p/M', 'E_t/M')
if strcmp(units,'au')
  xlabel('time (yr)')
  ylabel('Energy per unit mass (au/yr)^2')
else
  xlabel('time (s)')
  ylabel('Energy per unit mass (m/s)^2')
end
set(gca,'FontSize',15)
grid on

% Find and plot orbital length
l = integral(t,v);
l = l(end)/2;        % divided by 2 because we had two periods
K(1) = 1;              
ln(1) = 2*pi*a*K(1);
for i = 2:11
  K(i) = K(i-1) * (2*i-3)*(2*i-5)/(2*i-2)^2 * e^2;
  ln(i) = ln(i-1) + 2*pi*a*K(i);
end

figure(16)
n = 0:10;
plot( n,ln/unitLength,'o', [0,10],[l,l]/unitLength, ...
      'MarkerSize',10, 'LineWidth',1.5 )
xlabel('N')
if strcmp(units,'au')
  ylabel('Orbital length (au)')
else
  ylabel('Orbital length (m)')
end
set(gca,'FontSize',15)
grid on

end

