% Math 32a   Fall 2003  Richard Palais  Brandeis Univ.

%              Topics for Today
%
% 1)  Subscripting or Indexing
% 2)  Vectorizing code
% 3)  Functions 
% 4)  Relational Operators and Boolean Expressions 
% 5) Flow Control
%   a)   if  else
%   b)   for loops
%   c)   while loops
%   d)   switch statement
% 6)  2D and 3D Graphing: the plot, plot3, mesh, and surf commands
% 7)   M-Files
%   a) Script M-Files 
%   b) Function M-Files
%       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  To practice with Matlab on your own, choose
%  a topic and enter the lines, one by one, into
%  the Matlab Command Window, and make sure you
%  understand all the output you get back before
%  going on to the next topic.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Subscripting
%
a = [0.0: 0.1: 1]
a(1)
a(5)
a(end)
b = [0:5 ; 1:6; 2:7 ; 3:8 ; 4:9]
b(2:4, 4)
b(3, 2:end)
b(4,4) = 0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Stopwatch Functions
%
%  If you execute the command tic, Matlab starts 
% a computation time stopwatch. The next time you
% execute the command toc, Matlab reports how much
% computing time in seconds it has used.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%            Vectorizing Code
%
%  A matlab command is said to be "vectorized" if it can act
%  not only on scalars, but also on vectors and matrices 
%  (Often this just means it does to each element of an array 
%  what it would do to a scalar). But as the next two examples
%  show that's not always the case.
%
% Calcuate Factorial 150 two ways, with a loop and using vectorized product.
tic, factorial = 1; for i = 2:150, factorial = i*factorial; end, toc 
tic, factorial = prod(1:150); toc 
%
% A similar pair of computations for summing an arithmetic progression.
tic, SUM = 0; for i = 1:10000, SUM = i + SUM; end, toc 
tic, SUM = sum(1:10000); toc 
% Notice that most built-in Matlab functions are vectorized.
X = [1:10]
sin(X)
X^2
% Why the error message?
% The problem here is that for a matrix,  Matlab uses ^2 to mean 
% matrix squaring, which only makes sense for sauare matrices.
% For example:
X = [1,2;3,4]
% Vectorized squaring (i.e., squaring each element of an array
% is denoted by .^2.
X.^2
%  Similarly to multiply corresponding elements of two matrices 
% with the same number of rows and columns, use .*
Y = [4,3;2,1]
X.*Y
% Similarly for vecorized division, use ./
X./Y
% X/Y (without the dot) means matrix product of X with Y inverse.
X/Y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%            Inline  Functions 
%
%  Suppose we want to create the polynomial function 
%  f(x) = 3*x^2 + 5*x + 2
%  Create a string S that defines the polynomial and 
%  assign inline(S) to the name of the function:
f = inline('3*x^2 + 5*x + 2')
f(1)
f(2)
v = [0:5]
f(v)
%  We got an error because we didn't vectorize properly.
f = inline('3*x.^2 + 5*x + 2')
x = [-3: 0.1 : 1.5];
plot(x,f(x));
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Relational Operators and Boolean Expressions 
%
%  The basic relational operators are
%  > ,  < , >= , <=  with the obvious meanings, and 
%  == which means "equals" and ~= which means "not equal"
% If R is any one of these operators and x and y are 
% constants or variables then xRy is a boolean 
% expession, i.e., it eiher has the value 1 meaning 
% xRy is true, or it has the value o meaning xRy is false.
% DON'T CONFUSE (A == B) with A = B !! The latter assigns the 
% value of B to A, while (A == B) does not change the value
% of A and rather is an expression that has the value 1 if
% A is equal to B and otherwise has the value zero. 
% You can build up more complex boolean expressions by
% combining basic one with & (which means AND) and | (which
% means OR) and ~ (which as we have already seen means NOT).
% First check that the relational operators are vectorized:
A = 0:8 
B = 8-A
bool1 = A>4,
bool2 = (A==B)  
bool3 = bool1 & bool2
bool3 = bool1 | bool2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             Flow Control
%
%       if, else, and elseif
%
%  The general form of an if statement is:
%   
%   if booleanexpression 
%       command 1
%       command 2
%       ...
%   else
%       COMMAND1
%       COMMAND2
%        ...
%    end      
%
%  If booleanexpression evaluates to 1(true) the first 
%  group of commands is executed, otherwise the second
%  set of commands is executed. For example here is a
%  way to compute the maximum of a and b. To create a 
%  multi-path branch with more conditions, replace all 
%  but the final else by elseif.
 
max = 0;
a = 5;
b = 3;
if (a>b) 
  max = a;
else 
  max = b;
end
max

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       for loops
%
%  We have already seen several examples.
%  However ther is more than meets the eye.
%  The genral form is:
%
%  for x = Array
%    command 1
%    command 2
%    ...
%    end
%
%  First x is set equal to the first column of 
%  Array and all the commands executed, then x is 
%  set equal to the second column of Array and
%  the commands executed again and so on until
%  the last column of Array is reached.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       while loops
%
%  Here the general form is:
%
%   while booleanexpression
%    command 1
%    command 2
%    ...
%    end 
%
%  First, booleanexpression is evaluated, and if
%  it evaluates to 1 (true) the commands are 
%  executed, then booleanexpression is re-evaluated
%  and if it is true then the commands are evaluated
%  again, and so on until booleanexpression evaluates 
%  to 0 (false) at which point contol shifts to 
%  the first command after the end statement.
%  Of course something better happen durring execution 
%  of the commands that will eventually make 
%  booleanexpression false or we wre in an infinite loop.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        switch statement
%
%  The switch statement allows a program to switch 
%  its execution path to any one of a number of possible
%  different branches depending on the value of an expression
%  (the switch-expression).
%  It has the general form:
%
% switch <switch-expression>
%   case <case-expression> statement,statement,statement, ...;
%   case <case-expression> statement,statement,statement, ...;
%   ...
%  otherwise statement,statement,statement, ...;
%  end
%
%  First <case-expression> is evaluated. It should evaluate to
%  either a scalar or a string. Control then jumps to the 
%  statements following the first case-expression that 
%  matches the value of switch-expression, or if there is no 
%  match to the statements following otherwise (or if there is 
%  no otherwise statement, to the statement following end.
%  Example:
%
CurveType = 'Cubic';
switch CurveType
  case 'Linear'
      f = inline('a + b*x');
  case 'Parabola'
     f = inline('a + b*x + c*x^2');
  case 'Cubic'
     f = inline('a + b*x + c*x^2 + d* x^3');
  end
  f
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             Graphing
%   2D curves
%   A circle centerd at (a,b) with radius r is given parametrically by 
%   x = a + r*cos(t)  y =b +  r * sin(t)  with t in [0, 2*pi]
a = 1;
b = 2; 
r = 2;
t = [0: 0.1: 2*pi];
x = a + r* cos(t);
y = b + r* sin(t);
plot(x,y); axis equal
%%%%%%%%%%%%%%%
%  3D  Curves
%  Helix  x = cos(t); y = sin(t);   z = 0.2*t;
t = [0: 0.1: 6*pi];
x = cos(t);
y = sin(t);   
z = 0.2*t;
plot3(x,y,z);
rotate3d;      %This permits you to rotate the 3D object with the mouse.
%%%%%%%%%%%%%%%
%       Parametric Surfaces
%
%  In studying so-called polar spherical coordinates, you
%  probably learned that the point on the unit sphere
%  with longitude u and co-latitude v has the (x,y,z)
%  coordinates x = cos(u), y = sin(v), and z = cos(v), 
%  Here is the way to render the sphere in Matlab.
%  First, this is the command to create a two-dimensional grid of 
%  points (u,v) that will represent the longitude and co-latitude 
%  of points of the sphere. Let's divide both intervals [0, 2*pi]
%  and [0,pi] into 50 sub-intervals:
[u,v] = meshgrid([0 : 2*pi/50 : 2*pi],[0 : pi/50 : pi]);
%  Then we create the arrays of values of the components of the 
%  3D points (x,y,z) that have longitude u and co-latitude v:
x = cos(u).*sin(v);
y = sin(u).*sin(v);
z = cos(v);
%  Finally, the command mesh(x,y,z) takes the 3D grid that we have 
%  created and maps it into 3-space using "wireframe rendering".
%  This just means that each 2D gridlines is mapped to the 
%  3D segment joining the images of its endpoints.
mesh(x,y,z), axis equal
%  On the other hand, surf(x,y,z) renders the grid in "patch mode".
%  This means that each rectangle of the grid is filled in with a 
%  color that is determined by an algorithm we will discuss later.
surf(x,y,z), axis equal
%
% Exercise: Render the ellipsoid with semi-major axes a,b, and c.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Script M-Files
%
%  The first kind of M-File, called a Script File, is easy to 
%  explain. It is just a list of Matlab commands written 
%  just as you would write them into the Matlab command window. 
%  When you tyoe the name of the Script File (without the .m) 
%  into the command window, the effect is exactly the same as if 
%  you typed all the commands in the script File, one after another,
%  into the command window at the place. End of story. 
%  (Well, almost. you do have to make sure that the script file is in 
%  a directorythat is in matlab's search path in order for Matlab to 
%  be able to find it.)
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        Function M-Files
%
% The second kind of M-file, called a function M-File is a lot more
% sophisticated and more flexible, and it is the real key to what it
% takes to write powerful Matlab programs. They are much like subroutines
% in other programming languages that you may be familiar with. We will
% probably have a lot more to say about them in the weeks to come. 
% For an example, let's use the above switch statement to create a 
% function file called PlotCurve.m that will plot a line, parabola,  
% or cubic, depending on a string parameter the user enters.

 function PlotCurve(CurveType)

%  This function plots either a line, a parabola or a cubic curve depending 
%  on whether it is called with parameter 'Linear', 'Parabola', or 'Cubic'.

t = [-1.5: 0.05: 1.5];
switch CurveType
  case 'Linear'
      plot(t, 1 + 2*t); axis equal
  case 'Parabola'
      plot(t, 1 - 0.5*t + 2*t.^2); axis equal
  case 'Cubic'
      plot(t, 1 - 1.5*t + 2*t.^2 + 2.5* t.^3); axis equal
  end
