program DRSIVA c>> 2010-06-09 DRSIVA Krogh Used parameters for all dimenssions. c>> 2001-05-25 DRSIVA Krogh Minor change for making .f90 version. c>> 1996-06-14 DRSIVA Krogh Small change in output format c>> 1994-11-02 DRSIVA Krogh Changes to use M77CON c>> 1994-07-18 DRSIVA Krogh Last change. c--S replaces "?": DR?IVA, ?IVA, ?IVAF, ?IVAO c Sample driver for SIVA -- Set up to solve two second order equations. c integer INEQ, IFDIM, IKDIM, ITDIM, IYDIM parameter (INEQ=2,IFDIM=16*INEQ+1,IKDIM=6,ITDIM=4,IYDIM=4*INEQ) integer NEQ, KORD(IKDIM), IOPT(6) c--S Next line special: P=>D, X=>Q real TSPECS(ITDIM), Y(IYDIM), T, H, DELT, TFINAL integer NDIG real F(IFDIM), TOL equivalence (TSPECS(1), T), (TSPECS(2), H), (TSPECS(3), DELT), 1 (TSPECS(4), TFINAL) external SIVAO, SIVAF c++SP Default NDIG = 4 c++ Default NDIG = 10 c++ Substitute for NDIG below parameter (NDIG = 4 ) parameter (TOL = 10.E0 **(-NDIG)) c data NEQ, T, H, DELT, TFINAL / 1 2, 0.E0, 1.E0, 6.283185307179586477E0, 2.E1 /, 2 Y(1), Y(2), Y(3), Y(4) / 3 1.E0, 0.E0, 0.E0, 1.E0 / c c Set option for error control, local absolute error < TOL. data IOPT(1), IOPT(2), IOPT(3) / 1 16, 6, 3 /, c Group the system to be treated as a single unit, set tolerance value 2 KORD(6), F(3)/ 3 2, TOL/ c Set option for second order equations data IOPT(4), IOPT(5) / 1 17, 2 / c Flag end of options data IOPT(6) / 0/ c c Do the integration c KORD(1) = 0 100 continue call SIVA(TSPECS, Y, F, KORD, NEQ, SIVAF, SIVAO, 1 ITDIM, IYDIM, IFDIM, IKDIM, IOPT) if (KORD(1) .NE. 1) go to 100 stop end subroutine SIVAF(T, Y, F, KORD) c Sample derivative subroutine for use with SIVA c This evaluates derivatives for a simple two body problem. c integer KORD c--S Next line special: P=>D, X=>Q real T, Y(4), TP real F(2) c c Evaluate the derivatives c TP = Y(1)*Y(1) + Y(3)*Y(3) TP = 1.E0 / (TP * SQRT(TP)) F(1) = -Y(1) * TP F(2) = -Y(3) * TP return end subroutine SIVAO(TSPECS, Y, F, KORD) c Sample output subroutine for use with SIVA. c This subroutine gives output for a simple two body problem. c integer KORD c--S Next line special: P=>D, X=>Q real TSPECS(4), Y(4) real F(2) 1000 format (12X, 1 'RESULTS FOR A SIMPLE 2-BODY PROBLEM'// 2 8X, 'T', 13X, 'U/V', 11X, 'UP/VP', 9X, 'UPP/VPP') 1001 format (1P,SP,4E15.6 / 15X, 3E15.6/' ') c c Do the output c if (KORD .EQ. 1) then write (*, 1000) end if write (*, 1001) TSPECS(1), Y(1), Y(2), F(1), Y(3), Y(4), F(2) c return end