Matlab in Chemical Engineering at CMU

Contents

Colors, 3D Plotting, and Data Manipulation

Guest authors: Harrison Rose and Breanna Stillo

In this post, Harrison and Breanna present three-dimensional experimental data, and show how to plot the data, fit curves through the data, and plot surfaces. You will need to download mycmap.mat

function main

close all
clear all
clc

In order to generate plots with a consistent and hopefully attractive color scheme, we will generate these color presets. Matlab normalizes colors to [R G B] arrays with values between 0 and 1. Since it is easier to find RGB values from 0 to 255, however, we will simply normalize them ourselves. A good colorpicker is http://www.colorpicker.com.

pcol = [255,0,0]/255; % red
lcol = [135,14,179]/255; % a pinkish color

Raw Data

This is raw data from an experiment for three trials (a and b and c). The X and Y are independent variables, and we measured Z.

X_a = [8.3 8.3 8.3 8.3 8.3 8.3 8.3];
X_b = [11 11 11 11 11 11 11];
X_c = [14 14 14 14 14 14 14];

Y_a = [34 59 64 39 35 36 49];
Y_b = [39 32 27 61 52 57 65];
Y_c = [63 33 38 50 54 68 22];

Z_a = [-3.59833 7.62 0 4.233333333 -2.54 -0.635 7.209];
Z_b = [16.51 10.16 6.77 5.08 15.24 13.7 3.048];
Z_c = [36 20 28 37 40 32 10];

Plotting the raw data for all three trials:

As you can see, the raw data does not look like very much, and it is pretty hard to interperet what it could mean.

We do see, however, that since X_a is all of one value, X_b is all of another value, and X_c is all of a third, that the data lies entirely on three separate planes.

figure
hold on     % Use this to plot multiple series on a single figure
plot3(X_a,Y_a,Z_a,'.','Color',pcol)
plot3(X_b,Y_b,Z_b,'.','Color',pcol)
plot3(X_c,Y_c,Z_c,'.','Color',pcol)
hold off    % Use this to make sure the next plot command will not
% try to plot on this figure as well.
title('Raw Experimental Data for Trials A, B, and C')
xlabel('x Data')
ylabel('y Data')
zlabel('z Data')
grid on     % 3D data is easier to visualize with the grid. Normally
% the grid defaults to 'on' but using the 'hold on'
% command as we did above causes the grid to default to
% 'off'

A note on the view:

The command

 view(Az,El)

lets you view a 3D plot from the same viewpoint each time you run the code. To determine the best viewpoint to use, use the click the 'Rotate 3D' icon in the figure toolbar (it looks like a box with a counterclockwise arrow around it), and drag your plot around to view it from different angles. You will notice the text "Az: ## El: ##" appear in the lower left corner of the figure window. This stands for Azimuth and Elevation which represent a point in spherical coordinates from which to view the plot (the radius is fixed by the axes sizes). The command used here will always display the plot from azimuth = -39, and elevation = 10.

view(-39,10)

A closer look at the raw data:

figure
hold on
plot(Y_a,Z_a,'o','MarkerFaceColor',pcol,'MarkerEdgeColor','none')
title('Raw Data for Trial A, x = 8.3')
xlabel('y')
ylabel('z')
hold off

figure
hold on
plot(Y_b,Z_b,'o','MarkerFaceColor',pcol,'MarkerEdgeColor','none')
title('Raw Data for Trial B, x = 11')
xlabel('y')
ylabel('z')
hold off

figure
hold on
plot(Y_c,Z_c,'o','MarkerFaceColor',pcol,'MarkerEdgeColor','none')
title('Raw Data for Trial C, x = 14')
xlabel('y')
ylabel('z')
hold off

Fitting the raw data:

In this case, we expect the data to fit the shape of a binomial distribution, so we use the following fit function with three parameters:

    function v = mygauss(par, t)
        A  = par(1);
        mu = par(2);
        s  = par(3);

        v=A*exp(-(t-mu).^2./(2*s.^2));
    end

Fitting the data

res = 20;
Yfit=linspace(20,70,res);

