************************************************************************
************************************************************************
************************************************************************
      SUBROUTINE LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D)


C     This does the low level stuff for *NGSIZ, which is calculating the
C     angular size distance (and its derivative) between VD and VS.
C     This routine just does one leg of the integration; this is OK for
C     non-BOUNCE models; for BOUNCE models the second distance is put
C     together from separate calls to LOWLEV by *NGSIZ.  DEBUG is passed
C     to enable informational messages to be printed by LOWLEV, since
C     this would mean redundant code in *NGSIZ.

C     NSTEP determines the first guessed step size.  For the additional
C     distances in the BOUNCE case this must be much smaller for
C     numerical reasons.

C     YSTART are the starting values, which are replaced with the
C     results (needed when integrating further in the BOUNCE case).

C     Code written by Phillip Helbig, Hamburg Observatory, Sep 1995.


************************************************************************
************************************************************************
C     Declarations


************************************************************************
C     input variables from the calling routine (parameter list)


C     start and end values

      DOUBLE PRECISION VD, VS, NSTEP

      LOGICAL DEBUG



************************************************************************
C     output variables to the calling routine (parameter list)


C     the distance

      DOUBLE PRECISION D



************************************************************************
C     output to called routines (parameter list)


C     all called routines

      EXTERNAL ODEINT, ASDRHS, BSSTEP

C         CALL ODEINT(2,VD,VS,EPS,H1,0.0,BSSTEP,ASDRHS,
C     $               YSTART, NOK,NBAD)


C     2 is the dimension of YSTART
C     VD is the start value for the integration
C     VS is the  end  value for the integration
C     EPS is desired accuracy
C     H1 is a guessed first step size
C     the allowed minimum step size is 0.0
C     BSSTEP is a stepper routine
C     ASDRHS calculates the right hand sides of the differential equations

C     YSTART are the starting values; the end values replace these

C     constants in the parameter list

C     accuracy

      DOUBLE PRECISION EPS

      PARAMETER (EPS=1.0D-6)

C     start values, guessed first step size

      DOUBLE PRECISION YSTART(2), H1



************************************************************************
C     output from called routines (parameter list)


C     variables in the parameter list

C     NOK is the number of successful tries (not used in this version)
C     NBAD is the number of failed tries (not used in this version)

      INTEGER NOK, NBAD



************************************************************************
C     variables in common blocks


C     communication between here and BNGSIZ (q. v.)

      DOUBLE PRECISION YYSTRT(2)

      COMMON /LOWSIZ/ YYSTRT

C     communication between here and numerical integration stuff (q.v.)

      INTEGER MSTEP

      COMMON /MERROR/ MSTEP

C     save until the next call

      SAVE /LOWSIZ/, /MERROR/



************************************************************************
C     internally used quantities


C     variables

C     for the do loop

      INTEGER I



************************************************************************
************************************************************************
C     run-time initialisation


C     no error yet

      MSTEP = 0



************************************************************************
************************************************************************
C     calculations


C     for experimentally determined special cases

      IF( ( (VD. LE. 2.0D0) .OR. (VS .LT. (1.0D0 + EPS)) ) .AND.
     $    (ABS(VS-VD) .GE. 3.6D0)                              )  THEN

C        set the guessed first step size

         H1 = 2.0D0

C     in other cases

      ELSE

C     use a simply calculated guessed first step size

         H1 = (VS-VD)/NSTEP

      END IF

C     if D12 is practically = 0

      IF (ABS(H1) .LT. EPS/10.0D0) THEN

C         set it to zero to avoid numerical problems

          D = 0.0D0

C     otherwise calculate the angular size distance

      ELSE

C        call the solver

         CALL ODEINT(2,VD,VS,EPS,H1,0.0D0,BSSTEP,ASDRHS,
     $               YSTART, NOK,NBAD)

         IF (MSTEP .NE. 0) THEN

            IF (DEBUG) THEN

               IF (MSTEP .EQ. 1) THEN

                  PRINT*, 'ERROR IN BSSTEP'

               ELSE IF (MSTEP .EQ. 2) THEN

                  PRINT*, 'ERROR IN ODEINT: STEPSIZE TOO SMALL'

               ELSE IF (MSTEP .EQ. 3) THEN

                  PRINT*, 'ERROR IN ODEINT: TOO MANY STEPS'

               ELSE

                  PRINT*, 'CAN''T GET HERE'

               END IF

            END IF

            RETURN

         END IF

C        this is our result

         D = YSTART(1)

      END IF



************************************************************************
C    set the variables needed by BNGSIZ


      DO 100, I = 1, 2, 1

         YYSTRT(I) = YSTART(I)

100   CONTINUE



************************************************************************
************************************************************************
      END