% build matrix for pressure equation, based on the energy equation
% in fully conservative formulation
% with Mach uniform scaling 
% only the case ux>=

rhocent=0.5*([rholeft;rhonew]+[rhonew;rhoright]);
uxc=0.5*(uxnew(1:J)+uxnew(2:J+1));

% Apply a neumann bc at j=1

j=1;
  P(j,3)     =  Mref^2*(gamma-1)*( ...
               1/8*rhonew(j)*(-2*lambda*mxnew(j+1)/rhocent(j+1)^2 ...
               -2*lambda*mxnew(j)/(rhocent(j+1)*rhocent(j))) ...
	       -lambda^2/rhocent(j+1)*(gamma/(gamma-1)*pold(j)+1/2* ...
                rhonew(j)*uxc(j)^2)) ...
		-lambda^2/rhocent(j+1);
  P2ex =    	Mref^2*(gamma-1)*( ... 
               1/8*rhonew(j)*(2*lambda*mxnew(j+1)/(rhocent(j+1)*rhocent(j)) ...
                +2*lambda*mxnew(j)/rhocent(j)^2) ...
		-lambda^2/rhocent(j)*(gamma/(gamma-1)*pleft+1/2* ...
                rholeft*uxleft^2)-lambda*gamma/(gamma-1)*uxnew(j)) ...
		-lambda^2/rhocent(j);	
  P(j,2)  =  	 Mref^2*(gamma-1)*( ... 
               1/(gamma-1)+1/8*rhonew(j)*(2*lambda*mxnew(j+1)/rhocent(j+1)^2 ...
                -2*lambda*(mxnew(j+1)-mxnew(j))/(rhocent(j+1)*rhocent(j))-...
		2*lambda*mxnew(j)/rhocent(j)^2)+ ...
	       lambda^2*( ...
              (gamma/(gamma-1)*pold(j)+1/2*rhonew(j)*uxc(j)^2)/rhocent(j+1) ...
             +(gamma/(gamma-1)*pleft+1/2*rholeft*uxleft^2)/rhocent(j)) ...
             +lambda*gamma/(gamma-1)*uxnew(j+1)) ...
	     +lambda^2/rhocent(j)+lambda^2/rhocent(j+1)+P2ex;
for j=2:J-1
  P(j+1,3)=	Mref^2*(gamma-1)*( ...
               1/8*rhonew(j)*(-2*lambda*mxnew(j+1)/rhocent(j+1)^2 ...
               -2*lambda*mxnew(j)/(rhocent(j+1)*rhocent(j))) ...
	       -lambda^2/rhocent(j+1)*(gamma/(gamma-1)*pold(j)+1/2* ...
                rhonew(j)*uxc(j)^2)) ...
		-lambda^2/rhocent(j+1);	
  P(j,2)=	   Mref^2*(gamma-1)*( ... 
               1/(gamma-1)+1/8*rhonew(j)*(2*lambda*mxnew(j+1)/rhocent(j+1)^2 ...
                -2*lambda*(mxnew(j+1)-mxnew(j))/(rhocent(j+1)*rhocent(j))-...
		2*lambda*mxnew(j)/rhocent(j)^2)+ ...
	       lambda^2*( ...
        (gamma/(gamma-1)*pold(j)+1/2*rhonew(j)*uxc(j)^2)/rhocent(j+1) ...
       +(gamma/(gamma-1)*pold(j-1)+1/2*rhonew(j-1)*uxc(j-1)^2)/rhocent(j)) ...
          +lambda*gamma/(gamma-1)*uxnew(j+1)) ...
	  +lambda^2/rhocent(j)+lambda^2/rhocent(j+1);
  P(j-1,1)=    	Mref^2*(gamma-1)*( ... 
               1/8*rhonew(j)*(2*lambda*mxnew(j+1)/(rhocent(j+1)*rhocent(j)) ...
                +2*lambda*mxnew(j)/rhocent(j)^2) ...
		-lambda^2/rhocent(j)*(gamma/(gamma-1)*pold(j-1)+1/2* ...
                rhonew(j-1)*uxc(j-1)^2)-lambda*gamma/(gamma-1)*uxnew(j)) ...
		-lambda^2/rhocent(j);
end

% apply a Dirichlet bc at j=J

j=J;
  P(j,2)=1;
P2rhs=zeros(size(pold));
j=1;
P2rhs(j)  = Mref^2*(gamma-1)*(...
              1/8*rhoold(j)*(uxold(j)+uxold(j+1))^2 ...
             -1/8*rhonew(j)*(uxnew(j)+uxnew(j+1))^2-lambda*( ...
          gamma/(gamma-1)* ...
	 (pold(j)*mxnew(j+1)/rhocent(j+1)-pleft*mxnew(j)/rhocent(j))...
	 +1/2*(rhonew(j)*uxc(j)^2*mxnew(j+1)/rhocent(j+1) - ...
	         rholeft*uxleft^2*mxnew(j)/rhocent(j)))) ...
		 -lambda*(uxnew(j+1)-uxnew(j));

for j=2:J-1
  P2rhs(j)  = Mref^2*(gamma-1)*(...
              1/8*rhoold(j)*(uxold(j)+uxold(j+1))^2 ...
             -1/8*rhonew(j)*(uxnew(j)+uxnew(j+1))^2-lambda*( ...
          gamma/(gamma-1)* ...
	 (pold(j)*mxnew(j+1)/rhocent(j+1)-pold(j-1)*mxnew(j)/rhocent(j))...
	 +1/2*(rhonew(j)*uxc(j)^2*mxnew(j+1)/rhocent(j+1) - ...
	         rhonew(j-1)*uxc(j-1)^2*mxnew(j)/rhocent(j)))) ...
		 -lambda*(uxnew(j+1)-uxnew(j));
end 

% build pressure matrix

  Pmat=spdiags(P,[-1 0 1],J,J);
  
% solver pressure correction equation
  
  dp=(Pmat)\(P2rhs);


