! comments ! remove 6 spaces ! $-->& ! longer lines ! ! add F95 features ! ! new concept for variables in common blocks ! ! public release ! ! add stuff to INICOS ! ! 1 module to replace all old common blocks, USE ONLY? ! ZZREAL ! !revision history ! ! 15 September 2000 ! Since 18-JUN-1997 18:21:11.00, nothing has been done. ! I am moving from version 1.0 to version 2.1 in one go. ! 1.0 to 2.0 implies adding stuff from Rollin Thomas. ! 2.0 to 2.1 implies adding quantities other than the angular size ! distance. ! I have also taken the opportunity to correct some typos IN COMMENTS ! and make some changes for a cleaner compile. Technically, these were ! mistakes, but in practice, the behaviour of the program was not ! affected. Chief among these is the fact that I got rid of all COMMON ! blocks. Previously, the COMMON variables were intialised in a BLOCK ! DATA routine. This is not possible in a module, so I had commented ! that out but left the initialisation and the common blocks. I have ! now gotten rid of the common blocks, which required some name changes ! in INICOS (but the new version is actually more understandable). !*********************************************************************** !*********************************************************************** !*********************************************************************** MODULE HELBIG_ANGSIZ_OLD ! contains the routines for calculating the angular size distance IMPLICIT NONE PRIVATE PUBLIC :: INICOS, ANGSIZ, BNGSIZ, CHOICE, RISK, LOOPS!& ! ETAZ, MMM, POLEX, QQ, ASDRHS, BSSTEP, LOWLEV, ODEINT !*********************************************************************** !*********************************************************************** !*********************************************************************** ! BLOCK DATA COSANG !initialises common blocks used by INICOS, ANGSIZ and related !routines !Code written by Phillip Helbig, Hamburg Observatory, Sep 1995. !*********************************************************************** !*********************************************************************** !declaration of variables in common blocks !the usual cosmological variables DOUBLE PRECISION LAMBDA, OMEGA, ETA, K !COMMON /COSMOL/ LAMBDA, OMEGA, ETA, K !internal names for other parameter list variables and internal !stuff; see INICOS INTEGER :: ERROR, WMTYPE, LOOPS=0 DOUBLE PRECISION ZMAX LOGICAL NOTYET, VARETA, DEBUG, MAXZ, BOUNCE !COMMON /ANGINI/ ZMAX, ERROR, WMTYPE, ! $ NOTYET, VARETA, DEBUG, MAXZ, BOUNCE !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. INTEGER CHOICE !COMMON /WHICH/ CHOICE !to allow computation in the Eddington model, if desired LOGICAL RISK !COMMON /ARTHUR/ RISK !sometimes ANGSIZ could be called with a Z which is smaller than !ZMAX but large enough for numerical problems. DOUBLE PRECISION ZZREAL !COMMON /FIX/ ZZREAL !save !SAVE /COSMOL/, /ANGINI/, /WHICH/, /ARTHUR/, /FIX/ !new !communication between ANGSIZ and BNGSIZ (q. v.) DOUBLE PRECISION ZZ1, ZZ2, DD12 ! separate module? !COMMON /BNCSIZ/ ZZ1, ZZ2, DD12 !communication between ANGSIZ and ETAZ (q. v.) LOGICAL VALUE2 !COMMON /CNTRCT/ VALUE2 !communication between here and numerical integration stuff (q.v.) INTEGER MSTEP !COMMON /MERROR/ MSTEP !*********************************************************************** !*********************************************************************** !compile-time initialisation of common block variables DOUBLE PRECISION DUMMY, IDUMMY !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 this original implementation !in now 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 declare the COMMON !block ARTHUR in the calling routine and set RISK to .TRUE., in !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 /-3, IDUMMY, DUMMY/, $ NOTYET, VARETA, DEBUG, MAXZ, BOUNCE /5*.TRUE./ DATA CHOICE /0/ DATA RISK /.FALSE./ DATA ZZREAL /DUMMY/ !*********************************************************************** !*********************************************************************** ! END BLOCKDATA COSANG !*********************************************************************** !*********************************************************************** 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 README.TEX file. !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. ! Originally, I had used different names in the common block here, ! in order to use the REAL names for the internal calculations, which ! were also the names in the parameter list. With the variables now ! defined in the whole module, I have to use the same names everywhere, ! so I use the strange ones in the parameter list and in the internal ! calculations. In other words, I have merely reversed the roles. ! One might think that I could use the common block variables for ! the internal stuff, but this was avoided to avoid setting them unless ! it was determined that the input variables were all correct (this is ! probably overkill). ! ! This is somewhat confused by the fact that the input variables ! are single precision, so I needed different names for the corresponding ! internal ones. ! The names in the comments and in some of the (commented out) ! warnings etc are probably still inconsistent. This will be ! corrected in the rewrite to Fortran95. !*********************************************************************** !*********************************************************************** !declarations !*********************************************************************** !input variables from the calling routine (parameter list) !the cosmological parameters---input variables REAL LAMBDR, OMEGAR, ETAR !the cosmological parameters---internal counterparts DOUBLE PRECISION LAMNUL, OMENUL, ETAETA !customised ETA? display messages? LOGICAL ETAVAR, DEEBUG !*********************************************************************** !output variables to the calling routine (parameter list) !type of cosmological model, error message INTEGER TYPEWM, IERROR !maximum redshift !real output variable REAL ZMAXR !internal counterpart DOUBLE PRECISION ZZMAX !is there a maximum redshift? is this a bounce model? LOGICAL MMAXZ, BOUNCY !*********************************************************************** !output from called routines (parameter list) !this is the quantity Q(v)**2 (see Kayser et al., A&A, 1996) ! DOUBLE PRECISION QQ ! EXTERNAL QQ !*********************************************************************** !variables in common blocks !with different names than elsewhere since the same quantities !appear in a parameter list as well; additional variables are !explained ! note that the actual input variable names are different ! SUBROUTINE INICOS(LAMBDA,OMEGA,ETA,VARETA,DEBUG, !$ WMTYPE,MAXZ,ZMAX,BOUNCE,ERROR) !K gives us a measure for the curvature (see below) ! DOUBLE PRECISION LAMBDA, OMEGA, ETA, K ! COMMON /COSMOL/ LAMBDA, OMEGA, ETA, K ! INTEGER ERROR, WMTYPE ! DOUBLE PRECISION ZMAX !have we called anything yet? ! LOGICAL NOTYET, VARETA, DEBUG, MAXZ, BOUNCE ! COMMON /ANGINI/ ZZMAX, IERROR, TYPEWM, ! $ NOTYET, ETAVAR, DEEBUG, MMAXZ, BOUNCY !sometimes ANGSIZ could be called with a Z which is smaller than !ZMAX but large enough for numerical problems. ! DOUBLE PRECISION ZZREAL ! COMMON /FIX/ ZZREAL !save until the next call !in order to avoid recalculating cosmological quantities when this !is unnecessary ! SAVE /COSMOL/, /ANGINI/, /FIX/ !*********************************************************************** !internally used quantities !named constants !$\frac{\pi}{3}$, accuracy DOUBLE PRECISION PITHRD, EPS !The value of EPS should be about right for a 24 bit mantissa; !see the note concerning precision in ANGSIZ. PARAMETER (PITHRD = 1.0471975511965977D0, EPS = 1.0D-6) !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 .LT. 0.0) THEN ERROR = 1 IF (DEEBUG) THEN PRINT*, 'INICOS ERROR 1: OMEGA < 0.0; OMEGA = ', OMENUL END IF ELSE IF ( (.NOT. ETAVAR) .AND. $ ((ETAETA .LT. 0.0) .OR. (ETAETA .GT. 1.0)) ) THEN ERROR = 2 IF (DEEBUG) THEN PRINT*, 'INICOS ERROR 2: ETA OUT OF RANGE; ETA = ', ETAETA END IF ELSE ERROR = 0 END IF IERROR = ERROR !*********************************************************************** !abnormal exit IF (ERROR .GT. 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? IF ( (LAMNUL .NE. LAMBDA) .OR. (OMENUL .NE. OMEGA) .OR. $ NOTYET ) THEN ! this only has an effect the first time NOTYET = .FALSE. ! print*, 'notyet: ', notyet ! the sign of k is the sign of the curvature K = OMENUL + LAMNUL - 1.0D0 IF (LAMNUL .LT. 0.0) THEN MAXZ = .FALSE. BOUNCE = .FALSE. IF (K .LT. 0.0) THEN IF (OMENUL .GT. 0.0) THEN WMTYPE = 2 ELSE ! empty universe WMTYPE = 1 END IF ELSE IF (K .GT. 0.0) THEN WMTYPE = 4 ELSE ! k = 0 WMTYPE = 3 END IF ELSE IF (LAMNUL .EQ. 0.0) THEN MAXZ = .FALSE. BOUNCE = .FALSE. IF (K .LT. 0.0) THEN IF (OMENUL .GT. 0.0) THEN WMTYPE = 8 ELSE ! empty universe, `Milne-model' WMTYPE = 7 END IF ELSE IF (K .GT. 0.0) THEN ! `Misner-Thorne-Wheeler model' WMTYPE = 5 ELSE ! k = 0, Einstein-de Sitter model WMTYPE = 9 END IF ELSE ! LAMNUL > 0.0 IF (K .LT. 0.0) THEN MAXZ = .FALSE. BOUNCE = .FALSE. IF (OMENUL .GT. 0.0) THEN WMTYPE = 11 ELSE ! empty universe WMTYPE = 10 END IF ELSE IF (K .EQ. 0.0) THEN MAXZ = .FALSE. BOUNCE = .FALSE. IF (OMENUL .GT. 0.0) THEN WMTYPE = 13 ELSE ! empty universe, de Sitter model WMTYPE = 12 END IF ELSE ! k > 0 IF (OMENUL .EQ. 0.0) 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) .LT. EPS) THEN WMTYPE = 12 MAXZ = .FALSE. BOUNCE = .FALSE. ELSE ! normal code for this case WMTYPE = 19 MAXZ = .TRUE. BOUNCE = .TRUE. ZMAX = SQRT(LAMNUL/(LAMNUL - 1.0D0)) - 1.0D0 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 .GT. EPS) THEN ! Lemaitre model WMTYPE = 14 MAXZ = .FALSE. BOUNCE = .FALSE. ELSE IF (ABS(T) .LT. EPS) THEN IF (Q .GT. 0.0) THEN ! A1 model, Q .GT. 0.5 WMTYPE = 15 MAXZ = .FALSE. BOUNCE = .FALSE. ELSE ! A2 model, Q .LT. -1, Eddington model WMTYPE = 17 MAXZ = .TRUE. BOUNCE = .FALSE. ZMAX = (LAMNUL/SIGMA)**(1.0D0/3.0D0) - 1.0D0 END IF ELSE ! T < 0.0 IF (Q .GT. 0.0) THEN ! q .ge. 0.5 WMTYPE = 6 MAXZ = .FALSE. BOUNCE = .FALSE. ELSE ! q .le. -1.0, bounce model WMTYPE = 18 MAXZ = .TRUE. BOUNCE = .TRUE. ! calculate ZMAX with standard formula R = SQRT(K/(LAMNUL*3.0D0)) 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/v; we need z ! R_0/R = 1 + z ZMAX = 1.0D0/MAX(Y1,Y2,Y3) - 1.0D0 END IF END IF END IF END IF END IF LAMBDA = LAMNUL OMEGA = OMENUL !same result, better style LAMBDA = REAL(LAMNUL) OMEGA = REAL(OMENUL) END IF !*********************************************************************** !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. ! !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 by a !COMMON block 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 ! print*, "ZZREAL, ZMAX: ", zzreal, zmax IF (MAXZ) THEN NBAD = 0 ! sorry, standard FORTRAN77 has no DO WHILE, so here's a label ! which one can GOTO 100 CONTINUE TEMP = QQ(ZMAX + 1.0D0) IF (TEMP .LT. EPS) THEN ZMAX = ZMAX - EPS*10.0D0 NBAD = NBAD + 1 ! 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. IF (NBAD .EQ. 100) THEN IF (DEEBUG) THEN PRINT*, 'INICOS ERROR 3: ROUNDOFF ERROR, BAD ZMAX' ! Uncomment the next three lines for a more verbose ! error message, useful since some output variables are ! 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 ! This is just to force an exit without introducing an ! additional label. ZMAX = 0.0D0 ERROR = 3 * 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 *! above. * ** 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', ** $ TYPEWM, LAMBDA, OMEGA,MMAXZ,ZZMAX,TEMP,NBAD ** ** PRINT* * * END IF * END IF IERROR = ERROR GOTO 100 END IF END IF !*********************************************************************** !tidy things up !input variables ! SUBROUTINE INICOS(LAMBDA,OMEGA,ETA,VARETA,DEBUG, C! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ C! in !set the common block variables to those in the parameter list !LAMNUL and OMENUL have been set to LAMBDA and OMEGA above !if these have changed (and are valid values) ETA = ETAETA VARETA = ETAVAR DEBUG = DEEBUG ! COMMON /COSMOL/ LAMNUL, OMENUL, ETAETA, K ! COMMON /ANGINI/ ZZMAX, IERROR, TYPEWM, !$ NOTYET, ETAVAR, DEEBUG, MMAXZ, BOUNCY !output variables !$WMTYPE,MAXZ,ZMAX,BOUNCE,ERROR) C!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ C! out !set the output variables to the common block values, which if !necessary (new LAMBDA or OMEGA) have been recomputed TYPEWM = WMTYPE MMAXZ = MAXZ ! COMMON /FIX/ ZZREAL !use the real ZZMAX for returning to the calling routine, a dummy !value if it doesn't exist IF (MMAXZ) THEN ZZMAX = ZZREAL ELSE ZZMAX = 0.0D0 END IF ! print*, "ZZREAL, ZMAX: ", zzreal, zmax 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. !*********************************************************************** !*********************************************************************** !Declarations !*********************************************************************** !input variables from the calling routine (parameter list) !(1 + z) DOUBLE PRECISION V !*********************************************************************** !variables in common block !cosmological variables ! DOUBLE PRECISION LAMBDA, OMEGA, ETA, K ! COMMON /COSMOL/ LAMBDA, OMEGA, ETA, K !save until the next call ! SAVE /COSMOL/ !*********************************************************************** !*********************************************************************** !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. !Changes to version 2.0 (September 2000): !Where possible, we now use Rollin Thomas's code as the faster !default, but this can be overridden by (a) an optional input !parameter (on a case-by-case basis) or (b) by a call to the !new routine SELECT_METHOD. !*********************************************************************** !*********************************************************************** !Declarations !*********************************************************************** !modules ! USE ROLLIN_Z2DL !*********************************************************************** !input variables from the calling routine (parameter list) !redshifts !the real input variables REAL Z1R, Z2R !the internally used counterparts DOUBLE PRECISION Z1, Z2 !*********************************************************************** !output variables to the calling routine (parameter list) !error flag INTEGER IERROR !angular size distance (D14, D34 and D32 are the other distances !in bounce models) !the actual output variable REAL D12R !the internally used counterpart DOUBLE PRECISION D12 !*********************************************************************** !output to called routines (parameter list) !all called routines !does the actual integration ! CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) ! EXTERNAL LOWLEV !constants DOUBLE PRECISION NSTEP !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. PARAMETER (NSTEP=2.0D0) !variables !(1 + Z1) and (1 + Z2) DOUBLE PRECISION VD, VS !dependent variables (starting values) for the integration routine DOUBLE PRECISION YSTART(2) !*********************************************************************** !output from called routines (parameter list) ! CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) ! G = SQRT(QQ(VD)) !angular size distance and Q(v)**2 ! DOUBLE PRECISION D, QQ DOUBLE PRECISION D ! EXTERNAL QQ !*********************************************************************** !variables in common blocks !communication between here and INICOS (q. v.) ! INTEGER IERROR, WMTYPE !ppp is IERROR and ERROR here and in INICOS correct and consistent? ! DOUBLE PRECISION ZMAX ! LOGICAL NOTYET, DEBUG, VARETA, MAXZ, BOUNCE ! COMMON /ANGINI/ ZMAX, IERROR, WMTYPE, ! $ NOTYET, VARETA, DEBUG, MAXZ, BOUNCE !communication between here and BNGSIZ (q. v.) ! DOUBLE PRECISION ZZ1, ZZ2, DD12 ! COMMON /BNCSIZ/ ZZ1, ZZ2, DD12 !communication between here and ETAZ (q. v.) ! LOGICAL VALUE2 ! COMMON /CNTRCT/ VALUE2 !communication between here and numerical integration stuff (q.v.) ! INTEGER MSTEP ! COMMON /MERROR/ MSTEP !to allow computation in the Eddington model, if desired ! LOGICAL RISK ! COMMON /ARTHUR/ RISK !sometimes ANGSIZ could be called with a Z which is smaller than !ZMAX but large enough for numerical problems. ! DOUBLE PRECISION ZZREAL ! COMMON /FIX/ ZZREAL !save until the next call ! SAVE /ANGINI/, /BNCSIZ/, /CNTRCT/, /MERROR/, /ARTHUR/, /FIX/ ! SAVE /BNCSIZ/, /CNTRCT/, /MERROR/ !*********************************************************************** !internally used quantities !constants !just to check for overflow DOUBLE PRECISION TINY !ppp use this intrinsic elsewhere?????? ! DOUBLE PRECISION :: EPS=EPSILON(1.0D0) CHARACTER TEXT17*56 !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 is 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. PARAMETER (TINY=1.0D-38, $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! IERROR = -1 !*********************************************************************** ! abnormal return: no further calculation in this case RETURN !*********************************************************************** END IF !based on the cosmological model? !Changing the second condition to (IERROR .LE. 2) will suppress the !ERROR = 3 message here (though not in INICOS, of course). !Reporting is normally enabled, because though the result is not !necessarily useful in this case, the knowledgeable user might !wish to have this behaviour. IF ((ERROR .GE. 1) .AND. (ERROR .LE. 3)) THEN IF (DEBUG) THEN WRITE (*, '(T2,A,I1,A)') 'ERROR ', ERROR, ': INICOS ERROR' END IF !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 .EQ. 17) .AND. .NOT. RISK) THEN IERROR = 4 IF (DEBUG) THEN PRINT*, 'ERROR 4: COMPUTATION NOT RECOMMENDED IN WMTYPE 17' END IF !based on the redshifts? ELSE ! print*, "ZZREAL, ZMAX: ", zzreal, zmax IF (MAXZ) THEN IF (Z2 .GT. ZMAX) THEN ! print*, "ZZREAL, ZMAX: ", zzreal, zmax IF (Z2 .GT. ZZREAL) THEN ! print*, "ZZREAL, ZMAX: ", zzreal, zmax IERROR = 9 ! print*, "ZZREAL, ZMAX: ", zzreal, zmax IF (DEBUG) PRINT*, 'ERROR 9: Z2 > ZMAX' ! print*, "ZZREAL, ZMAX: ", zzreal, zmax ELSE IERROR = 10 IF (DEBUG) WRITE (*,'(2A)') 'Z2 ', TEXT17 END IF ELSE IF (Z1 .GT. ZMAX) THEN IF (Z1 .GT. ZZREAL) THEN IERROR = 7 IF (DEBUG) PRINT*, 'ERROR 7: Z1 > ZMAX' ELSE IERROR = 8 IF (DEBUG) WRITE (*,'(2A)') 'Z1 ', TEXT17 END IF END IF END IF IF (Z1 .LT. 0.0) THEN IERROR = 5 IF (DEBUG) THEN PRINT*, 'ERROR 5: Z1 < 0.0' END IF ELSE IF (Z2 .LT. 0.0) THEN IERROR = 6 IF (DEBUG) THEN PRINT*, 'ERROR 6: Z2 < 0.0' END IF END IF END IF !********************************************************************** !abnormal exit !Add .AND. (IERROR .NE. 3) if you really want it (see above). ! print*, "error, ierror: ", error, ierror IF (IERROR .NE. 0) THEN D12 = 0.0D0 !OK but better style is to do the assignment in any case 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. !determine which method is to be used !this section will be modified as new methods become available !determine where Rollin's code can be used ! need to add ALL other cases here as they become available; ! is there a better way---also for the whole process of selection? !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 .LE. TINY) THEN IERROR = 11 !subroutine for error handling? IF (DEBUG) THEN PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE' END IF 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 .NE. 0) THEN IERROR = 12 !********************************************************************** ! 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) !set the variables needed by BNGSIZ ZZ1 = Z1 ZZ2 = Z2 DD12 = 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. !*********************************************************************** !*********************************************************************** !Declarations !*********************************************************************** !output variables to the calling routine (parameter list) !error flag INTEGER IERROR !angular size distances---internal values DOUBLE PRECISION D14, D34, D32 !angular size distances--- real output variables REAL D14R, D34R, D32R !*********************************************************************** !output to called routines (parameter list) !all called routines !does the actual integration ! CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) ! EXTERNAL LOWLEV !constants DOUBLE PRECISION NSTEP !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 !smooth and smooth, which means values between checked values are !probably OK. PARAMETER (NSTEP=200.0D0) !variables !(1 + Z1) and (1 + Z2) DOUBLE PRECISION VD, VS !dependent variables (starting values) for the integration routine DOUBLE PRECISION YSTART(2) !*********************************************************************** !output from called routines (parameter list) ! CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) ! G = SQRT(QQ(VD)) !angular size distance and Q(v)**2 ! DOUBLE PRECISION D, QQ DOUBLE PRECISION D ! EXTERNAL QQ !*********************************************************************** !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 !communication between here and ANGSIZ (q. v.) DOUBLE PRECISION Z1, Z2, D12 COMMON /BNCSIZ/ Z1, Z2, D12 !communication between here and LOWLEV (q. v.) !YSTART(2) declared above COMMON /LOWSIZ/ YSTART !communication between the calling routine, here and ETAZ (q. v.) !CHOICE is needed here, and not only in the user's calling routine !and ETA(Z), since 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. ! INTEGER CHOICE ! COMMON /WHICH/ CHOICE !communication between here and ETAZ (q. v.) LOGICAL VALUE2 COMMON /CNTRCT/ VALUE2 !communication between here and numerical integration stuff (q.v.) INTEGER MSTEP COMMON /MERROR/ MSTEP !save until the next call ! SAVE /ANGINI/, /BNCSIZ/, /LOWSIZ/, /WHICH/, /CNTRCT/, /MERROR/ SAVE /BNCSIZ/, /LOWSIZ/, /CNTRCT/, /MERROR/ !*********************************************************************** !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 .LE. 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 .GE. 1) .AND. (ERROR .LE. 3)) THEN IF (DEBUG) THEN WRITE (*, '(T2,A,I1,A)') 'ERROR ', ERROR, ': INICOS ERROR' END IF END IF !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 !ANGSIZ OK? IF ( ( (ERROR .GE. 4) .AND. (ERROR .LE. 10) ) .OR. $ ( ERROR .EQ. 12 ) ) THEN IF (DEBUG) THEN WRITE (*, '(T2,A,I1,A)') 'ERROR ', ERROR, ': ANGSIZ ERROR' END IF END IF IERROR = ERROR !no need to check the redshifts, since these were checked by ANGSIZ !Add .AND. (ERROR .NE. 3) if you really want it (see ANGSIZ). IF ( (ERROR .NE. 0) .AND. (ERROR .NE. 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 .GE. 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 .LE. TINY) THEN IERROR = 11 IF (DEBUG) THEN PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE' END IF 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 .GE. Z1) and thus ((VS - VD) .GE. 0.0D0) END IF CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) IF (MSTEP .NE. 0) THEN IERROR = 12 IF (DEBUG) THEN PRINT*, 'ERROR 12: ERROR CALCULATING D14, EXPANDING PHASE' END IF !*********************************************************************** ! 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 .NE. 0) THEN IERROR = 12 IF (DEBUG) THEN PRINT*, 'ERROR 12: ERROR CALCULATING D14, CONTRACTING PHASE' END IF !*********************************************************************** ! 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 .LT. 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 .LE. TINY) THEN IERROR = 11 IF (DEBUG) THEN PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE' END IF 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 .NE. 0) THEN IERROR = 12 IF (DEBUG) THEN PRINT*, 'ERROR 12: ERROR CALCULATING D34, CONTRACTING PHASE' END IF !*********************************************************************** ! 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 .GE. 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 .LE. TINY) THEN IERROR = 11 IF (DEBUG) THEN PRINT*, 'ERROR 11: OVERFLOW ERROR POSSIBLE' END IF 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 .GE. Z1) and thus ((VS - VD) .GE. 0.0D0) END IF CALL LOWLEV(DEBUG,VD,VS,NSTEP, YSTART, D) IF (MSTEP .NE. 0) THEN IERROR = 12 IF (DEBUG) THEN PRINT*, 'ERROR 12: ERROR CALCULATING D32, CONTRACTING PHASE' END IF !*********************************************************************** ! 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 .NE. 0) THEN IERROR = 12 IF (DEBUG) THEN PRINT*, 'ERROR 12: ERROR CALCULATING D32, EXPANDING PHASE' END IF !*********************************************************************** ! 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. !*********************************************************************** !*********************************************************************** !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.0,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. LE. 2.0D0) .OR. (VS .LT. (1.0D0 + EPS)) ) .AND. $ (ABS(VS-VD) .GE. 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) .LT. 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 .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 ! 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 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. !*********************************************************************** !*********************************************************************** !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 100 I = 1, NVAR Y(I) = YINOUT(I) 100 CONTINUE !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 .GT. 0) XSAV = X - DXSAV*2.0 !at most NMAX steps DO 600 NSTP = 1, MAXSTP ! get the right-hand sides CALL DERIVS(X,Y, DYDX) DO 200 I = 1, NVAR ! scaling to monitor accuracy YSCAL(I) = ABS(Y(I)) + ABS(H*DYDX(I)) + TINY 200 CONTINUE ! 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 ! make sure the step size isn't too large IF ( ((X + H - X2)*(X + H - X1)) .GT. 0.0 ) H = X2 - X !*********************************************************************** ! abnormal return; set a flag and get out IF (MSTEP .EQ. 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 .EQ. H)THEN NOK = NOK + 1 ELSE NBAD = NBAD + 1 ENDIF ! finished? IF ( ((X - X2)*(X2 - X1)) .GE. 0.0 ) THEN ! replace input values with output values DO 400 I = 1, NVAR YINOUT(I) = Y(I) 400 CONTINUE IF (KMAX .NE. 0) THEN ! save final step KOUNT = KOUNT + 1 XP(KOUNT) = X DO 500 I = 1, NVAR YP(I,KOUNT) = Y(I) 500 CONTINUE 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)) .LT. HMIN ) THEN MSTEP = 2 RETURN END IF !*********************************************************************** ! set the next H value H = HNEXT 600 CONTINUE !*********************************************************************** ! 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. !*********************************************************************** !*********************************************************************** !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 100 I = 1, NVAR Y1(I) = YIN(I) Y2(I) = YIN(I) + H*DYDX(I) 100 CONTINUE !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 300 N = 2, NSTEP DO 200 I = 1, NVAR SWAP = Y1(I) + H2*YOUT(I) Y1(I) = Y2(I) Y2(I) = SWAP 200 CONTINUE ! increment X X = X + H ! derivative evaluation CALL DERIVS(X,Y2, YOUT) 300 CONTINUE !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 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. !*********************************************************************** !*********************************************************************** !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 100 J = 1, NV DY(J) = YEST(J) YOUT(J) = YEST(J) 100 CONTINUE !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 !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 ! 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 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. !*********************************************************************** !*********************************************************************** !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 .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 !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 ! 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. !*********************************************************************** !*********************************************************************** !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 100 I = 1, NV YSAV(I) = Y(I) DYSAV(I) = DYDX(I) 100 CONTINUE !GOTO 1 allows one to repeat the main part of the routine if the !step fails. 1 CONTINUE ! evaluate sequence of modified midpoint integrations DO 200 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 12 J = 1, NV ERRMAX = MAX(ERRMAX,ABS(YERR(J)/YSCAL(J))) 12 CONTINUE ! scale accuracy relative to tolerance ERRMAX = ERRMAX/EPS ! has the step converged? IF (ERRMAX .LT. ONE) THEN ! 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 !*********************************************************************** ! normal return RETURN ENDIF !*********************************************************************** !get this far if the step fails 200 CONTINUE !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. EQ. X) THEN MSTEP = 1 !*********************************************************************** ! abnormal return RETURN END IF !*********************************************************************** !repeat if one gets this far GOTO 1 !*********************************************************************** !*********************************************************************** 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. !*********************************************************************** !*********************************************************************** !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_OLD