% Dataset A
guesses=[20, 40, 20];
[pars residuals J]=nlinfit(Y_a, Z_a-min(Z_a), @mygauss, guesses);

A1=pars(1);
mu1=pars(2);
s1=pars(3);

Zfitfun_a=@(y) A1*exp(-(y-mu1).^2./(2*s1.^2))+min(Z_a);
% Note: We have to shift the dataset up to zero because a gaussian
% typically cannot go below the horizontal axis

% Generate a fit-line through the data

Zfit_a=Zfitfun_a(Yfit);

% Dataset B
guesses=[10, 25, 20];
[pars residuals J]=nlinfit(Y_b, Z_b, @mygauss, guesses);

A2=pars(1);
mu2=pars(2);
s2=pars(3);

Zfitfun_b=@(y) A2*exp(-(y-mu2).^2./(2*s2.^2));

% Generate a fit-line through the data
Zfit_b=Zfitfun_b(Yfit);

% Dataset c
guesses=[20, 60, 20];
[pars residuals J]=nlinfit(Y_c, Z_c, @mygauss, guesses);

A3=pars(1);
mu3=pars(2);
s3=pars(3);

Zfitfun_c=@(y) A3*exp(-(y-mu3).^2./(2*s3.^2));

% Generate a fit-line through the data
Zfit_c=Zfitfun_c(Yfit);

Plotting the fitted data:

% Trial A
figure
hold on
plot(Y_a,Z_a,'o','MarkerFaceColor',pcol,'MarkerEdgeColor','none')
title('Fitted Data for Trial A, x = 8.3')
xlabel('y')
ylabel('z')
plot(Yfit,Zfit_a,'Color',lcol,'LineWidth',2)
hold off

% Trial B
figure
hold on
plot(Y_b,Z_b,'o','MarkerFaceColor',pcol,'MarkerEdgeColor','none')
title('Fitted Data for Trial B, x = 11')
xlabel('y')
ylabel('z')
plot(Yfit, Zfit_b,'Color',lcol,'LineWidth',2)
hold off

% Trial C
figure
hold on
plot(Y_c,Z_c,'o','MarkerFaceColor',pcol,'MarkerEdgeColor','none')
title('Fitted Data for Trial C, x = 14')
xlabel('y')
ylabel('z')
plot(Yfit,Zfit_c,'Color',lcol,'LineWidth',2)
hold off

Generate a surface plot:

For every point along the fit-line for dataset A, connect it to the corresponding point along the fit-line for dataset B using a straight line. This linear interpolation is done automatically by the surf command. If more datasets were available, we could use a nonlinear fit to produce a more accurate surface plot, but for now, we will assume that the experiment is well-behaved, and that 11 and 8.3 are close enough together that we can use linear interpolation between them.

Interpolate along the fit lines to produce YZ points for each data series (X):

Fits through the fits

Yspace = linspace(25,65,res);
Xs = [8.3; 11; 14];
Xspan1 = linspace(8.3,14,res);
Xspan2 = linspace(8,14.3,res);
q = 0;
r = 0;
for i = 1:res
    Zs = [Zfitfun_a(Yspace(i)); Zfitfun_b(Yspace(i)); Zfitfun_c(Yspace(i))];
    Yspan = linspace(Yspace(i),Yspace(i),res);
    p(:,i) = polyfit(Xs,Zs,2);
    for j = 1:res
        Zfit_fit(i,j) = polyval(p(:,i),Xspan1(j));
    end
end

Plot a surface through the fits through the fits

figure
hold all
surf(Xspan1,Yspace,Zfit_fit,'EdgeColor','black');   % Generate the surface
%surf(Xsurf,Ysurf,Zsurf,'EdgeColor','black');   % Generate the surface

% plot leading and lagging lines (just for show!)
for i = 1:res
    Zs = [Zfitfun_a(Yspace(i)); Zfitfun_b(Yspace(i)); Zfitfun_c(Yspace(i))];
    Yspan = linspace(Yspace(i),Yspace(i),res);
    plot3(Xspan2,Yspan,polyval(p(:,i),Xspan2),'Color','black')
end

