Test your Fortran and cosmology skills!

Several years ago, I wrote some Fortran code which calculates the
angular-size distance (from which all other distances, except the
"distance" related to light-travel time, can be trivially derived) for a
very broad class of cosmological models (arbitrary lambda and Omega,
arbitrary smoothness parameter eta (which can even vary arbitrarily as a
function of redshift)):

http://www.astro.multivax.de:8000/ph...og/readme.html

The routines were originally written in Fortran77. I made a
quick-and-dirty conversion to fixed-form Fortran90 to allow them to be
USEd as a MODULE. From the point of view of using them as a black box,
both versions are probably good enough.

I've used them quite a bit myself as have other people (including Saul
Perlmutter's group, who used them for the paper which got them the 2011
Nobel Prize in physics). A few years ago, I found a bug (essentially a
typographical error) which probably no-one noticed since it affected
only a cosmological model which is a set of measure 0 (i.e. a curve in
the lambda-Omega plane) (had the buggy code actually been invoked, it
would have thrown another error).

I am now incorporating them into a much larger package, written in
Fortran95. An intermediate result will be a nicer Fortran95-version
with free-form source code and several changes to the code to make it
better, though the actual results produced will of course be the same.
This should be ready soon and I'll post it at the link above, which will
also contain a pointer to the new package when it is finished (hopefully
sometime this year). That will also be the last version of ANGSIZ,
since one can then use the new package which contains more
functionality.

As part of the rewrite, I noticed a real bug. I am certain that no-one
has been affected by it, since it can occur only in a very, very special
case. (I also noticed a typo which has no functional effect.) I'll
post the new version when it's ready (within a few days or at most
within a few weeks). However, it would be interesting if anyone can
find the bug or the typo (which are present in all versions at the web
page above). (Again, the circumstances required to invoke it are so
special that I'm sure no-one has been bitten by it.)

So, test your knowledge of Fortran and cosmology and let me know if you
find the bug or the typo!

On 6/2/12 8:22 AM, Phillip Helbig---undress to reply wrote:


So, test your knowledge of Fortran and cosmology and let me know if you
find the bug or the typo!

As a Pascal programmer, I find the Fortran idiosyncrasies almost
unfathomable in reproducing from code the original physical concepts.

Pascal

procedure 1
procedure 2
....
....
procedure n

main program

In any case
it would be helpful to have the equations
followed by a flow chart
and then the program.

In general, you appear to be quantifying distortions effects of 1+z

Richard D. Saam

On 6/11/2012 9:47 AM, Richard D. Saam wrote:
On 6/2/12 8:22 AM, Phillip Helbig---undress to reply wrote:


So, test your knowledge of Fortran and cosmology and let me know if you
find the bug or the typo!

As a Pascal programmer, I find the Fortran idiosyncrasies almost
unfathomable in reproducing from code the original physical concepts.

As an occasional Fortran programmer I do not find this
code an example of good programming style (in Fortran,
Pascal, or any language). What I downloaded was:

http://www.astro.multivax.de:8000/ph...ANGSIZ_OLD.FOR

Another problem is that it does not do anything! When
compiling is attempted, using
gfortran HELBIG_ANGSIZ_OLD.FOR
there is an error (no main program present). And with
the -c option it does compile but then there will not
be any routine called..

So the first thing Phillip has to do is to add a main
program and a sample input file as an installation test.
This should be directly downloadable together with the
program. Then we could test it, more or less like:

gfortran HELBIG_ANGSIZ_OLD.FOR HELBIG_test.f90

a.out HELBIG_input.data JB_output.data

Perhaps it was not the latest version I downloaded,
but I did not find it very easy to read:
1) Unclear names of variables (descriptive names
are extremely important for maintenance.)
2) Too many comments are only duplicating what should
already be clear from the lines of code (especially
if the previous point corrected)

Avoiding these things and using more powerful modern
Fortran syntax usually gives you 3 times shorter code
which is 3 times easier (and faster) to read.

Using another language than Fortran does not help, that
will at best give the same advantages as using modern
Fortran features instead of the old f77 subset (at best,
I write! Usually it doesn't.)

--
Jos

In article , "Richard D. Saam"
writes:

So, test your knowledge of Fortran and cosmology and let me know if you
find the bug or the typo!

As a Pascal programmer, I find the Fortran idiosyncrasies almost
unfathomable in reproducing from code the original physical concepts.

The original version is Fortran77. There is a quick-and-dirty
conversion to Fortran90 (which is not completely standard but probably
works with all compilers). As part of integrating this with something
else, I've produced a much nicer Fortran95 version which I will release
soon.

In any case
it would be helpful to have the equations
followed by a flow chart
and then the program.

That would be too easy. :-)

In article , Jos Bergervoet
writes:

As an occasional Fortran programmer I do not find this
code an example of good programming style (in Fortran,
Pascal, or any language).

As I noted in another post, the Fortran77 version has the limitations of
Fortran77. :-| Believe me, there is MUCH worse Fortran code out there.

What I downloaded was:

http://www.astro.multivax.de:8000/ph...ANGSIZ_OLD.FOR

Another problem is that it does not do anything! When
compiling is attempted, using
gfortran HELBIG_ANGSIZ_OLD.FOR
there is an error (no main program present). And with
the -c option it does compile but then there will not
be any routine called..

Right. You have to write the main program yourself.

So the first thing Phillip has to do is to add a main
program and a sample input file as an installation test.
This should be directly downloadable together with the
program. Then we could test it, more or less like:

gfortran HELBIG_ANGSIZ_OLD.FOR HELBIG_test.f90

a.out HELBIG_input.data JB_output.data

Soon I should have a nicer Fortran95 version and indeed will do
something like this since I am interested in results from different
platforms. I tested the original version on a large number of
platforms.

Perhaps it was not the latest version I downloaded,
but I did not find it very easy to read:
1) Unclear names of variables (descriptive names
are extremely important for maintenance.)

Fortran77. :-(

2) Too many comments are only duplicating what should
already be clear from the lines of code (especially
if the previous point corrected)

Better too many than too few comments!

Avoiding these things and using more powerful modern
Fortran syntax usually gives you 3 times shorter code
which is 3 times easier (and faster) to read.

Wait for the F95 version!

On 6/13/2012 1:00 PM, Phillip Helbig---undress to reply wrote:
..
As I noted in another post, the Fortran77 version has the limitations of
Fortran77. :-| Believe me, there is MUCH worse Fortran code out there.

In contrast to C, there never was an "obfuscated FORTRAN
contest" because it simply was too easy.. :^)

..
gfortran HELBIG_ANGSIZ_OLD.FOR HELBIG_test.f90

a.out HELBIG_input.data JB_output.data

Soon I should have a nicer Fortran95 version and indeed will do
something like this

Will wait..

2) Too many comments are only duplicating what should
already be clear from the lines of code (especially
if the previous point corrected)

Better too many than too few comments!

No, not in general. (Perhaps if all variables have
completely meaningless names, but you will avoid
that, hopefully..) Compare it with a physics
textbook: You really don't want a lengthy
description *in words* that the next equation is
equating some time derivative to the sum of the
spatial derivatives of three vector components
with indices cyclicly rotated, and then have that
followed by a short, crystal clear curl equation
for the Maxwell field!

What the code is doing should preferably be clear
form the lines of code. Comments should *only*
add the things we might not immediately realize.
(But of those comments there indeed often are
too few..)

--
Jos

Op zaterdag 2 juni 2012 15:22:29 UTC+2 schreef Phillip Helbig---undress to reply het volgende:
So, test your knowledge of Fortran and cosmology and let me know if you
find the bug or the typo!

In order to evaluate the subroutines I have "modified" the programs to be used in EXCEL i.e. Visual Basic.
For a copy and a description select:
http://users.telenet.be/nicvroom/World%20Model.htm

The subroutine INICOS is rather straightforward to be used and implemented.
For the programs ANGSIZ and BNGSIZ this is much more complicated.
Why are those programs written to be used in a cyclic universe (WMtype 17, 18 and 19) ? and not for WMtype = 13 ?
When I test them the results are always zero.

In the program INICOS I found some errors. See documentation.

Nicolaas Vroom
http://users.pandora.be/nicvroom/

On Mon, 25 Jun 12 07:06:30 GMT, Nicolaas Vroom wrote:
In order to evaluate the subroutines I have "modified" the programs to be used in EXCEL i.e. Visual Basic.
For a copy and a description http://users.telenet.be/nicvroom/World%20Model.htm

Wow, what an astounding job!

In the program INICOS I found some errors. See documentation.

You mean omega_k has the wrong sign, in which case the program works
correctly only for a flat universe. Aha, the *real* reason everybody
is so committed to a flat universe, because otherwise the darn
software -- does -- not -- work -- !!! :-)

On Monday, June 25, 2012 10:28:26 AM UTC+2, Eric Flesch wrote:
On Mon, 25 Jun 12 07:06:30 GMT, Nicolaas Vroom wrote:
In order to evaluate the subroutines I have "modified" the programs to be used in EXCEL i.e. Visual Basic.
For a copy and a description http://users.telenet.be/nicvroom/World%20Model.htm

Wow, what an astounding job!

That depents.
The biggest problem is to find the correct values to test the specific
WMtypes (Wold Model types). Specific for 18 and 19, because these two are
the ones (using Z1 and Z2) to test the subroutines ANGSIZ and BNGSIZ.
Specific for these two I need more intermediate results because the final
results are zero.

In the program INICOS I found some errors. See documentation.

You mean omega_k has the wrong sign, in which case the program works
correctly only for a flat universe. Aha, the *real* reason everybody
is so committed to a flat universe, because otherwise the darn
software -- does -- not -- work -- !!! :-)

I think you make the wrong conclusion.
The standard method is to use the equation:
Omega(Lambda) + Omega(M) + Omega(K) = 1
Phillip uses: Omega(Lambda) + Omega(M) - Omega(K) = 1
There is 'nothing' wrong with this except:
Standard meaning: when Omega(K) 0 then k = -1
In his case when Omega(K) 0 then k = -1

The biggest problem is the intricate relation
between the parameters WMtype, EPS and TINY.

For example IMO when you try Lambda = 2 Omega=0.5
you should get WMtype = 18 and ZMAX = 1
Instead you get WMtype = 17 and ZMAX = 0
For Lambda = 2 and Omega= 0.49999 ZMAX = 0.99486
For Lambda = 2 and Omega= 0.50001 ZMAX = 0

Nicolaas Vroom
http://users.pandora.be/nicvroom/

In article , Eric Flesch
writes:

You mean omega_k has the wrong sign, in which case the program works
correctly only for a flat universe. Aha, the *real* reason everybody
is so committed to a flat universe, because otherwise the darn
software -- does -- not -- work -- !!! :-)

I haven't checked yet, but I define K as LAMBDA + OMEGA - 1. I don't
call it OMEGA_K and it has the opposite sign of what some people call
OMEGA_K.