function plt2(pfix)
%
% close all; clear; clc;
% Vykresleni vysledku pro "shock-tube problem"
%
% vstupni data:
%
% pfix = 'ini' ==> pocatecni data (vytrorena programem meshfree.exe)
%         'f'   ==> vypoctena data (vytrorena programem meshfree.exe)
%
% 
% vzorove volani:
%
% plt2('f')
%
% pouziva:
%
% program e1rpex.exe, odpovidajici vstupy jsou e1rpex11.ini
%% soubor obsahujici souradnice a rychlosti (fiktivni_castice)
%fid=fopen([pfix '_xv.txt'],'rt');
fid=fopen('f_xv.txt','rt');
aa=fscanf(fid,'%i %i %g',3);
ntotal=aa(1); % celkovy pocet uzlu
dim=aa(2);    % dimenze
time=aa(3);   % cas
x=zeros(ntotal,dim);
vx=zeros(ntotal,dim);
for i=1:ntotal
  ii=fscanf(fid,'%i',1);
  for j=1:dim
    x(i,j)=fscanf(fid,'%g',1);
  end
  for j=1:dim
    vx(i,j)=fscanf(fid,'%g',1);
  end
end
fclose(fid);
%% soubor obsahujici souradnice a rychlosti (virtualni_castice)
%fid=fopen([pfix '_xv.txt'],'rt');
fid=fopen('f_xv_v.txt','rt');
aa=fscanf(fid,'%i %i %g',3);
ntotal=aa(1); % celkovy pocet uzlu
dim=aa(2);    % dimenze
time=aa(3);   % cas
x_v=zeros(ntotal,dim);
vx_v=zeros(ntotal,dim);
for i=1:ntotal
  ii=fscanf(fid,'%i',1);
  for j=1:dim
    x_v(i,j)=fscanf(fid,'%g',1);
  end
  for j=1:dim
    vx_v(i,j)=fscanf(fid,'%g',1);
  end
end
fclose(fid);
%% soubor obsahujici polohu fiktivnich castic
fid_fp=fopen('xv_vp.txt','rt');
bb=fscanf(fid_fp,'%i',1);
ntotal_vp=bb(1); % celkovy pocet uzlu
x_vp=zeros(ntotal_vp,dim); % dim = 2
vx_vp=zeros(ntotal,dim);
for i=1:ntotal_vp
  ii=fscanf(fid_fp,'%i',1);
  for j=1:dim
    x_vp(i,j)=fscanf(fid_fp,'%g',1);
  end
  for j=1:dim
    vx_vp(i,j)=fscanf(fid_fp,'%g',1);
  end
end
fclose(fid_fp);
%% soubor obsahujici polohu virtualnich castic
fid_vp=fopen('xv_vp_v.txt','rt');
bb=fscanf(fid_vp,'%i',1);
ntotal_vp=bb(1); % celkovy pocet uzlu
x_vp_v=zeros(ntotal_vp,dim); % dim = 2
vx_vp_v=zeros(ntotal,dim);
for i=1:ntotal_vp
  ii=fscanf(fid_vp,'%i',1);
  for j=1:dim
    x_vp_v(i,j)=fscanf(fid_vp,'%g',1);
  end
  for j=1:dim
    vx_vp_v(i,j)=fscanf(fid,'%g',1);
  end
end
fclose(fid_vp);
%% soubor obsahujici hmotnosti, hustotu, tlak a vnitrni energii (fikt. cast.)
% fid=fopen([pfix '_state.txt'],'rt');
fid=fopen('f_state.txt','rt');
mass=zeros(ntotal,1);
p=zeros(ntotal,1);
rho=zeros(ntotal,1);
u=zeros(ntotal,1);
for i=1:ntotal
  ii=fscanf(fid,'%i',1);
  mass(i)=fscanf(fid,'%g',1);
  rho(i)=fscanf(fid,'%g',1);
  p(i)=fscanf(fid,'%g',1);
  u(i)=fscanf(fid,'%g',1);
