************************************************************************ ************************************************************************ ************************************************************************ SUBROUTINE INICOS(LAMBDR,OMEGAR,ETAR,VARETA,DEBUG, C ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ C in C $ WMTYPE,MAXZ,ZMAXR,BOUNCE,ERROR) C ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ C out C checks the cosmological parameters, relays user wishes to ANGSIZ, C calculates the maximum possible redshift, indicates errors, is to C be called as initialisation before calling ANGSIZ, and gives some C indication as to the type of cosmological model. C LAMBDA and OMEGA are the usual $\lambda_{0}$ and $\Omega_{0}$, ETA C gives the fraction of homogeneously distributed matter. VARETA, C if .TRUE., indicates that ETA is not to be determined by the C (z-independent) value of ETA, but rather through the C REAL FUNCTION ETAZ(Z), which the user can modify as needed. C DEBUG, if .TRUE., enables the reporting of error messages. C C WMTYPE gives an indication as to the cosmological model, which can C be used in conjuction with the table in the README.TEX file. C BOUNCE tells whether or not the world model is a bounce model; in C these cases, ZMAX is finite and BOUNCE and MAXZ are .TRUE. For C the A2 models, BOUNCE is false, but ZMAX is also finite, and thus C MAXZ is .TRUE. (If ZMAX is analytically infinite, it is C computationally undefined, and MAXZ is .FALSE..) If an error C occurs, then ERROR is suitably defined, whether or not DEBUG is C .TRUE. C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declarations ************************************************************************ C input variables from the calling routine (parameter list) C the cosmological parameters---input variables REAL LAMBDR, OMEGAR, ETAR C the cosmological parameters---internal counterparts DOUBLE PRECISION LAMBDA, OMEGA, ETA C customised ETA? display messages? LOGICAL VARETA, DEBUG ************************************************************************ C output variables to the calling routine (parameter list) C type of cosmological model, error message INTEGER WMTYPE, ERROR C maximum redshift C real output variable REAL ZMAXR C internal counterpart DOUBLE PRECISION ZMAX C is there a maximum redshift? is this a bounce model? LOGICAL MAXZ, BOUNCE ************************************************************************ C output from called routines (parameter list) C this is the quantity Q(v)**2 (see Kayser et al., A&A, 1996) DOUBLE PRECISION QQ EXTERNAL QQ ************************************************************************ C variables in common blocks C with different names than elsewhere since the same quantities C appear in a parameter list as well; additional variables are C explained C SUBROUTINE INICOS(LAMBDA,OMEGA,ETA,VARETA,DEBUG, C $ WMTYPE,MAXZ,ZMAX,BOUNCE,ERROR) C K gives us a measure for the curvature (see below) DOUBLE PRECISION LAMNUL, OMENUL, ETAETA, K COMMON /COSMOL/ LAMNUL, OMENUL, ETAETA, K INTEGER IERROR, TYPEWM DOUBLE PRECISION ZZMAX C have we called anything yet? LOGICAL NOTYET, ETAVAR, DEEBUG, MMAXZ, BOUNCY COMMON /ANGINI/ ZZMAX, IERROR, TYPEWM, $ NOTYET, ETAVAR, DEEBUG, MMAXZ, BOUNCY 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 C in order to avoid recalculating cosmological quantities when this C is unnecessary SAVE /COSMOL/, /ANGINI/, /FIX/ ************************************************************************ C internally used quantities C named constants C $\frac{\pi}{3}$, accuracy DOUBLE PRECISION PITHRD, EPS C The value of EPS should be about right for a 24 bit mantissa; C see the note concerning precision in ANGSIZ. PARAMETER (PITHRD = 1.0471975511965977D0, EPS = 1.0D-6) C variables C accuracy monitor INTEGER NBAD C internal work variables (SIGMA and Q have the usual cosmological C meaning) DOUBLE PRECISION SIGMA, Q, T, R, PHI, Y1, Y2, Y3, TEMP ************************************************************************ ************************************************************************ C run-time initialisation C set the input variables LAMBDA = DBLE(LAMBDR) OMEGA = DBLE(OMEGAR) ETA = DBLE(ETAR) ************************************************************************ ************************************************************************ C calculations ************************************************************************ C check the input parameters IF (OMEGA .LT. 0.0) THEN IERROR = 1 IF (DEBUG) THEN PRINT*, 'INICOS ERROR 1: OMEGA < 0.0; OMEGA = ', OMEGA END IF ELSE IF ( (.NOT. VARETA) .AND. $ ((ETA .LT. 0.0) .OR. (ETA .GT. 1.0)) ) THEN IERROR = 2 IF (DEBUG) THEN PRINT*, 'INICOS ERROR 2: ETA OUT OF RANGE; ETA = ', ETA END IF ELSE IERROR = 0 END IF ERROR = IERROR ************************************************************************ C abnormal exit IF (ERROR .GT. 0) RETURN ************************************************************************ C calculate WMTYPE, BOUNCE, MAXZ and ZMAX, if necessary C has OMEGA or LAMBDA changed since the last call? .OR. C is this the first call? IF ( (LAMBDA .NE. LAMNUL) .OR. (OMEGA .NE. OMENUL) .OR. $ NOTYET ) THEN C this only has an effect the first time NOTYET = .FALSE. C the sign of k is the sign of the curvature K = OMEGA + LAMBDA - 1.0D0 IF (LAMBDA .LT. 0.0) THEN MMAXZ = .FALSE. BOUNCY = .FALSE. IF (K .LT. 0.0) THEN IF (OMEGA .GT. 0.0) THEN TYPEWM = 2 ELSE C empty universe TYPEWM = 1 END IF ELSE IF (K .GT. 0.0) THEN TYPEWM = 4 ELSE C k = 0 TYPEWM = 3 END IF ELSE IF (LAMBDA .EQ. 0.0) THEN MMAXZ = .FALSE. BOUNCY = .FALSE. IF (K .LT. 0.0) THEN IF (OMEGA .GT. 0.0) THEN TYPEWM = 8 ELSE C empty universe, `Milne-model' TYPEWM = 7 END IF ELSE IF (K .GT. 0.0) THEN C `Misner-Thorne-Wheeler model' TYPEWM = 5 ELSE C k = 0, Einstein-de Sitter model TYPEWM = 9 END IF ELSE C LAMBDA > 0.0 IF (K .LT. 0.0) THEN MMAXZ = .FALSE. BOUNCY = .FALSE. IF (OMEGA .GT. 0.0) THEN TYPEWM = 11 ELSE C empty universe TYPEWM = 10 END IF ELSE IF (K .EQ. 0.0) THEN MMAXZ = .FALSE. BOUNCY = .FALSE. IF (OMEGA .GT. 0.0) THEN TYPEWM = 13 ELSE C empty universe, de Sitter model TYPEWM = 12 END IF ELSE C k > 0 IF (OMEGA .EQ. 0.0) THEN C empty universe, Lanczos model C Near the de Sitter model (TYPEWM = 12) MMAXZ is C very large, so practically no value of z will be C rejected. However, this can lead to numerical C problems, since below QQ(MMAXZ) is taken to make sure C it is large enough. (Negative values are bad because C a square root is taken via ANGSIZ; small positive C values are bad because the reciprocal is used by C ANGSIZ and this can lead to overflow.) We therefore C declare the world model to be practically de Sitter C if it is practically (numerically) flat. In day to C day practice this should lead to no error. IF (ABS(K) .LT. EPS) THEN TYPEWM = 12 MMAXZ = .FALSE. BOUNCY = .FALSE. ELSE C normal code for this case TYPEWM = 19 MMAXZ = .TRUE. BOUNCY = .TRUE. ZZMAX = SQRT(LAMBDA/(LAMBDA - 1.0D0)) - 1.0D0 END IF ELSE C sigma and q are mathematically and computationally useful SIGMA = 0.5D0*OMEGA Q = SIGMA - LAMBDA C Stabell-Refsdal equation (10) T =-K+3.0D0*SIGMA**(2.0D0/3.0D0)*LAMBDA**(1.0D0/3.0D0) C T = 0 defines the A1 and A2 curves in Stabell and C Refsdal. On the VAX, if we set EPS = 0, there is no C problem; when T is analytically = 0 (for example, C LAMBDA = 2.0 and OMEGA = 0.5) then it is also C numerically = 0, and there is no conflict (at least C when LAMBDA and OMEGA are represented exactly, as in C the example above). Other machines, for instance IBM C RS6000, calculate a T of about 1.0E-15 in this case, C due probably to filling up the multiplicands with 32 C bits of zeros in the registers. Now, normally, this C would make no difference, since for the purposes of C determining the world model, we are interested in C whether or not a particular criterion is NUMERICALLY C fulfilled. However, when T is very small, one is C very close to the A curve, and this can give problems C in the calculation of the angular size distance, since C it can become arbitrarily large for z near ZMAX. C Therefore, it is a good idea to draw the A curve with C a blunt numerical pencil. However, the fudge factor C is very small, and should be of no consequence for C the practical cosmologist. (One might think that the C fudge factor should be increased to exclude models C which are ANALYTICALLY and numerically near the A C curve, since the angular size distance can become very C large in these cases as well. However, tests show C that machine precision tends to be coarse enough so C that if one is analytically off the curve, one is far C enough away so as to encounter no problems later on. C In fact, initially EPS/100000.0 was OK on the VAX/VMS C (naturally), ALPHA/VMS, SUN, Intel/Linux and IBM C RS/6000/AIX for the pwr and pwr2 architectures. C However, optimising for the ppc architecture required C a reduction to EPS/10.0. To be on the safe side, we C use EPS here---as in other places in this package--- C which is about what one would expect the tolerance of C a 24 bit mantissa to be.) IF (T .GT. EPS) THEN C Lemaitre model TYPEWM = 14 MMAXZ = .FALSE. BOUNCY = .FALSE. ELSE IF (ABS(T) .LT. EPS) THEN IF (Q .GT. 0.0) THEN C A1 model, Q .GT. 0.5 TYPEWM = 15 MMAXZ = .FALSE. BOUNCY = .FALSE. ELSE C A2 model, Q .LT. -1, Eddington model TYPEWM = 17 MMAXZ = .TRUE. BOUNCY = .FALSE. ZZMAX = (LAMBDA/SIGMA)**(1.0D0/3.0D0) - 1.0D0 END IF ELSE C T < 0.0 IF (Q .GT. 0.0) THEN C q .ge. 0.5 TYPEWM = 6 MMAXZ = .FALSE. BOUNCY = .FALSE. ELSE C q .le. -1.0, bounce model TYPEWM = 18 MMAXZ = .TRUE. BOUNCY = .TRUE. C calculate ZZMAX with standard formula R = SQRT(K/(LAMBDA*3.0D0)) PHI = ACOS(-OMEGA/(LAMBDA*2.0D0*R**3))/3.0D0 Y1 = 2.0D0*R*COS(PHI) Y2 = -2.0D0*R*COS(PHI - PITHRD) Y3 = -2.0D0*R*COS(PHI + PITHRD) C the result is R/R0 = 1/v; we need z C R_0/R = 1 + z ZZMAX = 1.0D0/MAX(Y1,Y2,Y3) - 1.0D0 END IF END IF END IF END IF END IF LAMNUL = LAMBDA OMENUL = OMEGA END IF ************************************************************************ C fixes for roundoff problems C If a maximum redshift exists, we must make sure that roundoff C errors in its calculation have not made it slightly too large, C thus making Q negative. If this happens, we decrease ZZMAX in C small steps until Q becomes positive. This will almost always C keep Q from becoming negative due to roundoff error at smaller z C values. Usually, no steps are needed. When they are, one or two C are usually enough. Very rarely does one need as many as ten. C C If desired, one can use the variable NBAD, which counts how often C a bad (too small) value for QQ is returned, to monitor this, by C printing it out or making it accessible to other routines by a C COMMON block or by returning it from this routine in the C parameter list. Of course, this requires changing the code. C C Numerically, a positive but too small value for QQ is almost as C bad as a negative one. Although no error will result when the C square root is taken (via ANGSIZ), since the reciprocal of QQ is C used, overflow can be a problem. Since the `typical' value of QQ C is of the order 1, artificially decreasing ZZMAX (at which QQ is C analytically = 0) so that QQ is a positive number of the order of C the machine accuracy instead of 0 presents no practical problems. C The value of EPS set above (1.0E-6) should be OK for a 24 bit C mantissa (which most machines have). C C We save the real value for returning to the calling routine, C since ZZMAX might be decreased a bit for internal work. ZZREAL = ZZMAX IF (MMAXZ) THEN NBAD = 0 C sorry, standard FORTRAN77 has no DO WHILE, so here's a label C which one can GOTO 100 CONTINUE TEMP = QQ(ZZMAX + 1.0D0) IF (TEMP .LT. EPS) THEN ZZMAX = ZZMAX - EPS*10.0D0 NBAD = NBAD + 1 C A limit of 100 iterations is set. Since EPS is of the order C of the machine precision, this means we've decreased ZMAX by C 1000 times this. For an EPS of 1.0E-6, this means about C 1.0E-3. If this happens, it means one is dangerously near C the A2 curve, and further reduction is pointless since, even C though one could eventually produce a safe ZMAX, it would be C far away from the real ZMAX, not in actual value (perhaps C just in the second or third decimal place) but rather in the C sense that the distance to this ZMAX can be appreciably C different from the distance to the real ZMAX, so that it C seems better to give an error and let the user know he is in C a dangerous area rather than give the impression that the C results can be trusted. In practice, one will hardly ever C be looking at a world model where this is the case and, even C if this is the case, usually many fewer than 100 iterations C solve the problem without producing any meaningful error. C Further, since a further reduction is undesirable and this C part of the routine is not computationally intensive (and, C as noted above, such world models will rarely be examined C in practice) it is more trouble than it is worth to use a C more efficient method, like interval bisection. C In any case, this is only needed (according to tests) for C world models with OMEGA < 1.0E-3 or so, so the whole C business of `correcting' ZMAX (as, indeed, the bounce models C as a whole) is is of very little practical importance. IF (NBAD .EQ. 100) THEN IF (DEBUG) THEN PRINT*, 'INICOS ERROR 3: ROUNDOFF ERROR, BAD ZMAX' C Uncomment the next three lines for a more verbose C error message, useful since some output variables are C not yet accessible to the calling routine. * * PRINT*, 'WMTYPE, LAMBDA, OMEGA, MAXZ, ZMAX, Q,NBAD', * $ TYPEWM, LAMBDA, OMEGA,MMAXZ,ZZMAX,TEMP,NBAD * * PRINT* END IF C This is just to force an exit without introducing an C additional label. ZZMAX = 0.0D0 IERROR = 3 * ELSE * *C One can change the code to monitor NBAD or print out a *C message, if desired. * * IF (DEBUG) THEN * *C Since this is really a warning, it shouldn't be *C printed as an error message, and it's not worth the *C trouble to differentiate between errors and warnings *C with and without debug. If one wants the verbosity, *C one can uncomment the next line. If something *C dangerous happens, a real error message is given; see *C above. * ** WRITE (*,'(T2,A,1PE6.1E1)') 'INICOS WARNING: Q < ',EPS * *C Uncomment the next three lines for a more verbose *C warning message, useful since some output variables are *C not yet accessible to the calling routine. ** ** PRINT*, 'WMTYPE, LAMBDA, OMEGA, MAXZ, ZMAX, Q,NBAD', ** $ TYPEWM, LAMBDA, OMEGA,MMAXZ,ZZMAX,TEMP,NBAD ** ** PRINT* * * END IF * END IF ERROR = IERROR GOTO 100 END IF END IF ************************************************************************ C tidy things up C input variables C SUBROUTINE INICOS(LAMBDA,OMEGA,ETA,VARETA,DEBUG, CC ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ CC in C set the common block variables to those in the parameter list C LAMNUL and OMENUL have been set to LAMBDA and OMEGA above C if these have changed (and are valid values) ETAETA = ETA ETAVAR = VARETA DEEBUG = DEBUG C COMMON /COSMOL/ LAMNUL, OMENUL, ETAETA, K C COMMON /ANGINI/ IERROR, TYPEWM, ZZMAX, C $ NOTYET, ETAVAR, DEEBUG, MMAXZ, BOUNCY C output variables C $WMTYPE,MAXZ,ZMAX,BOUNCE,ERROR) CC ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ CC out C set the output variables to the common block values, which if C necessary (new LAMBDA or OMEGA) have been recomputed WMTYPE = TYPEWM MAXZ = MMAXZ C COMMON /FIX/ ZZREAL C use the real ZZMAX for returning to the calling routine, a dummy C value if it doesn't exist IF (MAXZ) THEN ZMAX = ZZREAL ELSE ZMAX = 0.0D0 END IF BOUNCE = BOUNCY C ERROR has been set above ************************************************************************ C set the REAL output variables ZMAXR = REAL(ZMAX) ************************************************************************ ************************************************************************ END