% plot phase field configuration from MD++ simulations

%dirname = '../MD++.svn/runs/phasefield-0/';
%dirname = '../MD++.svn/runs/phasefield-multi-0/';
dirname = './phase_field/';
%filename = 's50.cn';
filename = 's75.cn';
%filename = 'inter0001.cn';
%filename = 'inter0006.cn';
%filename = 'inter0019.cn';

data = load([dirname,filename]);
NF = data(1); NX = data(2); NY = data(3); NZ = data(4);

x = [0:NX-1]; y = [0:NY-1]; z = [0:NZ-1];

for i=1:NF,
    tmp = reshape(data([5:NX*NY*NZ+4]+(i-1)*NX*NY*NZ), NX, NY, NZ);
    phi{i} = shiftdim(tmp,1);
end


if NF>=1
    figure(1);
    clf
    p_1 = patch(isosurface(x,y,z,phi{1},0.5));
    isonormals(x,y,z,phi{1},p_1)
    set(p_1,'FaceColor','cyan','EdgeColor','none');
    daspect([1 1 1])
    view(3); axis tight
    camlight 
    lighting gouraud
    xlim([0 NX]); ylim([0 NY]); zlim([0 NZ]);
    alpha(0.7);
end

if NF==2,
    figure(1);
    p_2 = patch(isosurface(x,y,z,phi{2},0.5));
    isonormals(x,y,z,phi{2},p_2)
    set(p_2,'FaceColor','yellow','EdgeColor','none');
    daspect([1 1 1])
    view(3); axis tight
    camlight 
    lighting gouraud
    xlim([0 NX]); ylim([0 NY]); zlim([0 NZ]);
    alpha(0.7);
end

if NF>=3
    p_3 = patch(isosurface(x,y,z,phi{3},0.5));
    isonormals(x,y,z,phi{3},p_3)
    set(p_3,'FaceColor','magenta','EdgeColor','none');
    daspect([1 1 1])
    view(3); axis tight
    camlight 
    lighting gouraud
    xlim([0 NX]); ylim([0 NY]); zlim([0 NZ]);
    alpha(0.7);
    
    figure(2)
    clf
    p_2 = patch(isosurface(x,y,z,phi{2},0.5));
    isonormals(x,y,z,phi{2},p_2)
    set(p_2,'FaceColor','yellow','EdgeColor','none');
    daspect([1 1 1])
    view(3); axis tight
    camlight 
    lighting gouraud
    xlim([0 NX]); ylim([0 NY]); zlim([0 NZ]);
    alpha(0.7);
    
    figure(3)
    clf
    p_1 = patch(isosurface(x,y,z,phi{1},0.5));
    isonormals(x,y,z,phi{1},p_1)
    set(p_1,'FaceColor','cyan','EdgeColor','none');
    p_2 = patch(isosurface(x,y,z,phi{2},0.5));
    isonormals(x,y,z,phi{2},p_2)
    set(p_2,'FaceColor','yellow','EdgeColor','none');
    daspect([1 1 1])
    view(3); axis tight
    camlight 
    lighting gouraud
    xlim([0 NX]); ylim([0 NY]); zlim([0 NZ]);
    alpha(0.7);
end


xlabel('x');ylabel('y');zlabel('z');  
title(sprintf('%s    %d x %d x %d',filename,NX,NY,NZ));