%     asg7.1, 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 ;
      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);
      Lz=1.0  ;
      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));
      load alf.dat
%     profiles
      t=[1. 3. 10. 30. 50.];
      [m,n]=size(alf);
      rd=alf(1,2:n);
      x=rw*rd(:);
      alpha=alf(2:m,1);
      bess=alf(2:m,2:n);
      p=zeros(length(t),length(x));
 for i=1:length(t)
      td=kap*t(i)/rw^2;
      for j=1:length(x)
        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
      semilogx(x,p)
      axis([.2 2000. -100. 1000.])
      title('Pressure Profile, t=1, 3, 10, 30, 50 days')
      xlabel('Distance, feet')
      ylabel('Pressure, psi')
      NX=input('NX of F.D.?  0 if none. ')
      if (NX~=0)
         hold on
         clear x p
         load prof1.dat
         x=prof1(1:NX,1);
         for i=1:5
            ist=1+NX*(i-1);
            ilst=NX*i;
            p=prof1(ist:ilst,2);
            semilogx(x,p,'kx')
         end
      end
      wait=input('Hit ENTER to continue');
      clf reset
%     Pressure history
      r=rw
      clear t p prof1
      t=logspace(-2,2,100);
      p=zeros(size(t));
  for i=1:length(t)
      td=kap*t(i)/r^2;
        sum=0.;
        nsum=0;
        term=1.0;
        while (abs(term) > 1.e-7*abs(sum)  & nsum < m-1)
          nsum=nsum+1;
          term=exp(-alpha(nsum)^2*td)*bess(nsum,1);
          sum=sum+term;
        end
        p(i)=pic+FK*((2/(red^2-1))*(1/4+td)...
               -redf...
               +pi*sum);
  end
      plot(t,p)
      axis([0 50 -100 600])
      title('Pressure History, p=pw')
      xlabel('Time, days')
      ylabel('Pressure, psi')
      hold on
      load prof1.dat
      r1=prof1(1,1);
      load prod1.dat
      tfd=prod1(:,1);
      pfd=prod1(:,4)-FK*log(rw/r1);
      plot(tfd,pfd,'kx')
      wait=input('Hit ENTER to continue');
      clf reset
      semilogx(t,p)
      axis([0.001 100. -100 600])
      title('Pressure History, p=pw')
      xlabel('Time, days')
      ylabel('Pressure, psi')
      hold on
      semilogx(tfd,pfd,'kx')
      wait=input('Hit ENTER to continue');
      
      