In 2012, I researched computational materials science/physics for Dr. Dimitrios Papaconstantopoulos. The following is my semester report. All source code is listed at the end of the report.
Category: Matlab
Matlab: Smooth Rotating Animation for Line Plots
I recently became stuck trying to create an animation which consists of a smooth rotation of a viewpoint around the Lorenz attractor. Methods I use for changing viewpoints with respect to surface objects were creating jerky, lagging animations when applied to my line plot.
This will work for any 3D line plot:
plot3(x,y,z);
axis vis3d
fps = 60; sec = 10;
vidObj = VideoWriter('plotrotation.avi');
vidObj.Quality = 100;
vidObj.FrameRate = fps;
open(vidObj);
for i=1:fps*sec
camorbit(0.9,-0.1);
writeVideo(vidObj,getframe(gcf));
end
close(vidObj);
This results in the following smooth animation:
Matlab: Lorenz Attractor
I’m a big fan of the Lorenz Attractor, which, when plotted, resembles the half open wings of a butterfly. This attractor was derived from a simplified model of convection in the earth’s atmosphere. One simple version of the Lorenz attractor is pictured below:
The Lorentz system is a set of ordinary differential equations notable for its chaotic solutions (see below). Here \(x\), \(y\) and \(z\) make up the system state, \(t\) is time, and \(\rho, \sigma, \beta\) are the system parameters.
The Lorentz attractor is a chaotic solution to this system found when \(\rho = 28, \sigma = 10, \beta = 8/3\).
The series does not form limit cycles nor does it ever reach a steady state.
We can calculate and render the aforementioned chaotic solution to this ODE as follows:
function loren3
clear;clf
global A B R
A = 10;
B = 8/3;
R = 28;
u0 = 100*(rand(3,1) - 0.5);
[t,u] = ode45(@lor2,[0,100],u0);
N = find(t>10); v = u(N,:);
x = v(:,1);
y = v(:,2);
z = v(:,3);
plot3(x,y,z);
view(158, 14)
function uprime = lor2(t,u)
global A B R
uprime = zeros(3,1);
uprime(1) = -A*u(1) + A*u(2);
uprime(2) = R*u(1) - u(2) - u(1)*u(3);
uprime(3) = -B*u(3) + u(1)*u(2);
This results in the figure:
Matlab: Coaxial Cylinders (Polar Coordinates)
Let’s say we want to create an aesthetically pleasing visualization of 2 coaxial cylinders.
To do this, we’ll be adjusting the lighting, proportions, and transparency of the figure.
The code to create this figure utilizes the coordinate transformation from Cartesian to 3D-Polar coordinates. Recall that the transformation between coordinate systems is as follows.
Instead of using Matlab’s built in cart2pol
method, we will manually convert each Cartesian (x, y, z) coordinate to it’s equivalent polar coordinate (r, phi, z) .
Side note:
I used phi above in the transformation equations and theta below in the code to make a point. The angular coordinate in the cylindrical (a.k.a polar) coordinate system is generally represented as either phi or theta.
These are equivalent, and the variable chosen is purely a matter of notation.
Y=-5:5;
theta=linspace(0,2*pi,40);
[Y,theta]=meshgrid(Y,theta);
r = 1.5
% calculate x and z
X=r*cos(theta);
Z=r*sin(theta);
hs = surf(X,Y,Z)
set(hs,'EdgeColor','None', ...
'FaceColor', [0.5 0.5 0.5], 'FaceLighting', 'phong');
alpha(0.7);
hold on
camlight right;
r = 0.5
% recalculate x and z
X=r*cos(theta);
Z=r*sin(theta);
hs = surf(X,Y,Z)
axis equal
set(hs,'EdgeColor','None', ...
'FaceColor', [0.5 0.5 0.5], 'FaceLighting', 'phong');
alpha(0.7);
axis off
camlight right;
lighting gouraud
view(140, 24)
% white background
set(gcf,'color','white')
Matlab: Create Mesh or Surface From Line Plot
This is a continuation of Matlab: Lorentz Attractor, however, these methods can be applied to any line plot or collection of points.
A slightly more aesthetically pleasing representation of the Lorentz Attractor can be achieved by adding axis off
. And altering the view’s azimuth and elevation: view(15, 48)
.
Now we’re talking. Let’s say I want to make a surface or mesh from this dandy line plot. Using surf
or mesh
will throw an error, since x, y, and z are all 1D vectors! Whatever shall we do!
Never fear, mathematicians will save the day.
Delaunay created a sweet method of triangulating points. If we treat this line plot as a collection of points, we can triangulate to find an approximate surface.
tri = delaunay(x,y);
plot(x,y,'.')
%determine amount of triangles
[r,c] = size(tri);
disp(r)
%plot with trisurf
h = trimesh(tri, x, y, z, 'FaceAlpha', 0.6);
alpha = 0.4
view(15, 48)
axis vis3d
axis off
l = light('Position',[-50 -15 29])
lighting phong
shading interp
What if we’d like a surface instead of the mesh? Then we’ll change trimesh
to trisurf
add transparency (alpha = 0.7
) and find:
Matlab: Rotating Sphere Animation
Let’s say you want to make a simple simulation of a sphere spinning in Matlab.
First, you set the pop-up window to have the title ‘Spheres,’ the window to have black background, and specify said window’s position and size.
format compact
set(gcf,'Menubar','none','Name','Spheres', ...
'NumberTitle','off','Position',[10,350,400,300], ...
'Color',[0 0 0]);
Matlab has a nice built in “sphere” surface function that you can invoke after you specify where you want it’s axes. Note that if you don’t include the
h = axes('Position',[0 0 1 1]);
command, the program will run the same except your view of the sphere will be farther away (smaller).
%Zooms in on sphere
h = axes('Position',[0 0 1 1]);
%Draws sphere
[X Y Z]=sphere(20);
C = Z^2 + Y^2; %creates color data to map onto sphere
hs = surf(X, Y, Z, C);
The sphere function creates a “curved” surface using quadrilaterals. The boundaries of these quadrilaterals can be hidden set(hs,'EdgeColor','None');
or shown set(hs,'EdgeColor',[0.5 0.5 0.5]);
%Shows wires
set(hs,'EdgeColor',[0.5 0.5 0.5]);
%Adjusts transparency
alpha('color');
alphamap('rampdown');
%Adjusts lighting
camlight right;
lighting phong
hidden off
axis off
axis equal
Stopping here, you have a static sphere.
Depending on whether you wish to just play the animation or if you wish to also save it as an ‘.avi’ file, you approach the animation in different ways. Either way, the animation looks like this:
To play the animation without saving it:
oh=findobj(gca,'type','surface');
%Spins about z axis.
for i = 1:36
axis off;
rotate(oh,[0 0 1],10);
M(i) = getframe(gca);
end
%Spins about y axis.
for i = 1:36
axis off;
rotate(oh,[0 1 0],10);
M(i) = getframe(gca);
end
%Spins about x axis.
for i = 1:36
axis off;
rotate(oh,[1 0 0],10);
M(i) = getframe(gca);
end
To save the movie as an ‘.avi’ for future use, you can modify the above code section slightly. You create a writer object and write the video to it frame by frame (using writeVideo) as part of each forloop.
%Create writerObj
writerObj = VideoWriter('sphere.avi');
open(writerObj)
%Animated movie of the rotation of the 3D globe
oh=findobj(gca,'type','surface');
%Spins about z axis.
for i = 1:36
axis off;
rotate(oh,[0 0 1],10);
M(i) = getframe(gca);
frame = getframe;
writeVideo(writerObj,frame);
end
%Spins about y axis.
for i = 1:36
axis off;
rotate(oh,[0 1 0],10);
M(i) = getframe(gca);
frame = getframe;
writeVideo(writerObj,frame);
end
%Spins about x axis.
for i = 1:36
axis off;
rotate(oh,[1 0 0],10);
M(i) = getframe(gca);
frame = getframe;
writeVideo(writerObj,frame);
end
close(writerObj);
Putting it all together:
format compact
set(gcf,'Menubar','none','Name','Spheres', ...
'NumberTitle','off','Position',[10,350,400,300], ...
'Color',[0 0 0]);
h = axes('Position',[0 0 1 1]);
%Draws sphere
[X Y Z]=sphere(20);
C = Z^2 + Y^2;
hs = surf(X, Y, Z, C);
%sets wireframe to visible
set(hs,'EdgeColor',[0.5 0.5 0.5]);
%Alters transparency of sphere
alpha('color');
alphamap('rampdown');
%Sets lighting of sphere
camlight right;
lighting phong
hidden off
axis equal
%Create writerObj
writerObj = VideoWriter('sphere.avi');
open(writerObj)
%Animated movie of the rotation of the 3D globe
oh=findobj(gca,'type','surface');
%Spins about z axis.
for i = 1:36
axis off;
rotate(oh,[0 0 1],10);
M(i) = getframe(gca);
frame = getframe;
writeVideo(writerObj,frame);
end
%Spins about y axis.
for i = 1:36
axis off;
rotate(oh,[0 1 0],10);
M(i) = getframe(gca);
frame = getframe;
writeVideo(writerObj,frame);
end
%Spins about x axis.
for i = 1:36
axis off;
rotate(oh,[1 0 0],10);
M(i) = getframe(gca);
frame = getframe;
writeVideo(writerObj,frame);
end
close(writerObj);
Matlab: 2 Methods of Generating a Rainbow
In my data visualization class, we had an assignment to create a “rainbow” (create and display 128 vertical stripes of color in one image, in RGB sequence). Something like this:
In Matlab, a colormap is an m-by-3 matrix of real numbers between 0.0 and 1.0. Each row is an RGB vector that defines one color. The [R,G,B] color attribute value for a dataset with attribute value ranging between 0 and 1.
My original solution is:
image(1:128);
colormap(jet(128))
print -dpng 'rainbow.png'
The command colormap(jet(128))
creates a rainbow colormap with 128 colors. Files in the color folder generate a number of colormaps, and each file (in this case, “jet”) accepts the colormap size as an argument.
While this solution is certainly concise and clear, it is not the solution the professor is looking for.
It turns out that the professor wants the map to be generated algorithmically (instead of using a built-in file from the color folder).
My second solution is:
image(1:128);
R(128,1) = 0;
G(128,1) = 0;
B(128,1) = 0;
dx = 0.8;
for f=0:(1/128):1
g = (6-2*dx)*f+dx;
index = int16(f*128 + 1);
R(index,1) = max(0,(3-abs(g-4)-abs(g-5))/2);
G(index,1) = max(0,(4-abs(g-2)-abs(g-4))/2);
B(index,1) = max(0,(3-abs(g-1)-abs(g-2))/2);
end
%concatenate arrays horizontally
map = horzcat(R, G, B);
%map colors onto image
colormap(map)
print -dpng 'rainbow.png'
This creates a different rainbow and satisfies the professor’s algorithmic requirement.