% Figure4p6cont.m
% Plot of Legendre & Legendre Figure 4.6
% Using filled contours
% rn sets the range of the x and y axes;
% Written by E. Gallagher, EEOS, 2/26/02
% Revised 7/8/2010
% Eugene.Gallagher@umb.edu
rn=4;
% 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]; % set the range for the plotting axes
% 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
contourf(X,Y,Z,10)
colormap(jet)
title('sx = sy = 1, r= 0')
V=caxis; % save the value of the coloring scheme
%axis(lims);
axis square
% 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
contourf(X,Y,Z,10)
title('sx = 1.43, sy = 1, r = 0')
caxis(V) % Use the colorscheme from the upper right graph, saved as V
%axis(lims);
axis square;
colorbar % Add the colorbar
% Now plot with a correlation of 4.6left with a correlation of -.745
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
contourf(X,Y,Z,10)
title('sx = 1.43, sy = 1, r = -.745')
caxis(V);% axis(lims);
axis square;
% 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)
contourf(X,Y,Z,10)
title('sx = sy = 1, r = -.745')
caxis(V);%axis(lims);
axis square;
figure(gcf) % display the figure
% print -depsc -tiff -r300 fig46c % saves the eps figure with a tiff preview
% print -depsc fig46c
% print -dtiff -r600 fig46c