MATLAB INTRODUCTION

Matlab is a commercial product that is used widely by students
and some (but not all) faculty and researchers at IGPP.  It
provides a "high-level" progamming environment for computing
with a special emphasis on matrix operations.  Many specialized
functions are "hard wired" into Matlab so there is no need to
use subroutine libraries for these tasks.  It is also very easy
to make plots of Matlab output--many plotting functions are
also built into the program.

Matlab is available for PCs, Macs, and UNIX workstations.
The student version of Matlab for the Mac currently (Sept. 2005)
sells for $99 at the UCSD bookstore.  We have a site license for 
Matlab on the Mac network at IGPP.  Mac laptops with MatLab
installed as part of our site license need to be connected to
the network for MatLab to work.

Matlab is more user friendly and easier to learn than 
standard languages like Fortran and C.  It is good for
many problems, but ultimately does not have the flexibility
of C or Fortran.  Complicated Matlab scripts will often
run much slower than the Fortran or C equivalent.

Matlab is great for making quick plots but achieving
publication quality plots can be frustrating (although
I am told this has improved in the most recent versions
of Matlab).

Matlab is well documented with online help (look to
the menu bar when you are running the program in the
mode where it opens a new window, i.e., not the -nojvm
mode).  There are also many books written about Matlab
that you may find useful if you become a serious user.
Here we will give a brief introduction to help you get
started and to demonstrate some of MatLabs capabilities. 

Matlab can be used interactively as a super hand calculator,
or, more powerfully, run using scripts (i.e., programs).


INTERACTIVE MODE

Matlab on the Macs is normally in the Applications/MATLAB7
folder.  The executable has a funny looking icon that says
"MATLAB 7.0" (if that is the version you have).  You may
want to drag this icon onto your dock to make it easier to
launch Matlab.

Alternatively, you can run Matlab from within a Unix shell
(e.g., X11 or Terminal for the Mac) by simply typing:

> matlab

If  you don't like the fancy window that will come up, enter

> matlab -nojvm

Here is what this looks like on my Mac:

shearer@katmai 18> matlab -nojvm

                                          < M A T L A B >
                              Copyright 1984-2004 The MathWorks, Inc.
                                     Version 7.0.0.19901 (R14)
                                           May 06, 2004

 
  To get started, type one of these: helpwin, helpdesk, or demo.
  For product information, visit www.mathworks.com.
 
>>        

At the >> prompt, one can directly enter commmands.  As an example,
this is one way to compute 2+2:

>> a = 2

a =

     2

>> b = 2

b =

     2

>> c = a + b

c =

     4

It is often annoying to have Matlab always print out the value
of each variable.  To avoid this, put a semicolon after the
commands:

>> a = 2;
>> b = 2;
>> c = a + b;
>> c

c =

     4

Only the final line produces output.  Semicolons can also be
used to string together more than one command on the same line:

>> a = 2; b = 2; c = a + b; c

c =

     4

However, perhaps because I am used to programming in Fortran,
I find this coding style somewhat annoying.  It allows one
to brag about performing tasks is a small number of "lines"
of code, but makes the code much harder to read.

Of course Matlab will also allow more complicated operations:

>> a = 2;
>> b = -12;
>> c = 16;
>> quad1 = (-b + sqrt(b^2 - 4*a*c)) / (2*a)

quad1 =

     4


Note that * multiplies, / divides, and ^ raises to a power
Parentheses are used to determine the order of operations.

Trig functions are also supported.  Like Fortran or C, the arguments
are in radians, not degrees:

>> y = sin (pi/6)

y =

    0.5000

Note that "pi" is hardwired to pi; you don't have to set this.

One of the interesting (confusing?) things about Matlab is
that variables are considered arrays, even if they are single
numbers (scalars).  To see what all the variables currently defined 
are, use the "whos" command:

>> whos
  Name        Size         Bytes  Class

  a           1x1              8  double array
  b           1x1              8  double array
  c           1x1              8  double array
  quad1       1x1              8  double array
  y           1x1              8  double array

Grand total is 5 elements using 40 bytes

To remove all variables from memory, use the "clear" command:

>> clear
>> whos
>> 

The power of Matlab comes from the fact that you can perform
operations directly on matrices.  Example:

>> a = [1,2,3; 4,2,0; 1,1,2]

a =

     1     2     3
     4     2     0
     1     1     2

>> b = inv(a)

