************************************************************************
************************************************************************
************************************************************************
      SUBROUTINE ANGSIZ(Z1R,Z2R, D12R,ERROR)
C                       ^^^^^^^  ^^^^^^^^^^
C                         in        out


C     calculates the angular size distance between two redshifts (see
C     Kayser et al., A&A, 1996).

C     Given two redshifts Z1 and Z2 this subroutine calculates the
C     angular size distance between these two redshifts, which is
C     returned as D12.  In some world models (bounce models) there are
C     two (or four) independent distances for the same redshift; ANGSIZ
C     returns as D12 the distance between two redshifts Z1 and Z2 such
C     that the universe is NOT in a state of CONTRACTION---neither at
C     the time of emission nor at the time of reception of the observed
C     radiation.  (The additional distances in bounce models are
C     calculated by the SUBROUTINE BNGSIZ.)
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     If the cosmological model is in error, the corresponding value for
C     ERROR is returned.  Otherwise, the redshifts are examined, to see
C     if the distance is defined.
C
C     Undefined distances (producing an error) are those where at least
C     one redshift is < 0 or > ZMAX, where ZMAX = \infty in big bang
C     models and a finite value in other models.
C
C     An error message is also returned if ANGSIZ is called before
C     calling INICOS.  Of course, if the input parameters to INICOS are
C     changed after a call, it is the responsibility of the user to call
C     INICOS again before calling ANGSIZ.
C
C     If indicated in the INICOS call, error messages, when required,
C     are printed to standard output.

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


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


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


C     redshifts

C     the real input variables

      REAL Z1R, Z2R

C     the internally used counterparts

      DOUBLE PRECISION Z1, Z2



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


C     error flag

      INTEGER ERROR

C     angular size distance (D14, D34 and D32 are the other distances
C     in bounce models)

C     the actual output variable

      REAL D12R

C     the internally used counterpart

      DOUBLE PRECISION D12



************************************************************************
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     This seems to work fine for ANGSIZ; the exact value isn't too
C     important, as in BNGSIZ.  This is about 25% faster than using
C     the BNGSIZ value of 200.

      PARAMETER (NSTEP=2.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 BNGSIZ (q. v.)

      DOUBLE PRECISION ZZ1, ZZ2, DD12

      COMMON /BNCSIZ/ ZZ1, ZZ2, DD12

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     to allow computation in the Eddington model, if desired

      LOGICAL RISK

      COMMON /ARTHUR/ RISK

C     sometimes ANGSIZ could be called with a Z which is smaller than
C     ZMAX but large enough for numerical problems.

      DOUBLE PRECISION ZZREAL

      COMMON /FIX/ ZZREAL

C     save until the next call

      SAVE /ANGINI/, /BNCSIZ/, /CNTRCT/, /MERROR/, /ARTHUR/, /FIX/



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


C     constants

C     just to check for overflow

      DOUBLE PRECISION TINY

      CHARACTER TEXT17*56

C     This is roughly the smallest REAL number for a 4 byte REAL.
C     While internally higher precision is used, the input and output
C     variables are REAL and thus this is a better choice.  'higher
C     precision' instead of the corresponding FORTRAN key words is used
C     here to enable the conversion to single precision via a global
C     search and replace, if so desired.  Memory requirements of these
C     routines is minimal; whether single or higher precision is faster
C     depends on the machine.  All internal quantities are higher
C     precision to avoid conversions.  Higher precision is not required
C     if one does not use BNGSIZ, but is used throughout because memory
C     requirements are minimal and it is even faster on some machines.
C     If single precision is desired, replacing the FORTRAN key words
C     for higher precision with REAL in all routines will do the trick
C     and can be done automatically.  (Additionally, `D0' should be
C     removed.

      PARAMETER (TINY=1.0D-38,
     $TEXT17='NOT REALLY LARGER THAN ZMAX BUT CLOSE ENOUGH FOR TROUBLE')

C     variables

C     just a temporary work variable

      DOUBLE PRECISION G



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


      Z1 = DBLE(Z1R)

      Z2 = DBLE(Z2R)



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


C     is a calculation of the angular size distance possible?

C     ERROR will be changed if something goes wrong

      ERROR = 0

C     has INICOS been called first?

      IF (NOTYET) THEN

         PRINT*, 'YOU MUST CALL INICOS FIRST!!!!!!!!!!!!!!!!!!!!!!!!!!!'

************************************************************************
C        abnormal return: no further calculation in this case

         RETURN

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

      END IF

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

C     Don't allow computation in the Eddington model, WMTYPE 17, unless
C     the user explicitly requests this.

      ELSE IF ((WMTYPE .EQ. 17) .AND. .NOT. RISK) THEN

         ERROR = 4

         IF (DEBUG) THEN

            PRINT*, 'ERROR 4: COMPUTATION NOT RECOMMENDED IN WMTYPE 17'

         END IF

C     based on the redshifts?

      ELSE

         IF (MAXZ) THEN

            IF      (Z2 .GT. ZMAX) THEN

               IF (Z2 .GT. ZZREAL) THEN

                  ERROR = 9

                  IF (DEBUG) PRINT*, 'ERROR 9: Z2 > ZMAX'

               ELSE

                  ERROR = 10

                  IF (DEBUG) WRITE (*,'(2A)') 'Z2 ', TEXT17

               END IF

            ELSE IF (Z1 .GT. ZMAX) THEN

               IF (Z1 .GT. ZZREAL) THEN

                  ERROR = 7

                  IF (DEBUG) PRINT*, 'ERROR 7: Z1 > ZMAX'

               ELSE

                  ERROR = 8

                  IF (DEBUG) WRITE (*,'(2A)') 'Z1 ', TEXT17

               END IF

            END IF

         END IF

         IF (Z1 .LT. 0.0) THEN

            ERROR = 5

            IF (DEBUG) THEN

               PRINT*, 'ERROR 5: Z1 < 0.0'

            END IF

         ELSE IF (Z2 .LT. 0.0) THEN

            ERROR = 6

            IF (DEBUG) THEN

               PRINT*, 'ERROR 6: Z2 < 0.0'

            END IF

         END IF

      END IF



***********************************************************************
C     abnormal exit


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

      IF (ERROR .NE. 0) THEN

         D12 = 0.0D0

         IERROR = ERROR

         RETURN

      END IF



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


***********************************************************************
C     the distance in the expanding universe


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

      VALUE2 = .FALSE.

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

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

         D12 = 0.0D0

         IERROR = ERROR

         RETURN

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

      END IF

      D12 = D

C     set the output variable

      D12R = REAL(D12)

C     set the variables needed by BNGSIZ

      ZZ1 = Z1

      ZZ2 = Z2

      DD12 = D12



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