%Plot the raw data:
plot3(X_a,Y_a,Z_a,'o','markersize',5,'MarkerFaceColor',pcol,'MarkerEdgeColor','none')
plot3(X_b,Y_b,Z_b,'o','markersize',5,'MarkerFaceColor',pcol,'MarkerEdgeColor','none')
plot3(X_c,Y_c,Z_c,'o','markersize',5,'MarkerFaceColor',pcol,'MarkerEdgeColor','none')

%Plot the fit lines:
plot3(linspace(8.3,8.3,res),Yfit,Zfit_a,'Color',lcol,'LineWidth',2)
plot3(linspace(11,11,res),Yfit,Zfit_b,'Color',lcol,'LineWidth',2)
plot3(linspace(14,14,res),Yfit,Zfit_c,'Color',lcol,'LineWidth',2)
view(-39,10)
alpha(.2)   % Make the plot transparent
title('z vs. x vs. y')
xlabel('x')
ylabel('y')
zlabel('z')
grid on

Add a line through the max

We must minimize the negative of each of our initial three fits to find the maximum of the fit.

f1 = @(y) -Zfitfun_a(y);
f2 = @(y) -Zfitfun_b(y);
f3 = @(y) -Zfitfun_c(y);

Ystar_a = fminbnd(f1,0,100);
Ystar_b = fminbnd(f2,0,100);
Ystar_c = fminbnd(f3,0,100);

Ystars = [Ystar_a; Ystar_b; Ystar_c];
Zstars = [Zfitfun_a(Ystar_a); Zfitfun_b(Ystar_b); Zfitfun_c(Ystar_c)];

hold on
plot3(Xs,Ystars,Zstars,'o','markersize',7,'MarkerFaceColor','white')

xy = polyfit(Xs,Ystars,2);
xz = polyfit(Xs,Zstars,2);

plot3(Xspan1,polyval(xy,Xspan1),polyval(xz,Xspan1),'Color','yellow','LineWidth',2)

A note on colormaps

To edit the color of the surface, you can apply a colormap. Matlab has several built-in colormaps (to see them, type 'doc colormap' in the Command Window. As you see here, however, we are using our own colormap, which has been stored in the file 'mycmap.mat'

To modify a colormap, type

 cmapeditor

in the Command Window or click 'Edit' 'Colormap...' in the figure window. For instructions on how to use the colormap editor, type 'doc colormapeditor' in the Command Window. If you have 'Immediate apply' checked, or you click 'Apply' the colormap will load onto the figure. To save a colormap, type the following into the Command Window:

 mycmap = get(gcf,'Colormap')

 save('mycmap')

s = load('mycmap.mat');
newmap = s.mycmap;
set(gcf,'Colormap',newmap)  % See corresponding 'Note on colormaps')

end

Summary

Matlab offers a lot of capability to analyze and present data in 3-dimensions.

% categories: plotting, data analysis

Published with MATLAB® 7.13

Contents

Interacting with graphs with context menus

John Kitchin

In Post 1489 we saw how to interact with a contour plot using mouse clicks. That is helpful, but limited to a small number of click types, e.g. left or right clicks, double clicks, etc... And it is not entirely intuitive which one to use if you don't know in advance. In Post 1513 , we increased the functionality by adding key presses to display different pieces of information. While that significantly increases what is possible, it is also difficult to remember what keys to press when there are a lot of them!

Today we examine how to make a context menu that you access by right-clicking on the graph. After you select a menu option, some code is run to do something helpful. We will examine an enthalpy-pressure diagram for steam, and set up functions so you can select one a few thermodynamic functions to display the relevant thermodynamic property under the mouse cursor.

See the video.

function main

The XSteam module makes it pretty easy to generate figures of the steam tables. Here we make the enthalpy-pressure graph.

clear all; close all
P = logspace(-2,3,500); %bar

temps = [1 50 100 150 200 300 400 500 600 700 800]; % degC
figure; hold all
for Ti=temps
    Hh = @(p) XSteam('h_PT',p,Ti);
    H = arrayfun(Hh,P);
    plot(H,P,'k-')
    % add a text label to each isotherm of the temperature
    text(H(end),P(end),sprintf(' %d ^{\\circ}C',Ti),'rotation',90);
