************************************************************************
************************************************************************
************************************************************************
      SUBROUTINE BNGSIZ(D14R,D34R,D32R,ERROR)
C                       ^^^^^^^^^^^^^^^^^
C                              out


C     calculates the additional angular size distances (see
C     Kayser et al., A&A, 1996) for bounce models.

C     Given two redshifts Z1 and Z2 this subroutine calculates the
C     angular size distances between these two redshifts, which are
C     returned as D14, D34 and D32.  In bounce models there are
C     two or four independent distances for the same redshift; BNGSIZ
C     returns as D14 the distance between two redshifts Z1 and Z2 such
C     that the universe is NOT in a state of CONTRACTION at Z1 but is
C     NOT in a state of EXPANSION at Z2.  For D32 we have the reverse
C     situation; for D34 EXPANSION NEITHER at Z1 NOR at Z2.
C
C     If ETA(Z) is the same in the expanding and contracting phases,
C     then there are only two independent distances: D12 as calculated
C     by ANGSIZ and D14.  In this case, D34 = D12 and D32 = D14.
C
C     The integer variable ERROR is 0 when no error is encountered, and
C     is otherwise set to a specific value for the corresponding error.
C     This can be useful in debugging.
C
C     BNGSIZ can only be called after calling ANGSIZ.  If it is called
C     for a non-bounce model, and error message is printed.

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


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


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


C     error flag

      INTEGER ERROR

C     angular size distances---internal values

      DOUBLE PRECISION D14, D34, D32

C     angular size distances--- real output variables

      REAL D14R, D34R, D32R


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


C     all called routines

C     does the actual integration

C      CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D)

      EXTERNAL LOWLEV

C     constants

      DOUBLE PRECISION NSTEP

C     One really needs a value this big to avoid numerical instability.
C     Nothing explodes or crashes with a too small value, but the result
C     can be incorrect.  Although it is difficult to independently check
C     the distances calculated by BNGSIZ, they're OK as far as I can
C     tell and this value for NSTEP also makes the distance continuous
C     smooth and smooth, which means values between checked values are
C     probably OK.

      PARAMETER (NSTEP=200.0D0)

C     variables

C     (1 + Z1) and (1 + Z2)

      DOUBLE PRECISION VD, VS

C     dependent variables (starting values) for the integration routine

      DOUBLE PRECISION YSTART(2)



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


C      CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D)
C      G = SQRT(QQ(VD))

C     angular size distance and Q(v)**2

      DOUBLE PRECISION D, QQ

      EXTERNAL QQ



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


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

      INTEGER IERROR, WMTYPE

      DOUBLE PRECISION ZMAX

      LOGICAL NOTYET, DEBUG, VARETA, MAXZ, BOUNCE

      COMMON /ANGINI/ ZMAX, IERROR, WMTYPE,
     $                NOTYET, VARETA, DEBUG, MAXZ, BOUNCE

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

      DOUBLE PRECISION Z1, Z2, D12

      COMMON /BNCSIZ/ Z1, Z2, D12

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

C     YSTART(2) declared above

      COMMON /LOWSIZ/ YSTART

C     communication between the calling routine, here and ETAZ (q. v.)

C     CHOICE is needed here, and not only in the user's calling routine
C     and ETA(Z), since a CHOICE value < 0 indicates that the dependence
C     of ETA on Z is different in the phases of expansion and
C     contraction in a bounce model.

      INTEGER CHOICE

      COMMON /WHICH/ CHOICE

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

      LOGICAL VALUE2

      COMMON /CNTRCT/ VALUE2

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

      INTEGER MSTEP

      COMMON /MERROR/ MSTEP

C     save until the next call

      SAVE /ANGINI/, /BNCSIZ/, /LOWSIZ/, /WHICH/, /CNTRCT/, /MERROR/



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


C     constants

C     just to check for overflow

      DOUBLE PRECISION TINY

C     This is roughly the smallest REAL number for a 4 byte REAL;
C     see the comment in ANGSIZ.

      PARAMETER (TINY=1.0D-38)

C     variables

C     just a temporary work variable

      DOUBLE PRECISION G



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


C     ERROR will be changed if something goes wrong

      ERROR = 0

C     is a calculation possible

C     based on the cosmological model?

