%     asg 12 1-D, 1-Phase pressure response in linear system
%     Constant flux in finite system; Carslaw and Jaeger
      clear all
      clf reset
      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
      plot(x,p)
      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,'rx')
         end
      end
   