%MATLAB SCRIPT FILE s322eg6.m  */
%
%    EXAMPLE 6 FOR STATISTICS 322
%    POSTED 9/12/96
%    EXPLORES GRAPHICS

ipar = 1 ;     %  1,2,3,4,5,6,7,8

format compact ;
disp([10 'MATLAB script s322eg6.m, with ipar = ' num2str(ipar) 10]) ;
          %  disp IS A USEFUL COMMAND, THAT I JUST DISCOVERED (BUT 
          %  WISH I HAD KNOWN ABOUT LONG AGO).  IT "DISPLAYS" THE 
          %  STRING WHICH IS ITS ARGUMENT (WITHOUT WRITING THE 
          %  "VARIABLE NAME" AS I DID PREVIOUSLY).


  %  MATLAB GRAPHICS CAN BE USED AT TWO LEVELS.
  %
  %  THE "HIGH LEVEL" IS INTENDED FOR VERY QUICK AND EASY USE.
  %  
  %  THE "LOW LEVEL" IS INTENDED TO BE VERY FLEXIBLE, BUT OFTEN TAKES 
  %  MUCH MORE EFFORT TO USE.
  %  
  %  YOU SHOULD START BY WORKING WITH THE HIGH LEVEL GRAPHICS, BUT BE
  %  AWARE THAT EVENTUALLY YOU WILL WANT TO DO THINGS THAT ARE JUST NOT
  %  AVAILABLE THERE.  WHEN THAT HAPPENS, KNOW THAT THERE IS MUCH MORE
  %  POSSIBLE (AT A COST OF MORE EFFORT) FROM THE LOW LEVEL GRAPHICS.
  %  WE WILL STUDY ONLY HIGH LEVEL GRAPHICS FOR A WHILE.

if ipar == 1 ;    %  THEN DO PART 1 OF "GRAPHICS"

  %  THE BASIC MATLAB COMMAND FOR 2-D PLOTTING IS:
  plot(rand(20,1)) ;
          %  USING A SINGLE VECTOR ARGUMENT SAYS: "USE THESE AS Y'S
          %  AND GENERATE X'S AS 1,2,..."

  %  NOTE THAT MATLAB CREATED A NEW WINDOW CALLED "Figure No. 1".  IT
  %  WILL PUT ALL FURTHER SPECIFIED STUFF IN THAT SAME WINDOW.  YOU CAN 
  %  MANIPULATE IT JUST YOU DO OTHER WINDOWS.  TRY MOVING IT AROUND,
  %  AND RESIZING IT.  IN PARTICULAR A LOCATION WHERE YOU CAN WORK 
  %  IN THE COMMAND WINDOW, WHILE LOOKING AT THE PLOT WINDOW IS 
  %  USEFUL.

  disp('any key to continue') ;
  pause ;
          %  TRY COMMENTING OUT THE PAUSE, AND NOTE THAT YOU NEVER SEE 
          %  THE FIRST PLOT, ONLY THE SECOND.  SO THE pause COMMAND IS
          %  IMPORTANT IS YOU WANT TO SEE SEVERAL PLOTS IN SUCCESSION.

  %  NOW TRY SOMETHING WITH BOTH X'S AND Y'S FED INTO "PLOT"
  %  NOTE THE ORDER IN WHICH THE POINTS ARE CONNECTED BY LINES.  
  plot(rand(20,1),rand(20,1)) ;
          %  I DON'T KNOW ABOUT THE UNIX SYSTEM, BUT ON MY PC I NEED TO
          %  CLICK ON THE FIGURE WINDOW TO SEE THE PLOT CHANGE (SOMETIMES).


