%     asg 7.1, 1-Phase pressure response in infinite system
%     Line source well in infinite 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;
      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));
%     profiles
      t=[1. 3. 10. 30.];
      x=logspace(log10(rw),log10(re));
      p=zeros(length(t),length(x));
 for i=1:length(t)
      for j=1:length(x)
      arg=x(j)^2/(4*kap*t(i));
      p(i,j)=pic-FK*0.5*ei(-arg)  ;
      end
end
      semilogx(x,p)
      axis([.1 1000. 0. 1000.])
      title('Pressure Profile, t=1, 3, 10, 30 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:4
            ist=1+NX*(i-1);
            ilst=NX*i;
            p(:)=prof1(ist:ilst,2);
            semilogx(x,p,'wx')
         end
      end
      wait=input('Hit ENTER to continue');
      clg
%     Pressure history
      clear t p
      t=logspace(-3,2,100);
      p=zeros(size(t));
  for i=1:length(t)
%     pressure at X(1)
      load prof1.dat
      r=prof1(1,1)       ;
      arg=r^2/(4*kap*t(i));
      p(i)=pic-FK*0.5*ei(-arg)  ;
  end
      plot(t,p)
      axis([0 50 -100 600])
      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,'wx')
      wait=input('Hit ENTER to continue');
      clg
      semilogx(t,p)
      axis([0.001 100. -100 600])
      title('Pressure History, x=x(1,1)')
      xlabel('Time, days')
      ylabel('Pressure, psi')
      hold on
      semilogx(tfd,pfd,'wx')
      wait=input('Hit ENTER to continue');
      clg
end
