c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
c
c =====================================================
subroutine ic(maxmx,meqn,mbc,mx,x,dx,q)
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
V = (Vj*(1.d0-cfact/gamma))
P = Pj*(1.d0+cfact)
write (6,*) 'D: ',D,1/V,-D*V,P
c
y(1)=0.d0
lambda=0.d0
c
c a=cfact*dsqrt(1.d0-lambda/clambda)
c P=Pj*(1.d0+a)
c xcnt = 0.d0
c open(unit=11, status='unknown', form='formatted',
c & file='p.dat')
c do i=0,20
c write (11,*) xcnt,P
c xcnt = xcnt + 10.d0
c enddo
c close (11)
c
if (x(mx) .lt. sloc) then
z=0.d0
zend=sloc-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. sloc ) 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.sloc ) then
z=sloc-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
c============================================================
subroutine fcn1(n,x,y,yprime)
c============================================================
implicit double precision (a-h,l,o-z)
include "cuser.i"
c
dimension y(n),yprime(n)
common /CJ/ Dj,D,Vj,Pj,clambda,cfact
c
a=cfact*dsqrt(1.d0-x/clambda)
P=Pj*(1.d0+a)
V=Vj*(1.d0-a/gamma)
DU=D*V
yprime(1)=DU*exp(E0/(P*V))/(1.d0-x)
return
end
c
c==========================================================
subroutine fcn2(n,x,y,yprime)
c==========================================================
implicit double precision (a-h,l,o-z)
include "cuser.i"
c
dimension y(n),yprime(n)
common /CJ/ Dj,D,Vj,Pj,clambda,cfact
c
a=cfact*dsqrt(1.d0-y(1)/clambda)
P=Pj*(1.d0+a)
V=Vj*(1.d0-a/gamma)
DU=D*V
yprime(1)=dhalbe*(1.d0-y(1))*exp(-E0/(P*V))/DU
c
return
end