elseif ipar == 2 ;    %  THEN DO PART 2 OF "GRAPHICS"

  %  HERE IS A SPIRAL, GENERATED AS A "PARAMETRIC CURVE"
  vt = 0:.02:3 ;
  vx = vt .* cos(2 * pi * vt) ;
  vy = vt .* sin(2 * pi * vt) ;
          %  NOTE THAT Y(X) NEED NOT BE A "FUNCTION", ANY 2 VECTORS
          %  CAN BE PLOTTED

  plot(vx,vy) ;

  isubpar = 0 ;     %  0,1,2,3,4
          %  DO SOME SUBPARTS, TO AVOID REPEATING ABOVE CODE

  if isubpar >= 1 ;
          %  THEN ADD A TITLE TO THE PLOT
    title(['Parametric Spiral, Using ' num2str(length(vx)) ' Points']) ;
          %  CAN PUT VALUES OF VARIABLES, ETC. INTO THE TITLE, BY
          %  PUTTING THEM INTO THE STRINGS.

    if isubpar >= 2 ;
          %  THEN ADD LABELS TO THE X AND Y AXES
      xlabel('Here is the X-axis') ;
      ylabel('and this one is the Y-axis') ;          
          %  VARIABLES CAN BE WORKED INTO THESE STRINGS IF DESIRED, TOO

      if isubpar >= 3 ;
          %  THEN WRITE SOME TEXT TO THE SCREEN
        text(-1,1.5,'We can add text') ;
          %  FIRST TWO VALUES ARE PLOT COORDINATES OF WHERE TO PUT IT

        if isubpar >= 4 ;
              %  THEN ADD LABELS TO EVERY 10TH POINT
          vxtrunc = vx(find(round(5*vt) == 5*vt)) ;
          vytrunc = vy(find(round(5*vt) == 5*vt)) ;
              %  TRUNCATED VERSIONS OF vx AND vy
          text(vxtrunc,vytrunc,'Here is one') ;
              %  SAME STRING GETS WRITTEN FOR EVERY POINT

        end ;
      end ;
    end ;
  end ;

  %  NOTE HOW THE ABOVE if BLOCK IS DIFFERENT FROM USING elseif.

  %  NOTE THAT ALL OF THESE THINGS ARE ADDED after THE GRAPH IS FIRST
  %  PLOTTED.  TRY COPYING THE title COMMAND TO BEFORE THE GRAPH, AND
  %  RUNNING THIS WITH isubpar = 0.  THE REASON THE TITLE DOES NOT SHOW
  %  IS THAT plot CLEARS AWAY EARLIER STUFF (ESSENTIALLY REINITIALIZES
  %  THE GRAPHICS WINDOW).

  %  BUT AFTER plot IS CALLED, YOU CAN DO ALL OF THIS STUFF, YOU CAN
  %  ALSO CHANGE IT BY TYPING NEW COMMANDS ON THE COMMAND LINE (NOT 
  %  JUST IN THE SCRIPT FILE).  TRY THIS BY ISSUING THE COMMAND:
  %          gtext('Here is an important point')
  %  ON THE COMMAND LINE, THEN CLICKING WITH THE MOUSE IN THE 
  %  GRAPHICS WINDOW.  NOTE THAT STRING GETS ADDED TO THE PLOT.
  %  ON MY PC, I (SOMETIMES) HAD TO CLICK ANOTHER WINDOW, THEN CLICK 
  %  BACK TO SEE THE RESULT.

  %  A SIMILARLY USEFUL WAY TO GET MATLAB SCRIPTS TO INTERACT WITH
  %  YOU IN GRAPHICS IS THE COMMAND:   ginput


