c ===================================================== subroutine exact(maxmx,meqn,mbc,mx,x,dx,t,q) c ===================================================== c c Copyright (C) 2002 Ralf Deiterding c Brandenburgische Universitaet Cottbus c implicit double precision (a-h,l,o-z) include "cuser.i" c dimension q(meqn,1-mbc:maxmx+mbc) dimension x(1-mbc:maxmx+mbc) c dimension y(1),c(24),w(1,9) common /CJ/ Dj,D,Vj,Pj,clambda,cfact external fcn1,fcn2 c Dj=dsqrt((gamma**2-1.d0)*q0/2.d0) & +dsqrt((gamma**2-1.d0)*q0/2.d0+gamma) D=dsqrt(f)*Dj Vj=gamma*(1.d0+D**2)/((gamma+1.d0)*D**2) Pj=(1.d0+D**2)/(gamma+1.d0) clambda=(D**2-gamma)**2/(2.d0*(gamma**2 & -1.d0)*q0*D**2) cfact=(D**2-gamma)/(1.d0+D**2) c n=1 tol=1.d-10 ind=1 nw=1 z=0.d0 zend=0.5d0 y(1)=0.d0 call dverk(n,fcn1,z,y,zend,tol,ind,c,nw,w) dhalbe=y(1) c y(1)=0.d0 lambda=0.d0 ssloc = sloc if (moving.eq.1) ssloc = sloc + D*t if (x(mx) .lt. ssloc) then z=0.d0 zend=ssloc-x(mx)-0.5d0*dx tol=1.d-8 ind=1 call dverk(n,fcn2,z,y,zend,tol,ind,c,nw,w) lambda = y(1) endif c mcount = 1 dmc = 1.d0/mcount do 150 i=mx,1,-1 if ( x(i) .gt. ssloc ) then q(1,i)=0.d0 q(2,i)=1.d0 q(3,i)=-D if (moving.eq.1) q(3,i) = 0.d0 q(4,i)=1.0d0/gamma1+q0+0.5d0*q(3,i)**2 else if ( lambda.lt.1.d0 .and. x(i)-0.5d0*dx.lt.ssloc ) then z=ssloc-x(i)+0.5d0*dx q(1,i) = 0.d0 q(2,i) = 0.d0 q(3,i) = 0.d0 q(4,i) = 0.d0 do m=1,mcount zend=z+dmc*dx tol=1.d-8 ind=1 call dverk(n,fcn2,z,y,zend,tol,ind,c,nw,w) lambda=y(1) a=cfact*dsqrt(1.d0-lambda/clambda) V=Vj*(1.d0-a/gamma) P=Pj*(1.d0+a) DU=-D*V if (moving.eq.1) DU = DU + D q(1,i)=q(1,i)+dmc*lambda/V q(2,i)=q(2,i)+dmc*(1.d0-lambda)/V q(3,i)=q(3,i)+dmc*DU/V q(4,i)=q(4,i)+dmc*(P/gamma1+(1.d0-lambda)/V*q0+ & 0.5d0*q(3,i)**2*V) z=zend enddo else a=cfact*dsqrt(1.d0-1.d0/clambda) V=Vj*(1.d0-a/gamma) P=Pj*(1.d0+a) DU=-D*V if (moving.eq.1) DU = DU + D q(1,i) = 1.d0/V q(2,i) = 0.d0 q(3,i) = DU/V q(4,i) = P/gamma1+q(2,i)*q0+0.5d0*q(3,i)**2*V endif 150 continue c return end c