-
Notifications
You must be signed in to change notification settings - Fork 41
/
plotdata.m
84 lines (79 loc) · 3.23 KB
/
plotdata.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
% A program to create a plot of the computed results
% from the 3D Fortran Navier-Stokes solver
clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
'defaultaxesfontweight','bold')
% Load data
% Get coordinates
tdata=load('./data/tdata.dat');
x=load('./data/xcoord.dat');
y=load('./data/ycoord.dat');
z=load('./data/zcoord.dat');
nplots = length(tdata);
Nx = length(x); Nt = length(tdata);
Ny = length(y); Nz = length(z);
% reshape coordinates to allow easy plotting
[yy,xx,zz]=meshgrid(x,y,z);
for i =1:nplots
%
% Open file and dataset using the default properties.
%
FILEX=['./data/omegax',num2str(9999999+i),'.datbin'];
FILEY=['./data/omegay',num2str(9999999+i),'.datbin'];
FILEZ=['./data/omegaz',num2str(9999999+i),'.datbin'];
FILEPIC=['./data/pic',num2str(9999999+i),'.jpg'];
fid=fopen(FILEX,'r');
[fname,mode,mformat]=fopen(fid);
omegax=fread(fid,Nx*Ny*Nz,'real*8');
omegax=reshape(omegax,Nx,Ny,Nz);
fclose(fid);
fid=fopen(FILEY,'r');
[fname,mode,mformat]=fopen(fid);
omegay=fread(fid,Nx*Ny*Nz,'real*8');
omegay=reshape(omegay,Nx,Ny,Nz);
fclose(fid);
fid=fopen(FILEZ,'r');
[fname,mode,mformat]=fopen(fid);
omegaz=fread(fid,Nx*Ny*Nz,'real*8');
omegaz=reshape(omegaz,Nx,Ny,Nz);
fclose(fid);
%
% Plot data on the screen.
%
omegatot=omegax.^2+omegay.^2+omegaz.^2;
figure(100); clf;
subplot(2,2,1); title(['omega x ',num2str(tdata(i))]);
p1 = patch(isosurface(xx,yy,zz,omegax,.0025),...
'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
p2 = patch(isocaps(xx,yy,zz,omegax,.0025),...
'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
isonormals(omegax,p1); lighting phong;
xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); colorbar;
subplot(2,2,2); title(['omega y ',num2str(tdata(i))]);
p1 = patch(isosurface(xx,yy,zz,omegay,.0025),...
'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
p2 = patch(isocaps(xx,yy,zz,omegay,.0025),...
'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
isonormals(omegay,p1); lighting phong;
xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); colorbar;
subplot(2,2,3); title(['omega z ',num2str(tdata(i))]);
p1 = patch(isosurface(xx,yy,zz,omegaz,.0025),...
'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
p2 = patch(isocaps(xx,yy,zz,omegaz,.0025),...
'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
isonormals(omegaz,p1); lighting phong;
xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); colorbar;
subplot(2,2,4); title(['|omega|^2 ',num2str(tdata(i))]);
p1 = patch(isosurface(xx,yy,zz,omegatot,.0025),...
'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
p2 = patch(isocaps(xx,yy,zz,omegatot,.0025),...
'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
isonormals(omegatot,p1); lighting phong;
xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
axis equal; axis square; view(3);
saveas(100,FILEPIC);
end