%
%   TRAJECT.M   Plotting and Function Evaluation with MATLAB    
%
%   This demo gives a simple example of the use of plots in MATLAB.      
%   It also shows how to use MATLAB to perform simple functional evaluations    
%   using vectorized operations. See the MATLAB HELP files for more description.
%                                                                               
%   This is an example of plotting orthogonal trajectories (ie. families of     
%   curves that intersect at right angles).  It uses the cartesian representation
%   of the curves.                                        
%       1.  family of parabolas   y = c*x^2                                     
%       2.  family of ellipses   x^2 + 2*y^2 = k = b^2                                
%
%   File prepared by J. R. White, UMass-Lowell (July 2003)
%

%
%   getting started
      clear all,  close all,   nfig = 0;
%                                                                               
%   define x-axis
      nx = 2000;  x = linspace(-5,5,nx);
%
%   create curves for original family  --  parabolas
      cc = [0.0 0.1 0.3 0.8 1.5];  ncc = length(cc);   % set constants
      y1 = zeros(nx,ncc);  y2 = zeros(nx,ncc);         % initialize function values 
      for n = 1:ncc
        c = cc(n);
        yy = c*x.^2;   y1(:,n) = yy';
        y2(:,n) = -y1(:,n);    % for negative constants
      end
%
%   create curves for orthogonal family  --  ellipses
      kk = [5 10 20 30];  nkk = length(kk);            % set constants
      z1 = zeros(nx,nkk);  z2 = zeros(nx,nkk);         % initialize function values 
      for n = 1:nkk
        k = kk(n);   zz = real(sqrt((k - x.^2)/2)); 
        for i = 1:nx 
          if zz(i) > 0,  z1(i,n) = zz(i);  else  z1(i,n) = NaN;   end
        end
        z2(:,n) = -z1(:,n);    % for negative roots
      end
%
%   Note:  If the argument of the 'sqrt' function is negative, the result will be a 
%   pure imaginary number.  The 'real' command is used in the above equations to extract 
%   only the real part.  If the combination of the constant k and the value of x is such
%   that the argument is negative, then z(x) simply takes on a value of NaN (not-a-number)
%   -- and the Matlab plot command simply ignores these.  Thus, we can be a little sloppy
%   here without impacting the desired plots.  Thus, in the above development, we do not 
%   need to adjust the range of x for each new coefficient. 
%
%   plot both families  (full x-y range)
      nfig = nfig+1;   figure(nfig)
      plot(x,y1,'LineWidth',2), hold on
      plot(x,y2,'LineWidth',2),  plot(x,z1,'LineWidth',2),  plot(x,z2,'LineWidth',2)
      axis('equal');  v1 = [-5 5 -5 5]; axis(v1);  grid;
      title('Traject:  Example of Orthogonal Trajectories')
      xlabel('x-axis'),ylabel('y-axis'),
%
%   plot both families  (first quadrant only)
      nfig = nfig+1;   figure(nfig)
      plot(x,y1,'LineWidth',2), hold on
      plot(x,z1,'LineWidth',2)
      axis('equal');  v2 = [0 5 0 5]; axis(v2);  grid;
      title('Traject:  Example of Orthogonal Trajectories')
      xlabel('x-axis'),ylabel('y-axis'),
%
%   end of problem