P5: Given Scripts/Functions

  function [x,y] = CirclePoints(d,alpha)
% Determines the equinox and solstice points for the eccentric model 
% with displacement d and rotation alpha. 
%
% x and y are row 4-vectors with the property that
%
%  (x(1),y(1)) = coordinate of the Vernal   Equinox  point
%  (x(2),y(2)) = coordinate of the Summer   Solstice point
%  (x(3),y(3)) = coordinate of the Autumnal Equinox  point
%  (x(4),y(4)) = coordinate of the Winter   Solstice point

c = cos(alpha); s = sin(alpha);
tau = sqrt(1-(s*d)^2);
xV = d*s^2 + c*tau; yV = s*(-d*c + tau);
xA = d*s^2 - c*tau; yA = s*(-d*c - tau);

c = cos(alpha+pi/2); s = sin(alpha+pi/2);
tau = sqrt(1-(s*d)^2);
xS = d*s^2 + c*tau; yS = s*(-d*c + tau);
xW = d*s^2 - c*tau; yW = s*(-d*c - tau);

x = [xV xS xA xW];
y = [yV yS yA yW];

  function DrawEccentric(d,alpha,sE,sA)
% Displays the eccentric model with displacement d and rotation
% alpha. 
% sE is a row 4-vector of estimated season lengths (in days).
% sA is a row 4-vector of actual season lengths (in days).
% Both sE and sA specify these lengths in Spring-Summer-Fall-Winter order.

% Plot Circle and coordinate axes
theta = linspace(0,2*pi);
plot(cos(theta),sin(theta),'k')
axis([-1.2 1.2 -1.2 1.2])
axis equal off
hold on
plot(linspace(-1.1,1.1),zeros(1,100),':k')
plot(zeros(1,100),linspace(-1.1,1.1),':k')
hold on

[x,y] = CirclePoints(d,alpha);

% Plot the equinox connector

xV = x(1); yV = y(1);
xA = x(3); yA = y(3);
plot([xV xA],[yV yA])
text(1.1*xV,1.1*yV,'V')
text(1.1*xA,1.1*yA,'A')
plot(xV,yV,'m*',xA,yA,'m*')

% Plot the solstice connector

xS = x(2); yS = y(2);
xW = x(4); yW = y(4);
plot([xS xW],[yS yW])
text(1.1*xS,1.1*yS,'S')
text(1.1*xW,1.1*yW,'W')
plot(xS,yS,'m*',xW,yW,'m*')

% Display the estimated and actual season lengths
text( 1.0,1.15,'\bfSpring'); text( 1.0, 1.0,pretty(sE(1))); text( 1.0,  .85,pretty(sA(1)))
text(-1.5,1.15,'\bfSummer'); text(-1.5, 1.0,pretty(sE(2))); text(-1.5,  .85,pretty(sA(2)))
text(-1.5,-.85,'\bfAutumn'); text(-1.5,-1.0,pretty(sE(3))); text(-1.5,-1.15,pretty(sA(3)))
text( 1.0,-.85,'\bfWinter'); text( 1.0,-1.0,pretty(sE(4))); text( 1.0,-1.15,pretty(sA(4)))

% Display the error and the parameter values
title( ['Error = ' pretty(max(abs(sE-sA)))])
text(-.8,-1.3,sprintf('d = %10.6f     alpha = %10.6f',d,alpha))

hold off

  function s = Pretty(d)
%
% s = Pretty(d) converts a decimal number assumed to specify time in days
% into a string that indicates the time in days, hours, and minutes.
%
% Input:
%   
%   d    A nonnegative real number (Earth Days)
%
% Output:
%
%   s    s string that reports d in days, hours, and minutes.
%        For example,  if 
%
%            days = 92 + 34*(1/24) + 21*(1/(24*60)), 
%
%        then s would be the string '92d 34m 21s'.
% 
%
% -------------------------------------------------------------------------

%  Convert days to numerical d-h-m form.

  days = floor(d);   
  hrs  = (d-days)*24;
  h    = floor(hrs);
  m    = round((hrs-h)*60);

% Adjust so that 0<=m<60 and 0<=h<24

  if m==60
     m = 0;
     h = h+1;
  end
  if h==24
     h = 0;
     days = days+1;
  end

% Convert d, h, and m to strings.

  sd = [num2str(days) 'd '];
  if days<=9
     sd = ['0' sd];
  end 
  if h<=9
     sh = ['0' num2str(h) 'h '];
  else
     sh = [num2str(h) 'h '];
  end
  if m<=9
     sm = ['0' num2str(m) 'm'];
  else
     sm = [num2str(m) 'm'];
  end

% Concatenate

  s = [sd sh sm];