elseif ipar == 3 ;    %  THEN DO PART 3 OF "GRAPHICS"
    
  %  TO GET MULTIPLE CURVES ON THE SAME PLOT, I.E. AN "OVERLAY",
  %  YOU CAN USE:

  %  1.  MATRIX ARGUMENTS TO plot:
  vx = (-1:.1:1)' ;
          %  RECALL ' MAKES THIS A COLUMN VECTOR
  my = [ones(length(vx),1) vx vx.^2 vx.^3] ;
          %  RECALL .^  VS.  ^  FOR EXPONENTIATION

  isubpar = 1 ;     %  1,2,3

  if isubpar == 1 ;
    plot(vx,my) ;
          %  ARGUMENTS OF PLOT CAN BE EITHER VECTORS OR MATRICES, AS
          %  LONG AS THEY "LINE UP".  HERE I CHOSE TO THINK OF
          %  "PLOTTING COLUMNS".

          %  YOU SHOULD LOOK CAREFULLY FOR THE CONSTANT FUNCTION IN
          %  THIS PLOT.


  %  OR ELSE YOU CAN USE:
  %  2.   MULTIPLE ARGUMENTS
  elseif isubpar == 2 ;
      linex = [0 0; 1 -1]' ;
      liney = [-1 1; 0 0]' ;
    plot(vx,my,'c',linex,liney,'m') ;
          %  AGAIN ARGUMENTS CAN BE EITHER VECTORS OR MATRICES
          %  HERE TWO "AXIS" LINES ARE DRAWN, THE MULTIPLE ARGUMENTS
          %  OF plot IS USEFUL, SINCE THE ARGUMENT MATRICES HAVE 
          %  DIFFERENT SIZES

          %  STRINGS AS EVERY THIRD ARGUMENT ALLOWS CONTROL OF COLORS
          %  c IS "CYAN" (LIGHT BLUE)
          %  m IS "MAGENTA" (PURPLE)
          %  w IS "WHITE"


  %  MULTIPLE ARGUMENTS ARE ALSO USEFUL TO GET DIFFERENT
  %  LINES TYPES AND COLORS FOR DIFFERENT CURVES:
  elseif isubpar == 3 ;
    plot(vx,my(:,1),'-r',vx,my(:,2),':g',vx,my(:,3),'-.b',...
                  vx,my(:,4),'--') ;
          %  RECALL THREE DOTS FOR CONTINUING A LONG LINE

          %  NOW THE STRINGS CONTROL BOTH LINE TYPES AND COLORS
          %  - is "solid line"
          %  : is "dotted line"
          %  -. is "dashed and dotted line"
          %  -- is "dashed line"


  %  CAN ANYBODY FIND A MATRIX WAY TO CONTROL LINE TYPES?
  %  I.E. A WAY TO CONTROL THE LINE TYPES IN THE COMMAND:
  %  plot(vx,my)

  end ;


elseif ipar == 4 ;    %  THEN DO PART 4 OF "GRAPHICS"

  %  HERE ARE SOME NESTED "CIRCLES":
  vt = (0:.1:1)' ;
  mx = cos(2 * pi * vt) * [.5 1 2 4] ;   
  my = sin(2 * pi * vt) * [.5 1 2 4] ;   
          %  THESE MATRIX MULTIPLICATIONS GIVES 4 COLUMN MATRICES

  plot(mx,my) ;

  isubpar = 0 ;     %  0,1,2,3,4

  %  TO MAKE A "SCATTERPLOT", WHERE THE SUCCESSIVE POINTS IN THE
  %  PLOT ARE NOT CONNECTED WITH LINE SEGMENTS, USE:
  if isubpar == 1 ;
    plot(mx,my,'*') ;
          %  * - PLOT A STAR, NOT CONNECTED BY LINES

  %  TO GET DIFFERENT SYMBOLS, CAN USE MULTIPLE ARGUMENTS OF PLOT
  %  (SOMETHING MORE CLEVER????)
  elseif isubpar == 2 ;
    plot(mx(:,1),my(:,1),'.',mx(:,2),my(:,2),'o',...
               mx(:,3),my(:,3),'x',mx(:,4),my(:,4),'+') ;
          %  . - PLOT A DOT, NOT CONNECTED BY LINES
          %  o - PLOT A CIRCLE, NOT CONNECTED BY LINES
          %  x - PLOT AN X, NOT CONNECTED BY LINES
          %  + - PLOT A PLUS SIGN, NOT CONNECTED BY LINES

  %  TO GET BOTH SYMBOLS, AND CONNECT THEM UP WITH LINES, CAN PLOT
  %  THE SAME CURVES TWICE WITHIN THE plot COMMAND
  elseif isubpar == 3 ;
    plot(mx,my,'ow',mx,my,'-r') ;

  %  ANOTHER APPROACH IS THE hold COMMAND:
  elseif isubpar == 4 ;
    plot(mx,my,'c') ;
    hold on ;
      plot(mx,my,'*g') ;
    hold off ;
          %  TRY COMMENTING OUT THE hold LINES, AND YOU GET ONLY THE
          %  SECOND OF THESE.  THE hold KEEPS THE PREVIOUS STUFF, AND
          %  YOU ADD THE NEW STUFF
          %  THIS IS USEFUL IN OTHER CONTEXTS AS WELL.

  end ;