end
fclose(fid);
%% soubor obsahujici hmotnosti, hustotu, tlak a vnitrni energii (virt. cast.)
% fid=fopen([pfix '_state.txt'],'rt');
fid=fopen('f_state_v.txt','rt');
mass_v=zeros(ntotal,1);
p_v=zeros(ntotal,1);
rho_v=zeros(ntotal,1);
u_v=zeros(ntotal,1);
for i=1:ntotal
  ii=fscanf(fid,'%i',1);
  mass_v(i)=fscanf(fid,'%g',1);
  rho_v(i)=fscanf(fid,'%g',1);
  p_v(i)=fscanf(fid,'%g',1);
  u_v(i)=fscanf(fid,'%g',1);
end
fclose(fid);
%% Presne reseni
sre=true; % ano
if  pfix=='ini', sre=false; end; % pro pocatecni data ne
sr=sre;
priklad=11;
priklad=sprintf('%02d',priklad);
while sr
  sr=0;
% cteni vstupniho souboru
  fid=fopen(['e1rpex',priklad,'.ini'],'rt');
  if fid < 0
    fprintf(sprintf('Warning: soubor %s neexistuje',['e1rpex',priklad,'.ini']));
    sr=0; srovnej=0;
    break;
  end;
  tline=fgetl(fid); DOMLEN=sscanf(tline,'%e',1);
  tline=fgetl(fid); DIAPHI=sscanf(tline,'%e',1);
  tline=fgetl(fid); CELLS=sscanf(tline,'%d',1);
  tline=fgetl(fid); GAMMA=sscanf(tline,'%e',1);
  tline=fgetl(fid); TIMEOU=sscanf(tline,'%e',1);
  tline=fgetl(fid); DL=sscanf(tline,'%e',1);
  tline=fgetl(fid); UL=sscanf(tline,'%e',1);
  tline=fgetl(fid); PL=sscanf(tline,'%e',1);
  tline=fgetl(fid); DR=sscanf(tline,'%e',1);
  tline=fgetl(fid); UR=sscanf(tline,'%e',1);
  tline=fgetl(fid); PR=sscanf(tline,'%e',1);
  tline=fgetl(fid); PSCALE=sscanf(tline,'%e',1);
  fclose(fid);
  
  fprintf('delka intervalu : %10.2f\n',DOMLEN);
  fprintf('pocet bunek     : %10d\n',CELLS);
  fprintf('cas             : %10.4f\n',TIMEOU);
  fprintf('hustota vlevo   : %10.4f    ',DL);fprintf('hustota vpravo  : %10.4f\n',DR);
  fprintf('rychlost vlevo  : %10.2f    ',UL);fprintf('rychlost vpravo : %10.2f\n',UR);
  fprintf('tlak vlevo      : %10.2f    ',PL);fprintf('tlak vpravo     : %10.2f\n',PR);
% vypocet
  fid=fopen('pr.txt','wt');
  fprintf(fid,'%2s',priklad);
  fclose(fid);  fprintf('E1RPEX is running ...\n');
  !e1rpex
  fprintf(['E1RPEX done\n',...
   '********************************************************\n'] );
  fid=fopen(['e1rpex',priklad,'.out'],'rt');
  if fid < 0
    error(sprintf('soubor %s neexistuje',['e1rpex',priklad,'.out']));
  end;
  XE=fscanf(fid,'%e',[5,CELLS]);
  fclose(fid);
end;  
%% priprava kresby
rx=0.1; ry=0.1; rd=0.1; rv=0.1; rp=0.1; ru=0.1;
xmin=min(x(:,1));xmax=max(x(:,1)); dx=xmax-xmin;
ymin=min(x(:,2));ymax=max(x(:,2)); dy=ymax-ymin;
xmin_vp=min(x_vp(:,1));xmax_vp=max(x_vp(:,1)); dx_vp=xmax-xmin;
close all;
%% castice
F1=figure(1);
clf(1);
%line([0 1e-3 1e-3 0 0],[0 0 1e-3 1e-3 0],'color','b','linewidth',2);
line(x_vp(:,1),x_vp(:,2),'linestyle','none',...
     'markeredgecolor','g', 'markerfacecolor','white',...
     'marker','o','markersize',2);
line(x(1:ntotal/2,1),x(1:ntotal/2,2),'linestyle','none',...
     'markeredgecolor','r', 'markerfacecolor','r',...
     'marker','o','markersize',2);
 line(x(ntotal/2+1:ntotal,1),x(ntotal/2+1:ntotal,2),'linestyle','none',...
     'markeredgecolor','b', 'markerfacecolor','w',...
     'marker','o','markersize',2);
 
