! when all clear need to compare with F77 for all ETA and VARETA ! order these points to do at the top ! does D change with ETAZ for ANGSIZ? ! problem with D32, related to turnaround? does it exist in old code? ! ONE usw. ! when is D negative? ! we don't set D to anything if e.g. Z2 > ZMAX ! set to DUMMY? also in other cases? ! small problem: ANGSIZ ERROR, then call BNGSIZ !THIS IS NOW SOVED ! need to really check ERROR and IERROR in ANGSIZ and BNGSIZ (and ! following routines) ! large problem: third distance is always 0 ! ! output unchanged for CORRECT calling sequence (INICOS error in ANGSIZ) ! IERROR and ERROR ! not specifically this code: 5 severities ! make sure the two lists correspond ! change EPS to 10*EPSILON and TINY to 3*TINY and compare results ! with extensive test program (including whole range of ETA stuff) ! check BNGSIZ by having CHOICE < 0 but same behaviour ! ! remove passing of module variables? ! search for block as in common block ! ONE, TWO, THREE etc ! change ' to "????? ! need to have the possibility for the user to specify his own ETAZ, ! perhaps with another name ! ! -2 big common block -1 many common blocks 0 keep common blocks ! 1 many modules for variables 2 just use module variables 3 all variables ! in separate module 4 like 3 but use only what is needed ! through 1 is clear not to do it 3 is like 2 but more complicated ! question is thus 2 or 4 will implement 2 for now ! ! muss ich Y = -Y for D32?!?!?!?!?!?!?! ! calling BNGSIZ before ANGSIZ is not caught! ! !! bug in TEXT17!!!!! 1_0-2 and 1_1-2 after all?! ! ! more named constants? perhaps wait for math library? !! comments !! remove 6 spaces !! $-->& !! operator symbols !! DO without CONTINUE !!indent comments (replace "! " with " !" etc) !!! case construct !! ERROR messages without block IF ! intent ! named constants (ONE) ! new concept for variables in common blocks ! remove ALL common blocks ! remove unneeded comments ! :: in declarations ! longer lines ! ! add F95 features ! ! ! public release ! ! add stuff to INICOS ! ! 1 module to replace all old common blocks, USE ONLY? ! ZZREAL ! !*********************************************************************** !*********************************************************************** !*********************************************************************** MODULE HELBIG_ANGSIZ ! contains the routines for calculating the angular size distance !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. ! revision history ! ! ! The next paragraph written 06 July 2007 summarises previous history: ! ! The original code was written by Phillip Helbig in September 1995, ! though this was essentially a cleanup of functionally identical code ! written for my own use during 1992 and 1993. The cleanup was done ! prior to making the code public. ! The original Fortran77 version was trivially converted to a Fortran90 ! version in June 1997. ! Between 18-JUN-1997 18:21:11.00 and writing this, nothing has been done ! except a bug fix to move 1_0-0 to 1_0-1 AND 1_1-0 to 1_1-1. ! ! ! July 2007: moving from 1_1-1 to 1_2-0 via the following changes, none ! of which should affect the output of the code (and don't on the ! machines I have tested it on): ! ! general cleanup: ! error messages with IF statement instead of block IF if possible ! more named constants (e.g. ONE=1.0D0) ! fixed- --> free-form source: ! o all comments now use "!" ! o unindented lines start in column 1 ! o "&" instead of "$" used for continuation ! o lines limited to 80 instead of 72 characters ! o declarations contain "::", whether necessary or not ! other Fortran90 features: ! o comparison symbols now used, e.g. ">=" instead of ".GE." ! o DO loops now have no CONTINUE statements but use END DO ! o the CASE construct is used where useful (which is nowhere) ! o arguments to routines now have the INTENT attribute ! o COMMON blocks replaced with MODULE variables (Note that this ! changes, at least in some cases, the interface between the ! user-supplied calling routine and the routines in this module, ! so this is a "dot release" instead of a "dash release". However, ! it is not a major release since no new functionality is provided.) ! o superfluous comments removed (many were related to the COMMON ! blocks or described the structure of the Fortran77 code since ! this was not always obvious from the code itself) ! Fortran95 features ! IMPLICIT NONE PRIVATE PUBLIC :: INICOS, ANGSIZ, BNGSIZ !procedures PUBLIC :: CHOICE, RISK, loops !variables !*********************************************************************** !*********************************************************************** !declaration of common module variables ! flags which can be set in the calling routine INTEGER :: CHOICE !value <0 means ETA(Z) depends on expansion or contraction LOGICAL :: RISK !to allow computation in the Eddington model, if desired !the usual cosmological variables set in INICOS DOUBLE PRECISION :: LAMBDA, OMEGA, ETA !helpful cosmological quantity DOUBLE PRECISION :: K !other input to INICOS LOGICAL :: VARETA, DEBUG ! things calculated by INICOS (and/or other routines) INTEGER :: ERROR, WMTYPE, loops=0 DOUBLE PRECISION :: ZMAX LOGICAL :: MAXZ, BOUNCE ! internal flag LOGICAL :: NOTYET !sometimes ANGSIZ could be called with a Z which is smaller than !ZMAX but large enough for numerical problems. DOUBLE PRECISION ZZREAL !the correct value for ZMAX !communication between ANGSIZ and BNGSIZ DOUBLE PRECISION :: Z1, Z2, D12 DOUBLE PRECISION, DIMENSION(2) :: YSTART !also used in call list to LOWLEV !communication between ANGSIZ and ETAZ LOGICAL VALUE2 !communication with numerical integration stuff INTEGER MSTEP !*********************************************************************** !*********************************************************************** !compile-time initialisation of common variables !Could be any values, really, in most cases, since most are set to !the values passed to INICOS or to values calculated by INICOS !before they are used (all of these but K are returned to the !calling routine with other names). Some strange values, however, !might help in troubleshooting, though the implementation of this module !in no way depends on this. The exceptions are a) NOTYET, which !must be initialised to .TRUE.---this value should also make !troubleshooting with the other LOGICAL values easier as well---, !b) CHOICE, which we initialize to 0, mainly to avoid a negative !value, which means that the dependence of ETA on Z is different !during the expanding and contracting phases of a bounce model !(if CHOICE is negative, 4 instead of 2 distances are calculated !by BOUNCE---this should only be done if the user specifically !requests it, by setting CHOICE to a negative value) and c) RISK, !which is initialised to .FALSE. The user can set RISK to .TRUE. in !order to make calculations in the Eddington model possible. ! !In Fortran77, the following statements were inside a BLOCKDATA routine, !since that was the only way to initialise COMMON-block variables. As !of Fortran90, I could do the initialisation with the declaration, and !have done so where it makes sense, but in these cases, I think it is !clearer that the discussion above precedes the initialisation below, !and also it makes since to use DATA when several variables are !initialised to the same value, even though one has to repeat the !variable name. INTEGER, PARAMETER :: IDUMMY=-222222 REAL, PARAMETER :: DUMMY=-2.2222222D-22 DATA CHOICE /0/, RISK /.FALSE./ DATA ERROR, WMTYPE /2*IDUMMY/, & LAMBDA, OMEGA, ETA, K, ZMAX, ZZREAL /6*DUMMY/, & VARETA, DEBUG, MAXZ, BOUNCE, NOTYET /5*.TRUE./ !*********************************************************************** !*********************************************************************** CONTAINS !*********************************************************************** !*********************************************************************** !*********************************************************************** SUBROUTINE INICOS(LAMBDR,OMEGAR,ETAR,ETAVAR,DEEBUG, & ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! in ! TYPEWM,MMAXZ,ZMAXR,BOUNCY,IERROR) ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! out !checks the cosmological parameters, relays user wishes to ANGSIZ, !calculates the maximum possible redshift, indicates errors, is to !be called as initialisation before calling ANGSIZ, and gives some !indication as to the type of cosmological model. !LAMBDA and OMEGA are the usual $\lambda_{0}$ and $\Omega_{0}$, ETA !gives the fraction of homogeneously distributed matter. VARETA, !if .TRUE., indicates that ETA is not to be determined by the !(z-independent) value of ETA, but rather through the !REAL FUNCTION ETAZ(Z), which the user can modify as needed. !DEBUG, if .TRUE., enables the reporting of error messages. ! !WMTYPE gives an indication as to the cosmological model, which can !be used in conjuction with the table in the User's Guide. !BOUNCE tells whether or not the world model is a bounce model; in !these cases, ZMAX is finite and BOUNCE and MAXZ are .TRUE. For !the A2 models, BOUNCE is false, but ZMAX is also finite, and thus !MAXZ is .TRUE. (If ZMAX is analytically infinite, it is !computationally undefined, and MAXZ is .FALSE..) If an error !occurs, then ERROR is suitably defined, whether or not DEBUG is !.TRUE. !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !A note on variable names. First, INICOS, ANGSIZ and BNGSIZ take !single-precision input parameters, but all internal calculations are !done in double precision: double precision internally for more !accuracy, but a single-precision interface so as not to suggest more !accuracy than actually exists. Thus, floating-point variables in the !parameter list have "R" (for REAL) at the end. With regard to !non-floating point variables such as VARETA and DEBUG, I originally !used these names in the parameter list and for internal calculation !in INICOS, then assigning them to variants such as ETAVAR and DEEBUG !in a COMMON block for communication with other routines. (In these !other routines, the variables could have the "proper" names since !association was via the name of the COMMON block, not via the names of !the variables.) This also applies to output variables such as WMTYPE !and BOUNCE. Thus, the "strange" names appeared only in a COMMON block !and only in INICOS. However, when using variables common to all !routines in the module, they should have the same name everywhere. !(It would be possible to put these variables in an additional module !and then rename them in the USE statement, thus keeping the naming !scheme the same in the Fortran77 and Fortran90 versions of the code, !but this is more trouble than it is worth, especially since getting rid !of the COMMON blocks makes the code more readable in any case.) Thus, !in the other routines I changed nothing, and here in INICOS now use the !"strange" names in the parameter list and for internal calculations, !and the common variables have the "normal" names. Although this means !that the "strange" names are a bit more visible, being used for !internal calculations, it also means that the common variables have the !same name when used by INICOS as when used by the other routines. (The !discuussion of the definition of the variables refers to the names of !the corresponding common module variables.) !*********************************************************************** !*********************************************************************** !declarations !*********************************************************************** !input variables from the calling routine (parameter list) !the cosmological parameters---input variables REAL, INTENT(IN) :: LAMBDR, OMEGAR, ETAR !the cosmological parameters---internal DOUBLE PRECISION counterparts DOUBLE PRECISION :: LAMNUL, OMENUL, ETAETA !customised ETA? display messages? LOGICAL, INTENT(IN) :: ETAVAR, DEEBUG !*********************************************************************** !output variables to the calling routine (parameter list) !type of cosmological model, error message INTEGER, INTENT(OUT) :: TYPEWM, IERROR !maximum redshift---real output variable REAL, INTENT(OUT) :: ZMAXR !maximum redshift---internal DOUBLE PRECISION counterpart DOUBLE PRECISION :: ZZMAX !is there a maximum redshift? is this a bounce model? LOGICAL, INTENT(OUT) :: MMAXZ, BOUNCY !*********************************************************************** !internally used quantities !named constants !$\frac{\pi}{3}$, accuracy, comparisons !The value of EPS should be about right for a 24 bit mantissa; !see the note concerning precision in ANGSIZ. DOUBLE PRECISION, PARAMETER :: PITHRD=1.0471975511965977D0, EPS = 1.0D-6, & ZERO = 0.0D0, ONE=1.0D0 !variables !accuracy monitor INTEGER :: NBAD !internal work variables (SIGMA and Q have the usual cosmological !meaning) DOUBLE PRECISION :: SIGMA, Q, T, R, PHI, Y1, Y2, Y3, TEMP !*********************************************************************** !*********************************************************************** !run-time initialisation !set the input variables LAMNUL = DBLE(LAMBDR) OMENUL = DBLE(OMEGAR) ETAETA = DBLE(ETAR) !*********************************************************************** !*********************************************************************** !calculations !*********************************************************************** !check the input parameters IF (OMENUL < ZERO) THEN ERROR = 1 IF (DEEBUG) PRINT*, 'INICOS ERROR 1: OMEGA < 0.0; OMEGA = ', OMENUL ELSE IF ((.NOT. ETAVAR) .AND. ((ETAETA < ZERO) .OR. (ETAETA > ONE))) THEN ERROR = 2 IF (DEEBUG) PRINT*, 'INICOS ERROR 2: ETA OUT OF RANGE; ETA = ', ETAETA ELSE ERROR = 0 END IF IERROR = ERROR !*********************************************************************** !abnormal exit IF (ERROR > 0) RETURN !*********************************************************************** !calculate WMTYPE, BOUNCE, MAXZ and ZMAX, if necessary !Has OMEGA or LAMBDA changed since the last call? .OR. is this the first call? ! !Note that no message is given if INICOS is changed without changing !the input parameters. While normally a mistake, there might be some !legitimate cases of this, such as during interactive use, with a !quick-and-dirty calling routine or when LAMBDA and OMEGA are not !directly specified in the calling routine. It thus seems better to !give no message at all. A better approach might be to have an INTEGER !variable DEBUG whose increasing value would increase the number of !messages output, e.g. none, success, informational, warning, error, !fatal error, where each category includes the previous. This doesn't !seem worth the trouble, especially for a rewrite. (Calling INICOS with !unchanged parameters would obviouly be a warning.) IF ( (LAMNUL /= LAMBDA) .OR. (OMENUL /= OMEGA) .OR. NOTYET) THEN !this only has an effect the first time NOTYET = .FALSE. !almost always executed but evaluation is more expensive !the sign of k is the sign of the curvature K = OMENUL + LAMNUL - ONE !determine WMTYPE, MAXZ and ZMAX IF (LAMNUL < ZERO) THEN MAXZ = .FALSE. BOUNCE = .FALSE. IF (K < ZERO) THEN IF (OMENUL > ZERO) THEN WMTYPE = 2 ELSE !empty universe WMTYPE = 1 END IF ELSE IF (K > ZERO) THEN WMTYPE = 4 ELSE !k = 0 WMTYPE = 3 END IF ELSE IF (LAMNUL == ZERO) THEN MAXZ = .FALSE. BOUNCE = .FALSE. IF (K < ZERO) THEN IF (OMENUL > ZERO) THEN WMTYPE = 8 ELSE !empty universe, `Milne-model' WMTYPE = 7 END IF ELSE IF (K > ZERO) THEN !`Misner-Thorne-Wheeler model' WMTYPE = 5 ELSE !k = 0, Einstein-de Sitter model WMTYPE = 9 END IF ELSE !LAMNUL > 0.0 IF (K < ZERO) THEN MAXZ = .FALSE. BOUNCE = .FALSE. IF (OMENUL > ZERO) THEN WMTYPE = 11 ELSE !empty universe WMTYPE = 10 END IF ELSE IF (K == ZERO) THEN MAXZ = .FALSE. BOUNCE = .FALSE. IF (OMENUL > ZERO) THEN WMTYPE = 13 ELSE !empty universe, de Sitter model WMTYPE = 12 END IF ELSE !k > 0 IF (OMENUL == ZERO) THEN !empty universe, Lanczos model !Near the de Sitter model (WMTYPE = 12) MAXZ is !very large, so practically no value of z will be !rejected. However, this can lead to numerical !problems, since below QQ(MAXZ) is taken to make sure !it is large enough. (Negative values are bad because !a square root is taken via ANGSIZ; small positive !values are bad because the reciprocal is used by !ANGSIZ and this can lead to overflow.) We therefore !declare the world model to be practically de Sitter !if it is practically (numerically) flat. In day to !day practice this should lead to no error. IF (ABS(K) < EPS) THEN WMTYPE = 12 MAXZ = .FALSE. BOUNCE = .FALSE. ELSE !normal code for this case WMTYPE = 19 MAXZ = .TRUE. BOUNCE = .TRUE. ZMAX = SQRT(LAMNUL/(LAMNUL - ONE)) - ONE END IF ELSE !sigma and q are mathematically and computationally useful SIGMA = 0.5D0*OMENUL Q = SIGMA - LAMNUL !Stabell-Refsdal equation (10) T = -K + 3.0D0*SIGMA**(2.0D0/3.0D0)*LAMNUL**(1.0D0/3.0D0) !T = 0 defines the A1 and A2 curves in Stabell and !Refsdal. On the VAX, if we set EPS = 0, there is no !problem; when T is analytically = 0 (for example, !LAMBDA = 2.0 and OMEGA = 0.5) then it is also !numerically = 0, and there is no conflict (at least !when LAMBDA and OMEGA are represented exactly, as in !the example above). Other machines, for instance IBM !RS6000, calculate a T of about 1.0E-15 in this case, !due probably to filling up the multiplicands with 32 !bits of zeros in the registers. Now, normally, this !would make no difference, since for the purposes of !determining the world model, we are interested in !whether or not a particular criterion is NUMERICALLY !fulfilled. However, when T is very small, one is !very close to the A curve, and this can give problems !in the calculation of the angular size distance, since !it can become arbitrarily large for z near ZMAX. !Therefore, it is a good idea to draw the A curve with !a blunt numerical pencil. However, the fudge factor !is very small, and should be of no consequence for !the practical cosmologist. (One might think that the !fudge factor should be increased to exclude models !which are ANALYTICALLY and numerically near the A !curve, since the angular size distance can become very !large in these cases as well. However, tests show !that machine precision tends to be coarse enough so !that if one is analytically off the curve, one is far !enough away so as to encounter no problems later on. !In fact, initially EPS/100000.0 was OK on the VAX/VMS !(naturally), ALPHA/VMS, SUN, Intel/Linux and IBM !RS/6000/AIX for the pwr and pwr2 architectures. !However, optimising for the ppc architecture required !a reduction to EPS/10.0. To be on the safe side, we !use EPS here---as in other places in this package--- !which is about what one would expect the tolerance of !a 24 bit mantissa to be.) IF (T > EPS) THEN !Lemaitre model WMTYPE = 14 MAXZ = .FALSE. BOUNCE = .FALSE. ELSE IF (ABS(T) < EPS) THEN IF (Q > ZERO) THEN !A1 model, Q > 0.5 WMTYPE = 15 MAXZ = .FALSE. BOUNCE = .FALSE. ELSE !A2 model, Q < -1, Eddington model WMTYPE = 17 MAXZ = .TRUE. BOUNCE = .FALSE. ZMAX = (LAMNUL/SIGMA)**(1.0D0/3.0D0) - ONE END IF ELSE !T < 0.0 IF (Q > ZERO) THEN !q >= 0.5 WMTYPE = 6 MAXZ = .FALSE. BOUNCE = .FALSE. ELSE !q <= -1.0, bounce model WMTYPE = 18 MAXZ = .TRUE. BOUNCE = .TRUE. !calculate ZMAX with standard formula R = SQRT(K/(LAMNUL*3.0D0)) !not the cosmological quantity PHI = ACOS(-OMENUL/(LAMNUL*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) !the result is R/R0 = 1/(1 + z); here R IS the cosmological R !we need z and R_0/R = 1 + z ZMAX = 1.0D0/MAX(Y1,Y2,Y3) - 1.0D0 END IF !Q END IF !T END IF !OMEGA END IF !K END IF !LAMBDA LAMBDA = LAMNUL OMEGA = OMENUL END IF !calculations !*********************************************************************** ! fixes for roundoff problems ! If a maximum redshift exists, we must make sure that roundoff ! errors in its calculation have not made it slightly too large, ! thus making Q negative. If this happens, we decrease ZMAX in ! small steps until Q becomes positive. This will almost always ! keep Q from becoming negative due to roundoff error at smaller z ! values. Usually, no steps are needed. When they are, one or two ! are usually enough. Very rarely does one need as many as ten. ! WMTYPE 17 with LAMBDA = 2.0 and OMEGA = 0.5 has ZMAX = 1. In this ! case, 82 iterations are needed. The redshift of 1.0 exists in the ! model, but the distance to it can be large. Thus, we need to make ! sure that the code behaves properly if one inputs 1.0 as Z2. ! ! If desired, one can use the variable NBAD, which counts how often ! a bad (too small) value for QQ is returned, to monitor this, by ! printing it out or making it accessible to other routines via including ! it in the PUBLIC statement or by returning it from this routine in the ! parameter list. Of course, this requires changing the code. ! ! Numerically, a positive but too small value for QQ is almost as ! bad as a negative one. Although no error will result when the ! square root is taken (via ANGSIZ), since the reciprocal of QQ is ! used, overflow can be a problem. Since the `typical' value of QQ ! is of the order 1, artificially decreasing ZMAX (at which QQ is ! analytically = 0) so that QQ is a positive number of the order of ! the machine accuracy instead of 0 presents no practical problems. ! The value of EPS set above (1.0E-6) should be OK for a 24 bit ! mantissa (which most machines have). ! ! We save the real value for returning to the calling routine, ! since ZMAX might be decreased a bit for internal work. ZZREAL = ZMAX IF_MAXZ: IF (MAXZ) THEN DO, NBAD = 0, 99, 1 TEMP = QQ(ZMAX + ONE) IF_EPS: IF (TEMP >= EPS) THEN EXIT ELSE !*!One can change the code to monitor NBAD or print out a !*!message, if desired. !* !*IF (DEEBUG) THEN !* !* !Since this is really a warning, it shouldn't be !* !printed as an error message, and it's not worth the !* !trouble to differentiate between errors and warnings !* !with and without debug. If one wants the verbosity, !* !one can uncomment the next line. If something !* !dangerous happens, a real error message is given; see !* !below. !* !* WRITE (*,'(T2,A,1PE6.1E1)') 'INICOS WARNING: Q < ',EPS !* !* !Uncomment the next three lines for a more verbose !* !warning message, useful since some output variables are !* !not (yet) accessible to the calling routine. !** !** PRINT*, 'WMTYPE, LAMBDA, OMEGA, MAXZ, ZMAX, Q, NBAD', & !** WMTYPE, LAMBDA, OMEGA, MAXZ, ZMAX, TEMP, NBAD !** PRINT* !* !*END IF ZMAX = ZMAX - EPS*10.0D0 !A limit of 100 iterations is set. Since EPS is of the order !of the machine precision, this means we've decreased ZMAX by !1000 times this. For an EPS of 1.0E-6, this means about !1.0E-3. If this happens, it means one is dangerously near !the A2 curve, and further reduction is pointless since, even !though one could eventually produce a safe ZMAX, it would be !far away from the real ZMAX, not in actual value (perhaps !just in the second or third decimal place) but rather in the !sense that the distance to this ZMAX can be appreciably !different from the distance to the real ZMAX, so that it !seems better to give an error and let the user know he is in !a dangerous area rather than give the impression that the !results can be trusted. In practice, one will hardly ever !be looking at a world model where this is the case and, even !if this is the case, usually many fewer than 100 iterations !solve the problem without producing any meaningful error. !Further, since a further reduction is undesirable and this !part of the routine is not computationally intensive (and, !as noted above, such world models will rarely be examined !in practice) it is more trouble than it is worth to use a !more efficient method, like interval bisection. ! !In any case, this is only needed (according to tests) for !world models with OMEGA < 1.0E-3 or so, so the whole !business of `correcting' ZMAX (as, indeed, the bounce models !as a whole) is is of very little practical importance. END IF IF_EPS END DO !we only get here if all 100 interations are done IF_NBAD: IF (NBAD == 100) THEN IF (DEEBUG) THEN PRINT*, 'INICOS ERROR 3: ROUNDOFF ERROR, BAD ZMAX' ! Uncomment the three lines beginning with "!*" for a more ! verbose error message, useful since some output variables ! are not accessible to the calling routine. ! !*PRINT*, 'WMTYPE, LAMBDA, OMEGA, MAXZ, ZMAX, Q, NBAD', & !* WMTYPE, LAMBDA, OMEGA, MAXZ, ZMAX, TEMP, NBAD !*PRINT* END IF ZMAX = 0.0D0 ERROR = 3 IERROR = ERROR END IF IF_NBAD END IF IF_MAXZ !*********************************************************************** !tidy things up !input variables !SUBROUTINE INICOS(LAMBDR,OMEGAR,ETAR,ETAVAR,DEEBUG, & !! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ !! in !set the common module variables to those in the parameter list !LAMBDA and OMEGA have been set to LAMNUL and OMENUL (which in turn had !been set to LAMBDAR and OMEGAR) above if LAMBDA or OMEGA has changed !(and both are valid values) ETA = ETAETA VARETA = ETAVAR DEBUG = DEEBUG !output variables ! TYPEWM,MMAXZ,ZMAXR,BOUNCY,IERROR) !! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ !! out !set the output variables to the values of the common module variable, !which if necessary (new LAMBDA or OMEGA) have been recomputed TYPEWM = WMTYPE MMAXZ = MAXZ !use the real ZZMAX for returning to the calling routine, a dummy !value if it doesn't exist IF (MAXZ) THEN ZZMAX = ZZREAL ELSE ZZMAX = 0.0D0 END IF BOUNCY = BOUNCE !ERROR has been set above !*********************************************************************** !set the REAL output variables ZMAXR = REAL(ZZMAX) !*********************************************************************** !*********************************************************************** END SUBROUTINE INICOS !*********************************************************************** !*********************************************************************** !*********************************************************************** DOUBLE PRECISION FUNCTION QQ(V) !Returns the square of Q(V); this is done for numerical purposes, !so that appropriate action (necessary because of roundoff error) !can be taken if QQ=Q**2 is negative. !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !Declarations !*********************************************************************** !input variables from the calling routine (parameter list) !(1 + z) DOUBLE PRECISION, INTENT(IN) :: V !*********************************************************************** !*********************************************************************** !calculations !Calculates the square of fuction Q(v) needed by the ANGSIZ routine !itself and by ASDRHS; in addition, INICOS uses QQ to make sure !that ZMAX isn't too large due to roundoff error. QQ = OMEGA*V**3 - K*V**2 + LAMBDA !*********************************************************************** !*********************************************************************** END FUNCTION QQ !*********************************************************************** !*********************************************************************** !*********************************************************************** SUBROUTINE ANGSIZ(Z1R,Z2R, D12R,IERROR) ! ^^^^^^^ ^^^^^^^^^^^ ! in out !calculates the angular size distance between two redshifts (see !Kayser et al., A&A, 1996). !Given two redshifts Z1 and Z2 this subroutine calculates the !angular size distance between these two redshifts, which is !returned as D12. In some world models (bounce models) there are !two (or four) independent distances for the same redshift; ANGSIZ !returns as D12 the distance between two redshifts Z1 and Z2 such !that the universe is NOT in a state of CONTRACTION---neither at !the time of emission nor at the time of reception of the observed !radiation. (The additional distances in bounce models are !calculated by the SUBROUTINE BNGSIZ.) ! !The integer variable ERROR is 0 when no error is encountered, and !is otherwise set to a specific value for the corresponding error. !This can be useful in debugging. ! !If the cosmological model is in error, the corresponding value for !ERROR is returned. Otherwise, the redshifts are examined, to see !if the distance is defined. ! !Undefined distances (producing an error) are those where at least !one redshift is < 0 or > ZMAX, where ZMAX = \infty in big bang !models and a finite value in other models. ! !An error message is also returned if ANGSIZ is called before !calling INICOS. Of course, if the input parameters to INICOS are !changed after a call, it is the responsibility of the user to call !INICOS again before calling ANGSIZ. ! !If indicated in the INICOS call, error messages, when required, !are printed to standard output. !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !Declarations !*********************************************************************** !input variables from the calling routine (parameter list) !redshifts !the real input variables REAL, INTENT(IN) :: Z1R, Z2R !*********************************************************************** !output variables to the calling routine (parameter list) !error flag INTEGER, INTENT(OUT) :: IERROR !angular size distance (D14, D34 and D32 are the other distances !in bounce models) !the actual output variable REAL, INTENT(OUT) :: D12R !*********************************************************************** !output to called routines (parameter list) !does the actual integration: CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) !constants !This seems to work fine for ANGSIZ; the exact value isn't too !important, as in BNGSIZ. This is about 25% faster than using !the BNGSIZ value of 200. DOUBLE PRECISION, PARAMETER :: NSTEP=2.0D0 !variables !(1 + Z1) and (1 + Z2) DOUBLE PRECISION :: VD, VS !*********************************************************************** !output from called routines (parameter list) !does the actual integration: CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) !angular size distance DOUBLE PRECISION :: D !*********************************************************************** !internally used quantities !constants !just to check for overflow DOUBLE PRECISION, PARAMETER :: TINY=1.0D-38 !This is roughly the smallest REAL number for a 4 byte REAL. !While internally higher precision is used, the input and output !variables are REAL and thus this is a better choice. 'higher !precision' instead of the corresponding FORTRAN key words is used !here to enable the conversion to single precision via a global !search and replace, if so desired. Memory requirements of these !routines are minimal; whether single or higher precision is faster !depends on the machine. All internal quantities are higher !precision to avoid conversions. Higher precision is not required !if one does not use BNGSIZ, but is used throughout because memory !requirements are minimal and it is even faster on some machines. !If single precision is desired, replacing the FORTRAN key words !for higher precision with REAL in all routines will do the trick !and can be done automatically. (Additionally, `D0' should be !removed. CHARACTER(LEN=56), PARAMETER :: TEXT17 = & 'NOT REALLY LARGER THAN ZMAX BUT CLOSE ENOUGH FOR TROUBLE' !variables !just a temporary work variable DOUBLE PRECISION :: G !*********************************************************************** !*********************************************************************** !run-time initialisation Z1 = DBLE(Z1R) Z2 = DBLE(Z2R) !*********************************************************************** !*********************************************************************** !calculations !is a calculation of the angular size distance possible? !ERROR will be changed if something goes wrong IERROR = 0 !has INICOS been called first? IF (NOTYET) THEN PRINT*, 'YOU MUST CALL INICOS FIRST!!!!!!!!!!!!!!!!!!!!!!!!!!!' !add ERROR = -1 to documentation and harmonise message style! IERROR = -1 ERROR = IERROR !*********************************************************************** !abnormal return: no further calculation in this case !*********************************************************************** RETURN END IF !based on the cosmological model? !Changing the second condition to (ERROR <= 2) will suppress !IERROR = 3 and the message here (though not in INICOS, of course). !Note that in contrast to 1 and 2, an ERROR of 3 does not necessarily !make further calculation meaningless, so one might want to change the !code accordingly. IF ((ERROR >= 1) .AND. (ERROR <= 3)) THEN IF (DEBUG) WRITE (*, '(T2,A,I1,A)') 'ERROR ', ERROR, ': INICOS ERROR' IERROR = ERROR !Don't allow computation in the Eddington model, WMTYPE 17, unless !the user explicitly requests this. !C might want to change this if this is possible with other methods ELSE IF ((WMTYPE == 17) .AND. .NOT. RISK) THEN IERROR = 4 IF (DEBUG) PRINT*, 'ERROR 4: COMPUTATION NOT RECOMMENDED IN WMTYPE 17' !based on the redshifts? ELSE ERROR = 0 !only reset ERROR if we get this far IF (MAXZ) THEN IF (Z2 > ZMAX) THEN IF (Z2 > ZZREAL) THEN IERROR = 9 IF (DEBUG) PRINT*, 'ERROR 9: Z2 > ZMAX' ELSE IERROR = 10 IF (DEBUG) WRITE (*,'(T2,2A)') 'Z2 ', TEXT17 END IF ELSE IF (Z1 > ZMAX) THEN IF (Z1 > ZZREAL) THEN IERROR = 7 IF (DEBUG) PRINT*, 'ERROR 7: Z1 > ZMAX' ELSE IERROR = 8 IF (DEBUG) WRITE (*,'(T2,2A)') 'Z1 ', TEXT17 END IF END IF END IF IF (Z1 < 0.0) THEN IERROR = 5 IF (DEBUG) PRINT*, 'ERROR 5: Z1 < 0.0' ELSE IF (Z2 < 0.0) THEN IERROR = 6 IF (DEBUG) PRINT*, 'ERROR 6: Z2 < 0.0' END IF END IF !********************************************************************** !abnormal exit !Note that only in the case of "hard errors" is ERROR set; otherwise, !it keeps the value it had in on entry. This is to allow the !possibility for calculation even if ERROR is /=0 on input. At present, !this is never the case, but if, for example, one wants to allow calculation !for ERROR = 3, one has only to remove the assignment to IERROR in that !case. Thus, we never set ERROR to 0 in this routine, not even for !an initialisation (since none is needed). For similar reasons, we use !IERROR and not ERROR as the test variable. In INICOS, ERROR is used !since this is more in keeping with the general usage of variables in !the parameter list and common module variables (or internal variables). !However, in INICOS, the error handling is always re-initialised and !ERROR and IERROR are always equal on exit. IF (IERROR /= 0) THEN D12 = 0.0D0 ERROR = IERROR RETURN END IF !********************************************************************** !if possible, calculate the angular size distance !********************************************************************** !the distance in the expanding universe !we are in the expanding phase now; set a flag in case ETAZ !depends on this VALUE2 = .FALSE. !some auxiliary variables for run time (not memory) optimisation VD = 1.0D0 + Z1 VS = 1.0D0 + Z2 !evaluate the stuff under the root (Q(V)) G = SQRT(QQ(VD)) IF (G <= TINY) THEN IERROR = 11 IF (DEBUG) PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE' END IF !initial values needed for the differential equation solver YSTART(1) = 0.0D0 !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 /= 0) THEN IERROR = 12 IF (DEBUG) PRINT*, 'ERROR 12: ERROR CALCULATING D12' !********************************************************************** !abnormal return: no further calculation if the !integration fails !********************************************************************** D12 = 0.0D0 ERROR = IERROR RETURN !********************************************************************** END IF D12 = D !set the output variable D12R = REAL(D12) !*********************************************************************** !*********************************************************************** END SUBROUTINE ANGSIZ !*********************************************************************** !*********************************************************************** !*********************************************************************** SUBROUTINE BNGSIZ(D14R,D34R,D32R,IERROR) ! ^^^^^^^^^^^^^^^^^^^^^ ! out !calculates the additional angular size distances (see !Kayser et al., A&A, 1996) for bounce models. !Given two redshifts Z1 and Z2 this subroutine calculates the !angular size distances between these two redshifts, which are !returned as D14, D34 and D32. In bounce models there are !two or four independent distances for the same redshift; BNGSIZ !returns as D14 the distance between two redshifts Z1 and Z2 such !that the universe is NOT in a state of CONTRACTION at Z1 but is !NOT in a state of EXPANSION at Z2. For D32 we have the reverse !situation; for D34 EXPANSION NEITHER at Z1 NOR at Z2. ! !If ETA(Z) is the same in the expanding and contracting phases, !then there are only two independent distances: D12 as calculated !by ANGSIZ and D14. In this case, D34 = D12 and D32 = D14. ! !The integer variable ERROR is 0 when no error is encountered, and !is otherwise set to a specific value for the corresponding error. !This can be useful in debugging. ! !BNGSIZ can only be called after calling ANGSIZ. If it is called !for a non-bounce model, and error message is printed. !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !Declarations !*********************************************************************** !output variables to the calling routine (parameter list) !error flag INTEGER, INTENT(OUT) :: IERROR !angular size distances--- real output variables REAL, INTENT(OUT) :: D14R, D34R, D32R !angular size distances---internal values DOUBLE PRECISION :: D14, D34, D32 !*********************************************************************** !output to called routines (parameter list) !does the actual integration: CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) !constants !One really needs a value this big to avoid numerical instability. !Nothing explodes or crashes with a too small value, but the result !can be incorrect. Although it is difficult to independently check !the distances calculated by BNGSIZ, they're OK as far as I can !tell and this value for NSTEP also makes the distance continuous !and smooth, which means values between checked values are probably OK. DOUBLE PRECISION, PARAMETER :: NSTEP=200.0D0 !variables !(1 + Z1) and (1 + Z2) DOUBLE PRECISION :: VD, VS !*********************************************************************** !output from called routines (parameter list) !does the actual integration: CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) DOUBLE PRECISION D !*********************************************************************** !internally used quantities !constants !just to check for overflow DOUBLE PRECISION TINY !This is roughly the smallest REAL number for a 4 byte REAL; !see the comment in ANGSIZ. PARAMETER (TINY=1.0D-38) !variables !just a temporary work variable DOUBLE PRECISION G !*********************************************************************** !*********************************************************************** !calculations !ERROR will be changed if something goes wrong IERROR = 0 !is a calculation possible !based on the cosmological model? !Changing the second condition to (IERROR <= 2) will suppress the !ERROR = 2 message here (though not in INICOS, of course). !Reporting is normally enabled, because the result is not !necessarily useful in this case, but the knowledgeable user might !wish to have this behaviour. IF ((ERROR >= 1) .AND. (ERROR <= 3)) THEN IF (DEBUG) WRITE (*, '(T2,A,I1,A)') 'ERROR ', ERROR, ': INICOS ERROR' END IF !is this a bounce model? IF (.NOT. BOUNCE) THEN IF (DEBUG) WRITE (*,'(T2,A,I2)') & 'ERROR 13: BNGSIZ INVALID FOR WMTYPE ', WMTYPE IERROR = 13 END IF !ANGSIZ OK? IF (((ERROR >= 4) .AND. (ERROR <= 10) ) .OR. (ERROR == 12)) THEN IF (DEBUG) WRITE (*, '(T2,A,I1,A)') 'ERROR ', ERROR, ': ANGSIZ ERROR' END IF IERROR = ERROR !no need to check the redshifts, since these were checked by ANGSIZ !Add .AND. (ERROR /= 3) if you really want it (see ANGSIZ). IF ( (ERROR /= 0) .AND. (ERROR /= 11) ) THEN !*********************************************************************** !abnormal return: no further calculation meaningful D14 = 0.0D0 D34 = 0.0D0 D32 = 0.0D0 ! add ERROR = -2 to documentation! IERROR = -2 RETURN !*********************************************************************** END IF !*********************************************************************** !if possible, calculate the angular size distance !calculate first to ZMAX... !we are in the expanding phase now; set a flag in case ETAZ !depends on this VALUE2 = .FALSE. IF (Z2 >= Z1) THEN !new initial values VD = 1.0D0 + Z2 VS = 1.0D0 + ZMAX !but the the initial values YSTART !are the result of the previous call to ANGSIZ ELSE !different new starting points, since we are integrating !starting at Z1 VD = 1.0D0 + Z1 VS = 1.0D0 + ZMAX !evaluate the stuff under the root (Q(V)) G = SQRT(QQ(VD)) IF (G <= TINY) THEN IERROR = 11 IF (DEBUG) PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE' END IF !and thus have the `normal' boundary conditions !initial values needed for the differential equation solver YSTART(1) = 0.0D0 !see Kayser et al., Eq. (36) for an explanation for the SIGN !YSTART(2) = SIGN((1.0D0/(VD*G)),(VS-VD)) YSTART(2) = 1.D00/(VD*G) !since (ZMAX >= Z1) and thus ((VS - VD) >= 0.0D0) END IF CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) IF (MSTEP /= 0) THEN IERROR = 12 IF (DEBUG) PRINT*, 'ERROR 12: ERROR CALCULATING D14, EXPANDING PHASE' !*********************************************************************** !abnormal return: no further calculation if the !integration fails D14 = 0.0D0 D34 = 0.0D0 D32 = 0.0D0 RETURN !*********************************************************************** END IF !...and then back to Z2 !we are in the contracting phase now; set a flag in case ETAZ !depends on this VALUE2 = .TRUE. VD = 1.0D0 + ZMAX VS = 1.0D0 + Z2 !call the solver for ZMAX--VS, the initial values YSTART !are the result of the previous call; the minus sign is needed !because we are reversing the integration direction YSTART(2) = -YSTART(2) CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) IF (MSTEP /= 0) THEN IERROR = 12 IF (DEBUG) PRINT*, 'ERROR 12: ERROR CALCULATING D14, CONTRACTING PHASE' !*********************************************************************** !abnormal return: no further calculation if the !integration fails D14 = 0.0D0 D34 = 0.0D0 D32 = 0.0D0 RETURN !*********************************************************************** END IF !this is our result for the second distance in bounce models D14 = D !*********************************************************************** !We're finished, unless ETA(Z) is different during contraction. !If ETA depends on Z only, that is, is the same in expanding !and contracting phases, then there are only two independent !distances, which we have already calculated. Otherwise, this !degeneracy vanishes, and we must calculate the other two. !Remember that a CHOICE value < 0 indicates that the dependence !of ETA on Z is different in the phases of expansion and !contraction in a bounce model. IF ( .NOT. ((CHOICE < 0) .AND. VARETA) ) THEN D34 = D12 D32 = D14 ELSE !*********************************************************************** !the third distance !we are in the contracting phase now; set a flag because ETAZ !depends on this VALUE2 = .TRUE. !some auxiliary variables for run time (not memory) optimisation VD = 1.0D0 + Z1 VS = 1.0D0 + Z2 !evaluate the stuff under the root (Q(V)) G = SQRT(QQ(VD)) IF (G <= TINY) THEN IERROR = 11 IF (DEBUG) PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE' END IF !initial values needed for the differential equation solver YSTART(1) = 0.0D0 !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 /= 0) THEN IERROR = 12 IF (DEBUG) PRINT*, 'ERROR 12: ERROR CALCULATING D34, CONTRACTING PHASE' !*********************************************************************** ! abnormal return: no further calculation if the ! integration fails D14 = 0.0D0 D34 = 0.0D0 D32 = 0.0D0 RETURN !*********************************************************************** END IF !this is our result for the third distance in bounce models D34 = D !*********************************************************************** !now for the fourth distance !calculate first to ZMAX... !we are in the contracting phase now; set a flag because ETAZ !depends on this VALUE2 = .TRUE. IF (Z2 >= Z1) THEN !new initial values VD = 1.0D0 + Z2 VS = 1.0D0 + ZMAX !but the the initial values YSTART !are the result of the previous call ELSE !different new starting points, since we are integrating !starting at Z1 VD = 1.0D0 + Z1 VS = 1.0D0 + ZMAX !evaluate the stuff under the root (Q(V)) G = SQRT(QQ(VD)) IF (G <= TINY) THEN IERROR = 11 IF (DEBUG) PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE' END IF !and thus have the `normal' boundary conditions !initial values needed for the differential equation solver YSTART(1) = 0.0D0 !see Kayser et al., Eq. (36) for an explanation for the SIGN !YSTART(2) = SIGN((1.0D0/(VD*G)),(VS-VD)) YSTART(2) = 1.0D0/(VD*G) !since (ZMAX >= Z1) and thus ((VS - VD) >= 0.0D0) END IF CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) IF (MSTEP /= 0) THEN IERROR = 12 IF (DEBUG) PRINT*, 'ERROR 12: ERROR CALCULATING D32, CONTRACTING PHASE' !*********************************************************************** !abnormal return: no further calculation if the !integration fails D14 = 0.0D0 D34 = 0.0D0 D32 = 0.0D0 RETURN !*********************************************************************** END IF !...and then back to Z2 !we are in the expanding phase now; set a flag because ETAZ !depends on this VALUE2 = .FALSE. VD = 1.0D0 + ZMAX VS = 1.0D0 + Z2 !call the solver for ZMAX--VS, the initial values YSTART !are the result of the previous call CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) IF (MSTEP /= 0) THEN IERROR = 12 IF (DEBUG) PRINT*, 'ERROR 12: ERROR CALCULATING D32, EXPANDING PHASE' !*********************************************************************** !abnormal return: no further calculation if the !integration fails D14 = 0.0D0 D34 = 0.0D0 D32 = 0.0D0 RETURN !*********************************************************************** END IF !this is our result for the fourth distance in bounce models D32 = D END IF !*********************************************************************** ! set the output variables D14R = REAL(D14) D34R = REAL(D34) D32R = REAL(D32) !*********************************************************************** !*********************************************************************** END SUBROUTINE BNGSIZ !*********************************************************************** !*********************************************************************** !*********************************************************************** SUBROUTINE LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) !This does the low level stuff for *NGSIZ, which is calculating the !angular size distance (and its derivative) between VD and VS. !This routine just does one leg of the integration; this is OK for !non-BOUNCE models; for BOUNCE models the second distance is put !together from separate calls to LOWLEV by *NGSIZ. DEBUG is passed !to enable informational messages to be printed by LOWLEV, since !this would mean redundant code in *NGSIZ. !NSTEP determines the first guessed step size. For the additional !distances in the BOUNCE case this must be much smaller for !numerical reasons. !YSTART are the starting values, which are replaced with the !results (needed when integrating further in the BOUNCE case). !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !Declarations !*********************************************************************** !input variables from the calling routine (parameter list) !start and end values DOUBLE PRECISION VD, VS, NSTEP LOGICAL DEBUG !*********************************************************************** !output variables to the calling routine (parameter list) !the distance DOUBLE PRECISION D !*********************************************************************** !output to called routines (parameter list) !all called routines ! EXTERNAL ODEINT, ASDRHS, BSSTEP ! EXTERNAL ASDRHS, BSSTEP ! CALL ODEINT(2,VD,VS,EPS,H1,0.0D0,BSSTEP,ASDRHS, !$ YSTART, NOK,NBAD) !2 is the dimension of YSTART !VD is the start value for the integration !VS is the end value for the integration !EPS is desired accuracy !H1 is a guessed first step size !the allowed minimum step size is 0.0 !BSSTEP is a stepper routine !ASDRHS calculates the right hand sides of the differential equations !YSTART are the starting values; the end values replace these !constants in the parameter list !accuracy DOUBLE PRECISION EPS PARAMETER (EPS=1.0D-6) !start values, guessed first step size DOUBLE PRECISION YSTART(2), H1 !*********************************************************************** !output from called routines (parameter list) !variables in the parameter list !NOK is the number of successful tries (not used in this version) !NBAD is the number of failed tries (not used in this version) INTEGER NOK, NBAD !*********************************************************************** !variables in common blocks !communication between here and BNGSIZ (q. v.) DOUBLE PRECISION YYSTRT(2) COMMON /LOWSIZ/ YYSTRT !communication between here and numerical integration stuff (q.v.) INTEGER MSTEP COMMON /MERROR/ MSTEP !save until the next call SAVE /LOWSIZ/, /MERROR/ !*********************************************************************** !internally used quantities !variables !for the do loop INTEGER I !*********************************************************************** !*********************************************************************** !run-time initialisation !no error yet MSTEP = 0 !*********************************************************************** !*********************************************************************** !calculations !for experimentally determined special cases IF (((VD <= 2.0D0) .OR. (VS < (1.0D0 + EPS))) .AND. (ABS(VS-VD) >= 3.6D0)) THEN ! set the guessed first step size H1 = 2.0D0 !in other cases ELSE !use a simply calculated guessed first step size H1 = (VS-VD)/NSTEP END IF !if D12 is practically = 0 IF (ABS(H1) < EPS/10.0D0) THEN ! set it to zero to avoid numerical problems D = 0.0D0 !otherwise calculate the angular size distance ELSE !call the solver CALL ODEINT(2,VD,VS,EPS,H1,0.0D0,BSSTEP,ASDRHS,YSTART, NOK,NBAD) IF (MSTEP /= 0) THEN IF (DEBUG) THEN IF (MSTEP == 1) THEN PRINT*, 'ERROR IN BSSTEP' ELSE IF (MSTEP == 2) THEN PRINT*, 'ERROR IN ODEINT: STEPSIZE TOO SMALL' ELSE IF (MSTEP == 3) THEN PRINT*, 'ERROR IN ODEINT: TOO MANY STEPS' ELSE PRINT*, 'CAN''T GET HERE' END IF END IF RETURN END IF !this is our result D = YSTART(1) END IF !*********************************************************************** ! set the variables needed by BNGSIZ DO, I = 1, 2, 1 YYSTRT(I) = YSTART(I) END DO !*********************************************************************** !*********************************************************************** END SUBROUTINE LOWLEV !*********************************************************************** !*********************************************************************** !*********************************************************************** SUBROUTINE ODEINT(NVAR,X1,X2,EPS,H1,HMIN,STEPR,DERIVS, & ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! in ! YINOUT, & ! ^^^^^^ ! inout ! NOK,NBAD) ! ^^^^^^^^ ! out !driver routine for integration of differential equations. !NVAR starting values are integrated from X1 to X2 with accuracy !EPS. H1 is a guessed first step size, HMIN the minimum allowed !step size, which can be zero. STEPR is the name of the stepper !routine to be used. DERIVS is a user-supplied routine for !calculating the right hand sides. ! !The starting values YINOUT are replaced with the results. ! !NOK and NBAD are the numer of good and bad steps taken. ! !Based on ideas from NUMERICAL RECIPES. !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !declarations !*********************************************************************** !named constants as array bounds INTEGER KMAXX, NMAX PARAMETER (KMAXX=200, NMAX=50) !*********************************************************************** !input variables in the parameter list !diemension of vector INTEGER NVAR !accuracy, guessed first step size, minimum step size, starting !and end point of integration, dependent variable DOUBLE PRECISION EPS, H1, HMIN, X1, X2, YINOUT(NVAR) !*********************************************************************** !output variables in the parameter list !number of good and bad steps for bookkeeping INTEGER NOK, NBAD !*********************************************************************** !output to called routines !all called routines !stepper routine, evaluation of right-hand sides EXTERNAL STEPR, DERIVS !variables !x an y for evaluation of right-hand sides, this result, step, !scaling vector for stepper routine DOUBLE PRECISION X, Y(NMAX), DYDX(NMAX), H, YSCAL(NMAX) !*********************************************************************** !output from called routines !number of successful steps and next step size DOUBLE PRECISION HDID, HNEXT !*********************************************************************** !variables in common blocks !bookkeeping (not used in this version, but put in in case someone !wants to debug or whatever) INTEGER KMAX, KOUNT DOUBLE PRECISION DXSAV, XP(KMAXX), YP(NMAX,KMAXX) COMMON /PATH/ DXSAV, XP, YP, KMAX, KOUNT !this is to keep track of (very rare) errors here and in related !routines INTEGER MSTEP COMMON /MERROR/ MSTEP !save until the next call SAVE /MERROR/ !*********************************************************************** !internally used quantities !named constants !maximum number of steps INTEGER MAXSTP !just for numerical stability; see below DOUBLE PRECISION TINY PARAMETER (MAXSTP=10000, TINY=1.0D-30) !variables !DO loops INTEGER I, NSTP !for bookkeeping DOUBLE PRECISION XSAV !*********************************************************************** !*********************************************************************** !run-time initialisation !reset stepsize error to zero whenever routine is called MSTEP = 0 !KMAX not initialised in DATA since it's in a common block; !BLOCK DATA not worth the effort KMAX = 0 !DXSAV isn't needed if KMAX = 0 !* DXSAV = 0 X = X1 H = SIGN(H1,X2-X1) NOK = 0 NBAD = 0 KOUNT = 0 !*********************************************************************** !*********************************************************************** !calculations !fill up work vector DO, I = 1, NVAR Y(I) = YINOUT(I) END DO !uncomment the next line if KMAX isn't 0; for KMAX = 0 neither !XSAV nor DXSAV is needed, so we don't need the next line !(This assures storage of the first step.) !IF (KMAX > 0) XSAV = X - DXSAV*2.0 !at most NMAX steps DO, NSTP = 1, MAXSTP !get the right-hand sides CALL DERIVS(X,Y, DYDX) DO, I = 1, NVAR !scaling to monitor accuracy YSCAL(I) = ABS(Y(I)) + ABS(H*DYDX(I)) + TINY END DO !storage for internal bookkeeping IF (KMAX > 0) THEN IF ( (ABS(X - XSAV)) > (ABS(DXSAV)) ) THEN IF( KOUNT < (KMAX - 1) ) THEN KOUNT = KOUNT + 1 XP(KOUNT) = X DO, I = 1, NVAR YP(I,KOUNT) = Y(I) END DO XSAV = X ENDIF ENDIF ENDIF !make sure the step size isn't too large IF ( ((X + H - X2)*(X + H - X1)) > 0.0D0 ) H = X2 - X !*********************************************************************** !abnormal return; set a flag and get out IF (MSTEP == 1) RETURN !This happens if there has been a failure in the stepper routine !on the previous call. !*********************************************************************** !call the stepper if everything's OK CALL STEPR(NVAR,DYDX,H,EPS,YSCAL,DERIVS, X,Y, HDID,HNEXT) IF (HDID == H)THEN NOK = NOK + 1 ELSE NBAD = NBAD + 1 ENDIF !finished? IF ( ((X - X2)*(X2 - X1)) >= 0.0D0 ) THEN !replace input values with output values DO, I = 1, NVAR YINOUT(I) = Y(I) END DO IF (KMAX /= 0) THEN !save final step KOUNT = KOUNT + 1 XP(KOUNT) = X DO, I = 1, NVAR YP(I,KOUNT) = Y(I) END DO ENDIF !*********************************************************************** !normal return; we have to get out now, since the error !conditions which follow this one are irrelevant if we're !successful, which means we've reached the next line RETURN ENDIF !*********************************************************************** !abnormal return (stepsize smaller than minimum); !set a flag and get out IF( (ABS(HNEXT)) < HMIN ) THEN MSTEP = 2 RETURN END IF !*********************************************************************** !set the next H value H = HNEXT END DO !*********************************************************************** !abnormal return (too many steps); set a flag and get out MSTEP = 3 RETURN !*********************************************************************** !*********************************************************************** END SUBROUTINE ODEINT !*********************************************************************** !*********************************************************************** !*********************************************************************** SUBROUTINE MMM(XIN,NVAR,YIN,DYDX,HTOT,NSTEP,DERIVS, YOUT) ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ^^^^ ! in out !implements the modified midpoint method for use in Bulirsch-Stoer !integration of differential equations. !Input at XIN are NV dependent variables YIN and their derivatives !DYDX as well as the total step to be taken, HTOT, and the number !of substeps NSTEP. DERIVS is an external routine for evaluating !the right hand sides. ! !The output is returned as YOUT. ! !Based on ideas from NUMERICAL RECIPES. !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !declarations !*********************************************************************** !named constants for array bounds INTEGER NMAX PARAMETER (NMAX=50) !*********************************************************************** !input variables in the parameter list !dimension of vector and number of substeps to be used INTEGER NVAR, NSTEP !independent variable, input vector and deriviative vector, total !step to be taken DOUBLE PRECISION XIN, YIN(NVAR), DYDX(NVAR), HTOT !*********************************************************************** !output variables in the parameter list !output vector DOUBLE PRECISION YOUT(NVAR) !*********************************************************************** !output to called routines !all called routines EXTERNAL DERIVS !variables !x and y for external evaluation routine DOUBLE PRECISION X, Y2(NMAX) !*********************************************************************** !internally used quantities !variables !for DO loops INTEGER I, N !increment, 2 temporaries, second vector variable (use see below) DOUBLE PRECISION H, H2, SWAP, Y1(NMAX) !*********************************************************************** !*********************************************************************** !calculations !set the stepsize H = HTOT/NSTEP !initialise the work vectors for the first step DO, I = 1, NVAR Y1(I) = YIN(I) Y2(I) = YIN(I) + H*DYDX(I) END DO !increment X X = XIN + H !temporarily store derivatives in YOUT CALL DERIVS(X,Y2, YOUT) !just to save some work in the DO loop H2 = 2.0D0*H !basic step DO, N = 2, NSTEP DO, I = 1, NVAR SWAP = Y1(I) + H2*YOUT(I) Y1(I) = Y2(I) Y2(I) = SWAP END DO !increment X X = X + H !derivative evaluation CALL DERIVS(X,Y2, YOUT) END DO !this is the last step DO, I = 1, NVAR YOUT(I) = 0.5D0*(Y1(I) + Y2(I) + H*YOUT(I)) END DO !*********************************************************************** !*********************************************************************** END SUBROUTINE MMM !*********************************************************************** !*********************************************************************** !*********************************************************************** SUBROUTINE POLEX(NV,XEST,YEST,NCALL,NUSE, YOUT,DY) ! ^^^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^ ! in out !extrapolates functions to X=0. !NV functions are evaluated at X = 0 by polynomial extrapolation !to a progressively smaller sequence of estimates XEST and the !function vectors YEST. This call is number NCALL. Use at most !NUSE estimates. ! !Extrapolated values are output as YOUT and their estimated error !as DY. ! !Based on ideas from NUMERICAL RECIPES. !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !declarations !*********************************************************************** !named constants as array bounds !maximum expected values of NCALL and NV INTEGER IMAX, NMAX PARAMETER (IMAX=13, NMAX=50) !*********************************************************************** !input variables from calling routine (parameter list) !number of functions, number of call, maximum number of previous !estimates to be used INTEGER NV, NCALL, NUSE !independent and dependent variables DOUBLE PRECISION XEST, YEST(NV) !*********************************************************************** !output variables to calling routine (parameter list) !output variable and estimated uncertainty DOUBLE PRECISION YOUT(NV), DY(NV) !*********************************************************************** !internally used quantities !variables !DO loops, calculated loop limit INTEGER J, K, MM !independent variable, table, work variables DOUBLE PRECISION X(IMAX), QCOL(NMAX,IMAX), D(NMAX), DELTA, F1, F2, Q !save until the next call SAVE QCOL, X !*********************************************************************** !*********************************************************************** !calculations X(NCALL) = XEST DO, J = 1, NV DY(J) = YEST(J) YOUT(J) = YEST(J) END DO !store first estimate in the first column IF(NCALL==1) THEN DO, J = 1, NV QCOL(J,1) = YEST(J) END DO ELSE !use at most NUSE previous estimates MM = MIN(NCALL,NUSE) DO, J = 1, NV D(J) = YEST(J) END DO DO, K = 1, (MM - 1) DELTA = 1.0D0/(X(NCALL - K) - XEST) F1 = XEST*DELTA F2 = X(NCALL - K)*DELTA ! propagate tableau one diagonal more DO, 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) END DO END DO DO, J = 1, NV QCOL(J,MM) = DY(J) END DO ENDIF !*********************************************************************** !*********************************************************************** END SUBROUTINE POLEX !*********************************************************************** !*********************************************************************** !*********************************************************************** DOUBLE PRECISION FUNCTION ETAZ(Z) !Calculates ETA as a function of Z. If desired, this value of ETA !can be used by ANGSIZ instead of a constant value. !This routine can be modified by the user to suit his or her needs. !Neither INICOS, ANGSIZ nor BNGSIZ need be modified. ! !Declaring the same COMMON block WHICH in the routine which calls !ANGSIZ makes it possible to investigate different behavior !of ETA(Z) without modifying the code. ! !This file in its original version is intended as an example. !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !declarations !*********************************************************************** !input variable from the calling routine DOUBLE PRECISION Z !*********************************************************************** !variables in common blocks !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 !This can also be declared in the routine which calls ANGSIZ, thus !making a decision about which z-dependence for ETA to use in this !calling routine. USE NEGATIVE CHOICE VALUES IF THE VALUE OF ETA !FOR THE SAME REDSHIFT CAN BE DIFFERENT IN THE CONTRACTING AND !EXPANDING PHASES OF A BOUNCE MODEL. (See VALUE2 below.) ! INTEGER CHOICE ! COMMON /WHICH/ CHOICE !LOGICAL VALUE2 !This can be used to find out if the integration is being done !in the contracting phase of a bounce model, in which case this !could lead to a different value for ETAZ than in the expanding !phase. ZMAX from the previous ANGINI COMMON block might be !useful in this case as well. !VALUE2 is .TRUE. during the contracting phase of a bounce model. !COMMON /CNTRCT/ VALUE2 !save until the next call ! SAVE /ANGINI/, /WHICH/, /CNTRCT/ !SAVE /CNTRCT/ !*********************************************************************** !internally used variable !just to keep things neat DOUBLE PRECISION Y !*********************************************************************** !*********************************************************************** !calculations Y = 1.0D0/(1.0D0 + Z) IF (CHOICE == 0) THEN ETAZ = 0.0D0 ELSE IF (CHOICE == -1) THEN ! Different ETA values before and after the bounce; ! the negative value of CHOICE is necessary for ANGSIZ, since ! some degeneracy is cancelled if ETAZ depends on VALUE2. IF (VALUE2) THEN ETAZ = 1.0D0 ELSE ! ETAZ = 0.0D0 ETAZ = 1.0D0 END IF ELSE IF (CHOICE == -2) THEN ! different ETA values before and after the bounce, continuous at ! ZMAX IF (VALUE2) THEN ETAZ = 1.0D0/(1.0D0 + ZMAX) ELSE ETAZ = Y END IF ELSE IF (CHOICE == 3) THEN ETAZ = 1.0D0 - Y ELSE IF (CHOICE == 4) THEN ETAZ = Y**2 ELSE IF (CHOICE == 5) THEN ETAZ = 1.0D0 - Y**2 ELSE IF (CHOICE == 6) THEN IF (Z <= 9) THEN ETAZ = LOG10(1.0D0 + Z) ELSE ETAZ = 1.0D0 END IF ELSE IF (CHOICE == 9) THEN ETAZ = SQRT(Y)*(1.0D0 - Y)**2 ELSE ETAZ = SQRT(1.0D0 - Y)*Y**2 END IF !take care of roundoff error IF (ETAZ > 1.0D0) THEN ETAZ = ETAZ - 0.0001D0 ELSE IF (ETAZ < 0.0D0) THEN ETAZ = ETAZ + 0.0001D0 END IF IF ( (ETAZ < 0.0D0) .OR. (ETAZ > 1.0D0) ) THEN !if programmed correctly, one can't get here, not even with !roundoff error PRINT*, 'SOMETHING''S WRONG!!!!!!' !just to make sure one is careful in the definition of ETAZ STOP END IF !*********************************************************************** !*********************************************************************** END FUNCTION ETAZ !********************************************************************** !*********************************************************************** !*********************************************************************** SUBROUTINE BSSTEP(NV,DYDX,HTRY,EPS,YSCAL,DERIVS, & ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! in ! X,Y, & ! ^^^ ! inout ! HDID,HNEXT) ! ^^^^^^^^^^ ! out ! USE HELBIG_ANGSIZ_OLD, ONLY: MMM, POLEX !Bulirsch-Stoer stepper routine for numerical integration of !differential equations. !NV dependent variables Y and their derivatives DYDX are input at X !along with the attempted stepsize HTRY, the required accuracy EPS !and the scaling vector YSCAL. DERIVS is a routine for evaluating !the right-hand sides. !X and Y are replaced with their output values. ! !HDID is the accomplished stepsize, HNEXT the next one. ! !Based on ideas from NUMERICAL RECIPES. (The Second Edition !discusses more involved techniques. They are not used here !because tests showed that for this application the difference in !computing time was negligible; thus I opted for the shorter and !clearer code.) !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !declarations !*********************************************************************** !named constants as array bounds INTEGER NMAX, IMAX PARAMETER (NMAX=10, IMAX=11) !*********************************************************************** !input variables from the calling routine (parameter list) !number of dependent variables INTEGER NV !derivatives, attempted stepsize, accuracy, scaling vector, !independent and dependent variables DOUBLE PRECISION DYDX(NV), HTRY, EPS, YSCAL(NV), X, Y(NV) !*********************************************************************** !output variables to the calling routine (parameter list) !accomplished stepsize, next stepsize DOUBLE PRECISION HDID, HNEXT !*********************************************************************** !output to called routines (parameter list) !all called routines !modified midpoint method, polynomial extrapolation, right-hand !sides ! EXTERNAL MMM, POLEX, DERIVS EXTERNAL DERIVS !named constants !number of previous estimates extrapolation should use INTEGER NUSE PARAMETER (NUSE=7) !variables !sequence of number of steps INTEGER NSEQ(IMAX) !x, y, stepsize and derivatives for modified mitpoint, x for !polynomial extrapolation DOUBLE PRECISION XSAV, YSAV(NMAX), DYSAV(NMAX), H, XEST !*********************************************************************** !output from called routines (parameter list) !output vector from polynomial extrapolation and its error DOUBLE PRECISION YSEQ(NMAX), YERR(NMAX) !*********************************************************************** !variables in common blocks !to let routines higher up now if something goes wrong INTEGER MSTEP COMMON /MERROR/ MSTEP !save until the next call SAVE /MERROR/ !*********************************************************************** !internally used quantities !named constants !ONE mustn't necessarily be `mathematically equivalent' to 1.0; !it's a parameter like the others---one could choose another value. !Has nothing to do with making it easier to code DOUBLE PRECISION !constants. DOUBLE PRECISION ONE, SHRINK, GROW PARAMETER (ONE=1.D0, SHRINK=.95D0, GROW=1.2D0) !variables !DO loops INTEGER I, J !for error checking DOUBLE PRECISION ERRMAX !*********************************************************************** !*********************************************************************** !initialisation of variables !internally used quantities DATA NSEQ /2,4,6,8,12,16,24,32,48,64,96/ !*********************************************************************** !*********************************************************************** !run-time initialisation MSTEP = 0 H = HTRY XSAV = X !*********************************************************************** !*********************************************************************** !calculations !save the starting values DO, I = 1, NV YSAV(I) = Y(I) DYSAV(I) = DYDX(I) END DO BIG_LOOP: DO !evaluate sequence of modified midpoint integrations DO, I = 1, IMAX !call the routine for modified midpoint method CALL MMM(XSAV,NV,YSAV,DYSAV,H,NSEQ(I),DERIVS, YSEQ) !series error is even, so squared XEST = (H/NSEQ(I))**2 !extrapolation with polynomials CALL POLEX(NV,XEST,YSEQ,I,NUSE, Y,YERR) !monitor local truncation error ERRMAX = 0.0D0 DO, J = 1, NV ERRMAX = MAX(ERRMAX,ABS(YERR(J)/YSCAL(J))) END DO !scale accuracy relative to tolerance ERRMAX = ERRMAX/EPS !has the step converged? IF (ERRMAX < ONE) THEN !set the return variables X = X + H HDID = H IF (I == NUSE) THEN HNEXT = H*SHRINK ELSE IF (I == NUSE - 1) THEN HNEXT = H*GROW ELSE HNEXT = (H*NSEQ(NUSE - 1))/NSEQ(I) ENDIF !*********************************************************************** ! normal return RETURN ENDIF !*********************************************************************** !get this far if the step fails END DO !Here it doesn't matter, with IMAX=11 and NUSE=7; !if one changes these values, one should think about whether !or not the integer division in the exponent is really wanted. !(In practice, one will probably never even get here.) !Anyway, this reduces the step size for another try if the step !fails. H = 0.25D0*H/2**((IMAX - NUSE)/2) !If H is too small, then set a flag and get out. IF (X + H == X) THEN MSTEP = 1 !*********************************************************************** ! abnormal return RETURN END IF !*********************************************************************** END DO BIG_LOOP !*********************************************************************** !*********************************************************************** END SUBROUTINE BSSTEP !*********************************************************************** !*********************************************************************** !*********************************************************************** SUBROUTINE ASDRHS(V,D, DDDV) ! ^^^ ^^^^ ! in out ! USE HELBIG_ANGSIZ_OLD, ONLY: QQ, ETAZ !calculates the Right Hand Side of the differential equation !for the Angular Size Distance !Input are the dependent variables V and D. ! !Output are the right-hand sides. !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !Major rewrite (nice Fortran95 style and features) July 2007. !*********************************************************************** !*********************************************************************** !declarations !*********************************************************************** !input variables from the calling routine (parameter list) !V = 1 + Z and is the basic independent variable, !D(1) is the distance and D(2) its derivative. DOUBLE PRECISION V, D(2) !*********************************************************************** !output variables to the calling routine (parameter list) !these are the result of this SUBROUTINE DOUBLE PRECISION DDDV(2) !*********************************************************************** !output variables from called routines (parameter list) ! G = SQRT(QQ(V)) ! IF (VARETA) ETA = ETAZ(V-1.0) !customised ETA, Q(v)**2 ! DOUBLE PRECISION ETAZ, QQ ! EXTERNAL ETAZ, QQ !*********************************************************************** !variables in common blocks !cosmological variables (see INICOS) ! DOUBLE PRECISION LAMBDA, OMEGA, ETA, K ! COMMON /COSMOL/ LAMBDA, OMEGA, ETA, K !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 ! save until the next call ! SAVE /COSMOL/, /ANGINI/ !*********************************************************************** !internally used quantities !variables !just to keep things neat DOUBLE PRECISION Q, G, GS !*********************************************************************** !*********************************************************************** !calculations Q = QQ(V) G = SQRT(Q) !and G's derivative GS = (3.0D0*OMEGA*V**2 - 2.0D0*K*V)/(2.0D0*G) !calculate the derivatives DDDV needed by the calling routine !D(1) and D(2) are D and D', so D' = D' (definition) DDDV(1) = D(2) !and D'' is calculated !if required replace ETA with the value given by ETAZ IF (VARETA) ETA = ETAZ(V-1.0D0) loops = loops + 1 DDDV(2) = -((3.0D0*OMEGA*ETA*V*D(1))/(2.0D0*Q)) - (2.0D0/V + GS/G)*D(2) !*********************************************************************** !*********************************************************************** END SUBROUTINE ASDRHS !*********************************************************************** !*********************************************************************** END MODULE HELBIG_ANGSIZ