c
c
program nonplanar
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c stability analysis for strong, slow, parallel, ideal MHD shocks
c
c Ref: Edelman, M.A. 1989, Astrofizika, 31 #2, 4071
c
c Version: Janueary 9, 2019
c
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
parameter (neq=11,i1=50,i2=50)
dimension y(neq),yprime(neq),yscal(neq)
common/steady/chib,chic,alpha,beta
common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
dimension q(4,2),qf(4,8),qdf(4,4),qnot(4)
dimension dr(100),di(100),br(100),bi(100)
c
data eps,max_step/1.0e-08,1000000/
c
c 1 2 3 4 5 6 7 2
c
c power law cooling function: lambda = C*rho**(beta-alpha)*P**alpha
c
write(6,*) 'Lambda = constant x rho^beta x temp^alpha'
write(6,*) ' '
write(6,*) 'Enter beta'
read(5,*) beta
write(6,*) 'Enter alpha, delta alpha, iter'
read(5,*) alpha,delta_alpha,max_alpha
c
write(6,*) 'Enter xcut--velocity cut-off at bottom of
& cooling layer'
read(5,*) xcut
xmin = - 0.01 * xcut
c
c chib: magnetic parameter = 0.5 * [U(mag)/U(flow)]
c
write(6,*) 'Enter magnetic parameter'
read(5,*) chib
c
if(chib.eq.0.0) then
nmin = 2
nmhd = 9
else
nmin = 4
nmhd = 11
end if
c
c perturbation wave vector (corrugation inverse wavelength)
c
write(6,*) 'Enter wave, delta kappa, iter '
read(5,*) wave,delta_wave,max_kappa
c
c eigenvalue guesses
c
write(6,*) 'Enter range for initial eigenvalues'
write(6,*) ' fr1,fr2,fi1,fi2'
write(6,*) ' byr1,byr2,byi1,byi2'
read(5,*) q(1,1),q(1,2),q(2,1),q(2,2)
read(5,*) q(3,1),q(3,2),q(4,1),q(4,2)
c
c max_alpha: number of cooling function temperature exponents to
c consider
c
do n_alpha = 1 , max_alpha
c
c calculate steady-state shock structure
c
call equilibrium
c
c max_kappa: number of corrugation wave vectors to consider
c
do n_kappa = 1 , max_kappa
c
c ======================================================================
c
c start iteration loop to find eigenvalues for given cooling function
c exponent and corrugation wave vector
c
c ======================================================================
c
111 continue
do 200 njacobi = 1 , 2*nmin+1
go to (1,2,3,4,5,6,7,8,9) njacobi
1 continue
fr = q(1,1)
fi = 0.5*(q(2,1)+q(2,2))
byr1 = 0.5*(q(3,1)+q(3,2))
byi1 = 0.5*(q(4,1)+q(4,2))
go to 10
2 continue
fr = q(1,2)
fi = 0.5*(q(2,1)+q(2,2))
byr1 = 0.5*(q(3,1)+q(3,2))
byi1 = 0.5*(q(4,1)+q(4,2))
go to 10
3 continue
fr = 0.5*(q(1,1)+q(1,2))
fi = q(2,1)
byr1 = 0.5*(q(3,1)+q(3,2))
byi1 = 0.5*(q(4,1)+q(4,2))
go to 10
4 continue
fr = 0.5*(q(1,1)+q(1,2))
fi = q(2,2)
byr1 = 0.5*(q(3,1)+q(3,2))
byi1 = 0.5*(q(4,1)+q(4,2))
go to 10
5 continue
fr = 0.5*(q(1,1)+q(1,2))
fi = 0.5*(q(2,1)+q(2,2))
byr1 = q(3,1)
byi1 = 0.5*(q(4,1)+q(4,2))
go to 10
6 continue
fr = 0.5*(q(1,1)+q(1,2))
fi = 0.5*(q(2,1)+q(2,2))
byr1 = q(3,2)
byi1 = 0.5*(q(4,1)+q(4,2))
go to 10
7 continue
fr = 0.5*(q(1,1)+q(1,2))
fi = 0.5*(q(2,1)+q(2,2))
byr1 = 0.5*(q(3,1)+q(3,2))
byi1 = q(4,1)
go to 10
8 continue
fr = 0.5*(q(1,1)+q(1,2))
fi = 0.5*(q(2,1)+q(2,2))
byr1 = 0.5*(q(3,1)+q(3,2))
byi1 = q(4,2)
go to 10
9 continue
fr = 0.5*(q(1,1)+q(1,2))
fi = 0.5*(q(2,1)+q(2,2))
byr1 = 0.5*(q(3,1)+q(3,2))
byi1 = 0.5*(q(4,1)+q(4,2))
go to 10
10 continue
c
c ======================================================================
c
c Start shock structure calculation for given frequency and transverse
c magnetic field vector, By
c
call jump_edelman(x,y,yprime,1)
htry = - 0.01
do j = 1 , max_step
call derivs(x,y,yprime)
do k = 1, nmhd
yscal(k) = abs(y(k))+abs(x*yprime(k))+1.0e-30
end do
call rkqc(y,yprime,nmhd,x,htry,eps,yscal,hdid,hnext)
htry = hnext
xmax = xcut - x
if(htry.lt.xmax) htry = xmax
if(x.le.xcut) go to 100
end do
c
c End of shock structure solution
c
c ======================================================================
c
100 continue
c
if (njacobi.eq.2*nmin+1) then
c
qnot(1) = - y(6)
qnot(2) = - y(7)
qnot(3) = - y(10)
qnot(4) = - y(11)
c
else
c
qf(1,njacobi) = y(6)
qf(2,njacobi) = y(7)
qf(3,njacobi) = y(10)
qf(4,njacobi) = y(11)
nprint=1
if(nprint.eq.1) then
write(6,*) ' '
write(6,20) alpha,wave,chic,chib
write(6,20) fr,fi,byr1,byi1,pzr,pzi
write(6,20) x,y(1),(qf(ii,njacobi),ii=1,4)
20 format(1x,1p8e10.2)
end if
c
end if
c
200 continue
c
do i = 1 , nmin
do j = 1 , nmin
qdf(i,j) = (qf(i,2*j)-qf(i,2*j-1))/(q(j,2)-q(j,1))
end do
end do
c
call gaussj(qdf,nmin,4,qnot,1)
on_key = 0
if(abs(qnot(1)).gt.1.0e-12) on_key = 1
if(abs(qnot(2)).gt.1.0e-12) on_key = 1
q_max = 0.0
do i = 1 , nmin
q(i,1) = 0.5 * (q(i,1)+q(i,2))
if(i.eq.1 .or. i.eq.2) then
ratio = abs(qnot(i)/q(i,1))
if(ratio.gt.q_max) q_max = ratio
end if
end do
q_max = amin1( q_max , 1.0e02 )
if (q_max.lt.0.2) then
q_max = 1.0
else
q_max = 0.2 / q_max
end if
do i = 1 , nmin
qnot(i) = q_max * qnot(i)
q(i,2) = q(i,1) + 2.0 * qnot(i)
end do
if(on_key.eq.1) go to 111
c
c ======================================================================
c
c end of iteration loop to find eigenvalues for given cooling function
c exponent and corrugation wave vector
c
c ======================================================================
c
cut=1.0e-20
if(abs(q(1,1)).le.cut) q(1,1)=cut
if(abs(q(2,1)).le.cut) q(2,1)=cut
if(abs(pzr).le.cut) pzr=cut
write(6,451) alpha,beta,wave,chic,chib
write(61,451) alpha,beta,wave,chib,
1 q(1,1),q(2,1),q(3,1),q(4,1)
write(62,451) alpha,beta,wave,chib,
1 q(1,1),q(2,1),pzr,pzi
451 format(1p8e11.3)
c
c ======================================================================
c
c Calculate and output shock structure for given frequency and By vector
c
go to 459
call jump_edelman(x,y,yprime,1)
write(63,452) x,(y(i),i=2,11)
452 format(1p12e10.2)
htry = - 0.01
do j = 1 , max_step
call derivs(x,y,yprime)
do k = 1, nmhd
yscal(k) = abs(y(k))+abs(x*yprime(k))+1.0e-30
end do
call rkqc(y,yprime,nmhd,x,htry,eps,yscal,hdid,hnext)
write(63,452) x,(y(i),i=2,11)
htry = hnext
xmax = xcut - x
if(htry.lt.xmax) htry = xmax
if(x.le.xcut) go to 300
end do
459 continue
c
c End of shock structure solution
c
c ======================================================================
c
300 continue
c
wave = wave + delta_wave
q(1,2) = 1.01 * q(1,1)
q(2,2) = 1.01 * q(2,1)
q(3,2) = 1.01 * q(3,1)
q(4,2) = 1.01 * q(4,1)
if ( n_kappa .eq. 1 ) then
dr( n_alpha ) = q(1,1)
di( n_alpha ) = q(2,1)
br( n_alpha ) = q(3,1)
bi( n_alpha ) = q(4,1)
end if
end do
wave = wave - max_kappa * delta_wave
alpha = alpha + delta_alpha
q(1,1) = dr( n_alpha )
q(1,2) = 1.01 * q(1,1)
q(2,1) = di(n_alpha)
q(2,2) = 1.01 * q(2,1)
q(3,1) = br(n_alpha)
q(3,2) = 1.01 * q(3,1)
q(4,1) = bi(n_alpha)
q(4,2) = 1.01 * q(4,1)
end do
end
c=============================================================
subroutine equilibrium
parameter (neq=11)
dimension y(neq),yprime(neq),yout(neq)
common/steady/chib,chic,alpha,beta
data max_step/100000/
chic=1000.0
100 continue
c
call jump_edelman(x,y,yprime,0)
c
do j = 1 , max_step
call derivs(x,y,yprime)
xstep=-0.01*abs(y(1)/yprime(1))
call rk4(y,yprime,1,x,xstep,yout)
y(1) = yout(1)
x = x + xstep
if(y(1).ge.-2.5e-08) go to 1000
end do
c
1000 continue
c
frac=(1.0-x)
if(abs(x).lt.1.0e-12) go to 2000
chic=chic*frac
go to 100
c
2000 continue
c
return
end
c====================================================================
function dv(x1)
common/steady/chib,chic,alpha,beta
data alpha1,beta1,ratio/1.0,1.0,1.0/
cool = chic*(-x1)**(alpha-beta)*(1.0+x1)**alpha
cool = cool+chic*ratio*(-x1)**(alpha1-beta1)*(1.0+x1)**alpha1
dv = -2.0*cool/(5.0+8.0*x1)
return
end
c====================================================================
SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT)
PARAMETER (neq=11)
DIMENSION Y(neq),DYDX(neq),YOUT(neq),Y1(neq),Y2(neq),Y3(neq)
C
HH=H/2.0
H6=H/6.0
XH=X+HH
C
DO 11 I=1,N
YOUT(I)=Y(I)+HH*DYDX(I)
11 CONTINUE
CALL DERIVS(XH,YOUT,Y1)
DO 12 I=1,N
YOUT(I)=Y(I)+HH*Y1(I)
12 CONTINUE
CALL DERIVS(XH,YOUT,Y2)
DO 13 I=1,N
YOUT(I)=Y(I)+H*Y2(I)
13 CONTINUE
CALL DERIVS(X+H,YOUT,Y3)
C
DO 14 I=1,N
YOUT(I)=Y(I)+H6*(DYDX(I)+2.0*(Y1(I)+Y2(I))+Y3(I))
14 CONTINUE
C
RETURN
END
c==================================================================
subroutine jump_edelman(x,y,yprime,n_order)
c
c n_order -- 0, steady state shock jump conditions
c n_order -- 1, perturbed shock jump conditions
c
parameter (neq=11)
real kzr,kzi
dimension y(neq),yprime(neq)
common/steady/chib,chic,alpha,beta
common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
complex zimag,zfreq,zkappa
complex deltab,zdby,zdrho,zdvy
complex zb1,zb2
c
zimag = cmplx(0.0,1.0)
zdby = cmplx(byr1,byi1)
zfreq = cmplx(fr,fi)
c
c zeroth order jump conditions
c
x = 1.00
y(1) = -0.25
if (n_order.eq.0) then
c
c zeroth order jump condition branch
c
return
c
else
c
c perturbed shock jump condition branch
c
c
c Require: kzr < 0 for outward propagation
c kzi > 0 for decaying wave at infinity
c
if (chib.le.0.5) then
b1r = 0.0
b2r = 0.0
b1i = 0.0
b2i = 0.0
byr1 = 0.0
byi1 = 0.0
pzr = 0.0
pzi = 0.0
else
call kappa(kzr,kzi)
write(6,*) kzr,kzi,' quadratic'
cc
call kappa_polynomial(kzr,kzi)
c
stop
ccc call kappa_edelman(kzr,kzi)
cc
pzr = kzr
pzi = kzi
zkappa = cmplx(kzr,kzi)
ckr = cos(kzr)*exp(-kzi)
cki = sin(kzr)*exp(-kzi)
bkr = (kzr*byr1+kzi*byi1)/(kzr**2+kzi**2)
bki = (kzr*byi1-kzi*byr1)/(kzr**2+kzi**2)
b1r = wave*bkr
b1i = wave*bki
zb1 = wave*(zdby/zkappa)
bkr = fi*pzr-fr*pzi
bki = -fr*pzr-fi*pzi
b2r = (byr1*bkr-byi1*bki)/(kzr**2+kzi**2)
b2i = (byi1*bkr+byr1*bki)/(kzr**2+kzi**2)
zb2 = zimag*zfreq*(zdby/zkappa)
end if
c
c
fin = 1.0/(fr**2+fi**2)
c
c
c y(2),y(3)--density perturbation
y(2) = - b1r
y(3) = - b1i
c y(4),y(5)--tangential v perturbation (y-component of velocity)
y10 = (+3.0*wave*fi*fin+4.0*(1.0-2.0*chib)*byr1)/(1.0-8.0*chib)
y11 = (+3.0*wave*fr*fin+4.0*(1.0-2.0*chib)*byi1)/(1.0-8.0*chib)
deltab = (+3.0*zimag*wave/zfreq+4.0*(1.0-2.0*chib)*zdby)
& /(1.0-8.0*chib)
y(4) = y10 + 4.0*b2r
y(5) = y11 + 4.0*b2i
c y(6),y(7)--perpendicular v perturbation (z-component of velocity)
y(6) = - 3.0
y(7) = 0.0
c y(8),y(9)--pressure perturbation
y(8) = + 2.0 - b1r
y(9) = - b1i
c y(10),y(11)--tangential B-field perturbation (y-component of B)
y(10) = y10
y(11) = y11
c
c
return
end if
end
c==================================================================
c=======================================================================
subroutine derivs(x,y,yprime)
c
c
c f(1): equilibrium momentum equation
c f(2),f(3): linearized continuity equation
c f(4),f(5): linearized y-component of momentum equation
c f(6),f(7): linearized z-component of momentum equation
c f(8),f(9): linearized energy equation
c f(10),f(11): linearized y-component of magnetic equation
c
c note--the z-component of the magnetic equation is algebraic and
c is used to eliminate the z-component of the magnetic field
c perturbation
c
parameter (neq=11)
dimension y(neq),yprime(neq),f(neq)
common/steady/chib,chic,alpha,beta
common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
c
data alpha1,beta1,ratio/1.0,1.0,1.0/
data gamma,gmin1/1.66666667,0.66666667/
c
c
frnew = fr
finew = fi
c
cc fr = finew
cc fi = -frnew
c
c
v0 = y(1)
rr = y(2)
ri = y(3)
yr = y(4)
yi = y(5)
zr = y(6)
zi = y(7)
pr = y(8)
pi = y(9)
byr = y(10)
byi = y(11)
fin = 1.0/(fr**2+fi**2)
bzr = v0*wave*fin*(fi*(yr-byr)-fr*(yi-byi))
bzi = v0*wave*fin*(fr*(yr-byr)+fi*(yi-byi))
p0 = 1.0+v0
dvdr = dv(v0)
dpdr = dvdr
cool = gmin1*chic*(p0)**alpha*(-v0)**(alpha-beta)
cool1 = gmin1*chic*ratio*(p0)**alpha1*(-v0)**(alpha1-beta1)
c
f(1) = dvdr
f(2) = -(rr*fr-ri*fi+x*dvdr/v0)+wave*v0*yi
f(3) = -(ri*fr+rr*fi) -wave*v0*yr
f(4) = -(yr*fr-yi*fi)-dvdr*yr-wave*p0*pi-wave*x*dpdr*fi*fin
1 +2.0*chib*wave*bzi
f(5) = -(yr*fi+yi*fr)-dvdr*yi+wave*p0*pr-wave*x*dpdr*fr*fin
1 -2.0*chib*wave*bzr
f(6) = -(zr*fr-zi*fi-x*dvdr/v0)-dvdr*rr+dpdr*pr
1 -2.0*dvdr*zr
f(7) = -(zr*fi+zi*fr)-dvdr*ri+dpdr*pi
1 -2.0*dvdr*zi
f(8) = -p0*(pr*fr-pi*fi-gamma*(rr*fr-ri*fi))
1 -cool*(x/v0+(alpha-1.0)*pr+(beta-alpha)*rr+fr*fin-zr)
1 -cool1*(x/v0+(alpha1-1.0)*pr+(beta1-alpha1)*rr+fr*fin-zr)
f(9) = -p0*(pr*fi+pi*fr-gamma*(rr*fi+ri*fr))
1 -cool*((alpha-1.0)*pi+(beta-alpha)*ri-fi*fin-zi)
1 -cool1*((alpha1-1.0)*pi+(beta1-alpha1)*ri-fi*fin-zi)
f(10) = -dvdr/v0*(yr-byr)+(byr*fr-byi*fi)/v0
f(11) = -dvdr/v0*(yi-byi)+(byi*fr+byr*fi)/v0
c
yprime(1) = f(1)
yprime(2) = (f(2)-f(6)-f(8)/v0)/(v0+gamma*p0)
yprime(3) = (f(3)-f(7)-f(9)/v0)/(v0+gamma*p0)
yprime(4) = (f(4)-2.0*chib*f(10))/(v0-2.0*chib)
yprime(5) = (f(5)-2.0*chib*f(11))/(v0-2.0*chib)
yprime(6) = f(2)/v0-yprime(2)
yprime(7) = f(3)/v0-yprime(3)
yprime(8) = f(8)/(v0*p0)+gamma*yprime(2)
yprime(9) = f(9)/(v0*p0)+gamma*yprime(3)
yprime(10) = yprime(4)-f(10)
yprime(11) = yprime(5)-f(11)
c
c
fr = frnew
fi = finew
c
c
return
end
c=====================================================================
SUBROUTINE GAUSSJ(A,N,NP,B,M)
PARAMETER (NMAX=50)
DIMENSION A(4,4),B(4),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX)
DO 11 J=1,N
IPIV(J)=0
11 CONTINUE
DO 22 I=1,N
BIG=0.
DO 13 J=1,N
IF(IPIV(J).NE.1)THEN
DO 12 K=1,N
IF (IPIV(K).EQ.0) THEN
IF (ABS(A(J,K)).GE.BIG)THEN
BIG=ABS(A(J,K))
IROW=J
ICOL=K
ENDIF
ELSE IF (IPIV(K).GT.1) THEN
PAUSE 'Singular matrix'
ENDIF
12 CONTINUE
ENDIF
13 CONTINUE
IPIV(ICOL)=IPIV(ICOL)+1
IF (IROW.NE.ICOL) THEN
DO 14 L=1,N
DUM=A(IROW,L)
A(IROW,L)=A(ICOL,L)
A(ICOL,L)=DUM
14 CONTINUE
DO 15 L=1,M
DUM=B(IROW)
B(IROW)=B(ICOL)
B(ICOL)=DUM
15 CONTINUE
ENDIF
INDXR(I)=IROW
INDXC(I)=ICOL
IF (A(ICOL,ICOL).EQ.0.) PAUSE 'Singular matrix.'
PIVINV=1./A(ICOL,ICOL)
A(ICOL,ICOL)=1.
DO 16 L=1,N
A(ICOL,L)=A(ICOL,L)*PIVINV
16 CONTINUE
DO 17 L=1,M
B(ICOL)=B(ICOL)*PIVINV
17 CONTINUE
DO 21 LL=1,N
IF(LL.NE.ICOL)THEN
DUM=A(LL,ICOL)
A(LL,ICOL)=0.
DO 18 L=1,N
A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18 CONTINUE
DO 19 L=1,M
B(LL)=B(LL)-B(ICOL)*DUM
19 CONTINUE
ENDIF
21 CONTINUE
22 CONTINUE
DO 24 L=N,1,-1
IF(INDXR(L).NE.INDXC(L))THEN
DO 23 K=1,N
DUM=A(K,INDXR(L))
A(K,INDXR(L))=A(K,INDXC(L))
A(K,INDXC(L))=DUM
23 CONTINUE
ENDIF
24 CONTINUE
RETURN
END
c=======================================================================
SUBROUTINE RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT)
PARAMETER (FCOR=.0666666667,
* ONE=1.,SAFETY=0.9,ERRCON=6.E-4,neq=11)
DIMENSION Y(neq),DYDX(neq),YSCAL(neq),YTEMP(neq),YSAV(neq)
1 ,DYSAV(neq)
PGROW=-0.20
PSHRNK=-0.25
XSAV=X
DO 11 I=1,N
YSAV(I)=Y(I)
DYSAV(I)=DYDX(I)
11 CONTINUE
H=HTRY
1 HH=0.5*H
CALL RK4(YSAV,DYSAV,N,XSAV,HH,YTEMP)
X=XSAV+HH
CALL DERIVS(X,YTEMP,DYDX)
CALL RK4(YTEMP,DYDX,N,X,HH,Y)
X=XSAV+H
IF(X.EQ.XSAV)PAUSE 'Stepsize not significant in RKQC.'
CALL RK4(YSAV,DYSAV,N,XSAV,H,YTEMP)
ERRMAX=0.
DO 12 I=1,N
YTEMP(I)=Y(I)-YTEMP(I)
ERRMAX=MAX(ERRMAX,ABS(YTEMP(I)/YSCAL(I)))
12 CONTINUE
ERRMAX=ERRMAX/EPS
IF(ERRMAX.GT.ONE) THEN
H=SAFETY*H*(ERRMAX**PSHRNK)
GOTO 1
ELSE
HDID=H
IF(ERRMAX.GT.ERRCON)THEN
HNEXT=SAFETY*H*(ERRMAX**PGROW)
ELSE
HNEXT=4.*H
ENDIF
ENDIF
DO 13 I=1,N
Y(I)=Y(I)+YTEMP(I)*FCOR
13 CONTINUE
RETURN
END
c========================================================================
subroutine kappa_edelman (kzr,kzi)
real kzr,kzi
common/steady/chib,chic,alpha,beta
common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
c
vpre = -1.0
wr = 2.0 * fi / ( wave * vpre )
wi = - 2.0 * fr / ( wave * vpre )
w1 = wr * wr - wi * wi
w2 = wr * wi
r1 = ( 2.0 * chib - 1.0 )
r2 = ( 2.0 * chib ) / ( 2.0 * chib - 1.0 ) * * 2
wrchi = 0.5 * wr / r1
wichi = 0.5 * wi / r1
b1 = ( 0.25 * w1 - r1 )
radical = sqrt ( b1 * b1 + 0.25 * w2 * w2 )
c
c real part of kappa
c
xlp = sqrt ( 0.5 * r2 * ( b1 + radical ) )
xrpm = 0.5 * ( wr / r1 ) - xlp
xrpp = 0.5 * ( wr / r1 ) + xlp
c
c imaginary part of kappa
c
xipm = 0.5 * ( wi / r1 ) + r2 * 0.25 * w2 / ( xrpm - wrchi )
xipp = 0.5 * ( wi / r1 ) + r2 * 0.25 * w2 / ( xrpp - wrchi )
c
c
xrpm = wave * xrpm
xrpp = wave * xrpp
c
c
xipm = wave * xipm
xipp = wave * xipp
c
c
write(6,90) fr,fi
write(6,110) xrpm,xipm,xrpp,xipp
90 format(1p2e12.4)
100 format(3h"-",1p2e12.3,4x,1p2e12.4)
110 format(3h"+",1p2e12.3,4x,1p2e12.4)
c
c check for outgoing and damping waves at + infinity
c
if(xrpm.ge.0.0) then
krpm = -1
else
krpm = 1
end if
if(xrpp.ge.0.0) then
krpp = -1
else
krpp = 1
end if
c
c
if(xipm.ge.0.0) then
kipm = 1
else
kipm = -1
end if
if(xipp.ge.0.0) then
kipp = 1
else
kipp = -1
end if
c
c
if (krpm.eq.1 .and. kipm.eq.1) then
kzr=xrpm
kzi=xipm
end if
if (krpp.eq.1 .and. kipp.eq.1) then
kzr=xrpp
kzi=xipp
end if
c
return
end
c========================================================================
subroutine kappa (kzr,kzi)
real kzr,kzi,kzrp,kzip,kzrm,kzim
common/steady/chib,chic,alpha,beta
common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
c
complex zfreq,zkappa
c
zfreq = cmplx(fr,fi)
c
r1 = ( 2.0 * chib - 1.0 ) * wave * * 2 + ( fr * * 2 - fi * * 2 )
r2 = 2.0 * fr * fi
phase = atan(r2/r1)
if ( r1 .le. 0.0 ) phase = phase + acos( - 1.0 )
ai = ( r1 * * 2 + r2 * * 2 ) * * 0.25 * cos( 0.5 * phase )
ar = - ( r1 * * 2 + r2 * * 2 ) * * 0.25 * sin( 0.5 * phase )
ar = sqrt( 2.0 * chib ) / ( 2.0 * chib - 1.0 ) * ar
ai = sqrt( 2.0 * chib ) / ( 2.0 * chib - 1.0 ) * ai
br = - fi / ( 2.0 * chib - 1.0 )
bi = fr / ( 2.0 * chib - 1.0 )
kzrm = br - ar
kzim = bi - ai
kzrp = br + ar
kzip = bi + ai
k_success = 0
if (kzrm.le.0.0) then
if (kzim.ge.0.0) then
write(6,*) 'good root for - sign'
kzr = kzrm
kzi = kzim
k_success = 1
end if
end if
if (kzrp.le.0.0) then
if (kzip.ge.0.0) then
write(6,*) 'good root for + sign'
kzr = kzrp
kzi = kzip
k_success = 1
end if
end if
if (k_success.eq.0) then
write(6,*) 'no good root for either sign'
kzr = 1.0
kzi = 1.0
end if
write(6,*) kzr,kzi
zkappa = cmplx(kzr,kzi)
return
end
c=======================================================================
subroutine kappa_polynomial(kzr,kzi)
c
c n_order -- 0, steady state shock jump conditions
c n_order -- 1, perturbed shock jump conditions
c
parameter (neq=11)
real kzr,kzi
dimension y(neq),yprime(neq)
common/steady/chib,chic,alpha,beta
common/eigen/fr,fi,wave,pzr,pzi,byr1,byi1
parameter(mroots=4,mroots1=mroots+1,maxm=101)
COMPLEX A(mroots1),ROOTS(mroots),AD(MAXM)
c
xmachsq=400.0e00
velocity=-1.0
gamma=(5.0/3.0)
fedr=fr
fedi=fi
fr= fedi
fi=-fedr
gmach=gamma/xmachsq
bsq=(wave*velocity)**2
fsqr = fr**2-fi**2
fsqi =+2.0*fr*fi
fcuber= fr*fsqr-fi*fsqi
fcubei= fr*fsqi+fi*fsqr
fquadr= fsqr**2-fsqi**2
fquadi=+2.0*fsqr*fsqi
c
c
c Dispersion Relation: fourth-order polynomial
c
c
c4r=(1.0-2.0*chib)*(1.0-gmach)
c4i=0.0
c3r=(2.0-2.0*chib-gmach)*2.0*fr
c3i=(2.0-2.0*chib-gmach)*2.0*fi
c2r=-bsq*2.0*chib-bsq*(1.0-2.0*chib)*gmach
& +(6.0-2.0*chib-gmach)*fsqr
c2i=(6.0-2.0*chib-gmach)*fsqi
c1r=(-2.0*bsq*(2.0*chib+gmach))*fr+4.0*fcuber
c1i=(-2.0*bsq*(2.0*chib+gmach))*fi+4.0*fcubei
c0r=fquadr-bsq*(2.0*chib+gmach)*fsqr
c0i=fquadi-bsq*(2.0*chib+gmach)*fsqi
c
c
a(5)=cmplx(c4r,c4i)
a(4)=cmplx(c3r,c3i)
a(3)=cmplx(c2r,c2i)
a(2)=cmplx(c1r,c1i)
a(1)=cmplx(c0r,c0i)
c
write(6,991) fr,fi,wave,chib
991 format(1p4e12.4)
c
c
call ZROOTS(A,mroots,ROOTS,POLISH)
c
print *,' '
do i=1,mroots
write(6,*) roots(i)/velocity
end do
print *,' '
c
c
c Mach number set to infinity
c
c
c4r=0.0
c4i=0.0
c3r=0.0
c3i=0.0
c2r=(1.0-2.0*chib)
c2i=0.0
c1r=2.0*fr
c1i=2.0*fi
c0r=fsqr-2.0*bsq*chib
c0i=fsqi
c
c
a(5)=(0.0,0.0)
a(4)=(0.0,0.0)
a(3)=cmplx(c2r,c2i)
a(2)=cmplx(c1r,c1i)
a(1)=cmplx(c0r,c0i)
c
write(6,991) ,fr,fi,wave,chib
c
c
call ZROOTS(A,2,ROOTS,POLISH)
c
print *,' '
do i=1,2
write(6,*) roots(i)/velocity
end do
print *,' '
c
c
fr=fedr
fi=fedi
return
end
SUBROUTINE ZROOTS(A,M,ROOTS,POLISH)
PARAMETER (EPS=1.E-8,MAXM=101)
parameter(mroots=4,mroots1=mroots+1)
COMPLEX A(mroots1),ROOTS(mroots),AD(MAXM),X,B,C
LOGICAL POLISH
DO 11 J=1,M+1
AD(J)=A(J)
11 CONTINUE
DO 13 J=M,1,-1
X=CMPLX(0.0,0.0)
CALL LAGUER(AD,J,X,EPS,.FALSE.)
IF(ABS(AIMAG(X)).LE.2.0*EPS**2*ABS(REAL(X))) X=CMPLX(REAL(X),0.)
ROOTS(J)=X
B=AD(J+1)
DO 12 JJ=J,1,-1
C=AD(JJ)
AD(JJ)=B
B=X*B+C
12 CONTINUE
13 CONTINUE
IF (POLISH) THEN
DO 14 J=1,M
CALL LAGUER(A,M,ROOTS(J),EPS,.TRUE.)
14 CONTINUE
ENDIF
DO 16 J=2,M
X=ROOTS(J)
DO 15 I=J-1,1,-1
IF(REAL(ROOTS(I)) .LE. REAL(X))GO TO 10
ROOTS(I+1)=ROOTS(I)
15 CONTINUE
I=0
10 ROOTS(I+1)=X
16 CONTINUE
RETURN
END
c=======================================================================
SUBROUTINE LAGUER(A,M,X,EPS,POLISH)
parameter(mroots=4,mroots1=mroots+1,maxm=101)
COMPLEX A(mroots1),X,DX,X1,B,D,F,G,H,SQ,GP,GM,G2,ZERO
LOGICAL POLISH
PARAMETER (ZERO=(0.0,0.0),EPSS=6.E-8,MAXIT=100)
DXOLD=CABS(X)
DO 12 ITER=1,MAXIT
B=A(M+1)
ERR=CABS(B)
D=ZERO
F=ZERO
ABX=CABS(X)
DO 11 J=M,1,-1
F=X*F+D
D=X*D+B
B=X*B+A(J)
ERR=CABS(B)+ABX*ERR
11 CONTINUE
ERR=EPSS*ERR
IF(CABS(B) .LE. ERR) THEN
DX=ZERO
RETURN
ELSE
G=D/B
G2=G*G
H=G2-2.0*F/B
SQ=CSQRT((M-1)*(M*H-G2))
GP=G+SQ
GM=G-SQ
IF(CABS(GP) .LT. CABS(GM)) GP=GM
DX=M/GP
ENDIF
X1=X-DX
IF(X.EQ.X1)RETURN
X=X1
CDX=CABS(DX)
IF(ITER.GT.6 .AND. CDX.GE.DXOLD) RETURN
DXOLD=CDX
IF(.NOT.POLISH)THEN
IF(CABS(DX) .LE. EPS*CABS(X)) RETURN
ENDIF
12 CONTINUE
PAUSE 'too many iterations'
RETURN
END
c=======================================================================