Munipack stable release

Yesterday, I finished work on a simple script to make of a regular release of Munipack. Geeks from Fedora Astronomy group included Munipack in their release. Unfortunately, they used a really historical version (it may come from past century). I think, it is due of may unconventional designing and releasing of Munipack (also Nightview) versions. New versions are generated once per day from cvs repository as nightly builds. Many of peoples are supposed that it is a development version only and it isn't suitable as a production release. It isn't true, but all peoples suppose it due to the naming scheme.

To suppress any confusion, the version issue is improved to a conventional naming scheme from now. The "stable" release will be also generated from cvs repository. The generation will less frequent (only once per month), the archive will named as munipack-0.4.?.tar.gz with incrementing only last digit. Also, the archive will generated if the following rules will be satisfied:

• the last changelog record is older then approximately 35 days
• the changelog records are different

I think, it will prevent incrementing of minor version without change of content and it will give me some time for developing and bug-fixing of a new code (you can still use of nightly builds in meantime).

Astronomical image processing guide (How to list of info keywords?)

A great advantage of the FITS format is support of any arbitrary information included into an image file. A place for the info is reserved in header of the FITS file. Every FITS header must contain a set of records (lines) containing certain mandatory keywords. The set may be followed by records with any other keywords.

The mandatory keywords are: SIMPLE, EXTEND and END. All others are optional. The records of header has the fixed structure:

1-8 - keyword
9 - = equal's sign
10-x - value (various lengths)
x -80 - comment (every text following slash '/')

All records are 80-bytes long and there is no separator between ones. The special keywords COMMENT and HISTORY has no defined a structure and introduce any text string (continuation lines are not allowed). Recent recommendations to contents of the records adds physical units to numerical values as a part of comment. Any set of keywords is not widely used. It means that various utilities may used different keyword for the same entity.

The most common keywords are included in following (artificial) example:

SIMPLE = T / file does conform to FITS standard
BITPIX = 16 / number of bits per data pixel
NAXIS = 2 / number of data axes
NAXIS1 = 382 / length of data axis 1
NAXIS2 = 255 / length of data axis 2
EXTEND = T / FITS dataset may contain extensions
EXPTIME = 60.000 / [s] Exposure time
DATE-OBS= '2008-06-06T01:02:41.390' / UTC of exposure start
FILTER = 'R ' / filter
OBJECT = 'Star' / Object name
COMMENT This file was written by XXX.
END