end

% add the saturated liquid line
hsat_l = arrayfun(@(p) XSteam('hL_p',p),P);
plot(hsat_l,P)

% add the saturated vapor line
hsat_v = arrayfun(@(p) XSteam('hV_p',p),P);
plot(hsat_v,P)

xlabel('Enthalpy (kJ/kg)')
ylabel('Pressure (bar)')
set(gca,'YScale','log') % this is the traditional view
% adjust axes position so the temperature labels aren't cutoff.
set(gca,'Position',[0.13 0.11 0.775 0.7])

Now we add a textbox to the graph where we will eventually print the thermodynamic properties. We also put a tooltip on the box, which will provide some help if you hold your mouse over it long enough.

textbox_handle = uicontrol('Style','text',...
    'String','Hover mouse here for help',...
     'Position', [0 0 200 25],...
     'TooltipString','Right-click to get a context menu to select a thermodynamic property under the mouse cursor.');

% we will create a simple cursor to show us what point the printed property
% corresponds too
cursor_handle = plot(0,0,'r+');

setup the context menu

hcmenu = uicontextmenu;

Define the context menu items

item1 = uimenu(hcmenu, 'Label', 'H', 'Callback', @hcb1);
item2 = uimenu(hcmenu, 'Label', 'S', 'Callback', @hcb2);
item3 = uimenu(hcmenu, 'Label', 'T', 'Callback', @hcb3);

Now we define the callback functions

    function hcb1(~,~)
        % called when H is selected from the menu
        current_point = get(gca,'Currentpoint');
        H = current_point(1,1);
        P = current_point(1,2);
        set(cursor_handle,'Xdata',H,'Ydata',P)

        % We read H directly from the graph
        str = sprintf('H   = %1.4g kJ/kg', H);
        set(textbox_handle, 'String',str,...
            'FontWeight','bold','FontSize',12)
    end

    function hcb2(~,~)
        % called when S is selected from the menu
        current_point = get(gca,'Currentpoint');
        H = current_point(1,1);
        P = current_point(1,2);
        set(cursor_handle,'Xdata',H,'Ydata',P)

        % To get S, we first need the temperature at the point
        % selected
        T = fzero(@(T) H - XSteam('h_PT',P,T),[1 800]);
        % then we compute the entropy from XSteam
        str = sprintf('S   = %1.4g kJ/kg/K',XSteam('s_PT',P,T));
        set(textbox_handle, 'String',str,...
            'FontWeight','bold','FontSize',12)
    end

    function hcb3(~,~)
        % called when T is selected from the menu
        current_point = get(gca,'Currentpoint');
        H = current_point(1,1);
        P = current_point(1,2);
        set(cursor_handle,'Xdata',H,'Ydata',P)

        % We have to compute the temperature that is consistent with
        % the H,P point selected
        T = fzero(@(T) H - XSteam('h_PT',P,T),[1 800]);
        str = sprintf('T   = %1.4g C',  T);
        set(textbox_handle, 'String',str,...
            'FontWeight','bold','FontSize',12)
    end

set the context menu on the current axes and lines

set(gca,'uicontextmenu',hcmenu)

hlines = findall(gca,'Type','line');
% Attach the context menu to each line
for line = 1:length(hlines)
    set(hlines(line),'uicontextmenu',hcmenu)
end

end

% categories: plotting
% tags: interactive, thermodynamics

Published with MATLAB® 7.13

Contents

Interacting with graphs with keypresses

John Kitchin

In Post 1489 we saw how to interact with a contour plot using mouse clicks. That is helpful, but limited to a small number of click types, e.g. left or right clicks, double clicks, etc... And it is not entirely intuitive which one to use if you don't know in advance. Today we examine how to make key presses while the mouse is over the graph do something useful. We will examine an enthalpy-pressure diagram for steam, and setup functions so you can press u, s, p, h, or t to display the relevant thermodynamic property under the mouse cursor.

