clear all;
close all;
%
T=load('times.dat')*1000; %milliseconds
ntimes=size(T);
ntimes=ntimes(1);
lastrow=ntimes;
firstrow=1;
nxpts=930;
ltot=9.298;
ldriver=2.08;
deltax=ltot/(nxpts-1);
X=[0:deltax:ltot]-ldriver;
%
% Window locations
%
xwintop=8.597-ldriver;
xwinbot=8.813-ldriver;
xwininttop=8.208-ldriver;
xwinintbot=8.328-ldriver;
%
% Declare 2D array size
%
xtfield=zeros(ntimes,nxpts);
intloc=zeros(ntimes,nxpts);
diaphloc=zeros(ntimes,nxpts);
%
% Load data from xt.exe
% The x-t data file can be changed from
%     pressure to density- or any other variable
%
dd=load('pressure.dat');
m1=load('moles1.dat');
m2=load('moles2.dat');
m3=load('moles3.dat');
%
% Position data for plotting
%
for row=firstrow:lastrow;
  startinc=(row-1)*nxpts+1;
  xtfield(row,:)=log10(dd(startinc:startinc+nxpts-1)+1)';
  xint=zeros(nxpts,1);
  d1=m1(startinc:startinc+nxpts-1);
  d2=m2(startinc:startinc+nxpts-1);
  d3=m3(startinc:startinc+nxpts-1);
  xint=d1./(d1+d2+d3);
  intloc(row,:)=xint';
  xdia=zeros(nxpts,1);
  xdia=d3./(d1+d2+d3);
  diaphloc(row,:)=xdia';
end;
%
%  Color Plot
%
figure;
[x,t]=meshgrid(X,T);
contourf(x,t,xtfield,64);shading flat;
hold;
ylabel('Time (ms)');
xlabel('x (m)');
title('He - air - CO_2, M=2.5');
contour(x,t,intloc,[0.5 0.5],'k');
contour(x,t,diaphloc,[0.5 0.5],'k');
xwin1x=[xwintop xwintop];
xwin1y=[T(1) T(ntimes)]*1000;
xwin2x=[xwinbot xwinbot];
xwin2y=[T(1) T(ntimes)]*1000;
plot(xwin1x,xwin1y,'w:',xwin2x,xwin2y,'w:');
xwin3x=[xwininttop xwininttop];
xwin3y=[T(1) T(ntimes)]*1000;
xwin4x=[xwinintbot xwinintbot];
xwin4y=[T(1) T(ntimes)]*1000;
plot(xwin3x,xwin3y,'w:',xwin4x,xwin4y,'w:');
%
%  Line Plot
%
figure;
contour(x,t,xtfield,20,'k');
hold;
ylabel('Time (ms)');
xlabel('x (m)');
title('He - air - CO_2, M=2.5');
contour(x,t,intloc,[0.5 0.5],'r');
contour(x,t,diaphloc,[0.5 0.5],'r');
xwin1x=[xwintop xwintop];
xwin1y=[T(1) T(ntimes)]*1000;
xwin2x=[xwinbot xwinbot];
xwin2y=[T(1) T(ntimes)]*1000;
plot(xwin1x,xwin1y,'b:',xwin2x,xwin2y,'b:');
xwin3x=[xwininttop xwininttop];
xwin3y=[T(1) T(ntimes)]*1000;
xwin4x=[xwinintbot xwinintbot];
xwin4y=[T(1) T(ntimes)]*1000;
plot(xwin3x,xwin3y,'b:',xwin4x,xwin4y,'b:');