We have more ways to handle with FITS header keywords. Practically, we will need to list of all keywords (as a variant of more or less command on unix's shell) or list of value of some specified record. Both situations can be easy coded. We introduce of utility FITSless which will print all keywords (full header) without any switches and a specified value when any keyword will presented, its value will be printed.

The values are usually different kinds (from computer point of view). For example, the observer's name will coded as a string, an exposure time will a real number, the date of start will a specially formatted string. To print the information, we will use only string representation, but in real code it will be better use of appropriate data type. Unfortunately, there is no simple way how to do it in Fortran, C and cFITSIO under recent versions of ones.

program FITSless

implicit none

integer :: status
! status ... FITS status (0=no error)

integer :: j,blocksize

character(len=666) :: name = 'image.fits'
! name .. fill with name of the image to open

character(len=666) :: keyword
! keyword to print the record

character(len=666) :: record
! header's record

character(len=666) :: value, comment
! record's value & comment

integer :: nhead, hpos
! number of records in header and current position

! get first command line parameter
call get_command_argument(1, keyword)

status = 0
call ftopen(25,name,0,blocksize,status)
if( status /= 0 ) stop 'File not found.'

! list specified keyword or full header
if( keyword /= ' ' ) then

write(*,*) 'List of a record specified by your keyword:'
write(*,*) '-------------------------------------------'

call ftgkys(25,keyword,value, comment, status)
! checking status of the operation
if( status == 0 ) then
write(*,*) trim(keyword),' = ',trim(value),' /',trim(comment)
else
write(*,*) 'Keyword "',trim(keyword),'" not found.'
end if

else

write(*,*) 'Full header list:'
write(*,*) '-----------------'

call ftghps(25,nhead,hpos,status)
do j = 1, nhead
call ftgrec(25,j,record,status)
write(*,*) trim(record)
enddo

end if

call ftclos(25,status)

end program FITSless

The code is self-explanatory, I recommend try this following experiments:

• try give as parameter an arbitrary text string or upper/lower case keywords
• process any printed information by some another text tool (sed, grep,..)
• play with listing of a record on given position in the header
• modify code to read an input file name by command-line parameters

A ridge of algorithms

I spend some time with the growth-curve method during last days. Its looks nice, but relative hard to control.

Some side product of the game is comparison of the original Stetson's algorithm for estimate of the arithmetical mean with new ones. To get mean, Stetson have introduced some rejection-type algorithm. I've no more detailed information about it - there is no reference in source code or a paper. Opposite with this, Munipack estimate of the mean by a method on base of minimisation of a general function. It minimises function in shape of the least-square near of minimum and suppress it to zero otherwise by method of the maximum like-hood. The chapter describing of a robust statistical methods in Numerical recipes is an ideal start point to get detailed view to this field. A source code of Munipack uses of those ideas, but it doesn't describe it. The new estimator has been included to Munipack approximately eight years ago.

I compared of the estimators on case of the ridge of magnitude on differences of following apertures. Both three-dimensional graphs shows number of the differences in cells specified by apertures and its value's ranges. The horizontal (x) axis takes magnitude differences, the horizontal (y) axis takes size of aperture and vertical (z) shows the number of values in that bin interval. The top side of the data cube shows projection of the number of values to a plane. Colors of the histogram maps number of points: dark - near zero, light - near maximum (hundred).

Muniphot's ridge

Stetson's ridge

Both images looks similar. Only one difference is a small systematic shift of the ridge to negative region on Stetson's graph. It is about 0.01 in magnitude, so I think, it is not important in real situations. Perhaps. It may means that the new method for the mean estimations works little bit better.

Ridges height strongly depends on aperture. Small values gives perfect estimation against apertures greater that fifteen (in my case) when the the histogram shows no peaks. The characteristics may be important for growth-curve fitting.

As test data for both graph, I used combined image of M67 as in previous post. However, graphs of ridges looks differently (gnuplot's pm3d feature has been used on graphs here). It is due of my mistake during data processing. The previous perfect ridge with strong peaks at large apertures comes from many zero's differences of bad magnitudes of 99.999. Keep smiling.

Astronomical image processing guide (How to save of an image to a FITS file?)

One of most common operation done on FITS files is creation of a new image. The creation process is straightforward. One is required to have an image (result of a mathematical algorithm) with known dimensions and that's all.

The FITS format supports up to seven-dimensional images without any strict limitation of its axis ranges. Practically, we've some limits due to physical properties of present computers. The most important is a size of the output image. The size is determined (approximately) as a product of image dimensions multiplied by number of bytes occupied by one pixel. For two-dimensional image of 100x100 pixels in both axis, represented by real numbers (4-byte, BITPIX=-32), we got size about 40 kilo-bytes. As we know from previous lecture, the images produced by an algorithm are usually saved as real data to preserve its numerical precision.

As example of a test image, I choose Bessel's function of zero kind J₀ which represents light's diffraction on cylindrical aperture. We can observe of square of Bessel's function in an ideal telescope as the image of a star (of course, we use only J₀ for simplicity, a star will looks differently). J₀ is non-standard Fortran function and may be not supported by all compilers (gfortran does it) so one is optional only and you can play with cosine under a strict-standard compiler.


An algorithm to save an image to FITS file is straightforward:

1. create image as an array
2. fill the array by an image
3. initiate a FITS
4. setup image parameters
5. save the data

program FITSsave

implicit none

integer :: status, bitpix, naxis, naxes(2)
! status ... FITS status (0=no error)
! naxis .. number of axes in the image (we set =2)
! naxes .. dimensions of the image (2-element array)

real, dimension(:,:), allocatable :: d
! data matrix

character(len=666) :: name = 'image.fits'
! name .. fill with name of the image to create

! aux
real :: x,y,r
integer :: i,j

! set dimensions of the new image
naxis = 2
naxes = (/ 100,100 /)

! set bipix of the image (try: 8,16,32 and -32)
bitpix = -32

! create a data storage (allocate memory) for the image
allocate(d(naxes(1),naxes(2)))

! fill image with values
do i = 1, naxes(1)
do j = 1, naxes(2)

! rectangular coordinates
! the left bottom pixel has 1,1
x = i
y = j

! distance from origin
r = sqrt(x**2 + y**2)

! set value
d(i,j) = cos(r/5.0)
!d(i,j) = besj0(r/5.0)
! uncomment for J0 (Bessel function of zero kind)
! J0 represents diffraction on cylindrical aperture
end do
end do

! save the data
status = 0
call ftinit(26,name,1,status)
call ftphps(26,bitpix,naxis,naxes,status)
call ftp2de(26,1,naxes(1),naxes(1),naxes(2),d,status)
call ftclos(26,status)

! free allocated memory
deallocate(d)

end program FITSsave

The code can be compiled and run by

host$ gfortran -Wall -o FITSsave FITSsave.f90 -L/usr/local/lib -lcfitsio
host$ ./FITSsave

the output file is named as image.fits and can be viewed by any FITS viewer, for example by ds9:

host$ ds9 image.fits


Notes.

Any handle with numerical operations in many computer languages may be little bit confusing. Fortran strictly distinguish between integer and real numbers. The notation 3/4 (both integers) products result 0 (reminder is forget), but 4/3 gives 1. Opposite with this, the notation 3.0/4.0 gives 0.75 (two significant places).

The function ftinit (the initial function for FITS file) can create only a new file. There is no way how to replace any existing file. This is simply a feature of cfitsio, not a bug. It means that you must remove the older file image.fits before run FITSsave again.

Astronomical image processing guide (How to list of values of a FITS image?)

Every FITS file is representation of an image usually created by an optical device. The image is quantized (sampled) to elementary cells called pixels. An information knows for any pixel are: a pixel coordinate (Cartesian x,y) and its captured optical flux (CCD's flux is a linear function of captured photons).

The pixels which represents an image, are rearranged and a captured flux is digitalised to a matrix. The matrix is saved to a FITS file by a defined algorithm. The pixels coordinates (integers) are arranged in that matrix by the way:

[M,1] [M,2] .. [M,N]
.. ..
[2,1] [2,2] .. [2,N]
[1,1] [1,2] .. [1,N]

The matrix represents an image of width of N pixels and M pixels of height. The origin and orientation is in usual mathematical fashion.

Every pixel is represented by a number. The kind of the number may be an integer and a real (real numbers are with fractional part). Raw images (pure product of a device) are usually represented by integer numbers in interval from 0 to 65535 (2^16) for CCD and from 0 to 4096 (2^12) for a digital camera. A processed images as result of mathematical operations are saved as a real numbers with floating point. (Arithmetical operations may reduce of its precision). The method (complicated on first sight) to store of data reflects many of astronomer's needs and save your disk space. The data representation is coded in parameter BITPIX in the fits header by the way (not all possibilities are included):

BITPIX bytes type range of values
8 1 integer 0 .. 255
16 2 integer 0 .. 65535
32 4 integer 0 .. 4294836225
-32 4 real -1e38 .. 1e38 (7 digits)

An operation on a image included in a FITS file is relative easy. Follow the instruction:

1. open of FITS
2. get of its size (width, height)
3. allocate memory for a matrix
4. read data
5. play with data

Fortran offers a very effective way for manipulation with matrixes. For a matrix D, we can select an i,j-element as D(i,j), a i-row D(i,:),i-column D(:,j) or submatrix D(1:10,50:60).

program FITSlist

implicit none

integer :: status, bitpix, naxis, naxes(2)
! status ... FITS status (0=no error)
! naxis .. number of axes in image (we require =2)
! naxes .. dimension of the image (2-element array)

integer :: i,j,blocksize,pcount,gcount,minvalue
logical :: extend, simple,anyf
! required by cFITSIO

real, dimension(:,:), allocatable :: d
! data matrix

character(len=666) :: name = 'image.fits'
! name .. fill with name of the image to open

status = 0
call ftopen(25,name,0,blocksize,status)
call ftghpr(25,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
allocate(d(naxes(1),naxes(2)))
call ftg2de(25,1,minvalue,naxes(1),naxes(1),naxes(2),d,anyf,status)
call ftclos(25,status)

! print value of a random pixel
write(*,*) '# d(1,1)=',d(1,1)

! print last tree values of first row
write(*,*) '# d(1,-10:)=',d(1,size(d,2)-3:)

! print of a submatrix with indexes
do i = 1,10
do j = 50,60
write(*,*) i,j,d(i,j)
end do
end do

end program FITSlist

The output of the code can be saved to a file by sequence of commands:

host$ gfortran -Wall -o FITSlist FITSlist.f90 -L/usr/local/lib -lcfitsio
host$ ./FITSlist > pixels

and easy plotted in a gnuplot with the command:

gnuplot> splot 'pixels'

How aureole shades stars

The photometry (magnitudes of stars) produced by Munipack is an aperture photometry, which means that a magnitude is determined only as a sum of a flux inside an artificial aperture. The aperture photometry is sufficient for a relative comparison of stars. It should be very useful in many applications like observing of variable objects (variable stars, asteroids, afterglows,..), but there are many of additional requirements (calibration of photometric system, precise photometry,...) in which the aperture photometry start to fail.

Briefly, an image of a star is theoretically a point, but Earth's atmosphere and telescope's optics spread the point to a blob with Gaussian central part and an outer parts showing an "aureole". The problem of the aperture photometry is our ignorance of outer radii which may include an important fraction of a total flux. If we will ignore the aureole's flux, any magnitudes will depends on atmospheric conditions (seeing spreads stars with dependence on its air mass or a night) and on magnitude of the magnitude (A golden rule of differential photometry is to compare objects with similar magnitude). It means that precise photometry requires more sophisticated approach.

Described problems bothers CCD's photometrics for a long time. They are described in articles of classics at 80'. Howell(1989) points to properties of the aperture photometry for differently bright stars. He shows that total flux of bright stars is determined relative precise, but any small error in determination of background for any faint star strongly affects total fluxes. The solution of the problem did published Stetson(1990) in that great article. He developed a method on base of cumulative distribution of magnitudes of stars (in astronomical dialect the grow-curve method). Please read the paper, I have no time to transcribe it here. The most important ones (from my point of view) are: The grow-curve method is complement to PSF method (PSF <=> precise relative magnitude, grow <=> precise calibrated magnitude). The model of cumulative profile of the star as an integral of PSF can be described by a sum of analytic functions (Gaussian, exponential and Moffat - Moffat is generalisation of Lorentz).

The grow-curve method for a large radii can be described by Moffat's function ~1/(1+r^2)^A (A=2 for Lorentz). We suppose that all stars on the image has the same cumulative (grow-curve) profile (integral of PSF modeled by Moffat). To construct of the mean grow-curve for all stars with no matter to its brightness, Stetson uses ratio of fluxes in et sequentia of apertures, eg. the following differences of aperture's magnitudes.

M 67

I try it on image of the M 67 open stellar cluster exposed on 2008-02-26 at R filter via 0.4m telescope at HaP MK Brno. The standard Munipack run has been done. About 250 stars has been found. The results included in graph in fashion of the Stetson's article are included.



I take out of these basic fact from the graphs:

• The scatter of the points is relative great over five magnitudes (not all are included). It is due to random errors, disturbance of magnitudes of overlap stars, wrong determination of background for faint stars.
• A pretty description of magnitudes displays frequency of magnitudes in an interval of difference magnitude and apertures. There is a spine visualising of most frequent values of the variables. I'd pleasure from behaviour of the spine at large radii where the values converges to zero. It means that the robust mean algorithm used in Munipack to determine of sky level is really good.
• The graph show differences about 0.1 magnitudes at radii usually applied to photometry. The function is not a profile itself. It means that the difference between a real and determined magnitude will a few tenth of magnitude. In other word, the calibration of Munipack's aperture magnitudes will depends on actual observation conditions and the calibrations derived from the magnitudes will contain an systematic error which will be depend on magnitude of a star.