C     Changing the second condition to (IERROR .LE. 2) will suppress the
C     ERROR = 2 message here (though not in INICOS, of course).
C     Reporting is normally enabled, because the result is not
C     necessarily useful in this case, but the knowledgeable user might
C     wish to have this behaviour.

      IF ((IERROR .GE. 1) .AND. (IERROR .LE. 3)) THEN

        IF (DEBUG) THEN

            WRITE (*, '(T2,A,I1,A)') 'ERROR ', IERROR, ': INICOS ERROR'

         END IF

      END IF

C     is this a bounce model?

      IF (.NOT. BOUNCE) THEN

         IF (DEBUG) THEN

            WRITE (*,'(T2,A,I2)')
     $            'ERROR 13: BNGSIZ INVALID FOR WMTYPE ', WMTYPE

         END IF

         IERROR = 13

      END IF

C     ANGSIZ OK?

      IF ( ( (IERROR .GE. 4) .AND. (IERROR .LE. 10) ) .OR.
     $     (                        IERROR .EQ. 12  )      ) THEN

        IF (DEBUG) THEN

            WRITE (*, '(T2,A,I1,A)') 'ERROR ', IERROR, ': ANGSIZ ERROR'

         END IF

      END IF

      ERROR = IERROR

C     no need to check the redshifts, since these were checked by ANGSIZ

C     Add .AND. (ERROR .NE. 3) if you really want it (see ANGSIZ).

      IF ( (IERROR .NE. 0) .AND. (IERROR .NE. 11) ) THEN

************************************************************************
C        abnormal return: no further calculation meaningful

         D14 = 0.0D0

         D34 = 0.0D0

         D32 = 0.0D0

         RETURN

************************************************************************

      END IF

************************************************************************
C     if possible, calculate the angular size distance


C     calculate first to ZMAX...

C     we are in the expanding phase now; set a flag in case ETAZ
C     depends on this

      VALUE2 = .FALSE.

      IF (Z2 .GE. Z1) THEN

C     new initial values

         VD = 1.0D0 + Z2

         VS = 1.0D0 + ZMAX

C        but the the initial values YSTART
C        are the result of the previous call to ANGSIZ

      ELSE

C        different new starting points, since we are integrating
C        starting at Z1

         VD = 1.0D0 + Z1

         VS = 1.0D0 + ZMAX

C        evaluate the stuff under the root (Q(V))

         G = SQRT(QQ(VD))

         IF (G .LE. TINY) THEN

            ERROR = 11

            IF (DEBUG) THEN

               PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE'

            END IF

         END IF

