C ***********************
C * GUDERLEY SOLUTION *
C**********************************************************************
c Version: 2007.01.18 (c) M.M.BASKO 27.01.1988
c-----------------------------------------------------------------------
c The SUBROUTINE GUDER(ISYM,GAMMA) from this package calculates all
c the key parameters of the Guderley solution for all
c 1 < \gamma < \gamma_1, when the solution is unique and passes
c through the saddle singular point P_3.
c To calculate the value of \gamma_1,
c call SUBROUTINE GAMM1(ISYM,G1,NITER).
c SUBROUTINE GUD(TIME,RAD,DEN0,RAD0,PRES0,U,RO,P) calculates the
c physical values of velocity U=U(TIME,RAD), density RO=RO(TIME,RAD)
c and pressure P=P(TIME,RAD) according to the Guderley solution.
c Before calling GUD, the SUBROUTINE GUDER must be called at least
c once.
c***********************************************************************
c***********************************************************************
PROGRAM RUN
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
LOGICAL LEXIST,LOPEN
CHARACTER FNAME*12,FSTAT*3
C DO ISYM=1,2
C CALL GAMM1(ISYM,G1,NITER)
C WRITE(0,*) 'ISYM=',ISYM,', NITER=',NITER,', gamma_1=',G1
C ENDDO
FNAME='p.dat'
FSTAT='OLD'
INQUIRE(FILE=FNAME,EXIST=LEXIST)
IF(.NOT.LEXIST) FSTAT='NEW'
OPEN(21,FILE=FNAME,STATUS=FSTAT,ACCESS='SEQUENTIAL'
*,FORM='FORMATTED')
GAMMA=5.D0/3.D0
ISYM=2
CALL GUDER(ISYM,GAMMA)
DEN0=1.D0
RAD0=0.01097
PRES0=1.D0
TIME =0.0160
DO J=1,401
RAD=10.D0**(-4.D0+.01D0*(J-1))
CALL GUD(TIME,RAD,DEN0,RAD0,PRES0,U,RO,P)
WRITE(21,9010) RAD,RO
ENDDO
WRITE(0,*) ' '
WRITE(0,*) '*** End of the Guderley run ***'
STOP
9010 FORMAT(1P2E12.4)
END
c***********************************************************************
SUBROUTINE GUDER(ISYM,GAMMA)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C ======================================================================
c This routine calculates the self-similarity index ALPHA=1/(1+BETA)
c and the coordinates of 6 fiducial points of the Guderley solution.
c ISYM = 1 -> cylinder;
c ISYM = 2 -> sphere;
c GAMMA is the adiabatic index.
c Called by:
c Calls :
c ======================================================================
DIMENSION Y(3),Y6(3)
COMMON/GUDR1/S,GAM,ALP,BET,CAP,S1,AL1,GM1,GM2
COMMON/GUDR2/V(6),Z(6),GL(6),CSIL(6),CADI
COMMON/GUDR3/EPS,QZV2,QCSV2,VINF,ZINF,GINF,QV1Z6,ZZER,GZER,GEXP6
EXTERNAL DZDV,DZSIDV,DVZGL,DVD1Z
C: check for gamma < gamma_1:
IF(GAMMA.LE.1.D0) GOTO 1
IF(ISYM.EQ.1 .AND. GAMMA.LT.1.9092084D0) GOTO 2
IF(ISYM.EQ.2 .AND. GAMMA.LT.1.8697678D0) GOTO 2
1 WRITE(*,9010) ISYM,GAMMA
WRITE(0,9010) ISYM,GAMMA
STOP
2 CONTINUE
C: initial assignments:
S=ISYM
S1=S+1.D0
GAM=GAMMA
GM1=GAM-1.D0
GM2=GAM+1.D0
c fiducial point #1 is the incident shock front (xi=1):
V(1)=2.D0/GM2
Z(1)=2.D0*GAM*GM1/GM2**2
GL(1)=LOG(GM2/GM1)
CSIL(1)=0.D0
EPS=1.D-12
NSTEP=10.D0+2.D0*SQRT(SQRT(1.D0/EPS))
EPS1=EPS+1.D-7
EPSV3=1.D-12
EPS6=1.D-12
C: integrate from the coalesced singular points "a" = "b" to V=V(1):
BEDU=1.D0/(2.D0+GAM+SQRT(8.D0*GAM))
BET=S*GAM*BEDU
AL1=1.D0+BET
ALP=1.D0/AL1
CAP=2.D0*BET/GAM
VA=.5D0*(1.D0+BEDU*(2.D0-GAM))
VA1=1.D0-VA
AQ=(S1*VA-CAP)/VA1
BQ=4.D0*VA+S-2.D0-BET*(1.D0+4.D0/GAM)
CQ=VA1*(4.D0+S*GM1+BET*(3.D0-GAM+2.D0/GAM)-2.D0*VA*(2.D0+S*GM1))
c initial direction of the integral curve:
Q=.5D0*(SQRT(BQ*BQ+4.D0*AQ*CQ)-BQ)/AQ
X=VA+EPS1
Y(1)=VA1**2+Q*EPS1
H=(V(1)-X)/NSTEP
CALL RUTIN(1,X,Y,DZDV,H,NSTEP)
IVA=1
IF(Y(1).LT.Z(1))IVA=-1
C IVA= 1 -> the Guderley solution passes through the lower point "a"
C IVA=-1 -> the Guderley solution passes through the upper point "b"
C: calculate the self-similarity index ALPHA: iterations by using
c the "chord" method:
HB=-.002D0*(1.D0+(2.D0-GAM)**2)/SQRT(2.D0*GAM)
FP=Y(1)-Z(1)
BEDP=BEDU
NITER=0
c: iterations for BED=BETA/(GAMMA*S):
10 BED=BEDP+HB
IF(BED.GE.BEDU) BED=.5D0*(BEDP+BEDU)
IF(BED.LE.-1.D0/(S*GAM)) BED=.5D0*(BEDP-1.D0/(S*GAM))
IF(IVA.EQ.1 .AND. BED.LE.0.D0) BED=.5D0*BEDP
20 NITER=NITER+1
BET=S*GAM*BED
AL1=1.D0+BET
ALP=1.D0/AL1
IF(NITER.GT.100) THEN
WRITE(*,9020) ALP,BET,Z(1),Y(1)
WRITE(0,9020) ALP,BET,Z(1),Y(1)
STOP
ENDIF
CAP=2.D0*BET/GAM
VA=1.D0+(2.D0-GAM)*BED
R=VA*VA-8.D0*BED
IF(R.LT.0.D0) R=0.D0
VA=.5D0*(VA+IVA*SQRT(R))
VA1=1.D0-VA
AQ=(S1*VA-CAP)/VA1
BQ=4.D0*VA+S-2.D0-BET*(1.D0+4.D0/GAM)
CQ=VA1*(4.D0+S*GM1+BET*(3.D0-GAM+2.D0/GAM)-2.D0*VA*(2.D0+S*GM1))
Q=.5D0*(SQRT(BQ*BQ+4.D0*AQ*CQ)-BQ)/AQ
X=VA+EPS1
Y(1)=VA1**2+Q*EPS1
H=(V(1)-X)/NSTEP
CALL RUTIN(1,X,Y,DZDV,H,NSTEP)
F=Y(1)-Z(1)
IF(ABS(F).LT.1.D-3*EPS.OR.ABS(HB).LT.1.D-15) GOTO 40
R=F/ABS(F)
IF(R*FP.GT.0.D0) GOTO 30
R=(BEDP*F-BED*FP)/(F-FP)
BEDP=BED
FP=F
HB=.3D0*ABS(HB)
IF(R.LT.BED) HB=-HB
BED=R
GOTO 20
30 HB=1.1D0*HB
IF(ABS(F).LT.ABS(FP)) GOTO 35
HB=-HB
GOTO 10
35 BEDP=BED
FP=F
GOTO 10
c: iterations for BED completed.
40 V(6)=CAP/S1
GEXP6=CAP/(1.D0-V(6))
QV1Z6=V(6)*(1.D0-V(6))**2*(AL1-V(6))/(3.D0+S-2.D0*V(6))
WRITE(0,9400) ISYM,GAM,ALP,NITER,F
WRITE(*,9410) ISYM,GAM,ALP,NITER,F
WRITE(*,9501) CSIL(1),V(1),Z(1),EXP(GL(1))
C: parameters of the fiducial point #2 (the saddle singular
c point "a"=P3):
c QZV2 = dZ/dV, QCSV2=dln(xi)/dV
V(2)=VA
Z(2)=VA1**2
QZV2=Q
QCSV2=-(2.D0+Q/VA1)/(Q*AQ+2.D0*VA+S-BET-CAP)
C: calculate ln(xi) at point #2:
X=VA+EPS1
Y(1)=Z(2)+QZV2*EPS1
Y(2)=QCSV2*EPS1
H=(V(1)-X)/NSTEP
CALL RUTIN(2,X,Y,DZSIDV,H,NSTEP)
CSIL(2)=-Y(2)
CADI=GAM*GL(1)+LOG(GM2*.5D0/GAM)
GL(2)=(CADI+LOG(Z(2))+2.D0*AL1*CSIL(2)+2.D0*BET*
&LOG(1.D0-V(2))/S1)/(GM1-2.D0*BET/S1)
C: parameters of the fiducial point #3 (near the node singular
c point P4 of V=Z=0):
X=CSIL(2)+EPS1
Y(1)=V(2)+EPS1/QCSV2
Y(2)=Z(2)+EPS1*QZV2/QCSV2
Y(3)=GL(2)+EPS1*(1.D0/QCSV2+S1*V(2))/(1.D0-V(2))
H=2.D0/NSTEP
NST3=0
110 NST3=NST3+1
CALL RUTIN(3,X,Y,DVZGL,H,1)
IF(Y(1).GT.EPSV3) GOTO 110
CSIL(3)=X
V(3)=Y(1)
Z(3)=Y(2)
GL(3)=Y(3)
GINF=EXP(GL(3))
ZINF=EXP(LOG(Z(3))+2.D0*AL1*CSIL(3))
VINF=EXP(LOG(V(3))+AL1*CSIL(3))
IF(ABS(V(3)).LT.1.D-6) THEN
RAB=-V(3)*(1.D0+.5D0*V(3)+V(3)**2/3.D0)
ELSE
RAB=LOG(1.D0-V(3))
ENDIF
F=(CADI+LOG(Z(3))+2.D0*AL1*CSIL(3)+2.D0*BET*RAB/S1)/
&(GM1-2.D0*BET/S1)-GL(3)
R1=EXP(GL(2))
WRITE(*,9502) CSIL(2),V(2),Z(2),R1,QZV2,QCSV2
WRITE(*,9503) CSIL(3),V(3),Z(3),GINF,VINF,ZINF,NST3,F
C-----------------------------------------------------------------------
C: REFLECTED SHOCK:
c Fiducial point #4 - immediately before the reflected shock (along the
c integration path);
c Fiducial point #5 - immediately after the reflected shock;
c Fiducial point #6 - at a distance of EPS6 along the coordinate
c \sigma=1/Z from the singular saddle point P6, where xi=0.
C: integrate until the Hugoniot image reaches V=V(6):
Y(1)=-Y(1)
H=-H
120 CALL RUTIN(3,X,Y,DVZGL,H,1)
R=1.D0-Y(1)
V5=1.D0-(GM1*R+2.D0*Y(2)/R)/GM2
IF(V5.GT.V(6)) GOTO 120
CS5=X
Z5=Y(2)+.5D0*GM1*(R*R-(1.D0-V5)**2)
H=-H
C: integrate from point #6 down to Z=Z5
X6=LOG(EPS6)
Z(6)=1.D0/EPS6
V(6)=V(6)+QV1Z6/Z(6)
Y6(1)=V(6)
Y6(2)=0.D0
Y6(3)=0.D0
c.......................................................................
c Y6(1)=V, Y6(2)=ln(xi), Y6(3)=ln(G).
c.......................................................................
H6=(LOG(Z(6))-LOG(Z5))/NSTEP
CALL RUTIN(3,X6,Y6,DVD1Z,H6,NSTEP)
F=Y6(1)-V5
IF(F.GE.0.D0) GOTO 128
R1=EXP(-X6)
WRITE(*,9511) Z5,V5,V(6),Y6(1),R1
9511 FORMAT ('After this integration, at Z=Z5=',1PE18.10/
&', we must have V5 < V(6) < Y6(1), whereas we have gotten'
&/'V5,V(6),Y6(1)=',1P3E18.10/'Z=EXP(-X6)=',1PE17.10)
STOP
C: integrate until the integral curve emerging from point #6 crosses
c the Hugoniot image:
128 H=0.1D0*H
NST5=0
130 NST5=NST5+1
CS5P=CS5
FP=F
IR1=0
140 CALL RUTIN(3,X,Y,DVZGL,H,1)
CS5=X
R=1.D0-Y(1)
V5=1.D0-(GM1*R+2.D0*Y(2)/R)/GM2
Z(5)=Y(2)+.5D0*GM1*(R*R-(1.D0-V5)**2)
H6=-LOG(Z(5))-X6
CALL RUTIN (3,X6,Y6,DVD1Z,H6,1)
F=Y6(1)-V5
IF (IR1.EQ.1) GOTO 150
IF (F.GT.0.D0) GOTO 130
c: crossing point:
CS5=(CS5P*F-CS5*FP)/(F-FP)
H=CS5-X
IR1=1
GOTO 140
150 CSIL(4)=X
CSIL(5)=X
V(4)=Y(1)
Z(4)=Y(2)
GL(4)=Y(3)
V(5)=Y6(1)
GL(5)=GL(4)+LOG((1.D0-V(4))/(1.D0-V(5)))
CSIL(6)=CSIL(5)-Y6(2)
GL(6)=GL(5)-Y6(3)
GZER=EXP(GL(6)-GEXP6*CSIL(6))
ZZER=Z(6)*EXP((2.D0+GEXP6)*CSIL(6))
G4= EXP(GL(4))
WRITE(*,9504) CSIL(4),V(4),Z(4),G4
G5=EXP(GL(5))
PI45=(G5*Z(5))/(G4*Z(4))
WRITE(*,9505) EXP(CSIL(5)),V(5),Z(5),G5,PI45,NST5,F
G6=EXP(GL(6))
WRITE(*,9506) CSIL(6),V(6),Z(6),G6,QV1Z6,ZZER,GZER,GEXP6
WRITE(*,9510)
RETURN
9010 FORMAT('STOP: either ISYM.NE.(1 or 2), or GAMMA > gamma_1'
&/'ISYM=',I4,' GAMMA=',1PE16.8)
9020 FORMAT('STOP in GUDER: too many iterations when calculating ALPHA'
&/'ALPHA,BETA=',1P2E18.10,' Z(1),Y(1)=',1P2E20.12)
9400 FORMAT(/'Key parameters of the Guderley solution for ',
&'ISYM=',I2,' GAMMA=',F14.12/'self-sim.index \alpha=',F14.12/
&'number of iter-s=',I4,', mismatch of Z(1)=',1PE10.2)
9410 FORMAT(/72('-')/'Key parameters of the Guderley solution for ',
&'ISYM=',I2,' GAMMA=',F14.12/'self-sim.index \alpha=',F14.12/
&'number of iter-s=',I4,', mismatch of Z(1)=',1PE10.2)
9501 FORMAT(/'Point #1 (behind the incident shock): ln(xi)=',1PE17.10/
&'V(1)=',1PE17.10,' Z(1)=',1PE17.10,' G(1)=',1PE17.10)
9502 FORMAT(/'Point #2 (the singular saddle point): ln(xi)=',1PE17.10/
&'V(2)=',1PE17.10,' Z(2)=',1PE17.10,' G(2)=',1PE17.10/
&'dZ/dV=',1PE17.10,' dln(xi)/dV=',1PE17.10)
9503 FORMAT(/'Point #3 (xi=infinity, V=Z=0): ln(xi)=',1PE17.10/
&'V(3)=',1PE17.10,' Z(3)=',1PE17.10,' G(3)=',1PE17.10/
&'V*xi**(1./alpha)=',1PE17.10,' Z*xi**(2./alpha)=',1PE17.10/
&'NST3=',I6,', mismatch with ln(G) from adiab.intgral=',1PE10.3)
9504 FORMAT(/'Point #4 (in front of the reflected shock): ln(xi)=',
&1PE17.10/'V(4)=',E17.10,' Z(4)=',1PE17.10,' G(4)=',1PE17.10)
9505 FORMAT(/'Point #5 (behind the reflected shock): xi=',1PE17.10/
&'V(5)=',1PE17.10,' Z(5)=',1PE17.10,' G(5)=',1PE17.10/
&'Amplitude of the reflected shock Pi_s=',1PE17.10/
&'NST5=',I6,', mismatch of V(5)=',1PE10.3)
9506 FORMAT(/'Point #6 near xi=0: ln(xi)=',1PE17.10/'V(6)=',1PE17.10,
&' Z(6)=',1PE17.10,' G(6)=',1PE17.10/'dV/d(1/Z)=',1PE17.10,
&' Z*xi**(2.+GEXP6)=',1PE17.10/'G/xi**GEXP6=',1PE17.10,' GEXP6='
&,1PE17.10)
9510 FORMAT(72('='))
END
c***********************************************************************
SUBROUTINE DZDV(V,Z,F)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c ======================================================================
c This routine calculates F(V,Z)=dZ/dV
c Called by:
c Calls :
c ======================================================================
COMMON/GUDR1/S,GAM,ALP,BET,CAP,S1,AL1,GM1,GM2
V1=1.D0-V
R=Z*(S1*V-CAP)-V*V1*(V1+BET)
F=Z*(GM1+(Z-V1*V1)*(2.D0*AL1-V*(2.D0+S1*GM1))/R)/V1
RETURN
END
c***********************************************************************
SUBROUTINE DZSIDV(V,Z,F)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c ======================================================================
c This routine calculates F_1(V,Z)=dZ/dV, F_2(V,Z)=dln(xi)/dV.
c Z(1)=Z, Z(2)=ln(xi).
c Called by: RUTIN
c Calls :
c ======================================================================
COMMON/GUDR1/S,GAM,ALP,BET,CAP,S1,AL1,GM1,GM2
DIMENSION Z(2),F(2)
V1=1.D0-V
R=Z(1)*(S1*V-CAP)-V*V1*(V1+BET)
R1=Z(1)-V1*V1
F(1)=Z(1)*(GM1+R1*(2.D0*AL1-V*(2.D0+S1*GM1))/R)/V1
F(2)=-R1/R
RETURN
END
c***********************************************************************
SUBROUTINE DVZGL(CSL,Y,F)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c ======================================================================
c This routine calculates F_1(V,Z)=dV/dln(xi), F_2(V,Z)=dZ/dln(xi),
c F_3(V,Z)=dln(G)/dln(xi).
c CSL=ln(xi), Y(1)=V, Y(2)=Z, Y(3)=ln(G).
c Called by: RUTIN
c Calls :
c ======================================================================
COMMON/GUDR1/S,GAM,ALP,BET,CAP,S1,AL1,GM1,GM2
DIMENSION Y(3),F(3)
V=Y(1)
V1=1.D0-V
Z=Y(2)
R=(Z*(S1*V-CAP)-V*V1*(V1+BET))/(V1*V1-Z)
F(1)=R
F(2)=GM1*(V1*R+V*(V1+BET))-Z*(2.D0+CAP/V1)
F(3)=(R+S1*V)/V1
RETURN
END
c***********************************************************************
SUBROUTINE DVD1Z(X,Y,F)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c ======================================================================
c This routine calculates F_1=dV/dln(1/Z), F_2=dln(xi)/dln(1/Z),
c F_3=dln(G)/dln(1/Z).
c X=ln(1/Z); Y(1)=V; Y(2)=ln(xi); Y(3)=ln(G)
c Called by: RUTIN
c Calls :
c ======================================================================
COMMON/GUDR1/S,GAM,ALP,BET,CAP,S1,AL1,GM1,GM2
DIMENSION Y(3),F(3)
V=Y(1)
V1=1.D0-V
SIG=EXP(X)
R=S1*V-CAP-SIG*V*V1*(BET+V1)
R1=1.D0-SIG*V1*V1
R2=-V1*R/(GM1*R+R1*(2.D0*AL1-V*(2.D0+GM1*S1)))
R1=-R1/R
F(1)=R2
F(2)=R1*R2
F(3)=R2*(1.D0+S1*V*R1)/V1
RETURN
END
c***********************************************************************
SUBROUTINE RUTIN(NEQ,X,Y,DYDX,H,NSTEP)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c ======================================================================
c This routine makes NSTEP steps of magnitude H along X by using
c the 4-th order Runge-Kutta scheme for a system dY/dX = F(X,Y).
c NEQ = number of equations;
c SUBROUTINE DYDX(X,Y,F) calculates F(X,Y).
c Called by: GUDER
c Calls :
c ======================================================================
DIMENSION A(5),Y(1)
DIMENSION Z(3),F(3),DY(3)
DATA A/6.D0,3.D0,3.D0,6.D0,6.D0/
DO 99 I=1,NSTEP
DO K=1,NEQ
Z(K)=Y(K)
DY(K)=0.D0
ENDDO
XF=X
DO J=1,4
CALL DYDX(XF,Z,F)
DO K=1,NEQ
DY(K)=DY(K)+H*F(K)/A(J)
Z(K)=Y(K)+A(J+1)*H*F(K)/6.D0
ENDDO
XF=X+A(J+1)*H/6.D0
ENDDO
DO K=1,NEQ
Y(K)=Y(K)+DY(K)
ENDDO
99 X=X+H
RETURN
END
C***********************************************************************
SUBROUTINE GUD(TIME,RAD,DEN0,RAD0,PRES0,U,RO,P)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c=======================================================================
c This routine calculates U=U(TIME,RAD), RO=RO(TIME,RAD),
c P=P(TIME,RAD) from the Guderley solution.
c DEN0 = initial density in front of the incident shock;
c RAD0 = initial radius of the incident shock;
c PRES0= initial pressure immediately behind the incident shock;
c One must call GUDER before the first call of GUD!
c=======================================================================
COMMON/GUDR1/S,GAM,ALP,BET,CAP,S1,AL1,GM1,GM2
COMMON/GUDR2/V(6),Z(6),GL(6),CSIL(6),CADI
COMMON/GUDR3/EPS,QZV2,QCSV2,VINF,ZINF,GINF,QV1Z6,ZZER,GZER,GEXP6
DIMENSION Y(3)
EXTERNAL DVZGL
R=RAD
IF(R.LT.1.D-16) R=1.D-16
EPS1=EPS+1.D-7
HCSIL=1.D0/(10.D0+SQRT(SQRT(1.D0/EPS)))
D1=SQRT(.5D0*GM2*PRES0/DEN0)
T=TIME-ALP*RAD0/D1
IF(ABS(T).LT.1.D-16) T=-1.D-16
5 CSL=LOG(R/RAD0)-ALP*LOG(D1*AL1*ABS(T)/RAD0)
IF(ABS(CSL).LT.11.D0) GOTO 7
IF(CSL.GT.1.D0) T=T*EXP(AL1*(CSL-10.D0))
IF(CSL.LT.-1.D0) R=R*EXP(-10.D0-CSL)
GOTO 5
7 IF (CSL.LT.CSIL(3)) GOTO 20
GG=GINF
VV=VINF*EXP(-AL1*CSL)
IF(T.GT.0.D0) VV=-VV
ZZ=ZINF*EXP(-2.D0*AL1*CSL)
GOTO 200
20 IF(T.GT.0.D0) GOTO 100
IF (CSL.LT.0.D0) GOTO 210
R1=CSL-CSIL(2)
IF(ABS(R1).GT.EPS1) GOTO 30
VV=V(2)+R1/QCSV2
ZZ=Z(2)+R1*QZV2/QCSV2
GG=EXP(GL(2)+R1*(1.D0/QCSV2+S1*V(2))/(1.D0-V(2)))
GOTO 200
30 NSTEP=1.D0+ABS(R1)/HCSIL
H=R1/NSTEP
DCSIL=EPS1*(R1/ABS(R1))
X=CSIL(2)+DCSIL
Y(1)=V(2)+DCSIL/QCSV2
Y(2)=Z(2)+DCSIL*QZV2/QCSV2
Y(3)=GL(2)+DCSIL*(1.D0/QCSV2+S1*V(2))/(1.D0-V(2))
GOTO 190
100 IF(CSL.GT.CSIL(6)) GOTO 110
GG=GZER*EXP(GEXP6*CSL)
ZZ=ZZER*EXP(-(2.D0+GEXP6)*CSL)
VV=CAP/S1+QV1Z6/ZZ
GOTO 200
110 J=5
IF(CSL.GE.CSIL(4)) J=4
R1=CSL-CSIL(J)
NSTEP=1.D0+ABS(R1)/HCSIL
H=R1/NSTEP
X=CSIL(J)
Y(1)=V(J)
Y(2)=Z(J)
Y(3)=GL(J)
190 CALL RUTIN(3,X,Y,DVZGL,H,NSTEP)
VV=Y(1)
ZZ=Y(2)
GG=EXP(Y(3))
200 R1=ALP*R/T
U=R1*VV
RO=DEN0*GG
P=R1*ZZ*R1*DEN0*GG/GAM
RETURN
210 U=0.D0
RO=DEN0
P=0.D0
RETURN
END
c***********************************************************************
SUBROUTINE GAMM1(ISYM,G1,NITER)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C ======================================================================
c This routine calculates the value of gamma = gamma_1, below which
c the Guderley solution is unique and passes through the canonical
c saddle singular point "a".
c ISYM = 1 -> cylinder;
c ISYM = 2 -> sphere;
c Called by:
c Calls : RUTIN
c ======================================================================
DIMENSION Y(3),Y6(3)
COMMON/GUDR1/S,GAM,ALP,BET,CAP,S1,AL1,GM1,GM2
COMMON/GUDR2/V(6),Z(6),GL(6),CSIL(6),CADI
COMMON/GUDR3/EPS,QZV2,QCSV2,VINF,ZINF,GINF,QVIZ6,ZZER,GZER,GEXP6
EXTERNAL DZDV
c: initial assignments:
S=ISYM
S1=S+1.D0
GAM=1.4d0
DGAM=2.D0/7.D0
NITER=0
IVA=1
100 NITER=NITER+1
IF(NITER.GT.100) THEN
WRITE(0,*) 'STOP in GAMSTA: NITER > 100'
STOP
ENDIF
GAM=GAM+DGAM*IVA
GM1=GAM-1.D0
GM2=GAM+1.D0
c point #1 is the incident shock front:
V(1)=2.D0/GM2
Z(1)=2.D0*GAM*GM1/GM2**2
GL(1)=LOG(GM2/GM1)
CSIL(1)=0.D0
EPS=1.D-10
NSTEP=10.D0+2.D0*SQRT(SQRT(1.D0/EPS))
EPS1=EPS+1.D-8
C: integrate from the coalesced singular points "a" = "b" to V=V(1):
BEDU=1.D0/(2.D0+GAM+SQRT(8.D0*GAM))
BET=S*GAM*BEDU
AL1=1.D0+BET
ALP=1.D0/AL1
CAP=2.D0*BET/GAM
VA=.5D0*(1.D0+BEDU*(2.D0-GAM))
VA1=1.D0-VA
AQ=(S1*VA-CAP)/VA1
BQ=4.D0*VA+S-2.D0-BET*(1.D0+4.D0/GAM)
CQ=VA1*(4.D0+S*GM1+BET*(3.D0-GAM+2.D0/GAM)-2.D0*VA*(2.D0+S*GM1))
c exit direction:
Q=.5D0*(SQRT(BQ*BQ+4.D0*AQ*CQ)-BQ)/AQ
X=VA+EPS1
Y(1)=VA1**2+Q*EPS1
H=(V(1)-X)/NSTEP
CALL RUTIN(1,X,Y,DZDV,H,NSTEP)
IVAP=IVA
IVA=1
IF(Y(1).LT.Z(1))IVA=-1
C IVA= 1 -> the sought-for solution passes through the lower point "a"
C IVA=-1 -> the sought-for solution passes through the upper point "b"
IF(IVA*IVAP.LT.0) DGAM=.3D0*DGAM
IF(DGAM.GT.3.D-12) GOTO 100
G1=GAM
RETURN
END