%title(sprintf('shock tube 2D, cas = %6.2f',time));
%axis([xmin-rx*dx,xmax+rx*dx,ymin-ry*dy,ymax+ry*dy]);
axis([xmin_vp-rx*dx_vp,xmax_vp+rx*dx_vp,ymin-ry*dy,ymax+ry*dy])
%axis([-1e-4 (1e-3+1e-4) -1e-4 (1e-3+1e-4)]);
%axis([-0.45 0.45 -0.05 0.05]);
%grid on;
x_label=xlabel('$x\ [m]$');
y_label=ylabel('$y\ [m]$');
set(x_label,'interpreter','latex')
set(y_label,'interpreter','latex')
 axis equal
% axis manual
%% rychlost
F2=figure(2);
clf(2);
q=quiver(x(:,1),x(:,2),vx(:,1),vx(:,2),1,'LineStyle','-',...
   'Color','b','Marker','none','ShowArrowHead','on');
line(x_vp(:,1),x_vp(:,2),'linestyle','none',...
     'markeredgecolor','black', 'markerfacecolor','black',...
     'marker','o','markersize',1);
% line([0 1e-3 1e-3 0 0],[0 0 1e-3 1e-3 0],'Color','b','LineWidth',2);
ii=20; idx2=(ii-1)*40+1:ii*40; 
jj=20; idx1=jj:40:39*40+jj;
% line(x(idx2,1),x(idx2,2),'linestyle','none',...
%      'markeredgecolor','r', 'markerfacecolor','r',...
%      'marker','o','markersize',3);  
% line(x(idx1,1),x(idx1,2),'linestyle','none',...
%      'markeredgecolor','g', 'markerfacecolor','r',...
%      'marker','o','markersize',3);  

%title(sprintf('rychlostni pole, cas = %6.2f',time));

axis([xmin-rx*dx,xmax+rx*dx,ymin-ry*dy,ymax+ry*dy]);
axis equal;
%axis([-1e-4 (1e-3+1e-4) -1e-4 (1e-3+1e-4)]);

%% rez horizontalni
F4=figure(4);
 clf(4); 
 
 %yi=ones(1,40)*(0.005);
 yi=ones(1,40)*(-0.002);
 %yi=ones(1,40)*0; %polovina
 %yi=ones(1,40)*0.01; % na okraji
 %yi=ones(1,160)*2e-3/160;
 %xi=2e-3/160:2e-3/80:4e-3-2e-3/160;
 xi=-0.3995:0.02:0.3995;  % toto zmenit
 zi=griddata(x(:,1),x(:,2),vx(:,1),xi,yi);
 zi_v=griddata(x_v(:,1),x_v(:,2),vx_v(:,1),xi,yi);
 
 dmin=min(zi);dmax=max(zi);
 if sre
   plot(XE(1,:)-DIAPHI,XE(3,:),'k-',xi,zi_v,'.b',xi,zi,'.r');
   dmin=min(dmin,min(XE(3,:))); dmax=max(dmax,max(XE(3,:)));
  % title(sprintf('Eulerova rovnice, rychlost, cas = %6.2f',TIMEOU));
 else
  plot(x,vx,'-b',x,vx,'.r');
 % title(sprintf('Eulerova rovnice,  rychlost, cas =  %8.4g',time));
end;
dd=dmax-dmin;
if dd~=0
  axis([-0.4,0.4,dmin-rv*dd,dmax+rv*dd]);
end
grid on;
x_label=xlabel('$x\ [m]$');
y_label=ylabel('$v\ [m/\!s]$');
set(x_label,'interpreter','latex')
set(y_label,'interpreter','latex')
leg = legend('presne reseni','SPH + virt. castice','SPH + fikt. castice','Location','Best');
 
 %plot(xi,zi,'r.');
 
