c
c
c =====================================================
subroutine von_neumann(E0,R0,rhoamb,pamb,dist,dim,
& gamma,toltheta,rho,p,u)
c =====================================================
c
implicit double precision (a-h,o-z)
c
C2=0.8510d0
if (dim.EQ.1.d0)then
C2=0.53864d0
end if
C1=(E0/C2/rhoamb)**(1.d0/(2.d0+dim))
time=(R0/C1)**((2.d0+dim)/2.d0)
tmp1 = C1*(time)**(2.d0/(2.d0+dim))
tmp2 = (gamma-1.d0)/(dim+2.d0*gamma-2.d0)
tmp3 = (gamma+1.d0)/(dim*gamma-dim+2.d0)
theta=1.
xwzlasttheta=0.
xwzlastx=0.
x=0.
iend=0
imax=1000
iblast=0
rho=rhoamb
p=pamb
u=0.d0
powx=-1.d0*((dim**2.d0+4.d0)*gamma**2.d0+
& (8.d0*dim-3.d0*dim**2.d0-4.d0)*gamma +
& (4.d0*dim**2.d0-8.d0*dim))/((dim+2.d0)*
& (2.d0*gamma+dim-2.)*(dim*gamma-dim+2.d0))
powr=((dim**2.d0+4.d0)*gamma**2.d0+
& (8.d0*dim-3.d0*dim**2.d0-4.d0)*gamma +
& (4.d0*dim**2.d0-8.d0*dim))/((2.d0-gamma)*
& (2.d0*gamma+dim-2.)*(dim*gamma-dim+2.d0))
powu=-1.d0*((dim**2.d0+4.d0)*gamma**2.d0+
& (8.d0*dim-3.d0*dim**2.d0-4.d0)*gamma +
& (4.d0*dim**2.d0-8.d0*dim))/((dim+2.d0)*
& (2.d0*gamma+dim-2.)*(dim*gamma-dim+2.d0))
powp=((dim**2.d0+4.d0)*gamma**2.d0+
& (8.d0*dim-3.d0*dim**2.d0-4.d0)*gamma +
& (4.d0*dim**2.d0-8.d0*dim))/((dim+2.d0)*
& (2.d0-gamma)*(dim*gamma-dim+2.d0))
base=1.
do i =0,imax
base=((dim*(2.d0-gamma)*theta+dim+2.d0*gamma-2.d0)/
& (3.d0*dim-dim*gamma+2.d0*gamma-2.d0))
x=tmp1*
& (theta)**(tmp2)*
& ((theta+1.d0)/2.d0)**(-2.d0/(dim+2.d0))*
& ((gamma+theta)/(gamma+1.d0))**
& (tmp3)*
& base**(powx)
if(i.eq.0.and.dist.gt.x)then
c we are out of the pressurized range
iblast=0
goto 10
endif
if(dabs(x-dist)/dist.lt.toltheta)then
c we found the position
iblast=1
goto 10
endif
if(i.eq.imax-1)then
c we did not found it
print*, ' von-neumann initialization at distance ', dist,
& ' failed \n'
endif
dxdtheta=tmp1*(tmp2)*(theta)**(tmp2-1.d0)*
& ((theta+1.d0)/2.d0)**(-2.d0/(dim+2.d0))*
& ((gamma+theta)/(gamma+1.d0))**(tmp3)*base**(powx)+
& tmp1*(theta)**(tmp2)*(-1.d0/(dim+2.d0))*
& ((theta+1.d0)/2.d0)**(-2.d0/(dim+2.d0)-1.d0)*
& ((gamma+theta)/(gamma+1.d0))**(tmp3)*base**(powx)+
& tmp1*(theta)**(tmp2)*
& ((theta+1.d0)/2.d0)**(-2.d0/(dim+2.d0))*
& (1.d0/(dim*gamma-dim+2.d0))*
& ((gamma+theta)/(gamma+1.d0))**(tmp3-1.d0)*base**(powx)+
& tmp1*(theta)**(tmp2)*
& ((theta+1.d0)/2.d0)**(-2.d0/(dim+2.d0))*
& ((gamma+theta)/(gamma+1.d0))**(tmp3)*
& dim*(2.d0-gamma)/(3.d0*dim-dim*gamma+2.d0*gamma-2.d0)*
& powx*base**(powx-1.d0)
xwzlasttheta=theta
xwzlastx=x
theta=xwzlasttheta+(dist-xwzlastx)/dxdtheta
if(theta.le.0.d0)then
theta=xwzlasttheta/2.d0
elseif(theta.ge.1.d0)then
theta=1.1d0*xwzlasttheta
if(theta.ge.1.d0)then
theta=0.99
endif
endif
end do
10 if(iblast.eq.1)then
rho=(gamma+1.d0)/(gamma-1.d0)*rhoamb*
& theta**((dim)/(dim+2.d0*gamma-2.d0))*
& ((gamma+theta)/(gamma+1.d0))**((-2.d0*(dim-1.d0))/
& (dim*gamma-dim+2.d0))*
& (base)**(powr)
u=4.d0*C1/(2.d0+dim)/(gamma+1.d0)*
& time**(-1.d0*dim/(dim+2.d0))*
& theta**((gamma-1.d0)/(dim+2.d0*gamma-2.d0))*
& ((theta+1.d0)/2.d0)**(dim/(dim+2.d0))*
& ((gamma+theta)/(gamma+1.d0))**
& (-1.d0*(dim-1.d0)*(gamma-1.d0)/(dim*gamma-dim+2.d0))*
& base**(powu)
p=8.d0/(dim+2.d0)**2.d0/(gamma+1.d0)*rhoamb*C1*C1*
& time**(-2.d0*dim/(dim+2.d0))*
& ((theta+1.d0)/2.d0)**(2.d0*dim/(dim+2.d0))*
& ((gamma+theta)/(gamma+1.d0))**
& (-2.d0*(dim-1.d0)*gamma/(dim*gamma-gamma+2.d0))*
& base**powp
endif
return
end