Galaxies in 3D

In this example you will read in a file containing data of around 22,000 galaxies, convert the data to a format that can be used by Maple, write it back to a file, and then read it in again to display the galaxies in 3-dimensional space depending on their distance from the Sun.

All the Maple sheets and source files are available at the bottom of this page.

This tutorial is not intended to be complete. You should read the online help pages on the various file I/O functions. A good starting point is ?fprintf for it includes links to all other related commands (or use the help browser /Programming/Input and Output/File Manipulation).

Note: Since a much larger galaxy catalogue is available on the Net, this document has been changed in April 1999. The resulting Maple plots now show much more detail than those published before on this page.

Table of contents

Introduction

This topic was inspired by an article in the August 1996 issue of ASTRONOMY magazine on 'A River In the Universe' by Adam Frank discussing findings that not all galaxies do not simply expand outward (driven by the forces of the Big Bang) but 'that a huge chunk of the universe appears to be heading off .. toward .. the direction of Virgo [and that there is] a vast current of matter, a bulk flow of galaxies superimposed on top of the universe's smooth expansion.'

The article mentioned above contains a plot showing that galaxies are not evenly distributed in space but that there are concentrations of galaxies in clusters forming walls billions of light years long. This plot also shows 'bubbles', vast amounts of nearly empty space with the galaxies situated on top of the bubbles.

Thanks to a catalogue of galaxies compiled and published in the Internet by John Huchra, it is possible to compute a 3dimensional plot of the galaxies with Maple using its file input and output commands. You will see `beams` not showing true structure but which are the result of exact surveys of small parts of the sky containing very dim galaxies (so-called `pencil-beam surveys`).

While the original file contains 38,909 galaxies plus a lot of data on each one, I have extracted only the necessary information (apart from the brightness data) to speed up computation. Also, all the galaxies where velocities were not available have been removed. Thus 22,746 galaxies are contained in the Maple input file.

Structure of the Input File

The Maple input file includes the following information:

File Structure
Data Format (length) Position
start		 Position
end
Name of galaxy Alphabetic (10) 1 10
Right ascension - hours Integer (2) 11 12
Right ascension - minutes Integer (2) 13 14
Right ascension - seconds Float (4.1) 15 18
Declination - sign Alphabetic (1) 19 19
Declination - degrees Integer (2) 20 21
Declination - minutes Integer (2) 22 23
Declination - seconds Integer (2) 24 25
Brightness (magnitude) Float (5.2) 26 30
Heliocentric velocity Integer (5) 31 35

The following chart shows the position of the data and the first four lines in the file: 

0        1         2         3
12345678901234567890123456789012345
-----------------------------------
0000+0810 0000 0.0+081000 0.0028663
0000-3614 0000 3.5-361454 0.0034560
0000-3612 0000 7.5-361250 0.0014881
0000-3618 0000 9.6-361808 0.0015356
For a Maple plot the names and magnitudes of the galaxies are not essential. The corresponding columns can be skipped easily using the substring function. Thus we need the data from columns

Reading in lines: readline

The built-in readline function is used to read in each line from a text file. The return of readline is a string. readline does not automatically read in the entire file, it only reads one line and in order to get the next line, it must be called again, preferrably in a loop. (For a function that reads in the entire file, download the utils package.) If it is called the first time, readline automatically opens the file; likewise, if it has reached the end of the file, it is automatically closed.

Syntax:

readline
Calling sequence: readline(filename);
Parameter: filename: the name of the file to be read line by line (a string)
Return: the line as a string, or 0 if end of file has been reached

We must split the resulting string into a sequence of different substrings, convert the sequence into a list by enclosing it in square brackets, and then convert each substring in the list to a numeric using parse combined with map. If the end of the file has been reached, readline returns the integer 0, so we must exit the loop reading and processing all the lines in the file using break.

So the core of the loop looks like this: 

> do
>    alt := readline(fdread);  # read in current line
>    if alt=0 then break fi; # if end of file, exit loop
>    a := map(parse,
>       [seq(substring(alt, i), i=[11 .. 12, 13 .. 14, 15 .. 18,
>        19 .. 21, 22 .. 23, 24 .. 25, 31 .. 35])]);
>       # convert data to a list of numerics
>    [... further conversion statements ...]
> od;

Converting the data

All the resulting numerics must be properly converted to create three-dimensional Cartesian coordinates.

The right ascension data can be converted first to degrees by 

> decimal := (hrs, min, sec) -> hrs + min/60 + sec/3600;
and then to radians through 
> dec * 15 * Pi / 180;
1/12 dec Pi

A full circle consists of either 360 degrees or 24 hours (360/24=15). So the x and y Cartesian coordinates of each galaxy are: 

> x := cos(dec/12*Pi);
x := cos(1/12 dec Pi) 
> y := sin(dec/12*Pi);
y := sin(1/12 dec Pi)

