Assignment 2 Scripts

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Problem A %%%%%%%%%%%%%%%%%%%%%%
% Script ShowBigDipper
  clear
% Big Dipper Evolution 

% The Big Dipper is a constellation consisting of
% seven stars which we index as follows:
%
%
%                 2 *      3*                  7 *
%       1 *                     4 *
%                               
%                                               6 *
%                                   5 * 

% The location of these stars varies with time. We have data
% for times t(1),...,t(5) where
  
  t  = [-50; -25; 0; 25; 50];  
  
% Think of t as (-50000BC, -25000BC, Now, 25000AD, 50000AD)
  
% The (x,y) coordinate of star j at time t(i) is specified by
% (x(i,j),y(i,j)) where

  x =[ 15.5   29.0   36.2   45.2   49.5   62.0  55.2;...
       16.0   28.2   34.9   44.6   49.2   61.5  56.5;...
       17.8   27.9   34.3   44.1   48.6   61.0  58.9;...
       18.5   27.5   33.9   43.8   48.5   60.5  61.2;...
       20.3   26.2   33.2   42.9   47.8   59.7  62.5];

  y = [ 0.0   5.0    5.1    5.5   0.0   6.7   14.8 ;...
       -0.4   5.3    4.3    5.2   0.0   4.6   14.3 ;...
        0.5   5.5    4.3    5.2   0.0   4.1   13.1 ;...
        0.4   5.4    4.3    5.7   0.0   4.4   12.7 ;...
       -1.2   5.5    4.4    6.2   0.0   5.0   13.2];

% Here is a 5-frame "movie" of the evolution...

   bigscreen          % An NCM function that makes the figure window big
   starsize = 18;     % Might want to play with this value.
   fps = 1;           % frames per second
   close all
   for k=1:5
      hold on
      fill([0  80  80  0  0],[-30 -30  50  50 -30],'k')
      plot(x(k,:),y(k,:),'r')
      plot(x(k,:),y(k,:),'*y','markersize',starsize)
      text(30,40,sprintf('Time = %6.0f',1000*t(k)),'color','y','fontsize',24)
      hold off
      axis([0 80 -30 50])
      axis equal off
      pause(2)                % Take a moment to look at the frame
      M(k) = getframe;
   end
   close all
   movie(M,3,1)               % Play the movie


% A2A
% Shows Big Dipper Evolution movie
  
  bigscreen
  fps = 1;  % frames per second
  r = 3;    % movie shown this number of times
  M = BigDipperMovie(linspace(-50,100,16)');
  shg
  movie(M,r,fps)
  
  
%%%%%%%%%%%%%%%%%%%%%  Problem B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   function ShowMultipleInterpN(n,m)
 % Illustrates MultipleInterpN
 % n and m are integers that satisfy 1 < = m < = n
 %
 % m-point polynomial interpolants of sin(x) across
 % [0,2pi] are displayed. n shouldn't be too much bigger
 % than 20. 
 
 x = linspace(0,2*pi,n)';
 y = sin(x);
 C = MultipleInterpN(x,y,m);
 x0 = linspace(0,2*pi)';
 y0 = sin(x0);
 [p,q] = size(C);
 for k=1:q
     xsub = x(k:k+m-1);
     pVals = HornerN(C(:,k),xsub,x0);
     subplot(2,1,1)
     plot(x0,y0,x0,pVals,xsub,sin(xsub),'o',x,sin(x),'*k')
     title(['Interpolation point indices = ' sprintf('%1d  ',k:k+m-1)])
     axis([0 2*pi -3 3])
     subplot(2,1,2)
     semilogy(x0,abs(pVals-y0)+eps)
     axis([0 2*pi 10^-6 10^1])
     title('Approximate Absolute Error')
     shg
     pause
 end
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% problem C %%%%%%%%%%%%%%%%%%%%

% Script A2C

% Check for correctness..

A = linspace(0,360,13)';
D = [ 408 89 -66 10 338 807 1238 1511 1583 1462 1183 804 408]';

Avals = linspace(0,360,200)'; 
F = newCSInterp(D(1:12));
Fvals = CSeval(F,360,Avals);
plot(Avals,Fvals,A,D,'o')
axis([-10 370 -200 1700])
set(gca,'xTick',linspace(0,360,13))
xlabel('Ascension (Degrees)')
ylabel('Declination (minutes)')

% Check efficiency

   nRepeat = 1;     % Number of instances that are timed         
% Make big enough so clock "granularity"  is not an issue...
 
   clc
   disp('  n     CSInterp Time/ newCSInterp time')
   disp('----------------------------------------')
   for n = 100:50:500
       f = randn(n,1);
       % Benchmark the original implementation on page SCMV p. 101
       tic 
       for k=1:nRepeat
           F = CSInterp(f);
       end
       t1 = toc;
       tic 
       for k=1:nRepeat
           F = newCSinterp(f);
       end
       t2 = toc;
       disp(sprintf('  %3d            %6.3f',n,t1/t2))
   end