% title(sprintf('rychlost, horizontalni rez, cas = %6.2f',time));
 
 F5=figure(5);
 clf(5);
 dmin=min(rho);dmax=max(rho);
  zi=griddata(x(:,1),x(:,2),rho(:,1),xi,yi);
  zi_v=griddata(x_v(:,1),x_v(:,2),rho_v(:,1),xi,yi);
  if sre
   plot(XE(1,:)-DIAPHI,XE(2,:),'k-',xi,zi_v,'.b',xi,zi,'.r');
   dmin=min(dmin,min(XE(2,:))); dmax=max(dmax,max(XE(2,:)));
   
  % title(sprintf('Eulerova rovnice, hustota, cas = %6.2f',TIMEOU));
   
  else
   plot(xi,zi,'r.');
 end;
 dd=dmax-dmin;
if dd~=0
  axis([xmin-rx*dx,xmax+rx*dx,dmin-rd*dd,dmax+rd*dd]);
  axis([-0.4,0.4,dmin-rd*dd,dmax+rd*dd]);
end
 grid on;
 x_label=xlabel('$x\ [m]$');
 y_label=ylabel('$\rho\ [kg/\!m^3] $');
 set(x_label,'interpreter','latex')
 set(y_label,'interpreter','latex')
 leg = legend('presne reseni','SPH + virt. castice','SPH + fikt. castice');
 
% title(sprintf('hustota, horizontalni rez, cas = %6.2f',time));
 
  F6=figure(6);
 clf(6); 
  zi=griddata(x(:,1),x(:,2),p(:,1),xi,yi);
  zi_v=griddata(x_v(:,1),x_v(:,2),p_v(:,1),xi,yi);
  
  
  dmin=min(zi);dmax=max(zi);
 if sre
   plot(XE(1,:)-DIAPHI,XE(4,:),'k-',xi,zi_v,'.b',xi,zi,'.r');
   dmin=min(dmin,min(XE(4,:))); dmax=max(dmax,max(XE(4,:)));
%   title(sprintf('Eulerova rovnice, tlak, cas = %6.2f',TIMEOU));
 else
  plot(x,p,'-b',x,p,'.r');
%  title(sprintf('Eulerova rovnice,  tlak, cas =  %8.4g',time));
end;
dd=dmax-dmin;
if dd~=0
  axis([xmin-rx*dx,xmax+rx*dx,dmin-rv*dd,dmax+rv*dd]);
  axis([-0.4,0.4,dmin-rv*dd,dmax+rv*dd]);
end
grid on;
x_label=xlabel('$x\ [m]$');
y_label=ylabel('$p\ [N/\!m^2]$');
set(x_label,'interpreter','latex')
set(y_label,'interpreter','latex')
leg = legend('presne reseni','SPH + virt. castice','SPH + fikt. castice');
 %plot(xi,zi,'r.');
 %axis square;
 
 % title(sprintf('tlak, horizontalni rez, cas = %6.2f',time));
 
   F7=figure(7);
 clf(7); 
  zi=griddata(x(:,1),x(:,2),u(:,1),xi,yi);
  zi_v=griddata(x_v(:,1),x_v(:,2),u_v(:,1),xi,yi); 
  
  dmin=min(zi);dmax=max(zi);
 if sre
   plot(XE(1,:)-DIAPHI,XE(5,:),'k-',xi,zi_v,'.b',xi,zi,'.r');
   dmin=min(dmin,min(XE(5,:))); dmax=max(dmax,max(XE(5,:)));
 %  title(sprintf('Eulerova rovnice, vnitrni energie, cas = %6.2f',TIMEOU));
 else
  plot(x,u,'-b',x,u,'.r');
 % title(sprintf('Eulerova rovnice,  vnitrni energie, cas =  %8.4g',time));
end;
dd=dmax-dmin;
if dd~=0
  axis([xmin-rx*dx,xmax+rx*dx,dmin-ru*dd,dmax+ru*dd]);
  axis([-0.4,0.4,dmin-ru*dd,dmax+ru*dd]);
end
grid on;
x_label=xlabel('$x\ [m]$');
y_label=ylabel('$e\ [J/\!kg]$');  
set(x_label,'interpreter','latex')
set(y_label,'interpreter','latex')
leg = legend('presne reseni','SPH + virt. castice','SPH + fikt. castice');
%set(leg,'interpreter','latex')
%plot(xi,zi,'r.');

 % title(sprintf('vnitrni energie, horizontalni rez, cas = %6.2f',time));