elseif ipar == 5 ;    %  THEN DO PART 5 OF "GRAPHICS"

  %  A DIFFERENT WAY TO DISPLAY SEVERAL PLOTS IN THE SAME WINDOW (AND
  %  PRINT THEM ON THE SAME PAGE) IS TO USE subplot

  %  HERE IS A "TORNADO":
  vt = 5:.02:25 ;
  vcenter = (vt - 5) .* (vt - 15) / 4 ;
  xt = vcenter + vt .* cos(2 * pi * vt) ;
  yt = vt .* sin(2 * pi * vt) ;
  zt = vt ;

  %  HERE IS A "DRAFTSMAN'S PLOT" OF THE EACH OF THE THREE 2 DIMENSIONAL
  %  PROJECTIONS:
  subplot(2,2,1) ; plot(xt,yt) ; xlabel('x') ; ylabel('y') ;
          %  subplot(2,2 SAYS "CREATE A 2X2 ARRAY OF PLOTS TO FILL",
          %  THE 1 SAYS "MAKE THE FIRST ONE THE CURRENT ONE"
  subplot(2,2,2) ; plot(zt,yt) ; xlabel('z') ; ylabel('y') ;
          %  NOW GO TO THE SECOND ONE, WHICH IS TO THE RIGHT
  subplot(2,2,3) ; plot(xt,zt) ; xlabel('x') ; ylabel('z') ;
          %  THIRD ONE IS LOWER LEFT    

  disp('Hit any key to continue') ;
  pause ;

  %  NOW LOOK AT A THREE DIMENSIONAL VERSION
  clf ;
          %  COMMAND "clEAR fIGURE", REMOVES OLD STUFF IN WINDOW.
          %  IF THIS IS MISSING, THEN THIS NEXT PLOT GOEW WHERE THE 
          %  LAST SUBPLOT WAS.  TRY IT BY COMMENTING OUT THIS PART.
          %  NOTE THAT THIS CLEARS THE AXES LABELS AS WELL.
  plot3(xt,yt,zt) ;
          %  PLOT THE TORNADO AS A CURVE IN 3-D
  

elseif ipar == 6 ;    %  THEN DO PART 6 OF "GRAPHICS"

  %  A VERY USEFUL GRAPHICAL DEVICE IS "CHANGE OF SCALE".  THIS HAS
  %  2 COMMON PURPOSES:
  isubpar = 1 ;     %  1,2

  %  1.  REVEALING HIDDEN STRUCTURE
  if isubpar == 1 ;
    vx = .001:.001:1 ;
    vy = 1 + cos(30*vx) + 1 ./ vx ;

    subplot(1,2,1) ; plot(vx,vy) ; ylabel('y') ;
    subplot(1,2,2) ; plot(vx,log10(vy)) ; ylabel('log10(y)') ;

    disp('Hit any key to replace the plot of log10(y) with a semilog plot') ;
    pause ;
    subplot(1,2,2) ; semilogy(vx,vy) ; ylabel('y') ;
          %  THESE ARE THE SAME PLOTS, BUT THE AXES OF THE SEMILOG HAVE
          %  TICK MARKS THAT MAKE THE LOG STRUCTURE MORE CLEAR.

    %  ANOTHER APPROACH TO THE PROBLEM OF "MISSING FEATURES" IS TO PLOT
    %  ONLY PART OF THE Y-AXIS:
    disp('Hit any key to replace the plot of log10(y) with a smaller scale') ;
    pause ;
    subplot(1,2,2) ; plot(vx,vy) ; ylabel('y') ; 
    axis([0 1 0 10]) ;
          %  ARGUMENTS ARE Xmin, Xmax, Ymin, Ymax

          %  NOTE THAT THIS TOGETHER WITH THE PLOT ON THE RIGHT SHOW THE
          %  FULL PICTURE OF THIS CURVE.


  %  2.  CLARIFYING RELATIONSHIPS
  elseif isubpar == 2 ;
    vx = 1:.1:10 ;
    vy = (10 * vx.^3) .* (1.3 - .6 * rand(1,length(vx))) ;

    subplot(2,2,1) ;
      plot(vx,vy,'.') ; 
      title('"Raw" data') ;
      text(-1,13500,'log(y) may reveal hidden structure') ;
          %  GOT THESE COORDINATES BY FIRST LOOKING AT THE PLOT
      text(-1,12000,'or show if this is exponential') ;
          %  SINCE LOG WILL THEN LOOK LINEAR
    subplot(2,2,2) ;
      semilogy(vx,vy,'.') ;
      title('semilogy shows log(data)') ;
      text(-1,10^4.5,'Shows not exponential, try loglog') ;
          %  NOTE 10^4.3 IS MORE UNDERSTANDABLE ON THIS PLOT
    subplot(2,2,3) ;
      loglog(vx,vy,'.') ;
      title('log(y) vs. log(x)') ;
      text(10^(-.1),10^4.5,'Shows log10(y) ~ a log10(x) + b') ;
      text(10^(-.1),10^4,'with a ~ 3,  b ~ 1 ') ;
      text(10^(-.1),10^3.5,'i.e. y ~ 10*x^3') ;
          %  NOW THE X COORDINATE NEEDS TO BE CHANGED TO MAKE THE
          %  TEXT APPEAR (OTHERWISE IT IS OFF THE PAGE)

  end ;

  %  ANOTHER WAY TO VIEW SEVERAL PLOTS SIMULTANEOUSLY IS WITH THE 
  %  figure COMMAND.  THIS ACTUALLY CREATES NEW WINDOWS WHICH CAN 
  %  THEN BE MOVED AROUND AS DESIRED.


elseif ipar == 7 ;    %  THEN DO PART 7 OF "GRAPHICS"

  %  HERE IS A QUICK LOOK AT SOME FANCY WAYS THAT MATLAB CAN LOOK AT 
  %  2D SURFACES IN 3D.

  clf ;
          %  CLEARS AWAY THE WINDOW STUFF FROM LAST PART

  vx = (0:.1:2)' ;
  vy = 0:.2:5 ;
          %  DIFFERENT VALUE SO YOU CAN SEE WHAT IS WHERE
          %  MATLAB'S SURFACE GRAPHICS WORK EASIEST IF THE x'S ARE
          %  ALONG COLUMNS, AND THE Y'S ALONG ROWS

  [mx,my] = meshgrid(vx,vy) ;
          %  THESE ARE MATRIX VERSIONS, WHERE mx HAS THE SAME COLUMN
          %  VECTOR DUPLICATED SO IT HAS AS MANY ROWS AS vy, AND VICE
          %  VERSA

  mz = 1 ./ (1 + (mx - 1).^2 + (my - 4).^2) ;
          %  A CAUCHY LIKE BUMP, CENTERED AT (1,4)
  mz = mz + exp(-((mx - 1).^2 + (my - 2).^2) / .1) ;
          %  A GAUSSIAN LIKE BUMP CENTERED AT (1,2) 

  mesh(vx,vy,mz) ;
          %  THIS IS A STANDARD "WIRE MESH PLOT" OF THE SURFACE
          %  THINK CAREFULLY ABOUT THE SYMMETRIES EXPECTED HERE,
          %  AND LOOK CAREFULLY AT THE x AND y AXES.

  disp('Hit any key to see a contour plot') ;
  pause ;
  contour(vx,vy,mz,20) ;

  disp('Hit any key to see a nice 3d contour plot') ;
  pause ;
  contour3(vx,vy,mz,20) ;

  disp('Hit any key to see a surface plot') ;
  pause ;
  surf(vx,vy,mz) ;
          %  THIS IS SIMILAR TO THE WIRE MESH, BUT HAS THE "PANELS
          %  FILLED IN".

  disp('Hit any key to see a colormap plot') ;
  pause ;
  pcolor(vx,vy,mz) ;
          %  THIS IS SIMILAR TO THE SURFACE PLOT, BUT IS "FLATTENED
          %  DOWN TO TWO DIMENSIONS".

  disp('Hit any key to see another colormap plot') ;
  pause ;
  pcolor(vx,vy,mz) ;
  shading flat ;
          %  THIS COLORMAP ELIMINATES THE "BOUNDARY LINES"

  disp('Hit any key for another color scheme') ;
  pause ;
  pcolor(vx,vy,mz) ;
  shading flat ;
  colormap(hot) ;
          %  CHECK OUT THIS DIFFERENT "COLOR PALETTE"

  disp('Hit any key to improve the resolution') ;
  pause ;
  vx = (0:.04:2)' ;
  vy = 0:.1:5 ;
  [mx,my] = meshgrid(vx,vy) ;
  mz = 1 ./ (1 + (mx - 1).^2 + (my - 4).^2) ;
  mz = mz + exp(-((mx - 1).^2 + (my - 2).^2) / .1) ;
  pcolor(vx,vy,mz) ;
  shading flat ;

  disp('Hit any key for an interesting "light source" plot') ;
  pause ;
  clf ;
  surfl(vx,vy,mz) ;
  shading flat ;
  colormap(gray) ;
          %  FOR THIS, THINK ABOUT A LIGHT SHINING ON THE SURFACE.
          %  AREAS AWAY FROM THE LIGHT ARE "SHADED".


elseif ipar == 8 ;    %  THEN DO PART 8 OF "GRAPHICS"

  %  PRINTING IS DONE WITH THE COMMAND print.

  plot(sin(2 * pi * (0:.02:3))) ;
     title('Plot from MATLAB script s322eg6.m') ;

  print -dps s322eg6.ps ;
          %  THIS COMMAND CAN BE ISSUED EITHER FROM THE SCRIPT FILE,
          %  OR FROM THE COMMAND LINE, WHEN YOU WANT TO PRINT WHAT YOU
          %  SEE ON THE SCREEN.

          %  -dps SAYS "ORDINARY POSTSCRIPT"
          %  s322eg6.ps IS THE FILE THIS IS PRINTED TOO
          %  DEPENDING ON YOUR SYSTEM, THERE ARE VARIOUS WAYS TO
          %  GET THIS TO PRINT DIRECTLY

          %  INSTEAD OF MY POSTING THIS FILE, I SUGGEST YOU GENERATE
          %  IT YOURSELF, AND TRY BOTH PRINTING IT, AND LOOKING AT IT
          %  WITH A VIEWER LIKE GHOSTVIEW.

          %  -dps CAN BE REPLACED (IF IT WORKS FOR YOU) BY:
          %      -dpsc COLOR POSTSCRIPT
          %      -dps2 "LEVEL 2" POSTSCRIPT
          %      (THERE ARE OTHER VARIATIONS, SEE THE MANUAL)

end ;