C        and thus have the `normal' boundary conditions

C        initial values needed for the differential equation solver

         YSTART(1) = 0.0D0

C        see Kayser et al., Eq. (36) for an explanation for the SIGN

C        YSTART(2) = SIGN((1.0D0/(VD*G)),(VS-VD))

         YSTART(2) = 1.D00/(VD*G)

C        since (ZMAX .GE. Z1) and thus ((VS - VD) .GE. 0.0D0)

      END IF

      CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D)

      IF (MSTEP .NE. 0) THEN

         ERROR = 12

         IF (DEBUG) THEN

            PRINT*, 'ERROR 12: ERROR CALCULATING D14, EXPANDING PHASE'

         END IF

************************************************************************
C        abnormal return: no further calculation if the
C        integration fails

         D14 = 0.0D0

         D34 = 0.0D0

         D32 = 0.0D0

         RETURN

************************************************************************

      END IF

C     ...and then back to Z2

C     we are in the contracting phase now; set a flag in case ETAZ
C     depends on this

      VALUE2 = .TRUE.

      VD = 1.0D0 + ZMAX

      VS = 1.0D0 + Z2

C     call the solver for ZMAX--VS, the initial values YSTART
C     are the result of the previous call; the minus sign is needed
C     because we are reversing the integration direction

      YSTART(2) = -YSTART(2)

      CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D)

      IF (MSTEP .NE. 0) THEN

         ERROR = 12

         IF (DEBUG) THEN

            PRINT*, 'ERROR 12: ERROR CALCULATING D14, CONTRACTING PHASE'

         END IF

************************************************************************
C        abnormal return: no further calculation if the
C        integration fails

         D14 = 0.0D0

         D34 = 0.0D0

         D32 = 0.0D0

         RETURN

************************************************************************

      END IF

C     this is our result for the second distance in bounce models

      D14 = D



************************************************************************
C     We're finished, unless ETA(Z) is different during contraction.
C     If ETA depends on Z only, that is, is the same in expanding
C     and contracting phases, then there are only two independent
C     distances, which we have already calculated.  Otherwise, this
C     degeneracy vanishes, and we must calculate the other two.


C     Remember that a CHOICE value < 0 indicates that the dependence
C     of ETA on Z is different in the phases of expansion and
C     contraction in a bounce model.

      IF ( .NOT. ((CHOICE .LT. 0) .AND. VARETA) ) THEN

         D34 = D12

         D32 = D14

      ELSE

************************************************************************
C        the third distance

C        we are in the contracting phase now; set a flag because ETAZ
C        depends on this

         VALUE2 = .TRUE.

C        some auxiliary variables for run time (not memory) optimisation

         VD = 1.0D0 + Z1

         VS = 1.0D0 + Z2

C        evaluate the stuff under the root (Q(V))

         G = SQRT(QQ(VD))

         IF (G .LE. TINY) THEN

            ERROR = 11

            IF (DEBUG) THEN

               PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE'

            END IF

         END IF

C        initial values needed for the differential equation solver

         YSTART(1) = 0.0D0

C        see Kayser et al., Eq. (36) for an explanation for the SIGN

         YSTART(2) = SIGN((1.0D0/(VD*G)),(VS-VD))

         CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D)

         IF (MSTEP .NE. 0) THEN

            ERROR = 12

         IF (DEBUG) THEN

            PRINT*, 'ERROR 12: ERROR CALCULATING D34, CONTRACTING PHASE'

         END IF

************************************************************************
C           abnormal return: no further calculation if the
C           integration fails

            D14 = 0.0D0

            D34 = 0.0D0

            D32 = 0.0D0

            RETURN

************************************************************************

         END IF

C        this is our result for the third distance in bounce models

         D34 = D

************************************************************************
C     now for the fourth distance

C     calculate first to ZMAX...

C     we are in the contracting phase now; set a flag because ETAZ
C     depends on this

      VALUE2 = .TRUE.

      IF (Z2 .GE. Z1) THEN

C        new initial values

         VD = 1.0D0 + Z2

         VS = 1.0D0 + ZMAX

C        but the the initial values YSTART
C        are the result of the previous call

      ELSE

C        different new starting points, since we are integrating
C        starting at Z1

         VD = 1.0D0 + Z1

         VS = 1.0D0 + ZMAX

C        evaluate the stuff under the root (Q(V))

         G = SQRT(QQ(VD))

         IF (G .LE. TINY) THEN

            ERROR = 11

            IF (DEBUG) THEN

               PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE'

            END IF

         END IF

C        and thus have the `normal' boundary conditions

C        initial values needed for the differential equation solver

         YSTART(1) = 0.0D0

C        see Kayser et al., Eq. (36) for an explanation for the SIGN

C        YSTART(2) = SIGN((1.0D0/(VD*G)),(VS-VD))

         YSTART(2) = 1.0D0/(VD*G)

C        since (ZMAX .GE. Z1) and thus ((VS - VD) .GE. 0.0D0)

      END IF

      CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D)

      IF (MSTEP .NE. 0) THEN

         ERROR = 12

         IF (DEBUG) THEN

            PRINT*, 'ERROR 12: ERROR CALCULATING D32, CONTRACTING PHASE'

         END IF

************************************************************************
C        abnormal return: no further calculation if the
C        integration fails

         D14 = 0.0D0

         D34 = 0.0D0

         D32 = 0.0D0

         RETURN

************************************************************************

      END IF

C     ...and then back to Z2

C     we are in the expanding phase now; set a flag because ETAZ
C     depends on this

      VALUE2 = .FALSE.

      VD = 1.0D0 + ZMAX

      VS = 1.0D0 + Z2

C     call the solver for ZMAX--VS, the initial values YSTART
C     are the result of the previous call

      CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D)

      IF (MSTEP .NE. 0) THEN

         ERROR = 12

         IF (DEBUG) THEN

            PRINT*, 'ERROR 12: ERROR CALCULATING D32, EXPANDING PHASE'

         END IF

************************************************************************
C        abnormal return: no further calculation if the
C        integration fails

         D14 = 0.0D0

         D34 = 0.0D0

         D32 = 0.0D0

         RETURN

************************************************************************

      END IF

C     this is our result for the fourth distance in bounce models

      D32 = D

      END IF



************************************************************************
C      set the output variables


      D14R = REAL(D14)

      D34R = REAL(D34)

      D32R = REAL(D32)



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