%     asg 6.4 & 6.5: 2-D, 1-Phase pressure response in finite system
%     Well in finite system, Matthews and Russell
%     Radial symmetry
      clear all
      clf reset
      perm=50. ;
      por=0.3 ;
      vis=1.0 ;
      c=2.e-5 ;
      Lx=1000.;
      Ly=1000.;
      Lz=1.0  ;
      rw=0.4  ;
      re=1128.38;
      red=re/rw;
      redf=(3*red^4-4*red^4*log(red)-2*red^2-1)...
            /(4*(red^2-1)^2);
      b=1.    ;
      q=-10.  ;
      q=4*q   ;
%     If rate is for one quadrant, then multiply by factor of 4
      pic=1000. ;
      kap=0.00633*perm/(por*c*vis);
      FK=887.2*q*vis/(perm*b*Lz)*(1/(2*pi));
%     profiles
      t=[1. 3. 10. 30. 50.];
      load alf.dat
      [m,n]=size(alf);
      rd=alf(1,2:n);
      r=zeros(size(rd));
      r(:)=rw*rd(:);
      alpha=alf(2:m,1);
      bess=alf(2:m,2:n);
      p=zeros(length(t),length(r));
 for i=1:length(t)
      td=kap*t(i)/rw^2;
      for j=1:length(r)
        sum=0.;
        nsum=0;
        term=1.0;
        while (abs(term) > 1.e-7*abs(sum)  & nsum < m)
          nsum=nsum+1;
          term=exp(-alpha(nsum)^2*td)*bess(nsum,j);
          sum=sum+term;
        end
        p(i,j)=pic+FK*((2/(red^2-1))*(rd(j)^2/4+td)...
               -(red^2*log(rd(j)))/(red^2-1)-redf...
               +pi*sum);
      end
end
      figure(1)
      clf reset
      semilogx(r,p)
      axis([10 2000. 0. 1000.])
      title('Pressure Profile, t=1, 3, 10, 30, 50 days')
      xlabel('Radial Distance, feet')
      ylabel('Pressure, psi')
      NX=input('NX of F.D.?  0 if none. ')
      if (NX>0)
         hold on
         clear r p
         load prof1.dat
         r=sqrt(2)*prof1(1:NX,1);   % Radial distance along diagonal
         for i=1:5
            ist=1+NX*(i-1);
            ilst=NX*i;
            p=diag(prof1(ist:ilst,2:(NX+1))); % Pressure along diagonal
            semilogx(r,p,'kx')
         end
      end
if(NX>0)
%     Contour plots of simulation
      figure(2)
      clf reset
      clear x p
      v=0:50:1000       ;
      tim=[1 3 10 30]   ;
      load prof1.dat
      x=prof1(1:NX,1);
      y=x               ;
      p=zeros(length(x),length(y));
      for i=1:4
         ist=1+NX*(i-1);
         ilst=NX*i;
         for j=1:length(y)
           p(:,j)=prof1(ist:ilst,j+1);
         end
         subplot(2,2,i)
         contour(x,y,p,v)
         title(['50 psi intv., t = '...
               ,int2str(tim(i)),' days'])
%        xlabel('Distance, feet')
%        ylabel('Distance, feet')
         axis('square')
         axis('equal')
         axis([0 1000 0 1000])
      end
      %     Pressure history
%     expression for quarter well in corner of grid block
      figure(3)
      clear t p
      clf reset
      t=logspace(-2,log10(50),100);
      pinf=zeros(size(t));
  for i=1:length(t)
      r=exp(-0.5)*(2*Lx/NX)/sqrt(pi);
      arg=r^2/(4*kap*t(i));
      pinf(i)=pic-FK*0.5*ei(-arg);   % Infinite reservoir
  end
      plot(t,pinf,'b-')
      hold on
%     Use image wells to represent a confined well in a square
      psq=zeros(size(t));
  for i=1:length(t)
      arg=(2*Lx)^2/(4*kap*t(i));
      arg2=(2*sqrt(2)*Lx)^2/(4*kap*t(i));
      arg3=(4*Lx)^2/(4*kap*t(i));
      arg4=(2*sqrt(5)*Lx)^2/(4*kap*t(i));
      arg5=(2*sqrt(5)*Lx)^2/(4*kap*t(i));
      arg6=(4*sqrt(2)*Lx)^2/(4*kap*t(i));
      psq(i)=pinf(i)-4*FK*0.5*ei(-arg)-4*FK*0.5*ei(-arg2)...
             -4*FK*0.5*ei(-arg3)-4*FK*0.5*ei(-arg4) ...
             -4*FK*0.5*ei(-arg5)-4*FK*0.5*ei(-arg6) ...
              ; % Well in square
  end
      plot(t,psq,'g-')
      % Finite radial reservoir
      clear t p rd prof1
      r=exp(-0.5)*(2*Lx/NX)/sqrt(pi);  
      rd=r/rw;
      t=logspace(-2,log10(50),100);
      p=zeros(size(t));
      j=14;
%     j should be between 13 and 14, to be exact
%     i.e., 151 < rd = 68/0.4 < 229
 for i=1:length(t)
      td=kap*t(i)/rw^2;
      sum=0.;
      nsum=0;
      term=1.0;
      while (abs(term) > 1.e-7*abs(sum)  & nsum < m)
        nsum=nsum+1;
        term=exp(-alpha(nsum)^2*td)*bess(nsum,j);
        sum=sum+term;
      end
      if (nsum < m)
        p(i)=pic+FK*((2/(red^2-1))*(rd^2/4+td)...
             -(red^2*log(rd))/(red^2-1)-redf...
             +pi*sum);
      else
        p(i)=pic
      end
end
      title('Pressure History, x=x(1,1)')
      xlabel('Time, days')
      ylabel('Pressure, psi')
      hold on
      load prod1.dat
      tfd=prod1(:,1);
      pfd=prod1(:,4);
      plot(tfd,pfd,'kx')
      plot(t,p,'r')
      
      figure(4) % Plot with log time
      clf reset
      semilogx(tfd,pfd,'kx')
      title('Pressure History, x=x(1,1)')
      xlabel('Time, days')
      ylabel('Pressure, psi')
      hold on
      semilogx(t,pinf,'b-')
      semilogx(t,psq,'g-')
      semilogx(t,p,'r-')
   	end