Subject: Making a FITS file

From: colbert at astro.umd.edu (Edward J.M. Colbert)
Date: 2 Aug 1993 13:43:47 GMT

Hi:

I would like to take a 3-dimensional float array (x,y,z)_i and make
a FITS image out of it with x and y on the pixel axes.  There must 
be a routine already written to do this.  Where can I get it?

Thanks,
Ed

---

From: pence at tetra.gsfc.nasa.gov (William Pence)
Date: Mon, 2 Aug 93 10:51:43 EDT

>be a routine already written to do this.  Where can I get it?

I'm not sure what you mean by "with x and y on the pixel axes", but
the following sample subroutine will take a 3-d array and write it out
as a 3-D FITS image (data cube) if that is what you wanted to do.
  
First, build a copy of the FITSIO subroutine library, available from the
/pub/fitsio/ directory in the anonymous ftp account on tetra.gsfc.nasa.gov
then link it will the following Fortran subroutine as a template for creating a
FITS file.  You will probably need to make some minor modifications to this
subroutine to meet your exact requirements.:

C-------------------------------------------------------------------------------
	subroutine make3d(array,nx,ny,nz)

C	Write the array of data out to a FITS file.

C	array = input floating point 3D array
C	nx = input size of the first axis of array
C	ny = input size of the second axis of array
C	nz = input size of the third axis of array

	integer nx,ny,nz
	real array(nx,ny,nz)
	integer status,blocksize,iunit,bitpix,naxis,naxes(3)
	integer pcount,gcount,group,fpixel,nelem
	logical simple,extend

C	initialize status, FITS blocksize, and the fortran unit number to use
        status=0
	blocksize=1
        iunit=21

C       Create the new FITS file called file3d.fits
        call ftinit(iunit, 'file3d.fits', blocksize, status)

C	Define datatype of FITS array as single precision floating point.
        bitpix=-32

C       Initialize values for required keywords
        simple=.true.
        naxis=3
        naxes(1)=nx
        naxes(2)=ny
	naxes(3)=nz
        pcount=0
        gcount=1
        extend=.true.

C       Write the required primary array keywords
        call ftphpr(iunit,simple,bitpix,naxis,naxes,pcount,gcount,
     &       extend,status)

C       Define the primary array structure (sets up internal FITSIO parameters)
        call ftpdef(iunit,bitpix,naxis,naxes,pcount,gcount,status)

C       Write the primary array of data (real*4 data array)
        group=1
        fpixel=1
	nelem=nx*ny*nz
        call ftppre(iunit,group,fpixel,nelem,array,status)

C	Close the FITS file
        call ftclos(iunit,status)

	end

C-------------------------------------------------------------------------

It would be useful to add more error checking to the above program
such as checking that the returned status values are = 0.  Also, one
would probably want to add some additional optional header keywords to
the FITS file.  This can be easily be done with other FITSIO subroutine calls.

-Bill Pence