[fitsbits] Adaptive mesh refinement

[William Thompson]

Folks:

This is a modified version of an email I sent earlier to the fitswcs mailing 
list.  I thought it probably made sense to open up the discussion to the wider 
FITS community.

I was asked how to store data from model calculations using adaptive mesh
refinement in FITS files.  To my mind, there are two possibilities, pixel lists 
or the -TAB projection.  Personally, I prefer the -TAB projection, for reasons 
listed below, but I think that one sentence in the definition makes it 
impracticle for this application.  Below I present my understanding of each of 
the options.  Please let me know if my understanding is incorrect or incomplete.

Scenario 1:  The -TAB projection.

The advantages of the -TAB projection are that it explicitly ties together the 
data array with the coordinate information, and that the coordinates are stored 
directly without needing any additional computations.

Suppose that we have a model with 1000 points in it, each of which has 
coordinates x(i), y(i), z(i).  My belief is that I could store this using the 
following keywords.  (Some keywords, like CUNITn, have been left out for brevity.)

NAXIS	= 1			/One data axis
NAXIS1	= 1000
WCSAXES	= 3			/Three coordinate axes

CTYPE1  = 'HEQX-TAB'		/Heliocentric Earth Equatorial coordinates
CRPIX1	= 1
CRVAL1	= 1
CDELT1	= 1
CTYPE2  = 'HEQY-TAB'
CRPIX2	= 1
CRVAL2	= 1
CDELT2	= 1
CTYPE3  = 'HEQZ-TAB'
CRPIX3	= 1
CRVAL3	= 1
CDELT3	= 1

PC1_1	= 1			/Use PC matrix to tie coordinates together.
PC1_2	= 0
PC1_3	= 0
PC2_1	= 1
PC2_2	= 1
PC2_3	= 0
PC3_1	= 1
PC3_2	= 0
PC3_3	= 1

PS1_0	= 'WCS-TAB'
PV1_1	= 1
PV1_2	= 1
PS1_1	= 'COORDS'
PV1_3	= 1

PS2_0	= 'WCS-TAB'
PV2_1	= 1
PV2_2	= 1
PS2_1	= 'COORDS'
PV2_3	= 2

PS3_0	= 'WCS-TAB'
PV3_1	= 1
PV3_2	= 1
PS3_1	= 'COORDS'
PV3_3	= 3

The coordinate array in the binary table should have the dimensions 
(3,1000,1,1).  However, the documentation disallows degenerate axes, so we'd 
have to substitute (3,1000,2,2), of which only 25% is actually used.  This is 
highly inefficient, and makes the scheme impracticle.  Is this a correct 
assessment?  Is anything wrong about the above header?


Scenario 2: Pixel lists

The most straightforward way to use pixel lists would seem to be to set 
TCRPXn=0, TCRVLn=0, TCDLTn=1, so that the coordinates can be stored directly in 
the table.  One would then have a binary table header such as the following

NAXIS   =                    2 /Binary table
NAXIS1  =                   16 /Number of bytes per row
NAXIS2  =                 1000 /Number of rows
PCOUNT  =                    0 /Random parameter count
GCOUNT  =                    1 /Group count
TFIELDS =                    4 /Number of columns

TFORM1  = '1E      '           /Real*4 (floating point)
TTYPE1  = 'DENSITY '           /Label for column 1

TFORM2  = '1E      '           /Real*4 (floating point)
TCTYP2  = 'HEQX    '           /
TCRPX2  =                    0 /
TCRVL2  =              0.00000 /
TCDLT2  =              1.00000 /

TFORM3  = '1E      '           /Real*4 (floating point)
TCTYP3  = 'HEQY    '           /
TCRPX3  =                    0 /
TCRVL3  =              0.00000 /
TCDLT3  =              1.00000 /

TFORM4  = '1E      '           /Real*4 (floating point)
TCTYP4  = 'HEQZ    '           /
TCRPX4  =                    0 /
TCRVL4  =              0.00000 /
TCDLT4  =              1.00000 /

Some care would have to be exercised when spherical coordinates are used.  To 
keep the philosophy of storing the coordinates directly, I think that one would 
have to treat the data as being in the plate carree (-CAR) projection.

Thank you,

Bill Thompson



-- 
William Thompson
NASA Goddard Space Flight Center
Code 612.1
Greenbelt, MD  20771
USA

301-286-2040
William.T.Thompson.1 at gsfc.nasa.gov