b =

   -0.6667    0.1667    1.0000
    1.3333    0.1667   -2.0000
   -0.3333   -0.1667    1.0000

>> c = a*b

c =

     1     0     0
     0     1     0
     0     0     1

Note how the a matrix is defined, using brackets
(don't try to use parenthesis!) and semicolons to
separate the different rows.  


EXAMPLE OF COMPUTING AND PLOTTING LEAST SQUARES FIT

Matrices need not be square.  For example, let's compute 
the least squares best-fitting line to the points (1,1), 
(2,1.9), (3,3.2),  and (4,4). First, let's input these points 
as x and y values:

>> x = [1,2,3,4]

x =

     1     2     3     4

>> y = [1,1.9,3.2,4]

y =

    1.0000    1.9000    3.2000    4.0000

Next, let's plot these points to see what they look like:

>> plot(x,y,'+')

Without the '+', the plot would contain a line drawn between
the points.  Note that the plot is rather awkward because
the plot limits are automatically set by the minimum and
maximum values of the data.  This can be fixed by using
one of the options from the axes command:

>> set(gca,'XLim',[0 5],'YLim',[0 5]), hold on
>> plot(x,y,'+')

but right away we see that Matlab is not as user friendly
as some other programs, such as plotxy, for making plots.

We want to find the best fitting line.  Let the equation 
of the line be y = a + bx, where a is the y-intercept
and b is the slope.  Then we can express the desired
fit to the data as the matrix equation:

| 1   |   | 1 1 |   | a |
| 1.9 | = | 1 2 | * | b |
| 3.2 |   | 1 3 |
| 4   |   | 1 4 |

where the goal is to adjust a and b so that the r.h.s.
best fits the l.h.s.  Note that the 2x4 matrix in this
case consists of a column vector with ones (that are multiplied
by a) and a column vector containing the x values of the
data points (which are multiplied by b, the slope).

This is a standard form for linear problems:  d=Gm where d
is a data vector with the y-values of the data, G is the
linear operator for predicting the data from the model,
and m is the desired model.  In this case the model is
simply the coefficients a and b of the line.  

We want d and m to be column vectors.  In our case, m will be 
a 2-vector that contains the y-intercept and slope of the best 
fitting line.  Matlab uses prime to indicate the transpose and 
thus we can enter:

>> G = [ [1,1,1,1]' x' ]

G =

     1     1
     1     2
     1     3
     1     4

>> d = y'

d =

    1.0000
    1.9000
    3.2000
    4.0000

Finally, we can use the standard formula for the least
squares solution to d=Gm. 
>> m = inv(G'*G)*G'*d

m =

   -0.0500
    1.0300

(Note to Geophysics Graduate students:  The formula for m
given above is the sort of thing that you are expected
to know how to derive before you take your departmental
exam.  See a linear algebra or statistics text for details)

The best-fitting line is of the form y= -0.05 + 1.03*x.
To plot this line, we compute the y values for x values
of 0.5 and 4.5:

>> x2 = [0.5,4.5]

x2 =

    0.5000    4.5000

>> d2 = [ 1,0.5; 1 4.5]*m

d2 =

    0.4650
    4.5850

>> plot(x2,d2)

Note that this works even though x2 is a row vector
and d2 is a column vector.

This problem can also be solved using the Matlab lscov
function:

>> m = lscov(G,d,eye(4))

m =

   -0.0500
    1.0300

where eye(4) is the Matlab name for the 4x4 identity
matrix.  When the covariance is known, this can be
replaced with the covariance matrix.


QUITTING MATLAB (VERY IMPORTANT)

Always exit the Unix Matlab program by entering quit, i.e.,

>> quit

DO NOT simply close the window or log out!!!!  This will
leave a Matlab job running and slow the machine down for
everyone else.  If this happens, type top or newtop to
find the job-id number and then say 'kill -9 <job-id>'

Note:  This was a problem when we ran MatLab on the Suns;
I don't actually know if it is also a problem with the
MatLab on the Macs.  To be safe, however, please make sure
that you quit Matlab before you logout.


SCRIPT FILE MODE

Matlab commands don't have to get very complicated before
it makes sense to use a script to save all of the commands.
This makes it much easier to make changes and to keep track
of what one is doing.

A Matlab script file is simply a series of Matlab commands.
The file name should have a ".m" suffix (just like C programs
end in ".c" and Fortran programs end in ".f")

For instance, our least squares fit example could be written
in the file linefit.m (using our favorite text editor) as:

% linefit.m is example of least squares fit to some points
clear             %always good to start with this
x = [1,2,3,4];
y = [1,1.9,3.2,4];
set(gca,'XLim',[0 5],'YLim',[0 5]), hold on
plot(x,y,'+')
G = [ [1,1,1,1]' x' ];
d = y';
m = inv(G'*G)*G'*d
x2 = [0.5,4.5];
d2 = [ 1,0.5; 1 4.5]*m;
plot(x2,d2)

Within the Matlab program one can then simply enter:

>> linefit

m =

   -0.0500
    1.0300

It is now, of course, much easier to make changes to the
data points or to the plot parameters without having to
type everything in again.  For example, here is how to add
axes labels and a "legend" to the plot:
                                                             (linefit2.m)
title('The best-fitting line to four points');
legend('The points','The best-fitting line',4);
xlabel('X coordinate');
ylabel('Y coordinate');


EXAMPLE OF FITTING QUADRATIC TO DATA

Our script can readily be modified to fit a quadratic or 
higher order polynomial to the data. For example, let's
find the best quadratic to the points (1,1.5), (2,1.9), (3,3.2),
and (4,4.0).  Note that these are identical to the points
given for the linefit.m example in the notes except that 
y1=1.5.  We thus define x and y as:
 
>> x = [1,2,3,4];
>> y = [1.5,1.9,3.2,4];

To fit y = a + bx + cx^2 we define G 

>> G = [ [1,1,1,1]' x' x'.^2 ];

Notice that the individual elements of a vector can be
squared using the .^ operator in Matlab.  We can then
obtain m as:

>> d = y';
>> m = inv(G'*G)*G'*d;

m =

   -0.1750
    1.1550
   -0.0250

Plotting this result is a little more complicated because
we need a large number of points to properly display a
curve.  Here is one way to do this:

>> x3 = 0.5:0.1:4.5;
>> n = length(x3);
>> G = [ones(1,n)' x3' x3'.^2 ];
>> d3 = G*m;
>> plot(x3,d3)

We first define x3 as a vector with values from 0.5 to 4.5
at increments of 0.1 (this is expressed in Matlab as 0.5:0.1:4.5).

We set n to the number of points in x3 using the Matlab length
function (41 in this case).

To avoid excessive typing, we use the Matlab "ones" function
to define the first column vector of G.  Related useful functions
are:
zeros -- a matrix filled with zeros
rand  -- a matrix filled with uniform random numbers (0 to 1)
randn -- a matrix filled with normally distributed random numbers
eye   -- identity matrix

The parabola does not fit the data much better!

Here is the complete script:
                                                         (quadfit.m)
% quadfit.m is example of quadratic fit to some points
clear; clf             %always good to start with this
x = [1,2,3,4];
y = [1.5,1.9,3.2,4];
set(gca,'XLim',[0 5],'YLim',[0 5]), hold on
plot(x,y,'+');
G = [ [1,1,1,1]' x' x'.^2 ];
d = y';
m = inv(G'*G)*G'*d
x3 = 0.5:0.1:4.5;
n = length(x3);
G = [ones(1,n)' x3' x3'.^2 ];
d3 = G*m;
plot(x3,d3)
title('The best-fitting quadratic to four points');
legend('The points','The best-fitting quadratic',4);
xlabel('X coordinate');
ylabel('Y coordinate');

Notice that we use "clf" following "clear" to erase the graphics
window.  Otherwise the new plot may land on top of an existing plot.

RANDOM NUMBER EXAMPLE

Here is another example--a program that makes a histogram
of some normally distributed random numbers:

% testran.m tests Matlab random numbers
clear                 %clear variables
clf                   %clear plot window
n = input('How many numbers? ');
y = randn(n,1);       %normally distributed random numbers
%hist(y)
x = -4:0.2:4;
hist(y,x)

In this case, notice how a number can be input using the
Matlab input function.  "randn" is a Matlab function that
returns normally distributed random numbers (mean=0,
standard deviation = 1).  The hist command produces
a histogram plot.  In this script, my first attempt
at the plot used a default binning spacing which I
did not like.  So I commented out this line by putting
a "%" sign at the beginning of the line.  In the second
hist command, I force the bin spacing to be at intervals
of 0.2.

We can compute the mean and standard deviation of these
numbers using the Matlab "mean" and "std" functions:

m = mean(y)
stdev = std(y)

As an exercise, we could go through the steps of computing
these ourselves by writing code that looks very similar to
BASIC or FORTRAN:

total=0;
for i=1:n
   total = total + y(i);
end
avg = total/n
var = 0;
for i=1:n
   var = var + (y(i) - avg)^2;
end
stdev = sqrt(var/(n-1)) 

True Matlab aficionados would sneer at the use of the for loops,
because they don't fully take advantage of the vector capabilities
of Matlab.  Rather, they would write something like:

total = sum(y);
avg = total/n
var = sum((y-avg).^2);
stdev = sqrt(var/(n-1))

Note that (y-avg) is a vector and thus ".^2" must be used 
to perform the squaring of its individual elements.  You 
will get an error if you use "^2" (try it!).  You can multiply
and divide individual array elements by using ".*" and "./"
and so on.

The complete script demonstrates all of these methods, as well
as showing how to make an xy-plot of sequential pairs of random numbers.
The "pause" function is used to wait until the user enters hits
return.


% testran.m tests Matlab random numbers
clear                 %clear variables
clf                   %clear plot window
n = input('How many numbers? ');
y = randn(n,1);       %normally distributed random numbers
%y
%hist(y,50)
x = -4:0.2:4;
hist(y,x)

pause

m = mean(y)
stdev = std(y)

pause

total=0;
for i=1:n
   total = total + y(i);
end
avg = total/n
var = 0;
for i=1:n
   var = var+(y(i) - avg)^2;
end
stdev = sqrt(var/(n-1))

pause

total = sum(y);
avg = total/n
var = sum((y - avg).^2);
stdev = sqrt(var/(n-1))
   
pause

x1 = rand(5000,1);
for i=1:4999
   x2(i+1)=x1(i); 
end
x2(1) = x1(5000);
n = length(x1)
n = length(x2)
plot(x1,x2,'.') 



SUMMARY OF SOME MATLAB COMMANDS

i:j:k    is regularly spaced vector i, i+j, i+2j, ..., k
A(:,j)   is jth column of A
A(i,:)   is ith row of A
a(j:k)   is a(j),a(j+1), .... , a(k)
A(:)     is vector with all elements of A, end-to-end

axes       places axes on figure (see manual for details)
axis       axis scaling and appearance 
break      exit for or while loop (as in C)
cd         change working directory
cla        clear current axis
clabel     contour plot elevation labels
clc        clear command window
clear      clear memory (remove variables)
clf        clear current figure (clf reset for total clear)
cond       matrix condition number
contour    make contour plot (see manual)
conv(a,b)  convolve vectors a and b
corrcoef   correlation coefficients (see manual)
cov        covariance matrix (see manual)
cputime    elapsed CPU time
cumprod    cumulative produce
cumsum     cumulative sum
deconv     deconvolution
del2       five-pint discrete Laplacian
delete     delete files and graphic objects
det        matrix determinant
diag       diagonal matrices and diagonals of matrix
diff       differences and approximate derivatives
dir        list files in current directory
disp       display text or matrix
echo       echo M-files during execution
eig        eigenvalues and eigenvectors (see manual)
else       used in if statements (as in C)
elseif     see above
end        terminate for, while, else statements (as in C)
eps        returns floating point accuracy
erf        error function
erfinv     inverse error function
exist      check if variable exists
expm       matrix exponential
eye(n)     n-by-n identity matrix
fclose     close file
fft        1-D fast Fourier transform (see manual)
fft2       2-D fast Fourier transform
figure     open new graphics window (see manual)
fill       fill 2-D polygons
find       find indices and values of nonzero elements
finite     detect infinities
fix        round to integer by going toward zero
fliplr     flip matrix left-right
flipud     flip matrix up-down
floor      round to integer by going toward -infinity
flops      counts floating point operations
fmin       minimize a function of one variable
fmins      minimize a function of several variables
fopen      open a file (similar to C)
for        for loop (similar to C)
format     control output display format
fplot      plot the graph of a function (see manual)
fprintf    write formatted data to file (see manual)
fread      read binary data from file
fscanf     read formatted data from file
fwrite     write binary data from matrix to file
fzero      zero of a function of one variable
ginput     graphical input using mouse
gradient   compute gradient of 2-D array
grid       grid lines for 2-D and 3-D plots
gtext      mouse placement of text on a graph
help       online documentation
hist       histogram plot
hold       hold the current graph so new stuff plots on top
if         conditionally execute statements (as in C)
ifft       inverse 1-D FFT
ifft2      inverse 2-D FFT
image      display image
imagesc    display image scaled to use full color map
input      request user input
interp1    1-D data interpolarion
interp2    2-D data interpolation
inv        matrix inverse
j          complex number i
length     length of a vector
line       create line object
load       retrieve Matlab variables from disk (see "save")
loglog     make log-log scale plot
lookfor    keyword search through help entries
lscov      least squares solution when covariance is known
lu         LU factorization of matrix
max        maximum elements of matrix
mean       mean value of vector or matrix
median     median value of vector or matrix
mesh       3-D mesh surface plot
min        minimum elements of matrix
nextpow2   next power of two (useful for FFTs)
nnl        solves non-negative least squares problem
ode23      solves ordinary differential equations
orient     hardcopy paper orientation
pi         3.1415927...
plot       linear 2-D plot
poly       polynomial with specified roots
polyfit    polynomial curve fitting
polyval    polynomial evaluation
print filename   saves fig. window as filename (Postscript)
qr         orthogonal-triangular matrix decomposition
quad       numerical evaluation of integrals
quit       terminate Matlab
quiver     plot vector field using little arrows
rand       random number
randn      normally distributed random number
randperm   random permutation
real       real part
roots      polynomial roots
rot90      rotate matrix 90 degrees
round      round to nearest integers
save       save workspace variables to disk
sort       sort elements in ascending order
sound      convert vector to sound
spline     cubic spline interpolation
std        standard deviation
surf       3-D shaded surface plot
svd        singular value decomposition
text       add text to current plot
title      graph title
trace      trace of matrix
unix       execute UNIX command from MATLAB
while      repeat statements (as in C)
who,whos   lists variable in memory
xlabel     label x-axis (also ylabel, zlabel)


PLOTTING IMAGES IN MATLAB

In geophysics, we often are interested in plotting images
such as bathymetry/topography, SAR data, etc.  Matlab provides
a way to make quick plots of these files.  For example, the
file siotopo.bin contains bathymetry data for the offshore
region near SIO.  The data are stored as a series of binary
numbers, in this case a 214 by 171 matrix.  Here is how to
display this file using Matlab:
                                              (plotsiotopo.m)
>> [fid,message] = fopen('siotopo.bin')
>> [z,count] = fread(fid,[171,214],'float');
>> fclose(fid);
>> imagesc(z);
>> print siotopo.ps;

We could have said "fid = fopen('siotopo.bin')" but the
message option provides more info if there is an error
openning the file.  Note that we need to know the dimensions
of the array and that the first 171 numbers are given in the
file first, then the next 171, etc.  The 'float' argument
indicates that the number are 4-byte reals.  In the case of
4-byte integers this would be 'int'

To make N-S the y-axis, the image must be rotated by 90 degrees:

>> z2=rot90(z,1);
>> imagesc(z2);

The units on the axes are just the elements (pixels) in the
image.  In the image, the x-spacing is 0.078 km (171 values), 
the y-spacing is 0.093 km (214 values).  We can add this
information into the plot as follows:

>> imagesc( [0:170]*0.078, [0:213]*0.093, z2); 

This works but the aspect ratio is wrong.  We can fix this by
entering:
                                           
>> axis equal;

Better, but now it has annoying white bars on the sides.
We can get rid of these and add a colorbar scale by entering:

>> axis tight;
>> colorbar;

There are also ways to edit the figures using the pulldown
menus in the figure window.  Of course, manual changes like
this negate many of the advantages of using scripts to make
the figures (easy repeatability, etc.).

Here is the complete script (plotsiotopo2.m):

%Matlab script to plot SIO topography data
clear; clf
[fid,message] = fopen('siotopo.bin')
[z,count] = fread(fid,[171,214],'float');
fclose(fid);
z2=rot90(z,1);
imagesc( [0:170]*0.078, [0:213]*0.093, z2);
axis equal;
axis tight;
colorbar;
print siotopo.ps

The final command makes a Postscript file of the plot, which
is handy for bringing the plot into other applications.

If we want a snazzy 3-D looking plot, use the surf command:

>> surf( [0:170]*0.078, [0:213]*0.093, z2);

(***oops, it's backwards!)

Nothing like a little vertical exaggeration to make things
look dramatic!  In my opinion, plots like this are fine
for talks but are annoying in papers.  Usually, map view plots
convey the information much more clearly.


MATLAB EXAMPLES

It is often easiest to learn a computer language by studying
and/or modifying examples of existing code.  Here are included
some examples from IGPP researchers.

Bob Parker provides 3 examples from his 229 class.  The file names
and their contents are listed below.  An example of running
them would be:

>> Ylm = sh(5,3);    %generate and plot spherical harmonic for l=5, m=3
>> turn              %rotates and turns the plot
>> Yl = rany(5);     %generates and plots sph. harm for l=5 with random m coef.

For more details, see Bob's website for this class at:
http://mahi.ucsd.edu/parker/SIO229/sio229.html

-------------------------------------------------------------------------
sh.m

function Ylm = sh(l, m);
%  function Ylm = sh(l, m);
%  Draw the equipotential surface of a multipole potential based on
%  Re Ylm, a spherical harmonic, degree l, order m
%  Also return numerical values evenly spaced in angle of Re Ylm.

n = 6;  %  Number of samples per half wavelength
N = max([30, l*6]);
theta = [0 : N]' *pi/N;
phi   =[-N : N]  *pi/N;

%  Create the surface of Re Ylm(theta,phi)/r^(l+1) = 1;
%  Solves for r on a theta, phi grid
gamma = 1/(l+1);
all = legendre(l, cos(theta));
if (l == 0) all=all'; end         % Compensate for error in legendre 
Ylm = all(m+1, :)' * cos(m*phi);
r = abs(Ylm) .^ gamma;

%  Convert to Cartesian coordinates
X =  r.* (sin(theta)*cos(phi)) ;
Y =  r.* (sin(theta)*sin(phi));
Z =  r.* (cos(theta)*ones(size(phi)));

%  Color according to size and sign of Ylm
C = Ylm;

surf(X, Y, Z, C)
axis equal

colormap hot
ti=['Surface  Y_l^m/r^{l+1}=1  with l=' int2str(l) ', m=' int2str(m)];
title(ti);
--------------------------------------------------------------------------

--------------------------------------------------------------------------
turn.m

%  turn
%  Rotates the viewing angle of the object in the current figure

dt=0.5;  % Pause time between frames in sec
az=-30;  % Initial azimuth, deg

pause(2) %  Gives time to move cursor into figure

% Rotate about horizontal axis
for elev = [-90 : 10 : 90] + 30
% disp([az, elev])
  view(az, elev)
  pause (dt)
end

% Rotate about vertical axis
elev=30;
for az = [0 : 10 : 360] + 30
% disp([az, elev])
  view(az, elev)
  pause (dt)
end
--------------------------------------------------------------------------

--------------------------------------------------------------------------
rany.m

function Yl = rany(l)
%  function Yl = rany(l)
%  Generates and contours a lat-long grid matrix of values of 
%  a single degree SH with random contributions of all orders.

%  Note this is still an eigenfunction of grad_1^2

%  Plots an equipotential surface of the associated multipole
%  potential


n = 6;  %  Number of samples per half wavelength
N = max([30, l*6]);
theta = [0 : N]' *pi/N;
phi   =[-N : N]  *pi/N;

%  Generate degree  l   Schmidt normalized functions (equal
%  energy for each order)
all = legendre(l, cos(theta), 'sch');

%  Random coeffs, normally distributed
a =  randn(1, l+1);
b =  randn(1, l+1);

%  Sum the different orders - contour as we go
Yl = a(1)* all(1,:)' * ones(size(phi));
for m = 1 : l
contour(Yl); pause(0.5)
  Yl = Yl + all(m+1, :)'  ...
         *(a(m+1)* cos(m*phi) + b(m+1)*sin(m*phi));
end
contour(Yl)
pause(2)

% Display equipotential surface
gamma=1/(1+l);
r = abs(Yl) .^ gamma;
%  Convert to Cartesian coordinates
X =  r.* (sin(theta)*cos(phi)) ;
Y =  r.* (sin(theta)*sin(phi));
Z =  r.* (cos(theta)*ones(size(phi)));

%  Color according to size and sign of Yl
C = Yl;

surf(X, Y, Z, C)
axis equal

colormap hot
ti=['Surface random harmonic function of degree l=' int2str(l)];
title(ti);
--------------------------------------------------------------------------