The z ('height') coordinate is given in degrees, so: 

> z := deg/180*Pi;
z := 1/180 deg Pi

In 1921 Edwin Hubble discovered that galaxies are moving away from each other, and also that the more distant the galaxies the faster their motion. He found the following distance equation including the famous Hubble constant and the radial or escape velocity of the objects.

The constant's value is still not exactly known. In 1932 Edwin Hubble estimated it at being 530 km/sec/Mpc, Mpc = Megaparsec, 1 Mpc = 3.0857 *1019 km. (1 Parsec are 3.2616 lightyears or 206,264.8 AU with 1 AU the distance between the Sun and the Earth.) A new estimate for the Hubble constant was published in ASTRONOMY magazine of August '96: According to two recent studies the Hubble constant is either between 68 and 78, or between 53 and 61. Knowing the exact value is important to evaluate whether the universe is expanding forever or after a halt is contracting again. One purpose of NASA's Hubble Space Telescope is to nail down the constant more exactly. 

> hubble_constant := 73*km/sec/Mpc;

> with(student):  # This statement is not necessary in Maple 6 and later

> equ := radialv*km/sec = hubble_constant * distance;

> isolate(equ, distance);

Thus the absolute distance is the equation above; its value is multiplied with the Cartesian x-, y-, and z-coordinates.

Opening a file for input: fopen

If you would like to create a new file, you may open it before with the command fopen. The function takes the name of the file to be opened, and the method of how the file will be opened. The method may either be READ or WRITE, in this case we use the latter specification.

Syntax:

fopen
Calling sequence: fopen(filename, method);
Parameter: filename: the name of the file to be opened
Return: the file descriptor (an integer)

As you see, fopen returns the number of the file it has opened. With all the other file I/O functions available in Maple, you may use this file descriptor instead of the complete filename.

Writing the new data to a text file: fprintf

The fprintf function stores the processed distance and Cartesian coordinates cx, cy, and cz of the galaxies to the file opened before. This function writes data to a file according to a specified format. The formatting specifications must be passed as a string, followed by the names of the variables which values are to be written. As with all the other file I/O functions in Maple, fprintf is very similar to the corrsponding I/0 command available in C, so the formatting string may include one or more %-placeholders. See the online help on fprintf for a list of all formatting specifications accepted by the function.

(Note: You do not need to open the file before using fprintf, if the file fprintf has to write data to has not been opened before, fprintf does this automatically for you.)

Syntax:

fprintf
Calling sequence: fprintf(filename or file descriptor, formatstring, sequence of variables);
Parameter: filename: the name of the file to be read line by line (a string)
Return: the line as a string

Since we would like to store floats to the new file, the formatting string includes %f placeholders. Each galaxy is stored in a separate line, so the formatting string must also includes the newline tag \n for line separation. Note that the name of the file to be created is stored in the variable outputfilename, whereas the file descriptor determined by fopen is assigned to the name fdwrite which is subsequently used by fprintf.

The statments: 

> outputfilename := "c:/rsinf.txt":
> fdwrite := fopen(outputfilename, WRITE); # open rsinfo for input
> do
>    alt := readline(fdread);  # read in current line
>    if alt=0 then break fi; # if end of file, exit loop
>    a := map(parse,
>       [seq(substring(alt, i), i=[11 .. 12, 13 .. 14, 15 .. 18,
>        19 .. 21, 22 .. 23, 24 .. 25, 31 .. 35])]);
>       # convert data to a list of numerics
>    [... further statements to evaluate distance and 
>         Cartesian coordinates cx, cy, cz ...]
     fprintf(fdwrite, "%f %f %f %f\n", distance, cx, cy, cz);
>       # write coordinates to separate lines
>    fi
> od;

Closing files: fclose

Finally, after all the data has been written, the new file has to be closed with fclose, so that it can be opened or read by another application or Maple I/O function.

Syntax:

fclose
Calling sequence: fclose(filename or file descriptor);
Parameter: filename: the name of the file to be closed (a string)
file descriptor: the integer returned by fopen
Return: NULL

The galaxydata procedure

The complete listing of all statements for Release 4 and 5 (Release 3 worksheet is included in the redshift archives available below): 
> restart:

