      PROGRAM REENTR
      IMPLICIT DOUBLEPRECISION(S)
C
C     ------------------------------------------------------------
C
C     Testexample for BVPSOL: Reentry Problem.
C     Ref. Stoer, Bulirsch: Einfuehrung in die Numerische
C     Mathematik 2
C     Springer Verlag 1978, P. 191 - 197
C
C*  Purpose           Testexample for code BVPSOL
C*  Version           1.1
C*  Revision          January 1991
C*  Latest Change     January 1991
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Standard Fortran 77 environment on PC's,
C                     workstations and hosts.
C
C     ------------------------------------------------------------
C
      INTEGER M,N,LRWK,LIWK
      PARAMETER (M=6,N=7,LRWK=4000,LIWK=6000)
      DOUBLE PRECISION X(M)
      DOUBLE PRECISION Y(N,M)
      DOUBLE PRECISION WORK(LRWK)
      INTEGER IW(LIWK)
      EXTERNAL F,R,BLDFX1
      INTEGER I,J,KPRINT,M1
      INTEGER IFCN,NSTEP,NACCPT,NREJCT,NDEC,NSOL,IOPT(5)
      DOUBLE PRECISION RK,BETAR,FIX
      DOUBLE PRECISION EPS,ZEIT
      COMMON /KONST/ RK,BETAR,FIX
      COMMON /DXSTAT/ IFCN,NSTEP,NACCPT,NREJCT,NDEC,NSOL
      REAL ETIME,STIME
      INTEGER IFAIL
C:    Begin
1     FORMAT(' Reentry  Forward  Direction  With Hamiltonian At T = 0')
      WRITE(6,1)
      NSTEP = 0
      IFCN = 0
      KPRINT = 0
C
      RK = 1.0D0/2.09D2
      BETAR = 890.34D0
      FIX = 8.1D0*4.0D0*DATAN(1.0D0)/180.0D0
      OPEN(15,FILE='main_bvpsol.data',STATUS='OLD')
      READ(15,*)(X(J),J=1,M)
      DO 2 J=1,M
        READ(15,*)(Y(I,J),I=1,4)
        READ(15,*)(Y(I,J),I=5,N)
2     CONTINUE
      CLOSE(15)
      M1 = M-1
C
      Y(2,1)=-FIX
      Y(3,1)=4.0D0/209.0D0
      Y(3,M)=2.5D0/209.0D0
      EPS = 1.0D-6
C
3     FORMAT(' BVPSOL with DIFEX1')
      WRITE(6,3)
C
C     MAXIT: Maximum number of iterations
      IOPT(1) = 30
C     NONLIN: Highly nonlinear problem
      IOPT(2) = 3
C     IBVPSL: 0 = use original BVPSOL code, 1 = use "BVPSOG" code
      IOPT(3) = 1
C     IPRINT: Maximum level printout
      IOPT(4) = 1
C     LUPRI: Print output unit
      IOPT(5) = 6
      CALL ZIBSEC(STIME,IFAIL)
      CALL BVPSOL (F,R,BLDFX1,N,M,X,Y,EPS,IOPT,INFO,LRWK,WORK,LIWK,IW)
      CALL ZIBSEC(ETIME,IFAIL)
      ZEIT = ETIME-STIME
4     FORMAT('0','IFCN:',I7,3X,'ZEIT:',F8.3)
      WRITE(6,4)IFCN,ZEIT
C
      END
      SUBROUTINE F(N,X,Y,D)
      IMPLICIT DOUBLEPRECISION(S)
      INTEGER N
      DOUBLE PRECISION X
      DOUBLE PRECISION Y(7),D(7)
C     system of ODE's for example Reentry
C:    End Parameter
      DOUBLE PRECISION BRXI,CDV,CDVV,CL,CLV,COSXI,COSY2,Q1,Q2,Q3,
     *Q4,RO,SINXI,SINY2,S2MRO,S2MROV,U1,U2,U3,WURZ,XI,Y1,Y2,Y3,Y4,
     *Y5,Y6,Y7,Y1K,Y12
      DOUBLE PRECISION RK,BETAR,FIX
      COMMON /KONST/ RK,BETAR,FIX
