
ncol=64;
pause on;
clg;
subplot(111);

disp(' ')
disp('  ----------------------- ')
disp(' ')
disp('  Basic snaphot plotting ')
disp(' ')
disp('  ----------------------- ')
disp(' ')
disp(' ')

file=input(' Give name of run : ');

hfile=sprintf('%s_hd',file);
eval( [' load ' hfile ]);
head=eval([hfile]);

nx=head(3)+1;
ny=head(4)+1;
nz=head(5)+1;

disp(' Which plane to visualise ? ')
disp(' 1) X-const ? ')
disp(' 2) Y-const ? ')
disp(' 3) Z-const ? ')
pans=input(' Which ? ');

if pans==1,
file1=sprintf('%s_Xxs',file);
file2=sprintf('%s_Xys',file);
file3=sprintf('%s_Xzs',file);
elseif pans==2,
file1=sprintf('%s_Yxs',file); 
file2=sprintf('%s_Yys',file);
file3=sprintf('%s_Yzs',file);
elseif pans==3,
file1=sprintf('%s_Zxs',file);  
file2=sprintf('%s_Zys',file); 
file3=sprintf('%s_Zzs',file); 
end

disp(' Loading snaphot files ...');
disp(' x ...');eval(['load -ascii ' file1]); snap1=eval([file1]);
disp(' y ...');eval(['load -ascii ' file2]); snap2=eval([file2]);
disp(' z ...');eval(['load -ascii ' file3]); snap3=eval([file3]);

if pans==1, ns=floor(length(snap1)/(ny*nz));n1=ny;n2=nz; end
if pans==2, ns=floor(length(snap1)/(nx*nz));n1=nx;n2=nz; end
if pans==3, ns=floor(length(snap1)/(nx*ny));n1=nx;n2=ny; end

if pans==1, filec1=sprintf('%s_yco',file);filec2=sprintf('%s_zco',file);end
if pans==2, filec1=sprintf('%s_xco',file);filec2=sprintf('%s_zco',file);end
if pans==3, filec1=sprintf('%s_xco',file);filec2=sprintf('%s_yco',file);end

eval(['load -ascii ' filec1]); x1=eval([filec1]);
eval(['load -ascii ' filec2]); x2=eval([filec2]);

in1=100;
in2=100;

new1=min(x1)+(0:in1-1)*(max(x1)-min(x1))/(in1-1);
new2=min(x2)+(0:in2-1)*(max(x2)-min(x2))/(in2-1);

disp(' ')
disp(sprintf(' %g snapshots will be displayed ... ',ns))
disp(' ')

disp(' Which colors for snapshots ')
disp(' 1) gray ')
disp(' 2) hsv')
disp(' 3) cool')
disp(' 4) hot')
disp(' 5) jet')
disp(' 6) france')
cans=input(' Which? ');

j0=input(' Give first snap to show : ');
jinc=input(' Give snap increment : ');

isnap=0;

for j=j0:jinc:ns;

isnap=isnap+1;

s1=reshape(snap1((j-1)*n2*n1+1:j*n2*n1),n2,n1);
s2=reshape(snap2((j-1)*n2*n1+1:j*n2*n1),n2,n1);
s3=reshape(snap3((j-1)*n2*n1+1:j*n2*n1),n2,n1);

is1=interp2(x1',x2,s1,new1',new2);
is2=interp2(x1',x2,s2,new1',new2);
is3=interp2(x1',x2,s3,new1',new2);

rmax=max([ max(max(is1)) max(max(is2)) max(max(is3)) ]);

if rmax ~=0,

is1=32+is1/rmax*32;
is2=32+is2/rmax*32;
is3=32+is3/rmax*32;

subplot(131),image(1:n1,1:n2,is1),axis('equal'),axis('off'),title(file1)
if cans==1, colormap('gray'); end
if cans==2, colormap('hsv'); end
if cans==3, colormap('cool'); end
if cans==4, colormap('hot'); end
if cans==5, colormap('jet'); end
if cans==6, mycolormap('france33'); end
subplot(132),image(1:n1,1:n2,is2),axis('equal'),axis('off'),title(file2) 
subplot(133),image(1:n1,1:n2,is3),axis('equal'),axis('off'),title(file3)
disp(' Press any key to continue ')
pause

end

end

