%     asg 12 1-D, 1-Phase pressure response in linear system
%     Constant flux in finite system; Carslaw and Jaeger
      clear all
      perm=50. ;
      por=0.3 ;
      vis=1.0 ;
      c=2.e-5 ;
      Lx=1000.;
      dy=100. ;
      dz=10.  ;
      b=1.    ;
      q=-10.  ;
      pic=1000. ;
      kap=0.00633*perm/(por*c*vis);
      F=5.615*q/(por*c*Lx*dy*dz);
      FK=887.2*q*vis/(dy*dz*perm*b);
%     profiles
      t=[1. 3. 10. 30.];
      x=linspace(0,1,10);
      x=Lx.*x;
      p=zeros(length(t),length(x));
 for i=1:length(t)
      pavg=pic+F*t(i);
      for j=1:length(x)
         steady=(3*x(j)^2-Lx^2)/(6*Lx^2);
         sum=0.;
         for n=1:40
            ex=-kap*(pi*n/Lx)^2*t(i);
            arg=pi*n*x(j)/Lx;
            sum=sum+(((-1)^n)/n^2)*exp(ex)*cos(arg);
         end
      p(i,j)=pavg+FK*Lx*(steady-(2*pi^(-2))*sum);
      end
 end
      figure(1)
      clf reset
      plot(x,p) % Analytical solution
      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);
            plot(x,p,'kx')  %Numerical solution
         end
      end
%     Pressure history
      figure(2)
      clf reset
      if NX>0 
          dx=Lx/NX;
      else
          dx=0;
      end
      clear t p
      t=linspace(0,50);
 for i=1:length(t)
      pavg=pic+F*t(i);
         steady=1/3;
         sum=0.;
         for n=1:40
            ex=-kap*(pi*n/Lx)^2*t(i);
            arg=pi*n;
            sum=sum+(((-1)^n)/n^2)*exp(ex)*cos(arg);
         end
      p(i)=pavg+FK*Lx*(steady-(2*pi^(-2))*sum);
 end
      plot(t,p) % Analytical solution
      title('Pressure History, x=Lx')
      xlabel('Time, days')
      ylabel('Pressure, psi')
      if (NX~=0)
         hold on
         load prod1.dat
         tfd=prod1(:,1);
         pfd=prod1(:,4);
         plot(tfd,pfd,'kx') % Numerical solution
      end
      figure(3)     % Plot square root of time
      clf reset
      sqt=zeros(length(t));
      sqt=sqrt(kap.*t.*Lx^(-2));
      plot(sqt,p)   % Analytical solution
      axis([0 1.5 600 1000])
      title('Pressure History, x=Lx')
      ylabel('Pressure, psi')
      xlabel('sqrt(kappa*t/Lx**2)')
      if (NX~=0)
         hold on
         sqtfd=sqrt(kap.*tfd.*Lx^(-2));
         plot(sqtfd,pfd,'kx')   % Numerical solution
      end