C:    Begin
      Y1 = Y(1)
      Y2 = Y(2)
      Y3 = Y(3)
      Y4 = Y(4)
      Y5 = Y(5)
      Y6 = Y(6)
      Y7 = Y(7)
      Y1K = 1.0D0/Y1
      Y12 = Y1*Y1
      U1 = 0.6D0*Y5
      U2 = 0.9D0*Y1*Y4
      U3 = -1.0D0/DSQRT(U1*U1+U2*U2)
      BRXI =-BETAR*Y3
      IF(BRXI.GT.150.0D0) BRXI = 150.0D0
      IF(BRXI.LT.-150.0D0) BRXI = -150.0D0
      RO = 2.704D-3*DEXP(BRXI)
      WURZ = DSQRT(RO)*Y12
      S2MRO = 2.66D4*RO
      S2MROV = S2MRO*Y1
      CDV = S2MROV*(1.174D0-0.9D0*U2*U3)
      CL = S2MRO*0.6D0*U1*U3
      CLV = CL*Y1
      CDVV = CDV*Y1
      XI = 1.0D0/(1.0D0+Y3)
      Q1 = Y1*RK
      Q2 = 3.2172D-4*XI
      Q3 = Q2*Y1K
      Q4 = Q1-Q3
      SINY2 = DSIN(Y2)
      COSY2 = DCOS(Y2)
      SINXI = SINY2*XI
      COSXI = COSY2*XI
      D(1)=-Y7*(CDVV+Q2*SINXI)
      D(2)=Y7*(CLV+COSXI*Q4)
      D(3)=Y7*Q1*SINY2
      D(4)=Y7*(-30.0D0*WURZ+2.0D0*Y4*CDV-Y5*(CL+COSXI*(RK+Q3*Y1K))
     *-Y6*RK*SINY2)
      D(5)=Y7*(Y4*Q2*COSXI+Y5*SINXI*Q4-Y6*COSY2*Q1)
      D(6)=Y7*(5.0D0*BETAR*Y1*WURZ-Y4*(BETAR*CDVV+2.0D0*Q2*SINXI*
     *XI)+Y5*(BETAR*CLV+COSXI*XI*(Q4-Q3)))
      D(7)=0.0D0
      RETURN
      END
      SUBROUTINE R(YA,YB,W)
      IMPLICIT DOUBLEPRECISION(S)
      DOUBLE PRECISION YA(7),YB(7),W(7)
C     Boundary conditions for example Reentry
C:    End Parameter
      DOUBLE PRECISION BRXIA,BRXIB,CDVVA,CDVVB,CLVA,CLVB,COSXIA,
     *COSXIB,Q1A,Q2A,Q1B,Q2B,ROA,ROB,SINA,SINB,SINXIA,SINXIB,
     *S2MROA,S2MROB,U1A,U2A,U3A,U1B,U2B,U3B,XIA,XIB,YA1,YA12,YA2,
     *YA3,YB1,YB12,YB2,YB3
      DOUBLE PRECISION RK,BETAR,FIX
      COMMON /KONST/ RK,BETAR,FIX
C:    Begin
      YA1 = YA(1)
      YA2 = YA(2)
      YA3 = YA(3)
      YB1 = YB(1)
      YB2 = YB(2)
      YB3 = YB(3)
      U1A = 0.6D0*YA(5)
      U2A = 0.9D0*YA(4)*YA1
      U3A = -1.0D0/DSQRT(U1A*U1A+U2A*U2A)
      U1B = 0.6D0*YB(5)
      U2B = 0.9D0*YB(4)*YB1
      U3B = -1.0D0/DSQRT(U1B*U1B+U2B*U2B)
      BRXIA = -890.34D0*YA3
      BRXIB = -890.34D0*YB3
      IF(BRXIA.GT.150.0D0) BRXIA = 150.0D0
      IF(BRXIA.LT.-150.0D0) BRXIA = -150.0D0
      IF(BRXIB.GT.150.0D0) BRXIB = 150.0D0
      IF(BRXIB.LT.-150.0D0) BRXIB = -150.0D0
      ROA = 2.704D-3*DEXP(BRXIA)
      ROB = 2.704D-3*DEXP(BRXIB)
      XIA = 1.0D0/(1.0D0+YA3)
      XIB = 1.0D0/(1.0D0+YB3)
      SINA = DSIN(YA2)
      COSXIA = DCOS(YA2)*XIA
      SINXIA = SINA*XIA
      YA12 = YA1*YA1
      S2MROA = ROA*2.66D4
      CDVVA = S2MROA*YA12*(1.174D0-0.9D0*U2A*U3A)
      CLVA = S2MROA*YA1*0.6D0*U1A*U3A
      Q1A = 3.2172D-4*XIA
      Q2A = YA1/2.09D2
      SINB = DSIN(YB2)
      COSXIB = DCOS(YB2)*XIB
      SINXIB = SINB*XIB
      YB12 = YB1*YB1
      S2MROB = ROB*2.66D4
      CDVVB = S2MROB*YB12*(1.174D0-0.9D0*U2B*U3B)
      CLVB = S2MROB*YB1*0.6D0*U1B*U3B
      Q1B = 3.2172D-4*XIB
      Q2B = YB1/2.09D2
      W(1)=YB(1)-0.27D0
      W(2)=YB(2)
      W(3)=YB(3)-2.5D0/209.0D0
      W(4)=YA(1)-0.36D0
      W(5)=YA(2)+FIX
      W(6)=YA(3)-4.0D0/2.09D2
C     Hamilton - Funktion
      W(7)=10.0D0*YA12*YA1*DSQRT(ROA)-YA(4)*(CDVVA+Q1A*SINXIA)+YA
     *(5)*(CLVA+COSXIA*(Q2A-Q1A/YA1))+YA(6)*Q2A*SINA
      RETURN
      END
