% Figure4p6.m
% Plot of Legendre & Legendre Figure 4.6
% rn sets the range of the x and y axes;
% Written by E. Gallagher, ECOS, 2/26/02
% Eugene.Gallagher@umb.edu
rn=3.5;
% set standard deviations, means and correlation:
s1=1;s2=1;mu1=0;mu2=0;rho=0;
% meshgrid is the built-in matlab function needed for all 3-d plots
[X,Y]=meshgrid(-rn:.2:rn, rn:-.2:-rn);
% Need to get the size of the X and Y arrays for repmat function
[r,c]=size(X);
% Set the limits for the graphs
lims=[-rn rn -rn rn 0 .2]; % set the range for the plotting axes
% da will be used to set the scaling of the x, y, and z axes
da=[25 25 1];
% v will be used to set the viewing angle
v=[-30 33];
% Plot 4.6 right
% Z is the value of the probability density function at X and Y
% Z uses the formula for the bivariate normal distribution shown on
% page 149 of Legendre & Legendre, 2nd edition.
Z=1/(2*pi*s1*s2*(1-rho)^2)*exp(-0.5*1/(1-rho^2)*...
    (((X-repmat(mu1,r,c))./repmat(s1,r,c)).^2 - 2*rho* ...
    (X-repmat(mu1,r,c))./repmat(s1,r,c).*...
    (Y-repmat(mu2,r,c))./repmat(s2,r,c)+...
    ((Y-repmat(mu2,r,c))./repmat(s2,r,c)).^2 ));
% subplot is used for multiple graphs, see help subplot for information
% I'll plot the 'tallest' 3-d graph first to save the colorbar scaling info
subplot(2,2,2) % plot this graph in the upper right position
surf(X,Y,Z,'FaceColor','interp',...
    'Edgecolor','none',...
    'Facelighting','phong')
colormap(hsv)
title('sx = sy = 1, r= 0')
V=caxis; % save the value of the coloring scheme
axis(lims)
daspect(da)
view(v)
camlight left
material shiny
% Plot Figure 4.6 (left)
s1=1.43; % change the standard deviation of X
Z=1/(2*pi*s1*s2*(1-rho)^2)*exp(-0.5*1/(1-rho^2)*...
    (((X-repmat(mu1,r,c))./repmat(s1,r,c)).^2 - 2*rho* ...
    (X-repmat(mu1,r,c))./repmat(s1,r,c).*...
    (Y-repmat(mu2,r,c))./repmat(s2,r,c)+...
    ((Y-repmat(mu2,r,c))./repmat(s2,r,c)).^2));
subplot(2,2,1) % Plot this graph in the 1st position
surf(X,Y,Z,'FaceColor','interp',...
    'Edgecolor','none',...
    'Facelighting','phong')
title('sx = 1.43, sy = 1, r = 0')
caxis(V) % Use the colorscheme from the upper right graph, saved as V
axis(lims);daspect(da);view(v);camlight left;material shiny
colorbar % Add the colorbar

% Now plot with a correlation of 4.6left with a correlation of -.745
lims=[-rn rn -rn rn 0 .05]; % reduce the z axis limit to 0.05
da=[60 60 1];
s1=1.43;rho=-.745; % change the correlation to -.745
Z=1/(2*pi*s1*s2*(1-rho)^2)*exp(-0.5*1/(1-rho^2)*...
    (((X-repmat(mu1,r,c))./repmat(s1,r,c)).^2 - 2*rho* ...
    (X-repmat(mu1,r,c))./repmat(s1,r,c).*...
    (Y-repmat(mu2,r,c))./repmat(s2,r,c)+...
    ((Y-repmat(mu2,r,c))./repmat(s2,r,c)).^2 ));
subplot(2,2,3) % plot in the lower left position
surf(X,Y,Z,'FaceColor','interp',...
    'Edgecolor','none',...
    'Facelighting','phong')
title('sx = 1.43, sy = 1, r = -.745')
caxis(V);axis(lims);daspect(da);view(v);camlight left;material shiny
% Now plot the lower right figure
s1=1;
Z=1/(2*pi*s1*s2*(1-rho)^2)*exp(-0.5*1/(1-rho^2)*...
    (((X-repmat(mu1,r,c))./repmat(s1,r,c)).^2 - 2*rho* ...
    (X-repmat(mu1,r,c))./repmat(s1,r,c).*...
    (Y-repmat(mu2,r,c))./repmat(s2,r,c)+...
    ((Y-repmat(mu2,r,c))./repmat(s2,r,c)).^2 ));
subplot(2,2,4)
surf(X,Y,Z,'FaceColor','interp',...
    'Edgecolor','none',...
    'Facelighting','phong')
title('sx = sy = 1, r = -.745')
caxis(V);axis(lims);daspect(da);view(v);camlight left;material shiny
figure(gcf) % display the figure
% print -depsc -tiff -r300 fig46 % saves the eps figure with a tiff preview
% print -depsc fig46
% print -dtiff -r600 fig46