ABCDEN/plotphabc.m at master · Birfy/ABCDEN · GitHub

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
84
85
86
87
88
fid=fopen('pha.dat','r');
a=fscanf(fid,'%f %f %f %f %f %f',[6,inf]);
a=a';

% Nx=24;
% Ny=32;
% Nz=32;

pha1(1:Nx,1:Ny,1:Nz)=0.0;
phb1(1:Nx,1:Ny,1:Nz)=0.0;
phc1(1:Nx,1:Ny,1:Nz)=0.0;


dx=0.3;
dy=0.3;
dz=0.3;

n=2;

[X,Y,Z]=meshgrid(dy:dy:Ny*dy*n,dx:dx:Nx*dx*n,dz:dz:Nz*dz*n);

for k=1:Nz
    for i=1:Nx
        for j=1:Ny
            pha1(i,j,k)=a((i-1)*Ny*Nz+(j-1)*Nz+k,1);
            phb1(i,j,k)=a((i-1)*Ny*Nz+(j-1)*Nz+k,2);
            phc1(i,j,k)=a((i-1)*Ny*Nz+(j-1)*Nz+k,3);
        end
    end
end


pha=zeros(Nx*n,Ny*n,Nz*n);
phb=zeros(Nx*n,Ny*n,Nz*n);
phc=zeros(Nx*n,Ny*n,Nz*n);
for l=1:n
    for m=1:n
        for o=1:n
            for k=1:Nz
                for i=1:Nx
                    for j=1:Ny
                        pha(i+Nx*(l-1),j+Ny*(m-1),k+Nz*(o-1))=pha1(i,j,k);
                        phb(i+Nx*(l-1),j+Ny*(m-1),k+Nz*(o-1))=phb1(i,j,k);
                        phc(i+Nx*(l-1),j+Ny*(m-1),k+Nz*(o-1))=phc1(i,j,k);
                    end
                end
            end
        end
    end
end

pa = patch(isosurface(X,Y,Z,pha,0.5));
patch(isocaps(X,Y,Z,pha,0.5),'facealpha',1,...
    'FaceColor','blue',...
    'EdgeColor','none',...
    'AmbientStrength',.2,...
    'SpecularStrength',.5,...
    'DiffuseStrength',.3);
isonormals(X,Y,Z,pha,pa)

set(pa,'FaceColor','blue','EdgeColor','none')

pb = patch(isosurface(X,Y,Z,phb,0.5));
patch(isocaps(X,Y,Z,phb,0.5),'facealpha',1,...
    'FaceColor','green',...
    'EdgeColor','none',...
    'AmbientStrength',.2,...
    'SpecularStrength',.5,...
    'DiffuseStrength',.3);
isonormals(X,Y,Z,phb,pb)
set(pb,'FaceColor','green','EdgeColor','none') %transparency

pc = patch(isosurface(X,Y,Z,phc,1));
patch(isocaps(X,Y,Z,phc,1),'facealpha',1,...
    'FaceColor','red',...
    'EdgeColor','none',...
    'AmbientStrength',.2,...
    'SpecularStrength',.5,...
    'DiffuseStrength',.3);
isonormals(X,Y,Z,phc,pc)
set(pc,'FaceColor','red','EdgeColor','none') %transparency

axis off
daspect([1 1 1]) %set the axial ratio
view([1 1 1]); axis tight
camlight(10,90);lighting phong
light;
light;