A5 Scripts
A5A
% Script A5A
% Examines ConicFit
bigscreen
% The unperturbed data...
x = [ 1.02 .95 .87 .77 .67 .56 .44 .30 .16 .01]';
y = [ 0.39 .32 .27 .22 .18 .15 .13 .12 .13 .15]';
n = 10;
% We'll contaminate with scaled versions of these noise vectors
rand('seed',0)
xrand = -1+2*rand(n,1);
yrand = -1+2*rand(n,1);
for noise = [0 .0001 .001 .01 .1]
[a,b,c,d,e] = ConicFit(x+noise*xrand,y+noise*yrand);
xmin = -1; xmax = 2 ; deltax = .02;
ymin = -1; ymax = 1.5 ; deltay = .02;
[X,Y] = meshgrid(xmin:deltax:xmax,ymin:deltay:ymax);
Z = a*X.^2 + b*X.*Y + c*Y.^2 + d*X + e*Y + 1;
contour(X,Y,Z,[0 0],'k')
hold on
plot(x,y,'*')
hold off
title(sprintf('a = %9.6f b = %9.6f c = %9.6f d = %9.6f e = %9.6f',a,b,c,d,e))
shg
pause
end
%%%%%%%%%% Sensitivity study %%%%%%%
% Examine how small changes in the data affect the outcome.
% Here are two data sets.
% One is best approximated with an ellipse, the other by a hyperbola.
xEllipse = [1.0144 0.9409 0.8736 0.7736 0.6787 0.5577 0.4404 0.3066 0.1507 0.0011]';
yEllipse = [0.3906 0.3234 0.2602 0.2177 0.1713 0.1483 0.1337 0.1218 0.1386 0.1569]';
xHyperb = [0.9638 0.8594 0.9058 0.8059 0.7569 0.5367 0.4439 0.3662 0.0669 -0.0793]';
yHyperb = [0.3959 0.3542 0.1715 0.1967 0.0934 0.1335 0.1674 0.1378 0.2161 0.2192]';
% Here is a measure of the separation of the two data sets
sep = norm(xEllipse-xHyperb,1)+norm(yEllipse-yHyperb,1);
while sep >.00001
% Look at midpoint of the datasets and approximate with a conic
xMid = (xEllipse+xHyperb)/2;
yMid = (yEllipse+yHyperb)/2;
[a,b,c,d,e] = ConicFit(xMid,yMid);
if c < 0
% A hyperbola renders the best fit at the midpoint
xHyperb = xMid; yHyperb = yMid;
else
% An ellipse renders the best fit at the midpoint
xEllipse = xMid; yEllipse = yMid;
end
% Plot the "endpoint" conics
[a,b,c,d,e] = ConicFit(xEllipse,yEllipse);
[X,Y] = meshgrid(-1:.02:2,-1:.02:1.5);
Z = a*X.^2 + b*X.*Y + c*Y.^2 + d*X + e*Y + 1;
contour(X,Y,Z,[0 0],'k')
hold on
text(-.9,-.6,sprintf('Ellipse: a = %9.6f b = %9.6f c = %9.6f d = %9.6f e = %9.6f',a,b,c,d,e))
plot(xEllipse,yEllipse,'*k')
[a,b,c,d,e] = ConicFit(xHyperb,yHyperb);
text(-.9,-.7,sprintf('Hyperbola: a = %9.6f b = %9.6f c = %9.6f d = %9.6f e = %9.6f',a,b,c,d,e))
Z = a*X.^2 + b*X.*Y + c*Y.^2 + d*X + e*Y + 1;
contour(X,Y,Z,[0 0],'r')
plot(xHyperb,yHyperb,'or')
hold off
sep = norm(xEllipse-xHyperb,1)+norm(yEllipse-yHyperb,1);
title(sprintf('Ellipse is black, Hyperbola is red, Data Separation = %5.2e',sep))
pause
end
A5C
% A5C % Examines Update close all % Using fmin which is soon to disappear ... warning off matlab:fmin:ObsoleteFunction % Season length data... sA = [92.7556 93.6507 89.8424 88.9938]; % Initial guess.. d = 0; alpha = 0; kMax = 6; for k = 1:kMax [dnext,alphanext] = Update(d,alpha); d = dnext; alpha = alphanext; sE = SeasonLength(d,alpha); DrawEccentric(d,alpha,sE,sA) pause end
CirclePoints
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];
DrawEccentric
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
Pretty
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];