Watch the video here: http://screencast.com/t/nxAiiBXr

function main

The XSteam module makes it pretty easy to generate figures of the steam tables. Here we make the enthalpy-pressure graph.

clear all; close all
P = logspace(-2,3,500); %bar

temps = [1 50 100 150 200 300 400 500 600 700 800]; % degC
figure; hold all
for Ti=temps
    Hh = @(p) XSteam('h_PT',p,Ti);
    H = arrayfun(Hh,P);
    plot(H,P,'k-')
    % add a text label to each isotherm of the temperature
    text(H(end),P(end),sprintf(' %d ^{\\circ}C',Ti),'rotation',90);
end

% add the saturated liquid line
hsat_l = arrayfun(@(p) XSteam('hL_p',p),P);
plot(hsat_l,P)

% add the saturated vapor line
hsat_v = arrayfun(@(p) XSteam('hV_p',p),P);
plot(hsat_v,P)

xlabel('Enthalpy (kJ/kg)')
ylabel('Pressure (bar)')
set(gca,'YScale','log') % this is the traditional view
% adjust axes position so the temperature labels aren't cutoff.
set(gca,'Position',[0.13 0.11 0.775 0.7])

Now we add a textbox to the graph where we will eventually print the thermodynamic properties. We also put a tooltip on the box, which will provide some help if you hold your mouse over it long enough.

textbox_handle = uicontrol('Style','text',...
    'String','',... % start out empty
     'Position', [0 0 200 25],...
     'TooltipString','Press u,s,p,h or t to get the thermodynamic property under the mouse cursor. Ctrl-q to close the figure.');

% we will create a simple cursor to show us what point the printed property
% corresponds too
cursor_handle = plot(0,0,'r+');

Call back functions

we need to define two callback functions that will get called when the mouse moves over the figure, and when a key is pressed.

    function mouse_motion(~,~)
        % since we don't use (gcbo,event) the input arguments, we
        % replace them by ~
        current_point = get(gca,'Currentpoint');
        H = current_point(1,1);
        P = current_point(1,2);
        if (P > 0.01 && P < 1000) && (H > 0.01 && H < 4500)
            set(cursor_handle,'XData',H,'YData',P)
        end
    end

    function keypress_callback(gcbo,event)
        % We get the point under the cursor.
        H = get(cursor_handle,'XData');
        P = get(cursor_handle,'YData');

        % solve  for T
        T = fzero(@(T) H - XSteam('h_PT',P,T),[1 800]);

        % now define a string for each type of key press.
        str = ''; % default string to display
        switch event.Key
            case 'u'
                % compute internal energy for the P,T under cursor
                str = sprintf('U   = %1.4g kJ/kg', XSteam('u_PT',P,T));
            case 's'
                % compute entropy for the P,T under cursor
                str = sprintf('S   = %1.4g kJ/kg/K',XSteam('s_PT',P,T));
            case 'p'
                str = sprintf('P   = %1.4g bar',P);
            case 'h'
                str = sprintf('H   = %1.4g kJ/kg', H);
            case 't'
                str = sprintf('T   = %1.4g C',  T);
            case 'uparrow'
                if strcmp(event.Modifier,'control')
                    set(cursor_handle,'YData', P + 10)
                else
                    set(cursor_handle,'YData', P + 100)
                end
            case 'downarrow'
                if strcmp(event.Modifier,'control')
                    set(cursor_handle,'YData', P - 10)
                else
                    set(cursor_handle,'YData', P - 100)
                end
            case 'leftarrow'
                if strcmp(event.Modifier,'control')
                    set(cursor_handle,'XData', H - 10)
                else
                    set(cursor_handle,'XData', H - 100)
                end
            case 'rightarrow'
                if strcmp(event.Modifier,'control')
                    set(cursor_handle,'XData', H + 10)
                else
                    set(cursor_handle,'XData', H + 100)
                end
            case 'q'
                % quit. We use a dialog box to make sure we really
                % want to quit, in case q was accidentally pressed.
                choice = questdlg('Are you sure you want to quit?','','yes','no','no');
                switch choice
                    case 'yes'
                        close
                        return
                    case 'no'
                        str = '';
                end
            otherwise
                str = '';
        end

        % and put the string in our text box
        set(textbox_handle, 'String',str,...
            'FontWeight','bold','FontSize',12)
    end

