Repatch (from linux) of OSX IRAF

1 files changed, 548 insertions, 0 deletions

diff --git a/pkg/plot/t_pradprof.x b/pkg/plot/t_pradprof.x
new file mode 100644
index 00000000..dbf7aae1
--- /dev/null
+++ b/pkg/plot/t_pradprof.x
@@ -0,0 +1,548 @@
+include <imhdr.h>
+include <gset.h>
+include	<math.h>
+
+# T_ PRADPROF -- Compute a radial profile using user specified coordinates
+# and plot or list the result. 
+
+procedure t_pradprof()
+
+int	images 			# the list of images
+real	xinit, yinit		# the initial guess for the profile center
+real	pradius			# the plotting radius
+real	paz1, paz2		# azimuth limits
+bool	center			# center the object before computing profile
+int	cboxsize		# the centering box width
+bool	list			# list output instead of plot output
+
+int	rboxsize, npts
+pointer	sp, imname, im, radius, azimuth, intensity
+real	xcntr, ycntr
+
+bool	clgetb()
+int	imtopenp(), imtgetim(), clgeti(), rp_radius()
+pointer	immap()
+real	clgetr()
+
+begin
+	# Allocate stack space.
+	call smark (sp)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+
+	# Get the radial profiling parameters. The width of the extraction
+	# box mut be odd.
+	images = imtopenp ("input")
+	xinit = clgetr ("xinit")
+	yinit = clgetr ("yinit")
+	pradius = clgetr ("radius")
+	paz1 = clgetr ("az1")
+	paz2 = clgetr ("az2")
+	rboxsize = 2 * (nint (pradius + 1.0)) + 1
+
+	# Get the centering parameters. The centering box must be odd.
+	center = clgetb ("center")
+	if (center) {
+	   cboxsize = clgeti ("cboxsize")
+	   if (mod (cboxsize, 2) == 0)
+	        cboxsize = cboxsize + 1
+	}
+
+	# List the radial profile instead of plotting it?
+ 	list = clgetb ("list")
+
+	# Allocate memory for vectors. 
+	call malloc (radius, rboxsize * rboxsize, TY_REAL)
+	call malloc (azimuth, rboxsize * rboxsize, TY_REAL)
+	call malloc (intensity, rboxsize * rboxsize, TY_REAL)
+
+	# Loop over all images
+	while (imtgetim (images, Memc[imname], SZ_FNAME) != EOF) {
+
+	    # Open the image.
+	    iferr (im = immap (Memc[imname], READ_ONLY, 0)) {
+		call eprintf ("Image %s not found\n")
+		call pargstr (Memc[imname])
+		next
+	    }
+
+	    # Find the star center, if center=yes.
+            if (center)
+	       call rp_cntr (im, xinit, yinit, cboxsize, xcntr, ycntr)
+	    else {
+               xcntr = xinit
+	       ycntr = yinit
+	    }
+
+	    # Get the radius and intensity vectors.
+	    npts = rp_radius (im, xcntr, ycntr, rboxsize, pradius, paz1, paz2,
+	        Memr[radius], Memr[azimuth], Memr[intensity])
+
+	    # Make list of the radial profile if list=yes or plot if list=no.
+   	    if (list) 
+		call rp_rlist (Memc[imname], xcntr, ycntr, paz1, paz2,
+		    Memr[radius], Memr[azimuth], Memr[intensity], npts)
+	    else 
+		call rp_rplot (Memc[imname], xcntr, ycntr, paz1, paz2,
+		    Memr[radius], Memr[azimuth], Memr[intensity], npts)
+
+	    call imunmap (im)
+	}
+
+	call mfree (radius, TY_REAL)
+	call mfree (intensity, TY_REAL)
+	call imtclose (images)
+	call sfree (sp)
+end
+
+
+# RP_CNTR -- Compute the star center using a simple 1D centroiding algorithm
+# on the x and y marginals, after thresholding at the mean. 
+
+procedure rp_cntr (im, xstart, ystart, boxsize, xcntr, ycntr)
+
+pointer	im			# pointer to the input image
+real	xstart, ystart		# starting coordinates
+int	boxsize			# width of the centering box
+real	xcntr, ycntr		# centered coordinates
+
+int	half_box, x1, x2, y1, y2
+int	ncols, nrows, nx, ny, try
+pointer	bufptr, sp, x_vect, y_vect
+real	xinit, yinit
+pointer	imgs2r()
+
+begin
+	# Initialize.
+	half_box = (boxsize - 1) / 2
+	xinit = xstart
+	yinit = ystart
+	ncols = IM_LEN (im, 1)
+	nrows = IM_LEN (im, 2)
+
+	try = 0
+	repeat {
+
+	    # Compute the extraction region.
+	    x1 = max (xinit - half_box, 1.0) + 0.5
+	    x2 = min (xinit + half_box, real (ncols)) + 0.5
+	    y1 = max (yinit - half_box, 1.0) + 0.5
+	    y2 = min (yinit + half_box, real (nrows)) + 0.5
+	    nx = x2 - x1 + 1
+	    ny = y2 - y1 + 1
+
+	    # Get the data.
+	    bufptr = imgs2r (im, x1, x2, y1, y2)
+
+	    call smark (sp)
+	    call salloc (x_vect, nx, TY_REAL)
+	    call salloc (y_vect, ny, TY_REAL)
+
+	    # Compute the marginals.
+	    call aclrr (Memr[x_vect], nx)
+	    call aclrr (Memr[y_vect], ny)
+	    call rp_rowsum (Memr[bufptr], Memr[x_vect], nx, ny)
+	    call rp_colsum (Memr[bufptr], Memr[y_vect], nx, ny)
+
+	    # Compute the centers.
+	    call rp_getcenter (Memr[x_vect], nx, xcntr)
+	    call rp_getcenter (Memr[y_vect], ny, ycntr)
+
+	    # Add in offsets to image coordinate system.
+	    xcntr = xcntr + x1
+	    ycntr = ycntr + y1
+
+	    call sfree (sp)
+
+	    # If the shifts are greater than a pixel to 1 more iteration.
+	    try = try + 1
+	    if (try == 1) {
+		if ((abs (xcntr-xinit) > 1.0) || (abs (ycntr-yinit) > 1.0)) {
+		    xinit = xcntr
+		    yinit = ycntr
+		}
+	    } else
+		break
+	}
+end
+
+
+# RP_RADIUS  -- Get the data and compute the radius and intensity vectors.
+
+int procedure rp_radius (im, xcntr, ycntr, rboxsize, pradius, paz1, paz2,
+	radius, azimuth, intensity)
+
+pointer	im			# pointer to the input image
+real	xcntr, ycntr		# the center of the extraction box
+int	rboxsize		# the width of the extraction box
+real	pradius			# the plotting radius
+real	paz1, paz2		# the azimuth limits
+real	radius[ARB]		# the output radius vector
+real	azimuth[ARB]		# the output azimuth vector
+real	intensity[ARB]		# the output intensity vector
+
+int	half_box, ncols, nrows, x1, x2, y1, y2, nx, ny, npts
+pointer	bufptr	
+real	xinit, yinit
+int	rp_vectors()
+pointer	imgs2r()
+
+begin
+	# Initialize.
+	half_box = (rboxsize - 1) / 2
+	xinit = xcntr
+	yinit = ycntr
+	ncols = IM_LEN(im,1)
+	nrows = IM_LEN(im,2)
+
+	# Get the data.
+	x1 = max (xinit - half_box, 1.0) + 0.5
+	x2 = min (xinit + half_box, real (ncols)) + 0.5
+	y1 = max (yinit - half_box, 1.0) + 0.5
+	y2 = min (yinit + half_box, real (nrows)) + 0.5
+	nx = x2 - x1 + 1
+	ny = y2 - y1 + 1
+	bufptr = imgs2r (im, x1, x2, y1, y2)
+
+	# Compute the radius and intensity vectors.
+	npts = rp_vectors (Memr[bufptr], nx, ny, x1, y1, xcntr, ycntr, 
+	    pradius, paz1, paz2, radius, azimuth, intensity) 
+
+	return (npts)
+end
+
+
+# RP_RLIST -- Print the intensity as a function of radial distance on the
+# standard output.
+
+procedure rp_rlist (imname,  xcntr, ycntr, paz1, paz2, radius, azimuth,
+	intensity, npts)
+
+char	imname[ARB]		# the name of the input image
+real	xcntr, ycntr		# the center of the radial profile
+real	paz1, paz2		# the azimuth limits
+real	radius[npts]		# the radius vector
+real	azimuth[npts]		# the azimuth vector
+real	intensity[npts]		# the intensity vector
+int	npts			# the number of points
+
+int	i
+
+begin
+	call printf ("# [%s]   xcntr:%7.2f   ycntr:%7.2f\n")
+	    call pargstr (imname)
+	    call pargr (xcntr)
+	    call pargr (ycntr)
+	call printf ("# az1:%7.2f   az2:%7.2f\n")
+	    call pargr (min(paz1,paz2))
+	    call pargr (max(paz1,paz2))
+
+	do i = 1, npts {
+	    call printf ("%7.2f   %g\n")
+		call pargr (radius[i])
+		call pargr (intensity[i])
+	}
+end
+
+
+# RP_RPLOT -- Plot intensity as a function of radial distance.
+
+procedure rp_rplot (imname, xcntr, ycntr, az1, az2,
+	radius, azimuth, intensity, npts)
+
+char	imname[ARB]
+int	npts
+real	xcntr, ycntr, az1, az2
+real	radius[npts], azimuth[npts], intensity[npts]
+
+char	device[SZ_LINE]
+int	mode
+pointer	gp	
+
+bool	clgetb()
+pointer	gopen()
+
+begin
+	call clgstr ("graphics", device, SZ_LINE)
+
+	if (clgetb("append"))
+	   mode = APPEND
+	else
+	   mode = NEW_FILE
+ 
+	gp = gopen (device, mode, STDGRAPH)
+	call rp_graph (gp, imname, xcntr, ycntr, az1, az2,
+	    mode, radius, intensity, npts) 
+	call gclose (gp)
+end
+
+
+# RP_VECTORS -- Compute the radius and intensity vectors.
+
+int procedure rp_vectors (a, nx, ny, x1, y1, xcntr, ycntr, pradius, paz1, paz2,
+	radius, azimuth, intensity)
+								
+real	a[nx,ny]		# the input data array
+int	nx, ny			# dimensions of the input array
+int	x1, y1			# lower left corner of input array
+real	xcntr, ycntr		# coordinates of center pixel
+real	pradius			# the plotting radius
+real	paz1, paz2		# the azimuth limits
+real	radius[ARB]		# the output radius vector
+real	azimuth[ARB]		# the output azimuth vector
+real	intensity[ARB]		# the output intensity vector
+
+int	i, j, npts
+real	dx, dy, az, pr2, az1, az2, r2, dy2
+
+begin
+	az1 = DEGTORAD (min (paz1, paz2))
+	az2 = DEGTORAD (max (paz1, paz2))
+	while (az1 < 0.) {
+	    az1 = az1 + TWOPI
+	    az2 = az2 + TWOPI
+	}
+	while (az1 > TWOPI) {
+	    az1 = az1 - TWOPI
+	    az2 = az2 - TWOPI
+	}
+	pr2 = pradius * pradius
+
+	npts = 0
+	do i = 1, ny {
+	    dy = (ycntr - y1 + 1 - i)
+	    dy2 = dy ** 2 
+	    do j = 1, nx {
+		dx = (xcntr - x1 + 1 - j)
+		r2 = dx ** 2 + dy2 
+		if (r2 > pr2)
+		    next
+		az = atan2 (dy, dx)
+		if (az < 0.)
+		    az = az + TWOPI
+		if (az < az1 || az > az2)
+		    next
+		npts = npts + 1
+		radius[npts] = sqrt (r2)
+		azimuth[npts] = RADTODEG (az)
+		intensity[npts] = a[j,i]
+	    }
+        }
+
+	return (npts)
+end
+
+
+# RP_ROWSUM -- Sum all the rows in a raster.
+
+procedure rp_rowsum (v, row, nx, ny)
+
+real	v[nx,ny]		# the input subraster
+real	row[ARB]		# the output summed row
+int	nx, ny			# the dimensions of the input subraster
+
+int	i, j
+
+begin
+	do i = 1, ny
+	    do j = 1, nx
+		row[j] = row[j] + v[j,i]
+end
+
+
+# RP_COLSUM -- Sum all the columns in a raster.
+
+procedure rp_colsum (v, col, nx, ny)
+
+real	v[nx,ny]		# the input subraster
+real	col[ARB]		# the output summed column
+int	nx, ny			# the dimensions of the input subraster
+
+int	i, j
+
+begin
+	do i = 1, ny
+	    do j = 1, nx
+		col[j] = col[j] + v[i,j]
+end
+
+
+# RP_GETCENTER -- Compute the centroid of an array.
+
+procedure rp_getcenter (v, nv, vc)
+
+real	v[ARB]			# the input array
+int	nv			# length of the input array
+real	vc			# the output centroid
+
+int	i
+real	sum1, sum2, sigma, cont
+
+begin
+	sum1 = 0.0
+	sum2 = 0.0
+
+	call aavgr (v, nv, cont, sigma)
+	do i = 1, nv
+	    if (v[i] > cont) {
+	        sum1 = sum1 + (i-1) * (v[i] - cont)
+	        sum2 = sum2 + (v[i] - cont)
+	    }
+
+	vc = sum1 / sum2
+end
+
+
+define	MTYPES	"|point|box|plus|cross|circle|hebar|vebar|hline|vline|diamond|"
+define	RP_GBUF	0.10	
+define	RP_SZTITLE	512
+define	DEF_IMTITLE	"imtitle"
+
+# RP_GRAPH -- Graph the radial profile.
+
+procedure rp_graph (gp, imname, xcntr, ycntr, az1, az2, mode, x, y, npts)
+
+pointer	gp		# GIO pointer
+char	imname[ARB]	# image name
+real	xcntr		# starting x coordinate 
+real	ycntr		# starting y coordinate 
+real	az1, az2	# azimuth limits
+int	mode		# Mode
+real	x[npts]		# X data
+real	y[npts]		# Y data
+int	npts		# Number of points
+
+int	i, marks[10], linepattern, patterns[4], clgeti(), btoi(), strdic()
+pointer	sp, marker, title, xlabel, ylabel
+real	x1, x2, y1, y2, wx1, wx2, wy1, wy2, vx1, vx2, vy1,vy2, temp, 
+	szmarker, clgetr()
+bool	clgetb(), streq()
+
+data	patterns/GL_SOLID, GL_DASHED, GL_DOTTED, GL_DOTDASH/
+data	marks/GM_POINT, GM_BOX, GM_PLUS, GM_CROSS, GM_CIRCLE, GM_HEBAR,
+	GM_VEBAR, GM_HLINE, GM_VLINE, GM_DIAMOND/
+
+begin
+	call smark (sp)
+	call salloc (marker, SZ_LINE, TY_CHAR)
+
+	# If a new graph setup all the axes and labeling options and then
+	# make the graph.
+
+	if (mode == NEW_FILE) {
+
+	    call gclear (gp)
+
+	    linepattern = 0
+
+	    x1 = clgetr ("wx1")
+	    x2 = clgetr ("wx2")
+	    y1 = clgetr ("wy1")
+	    y2 = clgetr ("wy2")
+
+	    if (IS_INDEF (x1) || IS_INDEF (x2))
+	        call gascale (gp, x, npts, 1)
+	    if (IS_INDEF (y1) || IS_INDEF (y2))
+	        call gascale (gp, y, npts, 2)
+
+	    call gswind (gp, x1, x2, y1, y2)
+	    call ggwind (gp, wx1, wx2, wy1, wy2)
+
+	    temp = wx2 - wx1
+	    if (IS_INDEF (x1))
+	        wx1 = wx1 - RP_GBUF * temp
+	    if (IS_INDEF (x2))
+	        wx2 = wx2 + RP_GBUF * temp
+
+	    temp = wy2 - wy1
+	    if (IS_INDEF (y1))
+	        wy1 = wy1 - RP_GBUF * temp
+	    if (IS_INDEF (y2))
+	        wy2 = wy2 + RP_GBUF * temp
+
+	    call gswind (gp, wx1, wx2, wy1, wy2)
+	    call gsetr (gp, G_ASPECT, 0.)
+	    call gseti (gp, G_ROUND, btoi (clgetb ("round")))
+
+	    if (clgetb("fill"))
+		call gsetr (gp, G_ASPECT, 0.0)
+	    else
+		call gsetr (gp, G_ASPECT, 1.0)
+
+	    i = GW_LINEAR
+	    if (clgetb ("logx"))
+		i = GW_LOG
+	    call gseti (gp, G_XTRAN, i)
+	    i = GW_LINEAR
+	    if (clgetb ("logy"))
+		i = GW_LOG
+	    call gseti (gp, G_YTRAN, i)
+
+	    # Set the view port
+	    vx1 = clgetr ("vx1")
+	    vx2 = clgetr ("vx2")
+	    vy1 = clgetr ("vy1")
+	    vy2 = clgetr ("vy2")
+	    call gsview (gp, vx1, vx2, vy1, vy2)
+
+	    if (clgetb ("box")) {
+
+	        # Get number of major and minor tick marks.
+	        call gseti (gp, G_XNMAJOR, clgeti ("majrx"))
+	        call gseti (gp, G_XNMINOR, clgeti ("minrx"))
+	        call gseti (gp, G_YNMAJOR, clgeti ("majry"))
+	        call gseti (gp, G_YNMINOR, clgeti ("minry"))
+
+	        # Label tick marks on axes?
+	        call gseti (gp, G_LABELTICKS,
+		    btoi (clgetb ("ticklabels")))
+
+	        # Fetch labels and plot title string. 
+		call salloc (title, RP_SZTITLE, TY_CHAR)
+		call salloc (xlabel, SZ_LINE, TY_CHAR)
+		call salloc (ylabel, SZ_LINE, TY_CHAR)
+
+		# Build system info string
+		call sysid (Memc[title], SZ_LINE)
+		call strcat ("\n", Memc[title], RP_SZTITLE)
+
+		# Build the title string
+		call clgstr ("title", Memc[marker], SZ_LINE)
+		if (streq (Memc[marker], DEF_IMTITLE)) {
+		   call sprintf (Memc[marker], SZ_LINE, 
+		       "Radial Plot of %s at [%0.2f,%0.2f] az=[%.1f,%.1f]")
+		       call pargstr (imname)
+		       call pargr (xcntr)
+		       call pargr (ycntr)
+		       call pargr (az1)
+		       call pargr (az2)
+		}
+		call strcat (Memc[marker], Memc[title], RP_SZTITLE)
+
+		call clgstr ("xlabel", Memc[xlabel], SZ_LINE)
+		call clgstr ("ylabel", Memc[ylabel], SZ_LINE)
+
+	        call glabax (gp, Memc[title], Memc[xlabel], Memc[ylabel])
+	    }
+	}
+
+	# Draw the data.
+	if (clgetb ("pointmode")) {
+	    call clgstr ("marker", Memc[marker], SZ_LINE)
+	    i = strdic (Memc[marker], Memc[marker], SZ_LINE, MTYPES)
+	    if (i == 0)
+		i = 2
+	    if (marks[i] == GM_POINT)
+		szmarker = 0.0
+	    else
+		szmarker = clgetr ("szmarker")
+	    call gpmark (gp, x, y, npts, marks[i], szmarker, szmarker)
+	}
+	else {
+	    linepattern = min (4, linepattern + 1)
+	    call gseti (gp, G_PLTYPE, patterns[linepattern])
+	    call gpline (gp, x, y, npts)
+	}
+	call gflush (gp)
+
+	call sfree (sp)
+end