ISAJET/utils/cern_lib/rkstp.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.9  2000/07/25 14:53:05  mclareni
   6 * Version 7.51 from author
   7 *
   8 *
   9 *#include "sys/CERNLIB_machine.h"
  10 #include "isajet/pilot.h"
  11 #if defined(CERNLIB_NOCERN)
  12 C-----------------------------------------------------------------------
  13       SUBROUTINE RKSTP(N,H,X,Y,SUB,W)
  14 C-----------------------------------------------------------------------
  15 C
  16 C     From CERN Program Library, routine D209, with error message for
  17 C     N.LT.1 replaced by STOP 99 to eliminate Kernlib error routine.
  18 C
  19       DIMENSION Y(N),W(N,3)
  20       LOGICAL MFLAG,RFLAG
  21       EXTERNAL SUB
  22 C
  23 C     ******************************************************************
  24 C
  25 C     THIS SUBROUTINE REPLACES X BY X+H AND ADVANCES THE SOLUTION OF THE
  26 C     SYSTEM OF DIFFERENTIAL EQUATIONS DY/DX=F(X,Y) FROM Y(X) TO Y(X+H)
  27 C     USING A FIFTH-ORDER RUNGE-KUTTA METHOD.
  28 C
  29 C     SUB IS THE NAME OF A SUBROUTINE SUB(X,Y,F) WHICH SETS THE VECTOR F
  30 C     TO THE DERIVATIVE AT X OF THE VECTOR Y.
  31 C
  32 C     W IS A WORKING-SPACE ARRAY, TREATED AS CONSISTING OF THREE CONSEC-
  33 C     UTIVE WORKING VECTORS OF LENGTH N.
  34 C
  35 C     ******************************************************************
  36 C
  37 C  START.
  38       IF (N.LT.1) STOP 99
  39       NLOCAL=N
  40       HLOCAL=H
  41       H2=0.5*HLOCAL
  42       H6=HLOCAL/6.
  43       XH=X+HLOCAL
  44       XH2=X+H2
  45       CALL SUB(X,Y,W(1,1))
  46       DO 1 J=1,NLOCAL
  47          W(J,2)=Y(J)+H2*W(J,1)
  48     1 CONTINUE
  49       CALL SUB(XH2,W(1,2),W(1,3))
  50       DO 2 J=1,NLOCAL
  51          W(J,1)=W(J,1)+2.*W(J,3)
  52          W(J,2)=Y(J)+H2*W(J,3)
  53     2 CONTINUE
  54       CALL SUB(XH2,W(1,2),W(1,3))
  55       DO 3 J=1,NLOCAL
  56          W(J,1)=W(J,1)+2.*W(J,3)
  57          W(J,2)=Y(J)+HLOCAL*W(J,3)
  58     3 CONTINUE
  59       CALL SUB(XH,W(1,2),W(1,3))
  60       DO 4 J=1,NLOCAL
  61          Y(J)=Y(J)+H6*(W(J,1)+W(J,3))
  62     4 CONTINUE
  63       X=XH
  64       RETURN
  65       END
  66 #endif