%     Contour plots of simulation & production history for 2-D IMPES
      clear all
      NX=input('NX = ?  ');
      ncases=input('ncases = ? ');
      M=input('[M(1), M(2) ...] = ?  ');
      LX=1000
      OOIP=42746.88
      SWR=0.2
      SOR=0.3
      MPV=OOIP*(1-SWR-SOR)/(1.-SWR)
      Q=20
      tcar=MPV/Q
      tim=[100 300 1000 3000]   ;
      vp=0:50:2000       ;
      vs=0:0.05:1        ;
     for ifile=1:ncases
      eval(['load prof' int2str(ifile) '.dat;']);
      x=eval(['prof' int2str(ifile) '(1:NX,1);']);
      y=x               ;
      p=zeros(length(x),length(y));
      figure(1)
      clf reset
      for i=1:4
         ist=1+NX*(i-1);
         ilst=NX*i;
         for j=1:length(y)
           eval(['p(:,j)=prof' int2str(ifile) '(ist:ilst,j+1);']);
         end
         subplot(2,2,i)
         contour(x,y,p,vp)
         title(['M = ',num2str(M(ifile)),', t = '...
               ,int2str(tim(i)),' days, ',num2str(tim(i)/tcar),' MPV'])
%        xlabel('Distance, feet')
%        ylabel('Distance, feet')
         axis('square')
         axis('equal')
         axis([0 1000 0 1000])
      end
      wait=input('Hit ENTER to continue');
     end
 
     for ifile=1:ncases
      figure(2)
      clf reset
       s=zeros(length(x),length(y));
      for i=1:4
         ist=1+NX*(i-1);
         ilst=NX*i;
         for j=1:length(y)
           eval(['s(:,j)=prof' int2str(ifile) '(ist:ilst,j+1+NX);']);
         end
         subplot(2,2,i)
         contour(x,y,s,vs)
         title(['M = ',num2str(M(ifile)),', t = '...
               ,int2str(tim(i)),' days ',num2str(tim(i)/tcar),' MPV'])
%        xlabel('Distance, feet')
%        ylabel('Distance, feet')
         axis('square')
         axis('equal')
         axis([0 1000 0 1000])
      end
      wait=input('Hit ENTER to continue');
    end

%      PLOT PRODUCTION HISTORY
      figure(3)
      clf reset
     for ifile=1:ncases
      clear prod   
      eval(['load prod' int2str(ifile) '.dat;'])
      eval(['prod=prod' int2str(ifile) ';']);
      time=prod(:,1)./tcar;
      Q1=prod(:,5);
      Q2=prod(:,6);
      cut=Q1./(Q1+Q2);
      plot(time,cut,'-')
      title(['NX= ',int2str(NX),' ,M = ' num2str(M)])
      xlabel('DIMENSIONLESS TIME')
      ylabel('WATER CUT')
      hold on
     end
      wait=input('Hit ENTER to continue');
      figure(4)
      clf reset
    for ifile=1:ncases
      clear prod   
      eval(['load prod' int2str(ifile) '.dat;'])
      eval(['prod=prod' int2str(ifile) ';']);
      time=prod(:,1)./tcar;
      er=-prod(:,7)./OOIP;
      plot(time,er,'-')
      title(['NX= ',int2str(NX),' ,M = ' num2str(M)])
      xlabel('DIMENSIONLESS TIME')
      ylabel('RECOVERY EFFICIENCY')
      hold on
    end
      wait=input('Hit ENTER to continue');