finally, attach the callback functions to the current figure.

set(gcf,'KeyPressFcn',@keypress_callback)
set(gcf,'WindowButtonMotionFcn',@mouse_motion)

end

% categories: plotting
% tags: interactive, thermodynamics
% post_id = 1513; %delete this line to force new post;
% permaLink = http://matlab.cheme.cmu.edu/2011/12/07/interacting-with-graphs-with-keypresses/;

Published with MATLAB® 7.13

Contents

Turkeyfy your plots

John Kitchin

Need to make a festive graph? Today we look at putting a background graphic in your plot.

close all; clear all

Boring plot for thanksgiving

After a plateful of turkey you won't be able to keep your eyes open looking at this graph!

x = linspace(0,400,10);
y = x.^2;
plot(x,y)
title('Happy Thanksgiving')
xlabel('x')
ylabel('y')

Spice it up!

Let's add a background image, use a festive plot color, and a festive font.

Read the figure in

[X,MAP] = imread('8happy-thanksgiving-day.gif');

Plot the figure

note we have to flip it over so it comes out right,and set the colormap. the image is scaled by the xy data we have above. The image will be upside down for now, and the y-axis will be going in the wrong direction. we will fix that later.

image([min(x) max(x)], [min(y) max(y)], flipud(X));
colormap(MAP)

Plot data over the image

Hint: use the listfonts command to see what fonts are available on your machine.

hold on
plot(x,y,'o-','color',cmu.colors('Crimson glory'),'linewidth',2,...
    'markerfacecolor',cmu.colors('Crimson glory'))

title('Happy Thanksgiving','FontName','Lucida Handwriting')
xlabel('x','FontName','Lucida Handwriting')
ylabel('y','FontName','Lucida Handwriting')

Finally fix the y-axis

this makes the graph look like what you want! Much more festive.

set(gca,'ydir','normal');

% categories: plotting
% tags: whimsical

Published with MATLAB® 7.13

Contents

3D plots of the steam tables

John Kitchin

For fun let's make some 3D plots of the steam tables. We will use the XSteam module to map out the pressure, entropy and temperature of steam, and plot the 3D mapping.

clear all; close all; clc

create vectors that we will calculate T for

P = logspace(-2,3); %bar
S = linspace(0,10);

create meshgrid

these are matrices that map out a plane of P,S space

[xdata,ydata] = meshgrid(P,S);
zdata = zeros(size(xdata));

calculate T for each (p,S) pair in the meshgrids

zdata = arrayfun(@(p,s) XSteam('T_ps',p,s),xdata,ydata);

create the 3D plot

we use a lighted surface plot here

surfl(xdata,ydata,zdata)
xlabel('P (bar)')
ylabel('S (kJ/kg/K)')
zlabel('T (K)')
zlim([0 400])

customize the look of the surface plot

this smooths the faces out

shading interp

colormap pink

Add the saturated vapor and liquid lines, and the critical point

hold all

Tsat = arrayfun(@(p) XSteam('Tsat_p',p),P);
sV = arrayfun(@(p) XSteam('sV_p',p),P);
sL = arrayfun(@(p) XSteam('sL_p',p),P);
plot3(P,sV,Tsat,'r-','linewidth',5)
plot3(P,sL,Tsat,'g-','linewidth',5)

tc = 373.946; % C
pc = 220.584525; % bar
sc =XSteam('sV_p',pc);
plot3(pc,sc,tc,'ro','markerfacecolor','b') % the critical point
legend '' 'sat. vapor' 'sat. liquid' 'critical point'

Warning: Legend not supported for patches with FaceColor 'interp' 

Summary

There are lots of options in making 3D plots. They look nice, and from the right perspective can help see how different properties are related. Compare this graph to the one in Post 1484 , where isobars had to be plotted in the 2d graph.

% categories: plotting
% tags: thermodynamics

Published with MATLAB® 7.12