> galaxydata := proc(input::{string, name}, output::{string, name},
>    maxdist::{integer, constant})
>
>    # by Alexander F. Walz, July 1997, May 1998, December 13, 1998,
>    # February 07, 1999
>    local a, d, ra_decimal, decl_decimal, cx, cy, cz, alt, dist,
>       dec, cartx, carty, cartz, counter, fdread, fdwrite, distance;
> 
>    # various functions
>    dist := radialv -> 1/73*radialv;   # distance of galaxies
>    dec := (hrs, min, sec) -> hrs + min/60 + sec/3600;
>      # conversion of sexagesimal values to decimal values
>    cartx := ra -> cos(ra/12*Pi);
>      # conversion of right ascension to cartesian x-coordinate
>    carty := ra -> sin(ra/12*Pi);
>      # conversion of right ascension to cartesian y-coordinate
>    cartz := deg -> sin(deg*Pi/180);
>      # conversion of right ascension to cartesian z-coordinate
> 
>    Digits := 7;  # only in order to reduce size of output file
>    counter := 0;  # number of galaxies within radius Mpc
>    fdread := convert(input, string);
>    fdwrite := fopen(convert(output, string), WRITE);
>      # open output file
>    do
>       alt := readline(fdread);  # read in current line
>       if alt=0 then break fi;
>       a := map(parse,
>          [seq(substring(alt, i), i=[11 .. 12, 13 .. 14, 15 .. 18,
>           19 .. 21, 22 .. 23, 24 .. 25, 31 .. 35])]);
>       d := dist(a[7]);
>       ra_decimal := dec(a[1], a[2], a[3]);
>       decl_decimal := 
>       `if`(a[4] < 0, -dec(-a[4], a[5], a[6]), dec(a[4], a[5], a[6]));
>       cx := evalf(d*cartx(ra_decimal));
>       cy := evalf(d*carty(ra_decimal));
>       cz := evalf(d*cartz(decl_decimal));
>       distance := evalhf(sqrt(cx^2+cy^2+cz^2));
>       if distance < maxdist then
>         counter := counter + 1;
>          fprintf(fdwrite, "%f %f %f %f\n", distance, cx, cy, cz);
>          # write coordinates to separate lines
>       fi
>    od;
>    fclose(fdwrite);  # close output file
>    counter
> end:
galaxydata converts all the original data (first argument) of the galaxies within the given range (third argument) and writes the new data to another file (second argument). It returns the number of galaxies in the new file.

Syntax: 

> galaxydata(
>    pathname of source file, 
>    pathname of target file, 
>    distance to Sun in Mpc);

> galaxydata(`c:/data/redshift.txt`, `c:/data/rsinf.txt`, infinity);
14729

Reading in numerical data: readdata

The readdata is used to read in numerical data, i.e. integers and floats, from a file with a given format. The values must be separated with whitespaces, blanks or tabulators (\t). It returns a list of values if the file consists of only one column, or a list of lists otherwise. By default, readdata returns the data as floats, but you may pass the name integer if there are only integers in the file.

Syntax:

readdata
Calling sequence: readdata(filename or file descriptor,
number of columns, 'integer');
Parameter: filename: the name of the file to be closed (a string)
file descriptor: the integer returned by fopen
n: number of columns, an integer (default is 1)
'integer': (optional) the name
Return: a list or list of lists of numerics

The following statements read in the file created by galaxydata, assigns the numerical data (a list of lists) to the name d, applies the map and the user-defined function f to that list in order to extract only those galaxies that are within a given distance, and finally creates a plot. 

> restart:

> d := readdata(`c:/data/rsinf.txt`, float, 4):

> f := (x, parsec) -> `if`(x[1] < parsec, x[2 .. 4], NULL):

> g := map(f, d, 150):

> plots[pointplot3d](g,
>    symbol=POINT, color=navy, orientation=[0, 0],
>    axes=FRAME, scaling=constrained);
The plot of the galaxies within 150 Mpc of our Galaxy:

The plot computed by Margaret Geller and John Huchra published in ASTRONOMY 08/96, p. 46:

A plot in 3D clearly shows the pencil-beam surveys: 

> str := plots[pointplot3d](g, symbol=POINT, color=navy, 
>    orientation=[-24, 75], axes=BOX, scaling=constrained):

> str;

Creating VRML files

Users of Maple V R5 and later can also generate a VRML file of the above PLOT3D structure with the plottools package in order to examine it in realtime in a VRML viewer, such as Cosmo Player shipped with the Netscape Communicator 4. 

> with(plottools, vrml): 

> currentdir(`c:/data`); # (vrml does not accept subdirectories)

> vrml(str, `redshift.wrl`, emissive_color=COLOR(RGB,1,1,1));

Files for download

The following archives contain all of the above decribed procedures and the extracted galaxy catalogue file for your own exploration with Maple V Release 3 through Maple 7.

In addition you may also download the precomputed rsinf.txt file generated by the galaxydata procedure to that you may not create the file yourself (which takes a lot of time) and also the above mentioned redshift.wrl VRML file containing all galaxies.


back

MAPLE UNIVERSE GALAXIES 2.0.1 current as of December 25, 2001
Author: Alexander F. Walz, alexander.f.walz@t-online.de
Original file location: http://www.math.utsa.edu/mirrors/maple/mpluniv1.htm