function [R]=Quasi_1D_exit_bc(delt,delx,Cv,a,To,po,gamma,rhouI,rhouI_1,rhoI,rhoI_1,sI,sI_1) gamma_fac=(gamma-1)/(gamma+1) %Used coefficient uI=rhouI/rhoI; uI_1=rhouI_1/rhoI_1; %Get dependent variables %Temperature TI=To*(1-gamma_fac*uI^2/a^2); TI_1=To*(1-gamma_fac*uI_1^2/a^2); %Pressure pI=po*(1+gamma_fac*uI^2/a^2)^(gamma/(gamma-1)); pI_1=po*(1+gamma_fac*uI_1^2/a^2)^(gamma/(gamma-1)); %Speed of Sound cI=sqrt(gamma*pI/rhoI); cI_1=sqrt(gamma*pI_1/rhoI_1); %Average the variables u=(uI+uI_1)/2; p=(pI+pI_1)/2; c=(cI+cI_1)/2; rho=(rhoI+rhoI_1)/2; %Eigenvalues lam_1=u*delt/delx; lam_2=(u+c)*delt/delx; lam_3=(u-c)*delt/delx; %Residuals R1=(-lam_1)/(1+lam_1)*(rhoI-rhoI_1-(1/c^2)*(pI-pI_1)); R2=(-lam_2)/(1+lam_2)*(pI-pI_1+rho*c*(uI-uI_1))-... gamma*delt*p*u/(1+lam_2)*((sI-sI-1)/(sI*delx)); R3=(-lam_3)/(1+lam_3)*(pI-pI_1-rho*c*(uI-uI_1))-... gamma*delt*p*u/(1+lam_3)*((sI-sI-1)/(sI*delx)); %Use Residuals to obtain differences if (u/c) < 1 delp=(R2+R3)/2; else delp=0; end delrho=R1+delp/c^2; delu=(R2-delp)/(rho*c); pI=pI+delp; uI=uI+delu; rhoI=rhoI+delrho; TI=pI/((gamma-1)*Cv*rhoI); R=[R1,R2,R3];