************************************************************************ ************************************************************************ ************************************************************************ BLOCK DATA COSANG C initialises common blocks used by INICOS, ANGSIZ and related C routines C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declaration of variables in common blocks C the usual cosmological variables DOUBLE PRECISION LAMBDA, OMEGA, ETA, K COMMON /COSMOL/ LAMBDA, OMEGA, ETA, K C internal names for other parameter list variables and internal C stuff; see INICOS INTEGER ERROR, WMTYPE DOUBLE PRECISION ZMAX LOGICAL NOTYET, VARETA, DEBUG, MAXZ, BOUNCE COMMON /ANGINI/ ZMAX, ERROR, WMTYPE, $ NOTYET, VARETA, DEBUG, MAXZ, BOUNCE C a CHOICE value < 0 indicates that the dependence of ETA on Z is C different in the phases of expansion and contraction in a bounce C model. INTEGER CHOICE COMMON /WHICH/ CHOICE 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 SAVE /COSMOL/, /ANGINI/, /WHICH/, /ARTHUR/, /FIX/ ************************************************************************ ************************************************************************ C compile-time initialisation of common block variables DOUBLE PRECISION DUMMY, IDUMMY C Could be any values, really, in most cases, since most are set to C the values passed to INICOS or to values calculated by INICOS C before they are used (all of these but K are returned to the C calling routine with other names). Some strange values, however, C might help in troubleshooting, though this original implementation C in now way depends on this. The exceptions are a) NOTYET, which C must be initialised to .TRUE.---this value should also make C troubleshooting with the other LOGICAL values easier as well--- C b) CHOICE, which we initialize to 0, mainly to avoid a negative C value, which means that the dependence of ETA on Z is different C during the expanding and contracting phases of a bounce model C (if CHOICE is negative, 4 instead of 2 distances are calculated C by BOUNCE---this should only be done if the user specifically C requests it, by setting CHOICE to a negative value) and c) RISK, C which is initialised to .FALSE. The user can declare the COMMON C block ARTHUR in the calling routine and set RISK to .TRUE., in C order to make calculations in the Eddington model. PARAMETER (DUMMY=-3.473954D-11, IDUMMY=-2222) DATA LAMBDA, OMEGA, ETA, K /4*DUMMY/ DATA ERROR, WMTYPE, ZMAX /-1, IDUMMY, DUMMY/, $ NOTYET, VARETA, DEBUG, MAXZ, BOUNCE /5*.TRUE./ DATA CHOICE /0/ DATA RISK /.FALSE./ DATA ZZREAL /DUMMY/ ************************************************************************ ************************************************************************ END ************************************************************************ ************************************************************************ ************************************************************************ 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 ************************************************************************ ************************************************************************ ************************************************************************ DOUBLE PRECISION FUNCTION QQ(V) C Returns the square of Q(V); this is done for numerical purposes, C so that appropriate action (necessary because of roundoff error) C can be taken if QQ=Q**2 is negative. C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C Declarations ************************************************************************ C input variables from the calling routine (parameter list) C (1 + z) DOUBLE PRECISION V ************************************************************************ C variables in common block C cosmological variables DOUBLE PRECISION LAMBDA, OMEGA, ETA, K COMMON /COSMOL/ LAMBDA, OMEGA, ETA, K C save until the next call SAVE /COSMOL/ ************************************************************************ ************************************************************************ C calculations C Calculates the square of fuction Q(v) needed by the ANGSIZ routine C itself and by ASDRHS; in addition, INICOS uses QQ to make sure C that ZMAX isn't too large due to roundoff error. QQ = OMEGA*V**3 - K*V**2 + LAMBDA ************************************************************************ ************************************************************************ END ************************************************************************ ************************************************************************ ************************************************************************ 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 ************************************************************************ ************************************************************************ ************************************************************************ 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 ************************************************************************ ************************************************************************ ************************************************************************ 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 ************************************************************************ ************************************************************************ ************************************************************************ SUBROUTINE ODEINT(NVAR,X1,X2,EPS,H1,HMIN,STEPR,DERIVS, C ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ C in C $ YINOUT, C ^^^^^^ C inout C $ NOK,NBAD) C ^^^^^^^^ C out C driver routine for integration of differential equations. C NVAR starting values are integrated from X1 to X2 with accuracy C EPS. H1 is a guessed first step size, HMIN the minimum allowed C step size, which can be zero. STEPR is the name of the stepper C routine to be used. DERIVS is a user-supplied routine for C calculating the right hand sides. C C The starting values YINOUT are replaced with the results. C C NOK and NBAD are the numer of good and bad steps taken. C C Based on ideas from NUMERICAL RECIPES. C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declarations ************************************************************************ C named constants as array bounds INTEGER KMAXX, NMAX PARAMETER (KMAXX=200, NMAX=50) ************************************************************************ C input variables in the parameter list C diemension of vector INTEGER NVAR C accuracy, guessed first step size, minimum step size, starting C and end point of integration, dependent variable DOUBLE PRECISION EPS, H1, HMIN, X1, X2, YINOUT(NVAR) ************************************************************************ C output variables in the parameter list C number of good and bad steps for bookkeeping INTEGER NOK, NBAD ************************************************************************ C output to called routines C all called routines C stepper routine, evaluation of right-hand sides EXTERNAL STEPR, DERIVS C variables C x an y for evaluation of right-hand sides, this result, step, C scaling vector for stepper routine DOUBLE PRECISION X, Y(NMAX), DYDX(NMAX), H, YSCAL(NMAX) ************************************************************************ C output from called routines C number of successful steps and next step size DOUBLE PRECISION HDID, HNEXT ************************************************************************ C variables in common blocks C bookkeeping (not used in this version, but put in in case someone C wants to debug or whatever) INTEGER KMAX, KOUNT DOUBLE PRECISION DXSAV, XP(KMAXX), YP(NMAX,KMAXX) COMMON /PATH/ DXSAV, XP, YP, KMAX, KOUNT C this is to keep track of (very rare) errors here and in related C routines INTEGER MSTEP COMMON /MERROR/ MSTEP C save until the next call SAVE /MERROR/ ************************************************************************ C internally used quantities C named constants C maximum number of steps INTEGER MAXSTP C just for numerical stability; see below DOUBLE PRECISION TINY PARAMETER (MAXSTP=10000, TINY=1.0D-30) C variables C DO loops INTEGER I, NSTP C for bookkeeping DOUBLE PRECISION XSAV ************************************************************************ ************************************************************************ C run-time initialisation C reset stepsize error to zero whenever routine is called MSTEP = 0 C KMAX not initialised in DATA since it's in a common block; C BLOCK DATA not worth the effort KMAX = 0 C DXSAV isn't needed if KMAX = 0 * DXSAV = 0 X = X1 H = SIGN(H1,X2-X1) NOK = 0 NBAD = 0 KOUNT = 0 ************************************************************************ ************************************************************************ C calculations C fill up work vector DO 100 I = 1, NVAR Y(I) = YINOUT(I) 100 CONTINUE C uncomment the next line if KMAX isn't 0; for KMAX = 0 neither C XSAV nor DXSAV is needed, so we don't need the next line C (This assures storage of the first step.) C IF (KMAX .GT. 0) XSAV = X - DXSAV*2.0 C at most NMAX steps DO 600 NSTP = 1, MAXSTP C get the right-hand sides CALL DERIVS(X,Y, DYDX) DO 200 I = 1, NVAR C scaling to monitor accuracy YSCAL(I) = ABS(Y(I)) + ABS(H*DYDX(I)) + TINY 200 CONTINUE C storage for internal bookkeeping IF (KMAX .GT. 0) THEN IF ( (ABS(X - XSAV)) .GT. (ABS(DXSAV)) ) THEN IF( KOUNT .LT. (KMAX - 1) ) THEN KOUNT = KOUNT + 1 XP(KOUNT) = X DO 300 I = 1, NVAR YP(I,KOUNT) = Y(I) 300 CONTINUE XSAV = X ENDIF ENDIF ENDIF C make sure the step size isn't too large IF ( ((X + H - X2)*(X + H - X1)) .GT. 0.0 ) H = X2 - X ************************************************************************ C abnormal return; set a flag and get out IF (MSTEP .EQ. 1) RETURN C This happens if there has been a failure in the stepper routine C on the previous call. ************************************************************************ C call the stepper if everything's OK CALL STEPR(NVAR,DYDX,H,EPS,YSCAL,DERIVS, X,Y, HDID,HNEXT) IF (HDID .EQ. H)THEN NOK = NOK + 1 ELSE NBAD = NBAD + 1 ENDIF C finished? IF ( ((X - X2)*(X2 - X1)) .GE. 0.0 ) THEN C replace input values with output values DO 400 I = 1, NVAR YINOUT(I) = Y(I) 400 CONTINUE IF (KMAX .NE. 0) THEN C save final step KOUNT = KOUNT + 1 XP(KOUNT) = X DO 500 I = 1, NVAR YP(I,KOUNT) = Y(I) 500 CONTINUE ENDIF ************************************************************************ C normal return; we have to get out now, since the error C conditions which follow this one are irrelevant if we're C successful, which means we've reached the next line RETURN ENDIF ************************************************************************ C abnormal return (stepsize smaller than minimum); C set a flag and get out IF( (ABS(HNEXT)) .LT. HMIN ) THEN MSTEP = 2 RETURN END IF ************************************************************************ C set the next H value H = HNEXT 600 CONTINUE ************************************************************************ C abnormal return (too many steps); set a flag and get out MSTEP = 3 RETURN ************************************************************************ ************************************************************************ END ************************************************************************ ************************************************************************ ************************************************************************ SUBROUTINE BSSTEP(NV,DYDX,HTRY,EPS,YSCAL,DERIVS, C ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ C in C $ X,Y, C ^^^ C inout C $ HDID,HNEXT) C ^^^^^^^^^^ C out C Bulirsch-Stoer stepper routine for numerical integration of C differential equations. C NV dependent variables Y and their derivatives DYDX are input at X C along with the attempted stepsize HTRY, the required accuracy EPS C and the scaling vector YSCAL. DERIVS is a routine for evaluating C the right-hand sides. C X and Y are replaced with their output values. C C HDID is the accomplished stepsize, HNEXT the next one. C C Based on ideas from NUMERICAL RECIPES. (The Second Edition C discusses more involved techniques. They are not used here C because tests showed that for this application the difference in C computing time was negligible; thus I opted for the shorter and C clearer code.) C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declarations ************************************************************************ C named constants as array bounds INTEGER NMAX, IMAX PARAMETER (NMAX=10, IMAX=11) ************************************************************************ C input variables from the calling routine (parameter list) C number of dependent variables INTEGER NV C derivatives, attempted stepsize, accuracy, scaling vector, C independent and dependent variables DOUBLE PRECISION DYDX(NV), HTRY, EPS, YSCAL(NV), X, Y(NV) ************************************************************************ C output variables to the calling routine (parameter list) C accomplished stepsize, next stepsize DOUBLE PRECISION HDID, HNEXT ************************************************************************ C output to called routines (parameter list) C all called routines C modified midpoint method, polynomial extrapolation, right-hand C sides EXTERNAL MMM, POLEX, DERIVS C named constants C number of previous estimates extrapolation should use INTEGER NUSE PARAMETER (NUSE=7) C variables C sequence of number of steps INTEGER NSEQ(IMAX) C x, y, stepsize and derivatives for modified mitpoint, x for C polynomial extrapolation DOUBLE PRECISION XSAV, YSAV(NMAX), DYSAV(NMAX), H, XEST ************************************************************************ C output from called routines (parameter list) C output vector from polynomial extrapolation and its error DOUBLE PRECISION YSEQ(NMAX), YERR(NMAX) ************************************************************************ C variables in common blocks C to let routines higher up now if something goes wrong INTEGER MSTEP COMMON /MERROR/ MSTEP C save until the next call SAVE /MERROR/ ************************************************************************ C internally used quantities C named constants C ONE mustn't necessarily be `mathematically equivalent' to 1.0; C it's a parameter like the others---one could choose another value. C Has nothing to do with making it easier to code DOUBLE PRECISION C constants. DOUBLE PRECISION ONE, SHRINK, GROW PARAMETER (ONE=1.D0, SHRINK=.95D0, GROW=1.2D0) C variables C DO loops INTEGER I, J C for error checking DOUBLE PRECISION ERRMAX ************************************************************************ ************************************************************************ C initialisation of variables C internally used quantities DATA NSEQ /2,4,6,8,12,16,24,32,48,64,96/ ************************************************************************ ************************************************************************ C run-time initialisation MSTEP = 0 H = HTRY XSAV = X ************************************************************************ ************************************************************************ C calculations C save the starting values DO 100 I = 1, NV YSAV(I) = Y(I) DYSAV(I) = DYDX(I) 100 CONTINUE C GOTO 1 allows one to repeat the main part of the routine if the C step fails. 1 CONTINUE C evaluate sequence of modified midpoint integrations DO 200 I = 1, IMAX C call the routine for modified midpoint method CALL MMM(XSAV,NV,YSAV,DYSAV,H,NSEQ(I),DERIVS, YSEQ) C series error is even, so squared XEST = (H/NSEQ(I))**2 C extrapolation with polynomials CALL POLEX(NV,XEST,YSEQ,I,NUSE, Y,YERR) C monitor local truncation error ERRMAX = 0.0D0 DO 12 J = 1, NV ERRMAX = MAX(ERRMAX,ABS(YERR(J)/YSCAL(J))) 12 CONTINUE C scale accuracy relative to tolerance ERRMAX = ERRMAX/EPS C has the step converged? IF (ERRMAX .LT. ONE) THEN C set the return variables X = X + H HDID = H IF (I. EQ. NUSE) THEN HNEXT = H*SHRINK ELSE IF (I .EQ. NUSE - 1) THEN HNEXT = H*GROW ELSE HNEXT = (H*NSEQ(NUSE - 1))/NSEQ(I) ENDIF ************************************************************************ C normal return RETURN ENDIF ************************************************************************ C get this far if the step fails 200 CONTINUE C Here it doesn't matter, with IMAX=11 and NUSE=7; C if one changes these values, one should think about whether C or not the integer division in the exponent is really wanted. C (In practice, one will probably never even get here.) C Anyway, this reduces the step size for another try if the step C fails. H = 0.25D0*H/2**((IMAX - NUSE)/2) C If H is too small, then set a flag and get out. IF (X + H. EQ. X) THEN MSTEP = 1 ************************************************************************ C abnormal return RETURN END IF ************************************************************************ C repeat if one gets this far GOTO 1 ************************************************************************ ************************************************************************ END ************************************************************************ ************************************************************************ ************************************************************************ SUBROUTINE MMM(XIN,NVAR,YIN,DYDX,HTOT,NSTEP,DERIVS, YOUT) C ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ^^^^ C in out C implements the modified midpoint method for use in Bulirsch-Stoer C integration of differential equations. C Input at XIN are NV dependent variables YIN and their derivatives C DYDX as well as the total step to be taken, HTOT, and the number C of substeps NSTEP. DERIVS is an external routine for evaluating C the right hand sides. C C The output is returned as YOUT. C C Based on ideas from NUMERICAL RECIPES. C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declarations ************************************************************************ C named constants for array bounds INTEGER NMAX PARAMETER (NMAX=50) ************************************************************************ C input variables in the parameter list C dimension of vector and number of substeps to be used INTEGER NVAR, NSTEP C independent variable, input vector and deriviative vector, total C step to be taken DOUBLE PRECISION XIN, YIN(NVAR), DYDX(NVAR), HTOT ************************************************************************ C output variables in the parameter list C output vector DOUBLE PRECISION YOUT(NVAR) ************************************************************************ C output to called routines C all called routines EXTERNAL DERIVS C variables C x and y for external evaluation routine DOUBLE PRECISION X, Y2(NMAX) ************************************************************************ C internally used quantities C variables C for DO loops INTEGER I, N C increment, 2 temporaries, second vector variable (use see below) DOUBLE PRECISION H, H2, SWAP, Y1(NMAX) ************************************************************************ ************************************************************************ C calculations C set the stepsize H = HTOT/NSTEP C initialise the work vectors for the first step DO 100 I = 1, NVAR Y1(I) = YIN(I) Y2(I) = YIN(I) + H*DYDX(I) 100 CONTINUE C increment X X = XIN + H C temporarily store derivatives in YOUT CALL DERIVS(X,Y2, YOUT) C just to save some work in the DO loop H2 = 2.0D0*H C basic step DO 300 N = 2, NSTEP DO 200 I = 1, NVAR SWAP = Y1(I) + H2*YOUT(I) Y1(I) = Y2(I) Y2(I) = SWAP 200 CONTINUE C increment X X = X + H C derivative evaluation CALL DERIVS(X,Y2, YOUT) 300 CONTINUE C this is the last step DO 400 I = 1, NVAR YOUT(I) = 0.5D0*(Y1(I) + Y2(I) + H*YOUT(I)) 400 CONTINUE ************************************************************************ ************************************************************************ END ************************************************************************ ************************************************************************ ************************************************************************ SUBROUTINE POLEX(NV,XEST,YEST,NCALL,NUSE, YOUT,DY) C ^^^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^ C in out C extrapolates functions to X=0. C NV functions are evaluated at X = 0 by polynomial extrapolation C to a progressively smaller sequence of estimates XEST and the C function vectors YEST. This call is number NCALL. Use at most C NUSE estimates. C C Extrapolated values are output as YOUT and their estimated error C as DY. C C Based on ideas from NUMERICAL RECIPES. C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declarations ************************************************************************ C named constants as array bounds C maximum expected values of NCALL and NV INTEGER IMAX, NMAX PARAMETER (IMAX=13, NMAX=50) ************************************************************************ C input variables from calling routine (parameter list) C number of functions, number of call, maximum number of previous C estimates to be used INTEGER NV, NCALL, NUSE C independent and dependent variables DOUBLE PRECISION XEST, YEST(NV) ************************************************************************ C output variables to calling routine (parameter list) C output variable and estimated uncertainty DOUBLE PRECISION YOUT(NV), DY(NV) ************************************************************************ C internally used quantities C variables C DO loops, calculated loop limit INTEGER J, K, MM C independent variable, table, work variables DOUBLE PRECISION X(IMAX), QCOL(NMAX,IMAX), D(NMAX), $ DELTA, F1, F2, Q C save until the next call SAVE QCOL, X ************************************************************************ ************************************************************************ C calculations X(NCALL) = XEST DO 100 J = 1, NV DY(J) = YEST(J) YOUT(J) = YEST(J) 100 CONTINUE C store first estimate in the first column IF(NCALL.EQ.1) THEN DO 200 J = 1, NV QCOL(J,1) = YEST(J) 200 CONTINUE ELSE C use at most NUSE previous estimates MM = MIN(NCALL,NUSE) DO 300 J = 1, NV D(J) = YEST(J) 300 CONTINUE DO 500 K = 1, (MM - 1) DELTA = 1.0D0/(X(NCALL - K) - XEST) F1 = XEST*DELTA F2 = X(NCALL - K)*DELTA C propagate tableau one diagonal more DO 400 J = 1, NV Q = QCOL(J, K) QCOL(J,K) = DY(J) DELTA = D(J) - Q DY(J) = F1*DELTA D(J) = F2*DELTA YOUT(J) = YOUT(J) + DY(J) 400 CONTINUE 500 CONTINUE DO 600 J = 1, NV QCOL(J,MM) = DY(J) 600 CONTINUE ENDIF ************************************************************************ ************************************************************************ END ************************************************************************ ************************************************************************ ************************************************************************ SUBROUTINE ASDRHS(V,D, DDDV) C ^^^ ^^^^ C in out C calculates the Right Hand Side of the differential equation C for the Angular Size Distance C Input are the dependent variables V and D. C C Output are the right-hand sides. C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declarations ************************************************************************ C input variables from the calling routine (parameter list) C V = 1 + Z and is the basic independent variable, C D(1) is the distance and D(2) its derivative. DOUBLE PRECISION V, D(2) ************************************************************************ C output variables to the calling routine (parameter list) C these are the result of this SUBROUTINE DOUBLE PRECISION DDDV(2) ************************************************************************ C output variables from called routines (parameter list) C G = SQRT(QQ(V)) C IF (VARETA) ETA = ETAZ(V-1.0) C customised ETA, Q(v)**2 DOUBLE PRECISION ETAZ, QQ EXTERNAL ETAZ, QQ ************************************************************************ C variables in common blocks C cosmological variables (see INICOS) DOUBLE PRECISION LAMBDA, OMEGA, ETA, K COMMON /COSMOL/ LAMBDA, OMEGA, ETA, K 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 save until the next call SAVE /COSMOL/, /ANGINI/ ************************************************************************ C internally used quantities C variables C just to keep things neat DOUBLE PRECISION Q, G, GS ************************************************************************ ************************************************************************ C calculations Q = QQ(V) G = SQRT(Q) C and G's derivative GS = (3.0D0*OMEGA*V**2 - 2.0D0*K*V)/(2.0D0*G) C calculate the derivatives DDDV needed by the calling routine C D(1) and D(2) are D and D', so D' = D' (definition) DDDV(1) = D(2) C and D'' is calculated C if required replace ETA with the value given by ETAZ IF (VARETA) ETA = ETAZ(V-1.0D0) DDDV(2) = -((3.0D0*OMEGA*ETA*V*D(1))/(2.0D0*Q)) - $ (2.0D0/V + GS/G)*D(2) ************************************************************************ ************************************************************************ END ************************************************************************ ************************************************************************ ************************************************************************ DOUBLE PRECISION FUNCTION ETAZ(Z) C Calculates ETA as a function of Z. If desired, this value of ETA C can be used by ANGSIZ instead of a constant value. C This routine can be modified by the user to suit his or her needs. C Neither INICOS, ANGSIZ nor BNGSIZ need be modified. C C Declaring the same COMMON block WHICH in the routine which calls C ANGSIZ makes it possible to investigate different behavior C of ETA(Z) without modifying the code. C C This file in its original version is intended as an example. C Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. ************************************************************************ ************************************************************************ C declarations ************************************************************************ C input variable from the calling routine DOUBLE PRECISION Z ************************************************************************ 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 This can also be declared in the routine which calls ANGSIZ, thus C making a decision about which z-dependence for ETA to use in this C calling routine. USE NEGATIVE CHOICE VALUES IF THE VALUE OF ETA C FOR THE SAME REDSHIFT CAN BE DIFFERENT IN THE CONTRACTING AND C EXPANDING PHASES OF A BOUNCE MODEL. (See VALUE2 below.) INTEGER CHOICE COMMON /WHICH/ CHOICE LOGICAL VALUE2 C This can be used to find out if the integration is being done C in the contracting phase of a bounce model, in which case this C could lead to a different value for ETAZ than in the expanding C phase. ZMAX from the previous ANGINI COMMON block might be C useful in this case as well. C VALUE2 is .TRUE. during the contracting phase of a bounce model. COMMON /CNTRCT/ VALUE2 C save until the next call SAVE /ANGINI/, /WHICH/, /CNTRCT/ ************************************************************************ C internally used variable C just to keep things neat DOUBLE PRECISION Y ************************************************************************ ************************************************************************ C calculations Y = 1.0D0/(1.0D0 + Z) IF (CHOICE .EQ. 0) THEN ETAZ = 0.0D0 C Different ETA values before and after the bounce; C the negative value of CHOICE is necessary for ANGSIZ, since C some degeneracy is cancelled if ETAZ depends on VALUE2. ELSE IF (CHOICE .EQ. -1) THEN IF (VALUE2) THEN ETAZ = 1.0D0 ELSE ETAZ = 0.0D0 END IF C different ETA values before and after the bounce, continuous at C ZMAX ELSE IF (CHOICE .EQ. -2) THEN IF (VALUE2) THEN ETAZ = 1.0D0/(1.0D0 + ZMAX) ELSE ETAZ = Y END IF ELSE IF (CHOICE .EQ. 3) THEN ETAZ = 1.0D0 - Y ELSE IF (CHOICE .EQ. 4) THEN ETAZ = Y**2 ELSE IF (CHOICE .EQ. 5) THEN ETAZ = 1.0D0 - Y**2 ELSE IF (CHOICE .EQ. 6) THEN IF (Z .LE. 9) THEN ETAZ = LOG10(1.0D0 + Z) ELSE ETAZ = 1.0D0 END IF ELSE IF (CHOICE .EQ. 9) THEN ETAZ = SQRT(Y)*(1.0D0 - Y)**2 ELSE ETAZ = SQRT(1.0D0 - Y)*Y**2 END IF C take care of roundoff error IF (ETAZ .GT. 1.0D0) THEN ETAZ = ETAZ - 0.0001D0 ELSE IF (ETAZ .LT. 0.0D0) THEN ETAZ = ETAZ + 0.0001D0 END IF IF ( (ETAZ .LT. 0.0D0) .OR. (ETAZ .GT. 1.0D0) ) THEN C if programmed correctly, one can't get here, not even with C roundoff error PRINT*, 'SOMETHING''S WRONG!!!!!!' C just to make sure one is careful in the definition of ETAZ STOP END IF ************************************************************************ ************************************************************************ END