P5: Given Scripts/Functions
Circle Pointsfunction [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 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];