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