# Copyright restrictions apply - see stsdas$copyright.stsdas 
# 
include <tbset.h>

include	"elcolnames.h"

define	PI		3.141592654
define  PI2             (PI / 2.)
define	DPI		(2. * PI)
define	EPSILON		1.e-8

define	MAX_NVERTICES	90		# BEWARE ! This number must be equal
					# or less than the maximum number of
					# polygon vertices defined in APPHOT !

# T_ELAPERT -- Generates elliptical apertures for use in APPHOT 
# polygonal aperture photometry. Apertures have the form of x,y pair
# lists. Input is a SDAS table whith isophote parameters.

procedure t_elapert ()

pointer	tp, cptr			# table pointers
char	table[SZ_FNAME]			# table whith isophotal data
char	polygons[SZ_FNAME]		# vertices are stored here
char	coords[SZ_FNAME]		# centers are stored here

int	n, row, maxrow
int	outp, outc			# output file pointers
real	r, fi				# polar coordinates
real	x, y				# image coordinates
real	a0				# semi-major axis length
real	x0, y0				# center of isophote on frame
real	eps0, teta0			# elipticity, pos. angle
real	fistep

pointer	tbtopn()
int	tbpsta(), open(), strlen()

begin
	call clgstr ("table",  table,  SZ_FNAME)
	call clgstr ("polygons", polygons, SZ_FNAME)
	call clgstr ("coords", coords, SZ_FNAME)

	# Open output file(s) and isophote table.
	outp = open (polygons, NEW_FILE, TEXT_FILE)
	if (strlen(coords) > 0)
	    outc = open (coords, NEW_FILE, TEXT_FILE)
	tp = tbtopn (table, READ_ONLY, 0)
	maxrow = tbpsta (tp, TBL_NROWS)

	# Scan table.
	do row = 1, maxrow {

	    # Get row data.
	    call tbcfnd (tp, ES_CA, cptr, 1)
	    call tbegtr (tp, cptr, row, a0)

	    if (a0 > 0.) {
	        call tbcfnd (tp, ES_CEPS, cptr, 1)
	        call tbegtr (tp, cptr, row, eps0)
	        call tbcfnd (tp, ES_CTETA, cptr, 1)
	        call tbegtr (tp, cptr, row, teta0)
	        teta0 = teta0 / 180. * PI + PI2
	        call tbcfnd (tp, ES_CX0, cptr, 1)
	        call tbegtr (tp, cptr, row, x0)
	        call tbcfnd (tp, ES_CY0, cptr, 1)
	        call tbegtr (tp, cptr, row, y0)

	        # Sets fistep in order not to exceed the maximun number
	        # of polygon vertices allowable in APPHOT.
	        fi = 0.
	        r  = a0
	        n  = 0
	        while (fi < DPI) {
	            # Step over ellipse.
	            fi = fi + 1./ r
	            r = a0 * (1. - eps0) / sqrt (((1. - eps0) * cos (fi))**2 + 
		        (sin (fi))**2)
	            n = n + 1
	        }
	        if (n > MAX_NVERTICES)
	            fistep = real (n) / real (MAX_NVERTICES)
	        else
	            fistep = 1.

	        # Generate aperture.
	        fi = 0.
	        r  = a0
	        while (fi < DPI) {
	            # Get image coordinates of (r, fi) pixel
	            x = r * cos (fi + teta0) + x0
	            y = r * sin (fi + teta0) + y0
	            # Print
	            call fprintf (outp, "%g  %g\n")
	                call pargr (x)
	                call pargr (y)
	            # Step over ellipse.
	            fi = fi + fistep/ r
	            r = a0 * (1. - eps0) / sqrt (((1. - eps0) * cos (fi))**2 + 
		        (sin (fi))**2)
	        }
	        # Mark end of aperture.
	        call fprintf (outp, ";\n")
	        # Write center.
	        if (strlen(coords) > 0) {
	            call fprintf (outc, "%g  %g\n;\n")
	                call pargr (x0)
	                call pargr (y0)
	        }
	    }
	}

	call close (outp)
	if (strlen(coords) > 0)
	    call close (outc)
	call tbtclo (tp)
end