Repatch (from linux) of OSX IRAF

1 files changed, 563 insertions, 0 deletions

diff --git a/noao/imred/vtel/rmap.x b/noao/imred/vtel/rmap.x
new file mode 100644
index 00000000..313b03db
--- /dev/null
+++ b/noao/imred/vtel/rmap.x
@@ -0,0 +1,563 @@
+include <mach.h>
+include	<imhdr.h>
+include	"vt.h"
+include "numeric.h"
+
+define	LEN_HISTO	1025
+define	SPACING		1
+
+# RMAP -- Project a full disk solar image [2048x2048] into a square
+# image [180x180] such that lines of latitude and longitude are
+# perpendicular straight lines.
+
+procedure t_rmap()
+
+char	inputimage[SZ_FNAME]	# input image
+char	outputimage[SZ_FNAME]	# output data image
+char	outweight[SZ_FNAME]	# output weight image
+char	outabs[SZ_FNAME]	# output absolute value image
+char	histoname[SZ_FNAME]	# output histogram name
+
+real	bzero			# latitude of sub-earth point
+real	el[LEN_ELSTRUCT]	# ellipse parameters data structure
+pointer	inputim, outputim, outw, outa, sp
+pointer inim_subras_ptr, outim_subras_ptr
+pointer	outwei_subras_ptr, outav_subras_ptr
+int	outputrow, wvlngth
+int	inim_subras_bottom
+double	meanf, meanaf, zcm, muzero
+int	numpix
+bool	skip, helium
+real	tempr
+
+real	imgetr()
+int	imgeti()
+pointer	immap()
+pointer imgs2s(), imps2s(), imps2i()
+double	rmap_mode()
+errchk	immap, imgs2s, imps2s, imps2i, checkimmem, rowmap
+
+begin
+	# Get parameters from the cl.
+
+	# Image names.
+	call clgstr ("inputimage", inputimage, SZ_FNAME)
+	call clgstr ("outputimage", outputimage, SZ_FNAME)
+	call clgstr ("outweight", outweight, SZ_FNAME)
+	call clgstr ("outabs", outabs, SZ_FNAME)
+	call clgstr ("histoname", histoname, SZ_FNAME)
+
+	# Open images.
+	inputim = immap (inputimage, READ_ONLY, 0)
+	wvlngth = imgeti (inputim, "wv_lngth")
+	helium = false
+	if (wvlngth == 10830)
+	    helium = true
+	outputim = immap (outputimage, NEW_COPY, inputim)
+	outw = immap (outweight, NEW_COPY, inputim)
+	if (!helium)
+	    outa = immap (outabs, NEW_COPY, inputim)
+
+	# Compute mode estimate from the input image.
+	muzero = rmap_mode (inputim, histoname, helium)
+
+	# Define some parameters for output images.
+	IM_LEN(outputim, 1) = DIM_SQUAREIM
+	IM_LEN(outputim, 2) = DIM_SQUAREIM
+	IM_PIXTYPE(outputim) = TY_INT
+
+	IM_LEN(outw, 1) = DIM_SQUAREIM
+	IM_LEN(outw, 2) = DIM_SQUAREIM
+
+	if (!helium) {
+	    IM_LEN(outa, 1) = DIM_SQUAREIM
+	    IM_LEN(outa, 2) = DIM_SQUAREIM
+	    IM_PIXTYPE(outa) = TY_INT
+	}
+
+	# Get latitude of sub-earth point from input image header.
+	bzero = imgetr (inputim, "B_ZERO")
+
+	# Ellipse parameters.
+	E_XCENTER(el) = imgetr (inputim, "E_XCEN")
+	E_YCENTER(el) = imgetr (inputim, "E_YCEN")
+	E_XSEMIDIAMETER(el) = imgetr (inputim, "E_XSMD")
+	E_YSEMIDIAMETER(el) = imgetr (inputim, "E_YSMD")
+
+	# Remove the elipse parameters from the header records of the
+	# output images
+
+	call imdelf (outputim, "E_XCEN")
+	call imdelf (outputim, "E_YCEN")
+	call imdelf (outputim, "E_XSMD")
+	call imdelf (outputim, "E_YSMD")
+
+	call imdelf (outw, "E_XCEN")
+	call imdelf (outw, "E_YCEN")
+	call imdelf (outw, "E_XSMD")
+	call imdelf (outw, "E_YSMD")
+	call imaddb (outw, "WEIGHTS", YES)
+
+	if (!helium) {
+	    call imdelf (outa, "E_XCEN")
+	    call imdelf (outa, "E_YCEN")
+	    call imdelf (outa, "E_XSMD")
+	    call imdelf (outa, "E_YSMD")
+	    call imaddb (outa, "ABS_VALU", YES)
+	}
+
+	# Set the variable that keeps track of where in the input image the
+	# bottom of the subraster is, map in the initial subraster.
+
+	inim_subras_bottom = 1
+	inim_subras_ptr = imgs2s (inputim, 1, DIM_VTFD, inim_subras_bottom,
+	    inim_subras_bottom+DIM_IN_RAS-1)
+
+	# Map the outputimages into memory.
+	outim_subras_ptr = imps2i (outputim, 1, DIM_SQUAREIM, 1, DIM_SQUAREIM)
+	outwei_subras_ptr = imps2s (outw, 1, DIM_SQUAREIM, 1, DIM_SQUAREIM)
+	if (!helium)
+	    outav_subras_ptr = imps2i (outa, 1, DIM_SQUAREIM, 1, DIM_SQUAREIM)
+	else {
+	    call smark (sp)
+	    call salloc (outav_subras_ptr, DIM_SQUAREIM*DIM_SQUAREIM, TY_INT)
+	}
+
+	# Initialize meanf, meanaf, numpix.
+	meanf = 0.0
+	meanaf = 0.0
+	numpix = 0
+
+	# Map the input image into the output image by output image rows.
+	do outputrow = 1, DIM_SQUAREIM {
+
+	    # Check the current input subraster to see if it covers
+	    # the next output row to be mapped and map in a new subraster
+	    # if necessary.
+
+	    call checkimmem (inim_subras_bottom, bzero, inputim, outputrow,
+	        inim_subras_ptr, el, skip)
+
+	    # If checkimmem returns skip = true then this row is not contained
+	    # in the input image so fill it with zeros and skip it.
+
+	    if (skip) {
+		# Fill the empty row with zeros.
+		call emptyrow (outputrow, Memi[outim_subras_ptr],
+		    Mems[outwei_subras_ptr], Memi[outav_subras_ptr])
+		next
+	    }
+
+	    # Map this pixel row.
+	    call rowmap (inim_subras_bottom, Mems[inim_subras_ptr], bzero,
+		outputrow, Memi[outim_subras_ptr], Mems[outwei_subras_ptr],
+		Memi[outav_subras_ptr], el, muzero, meanf, meanaf, numpix,
+		helium)
+	}
+
+	# Put the mean field, the number of pixels, the zero corrected mean
+	# absolute field, the mode estimate, the zero corrected mean field,
+	# and the standard deviation in the output image header.
+
+	meanaf = meanaf/double(numpix)
+	meanf = meanf/double(numpix)
+	zcm = meanf - muzero
+	tempr = real(meanf)
+	call imaddr (outputim, "MEAN_FLD", tempr)
+	call imaddi (outputim, "NUMPIX", numpix)
+	if (!helium) {
+	    tempr = real(meanaf)
+	    call imaddr (outputim, "MEANAFLD", tempr)
+	    tempr = real(muzero)
+	    call imaddr (outputim, "MUZERO", tempr)
+	    tempr = real(zcm)
+	    call imaddr (outputim, "ZCM", tempr)
+	}
+
+	# Close images.
+	call imunmap (inputim)
+	call imunmap (outputim)
+	call imunmap (outw)
+	if (!helium)
+	    call imunmap (outa)
+	if (helium)
+	    call sfree (sp)
+end
+
+
+# CHECKIMMEM -- Check this row to see if the input subraster in memory
+# covers it and if it doesn't, map in a new subraster.
+
+procedure checkimmem (inim_subras_bottom, bzero, inputim, outputrow,
+	inim_subras_ptr, el, skip)
+
+int	inim_subras_bottom	# current bottom of the loaded subraster
+real	bzero			# latitude of sub-earth point for this image
+pointer	inputim			# pointer to input image
+int	outputrow		# which output row to map
+pointer	inim_subras_ptr		# input image subraster pointer
+real	el[LEN_ELSTRUCT]	# ellipse parameters data structure
+bool	skip			# returned flag saying to skip this line
+
+real	x ,y
+int	ymax, ymin
+real	uplat, downlat, lminusl0, latitude
+pointer	imgs2s()
+errchk	imgs2s
+
+begin
+	skip = false
+
+	# Find values for the latitudes of the upper and lower edges of this
+	# pixel row.
+
+	uplat = 180./3.1415926*asin(float(outputrow - 90)/90.)
+	downlat = 180./3.1415926*asin(float(outputrow - 91)/90.)
+
+	# Check to see if this row is either completely off the image or
+	# partially off the image.  If it is off the image then return
+	# skip = true.  If it is partially off the image then truncate
+	# the appropriate boundary latitude at the image boundary.
+
+	if (bzero > 0) {
+	    if ( downlat < (-90 + bzero) && uplat < (-90 + bzero)) {
+
+		# This row is not on the image.
+		skip = true
+		return
+	    }
+	    if (downlat < (-90 + bzero))
+		downlat = -90 + bzero
+	} else {
+	    if ( downlat > (90 - bzero) && uplat > (90 - bzero)) {
+
+		# This row is not on the image.
+		skip = true
+		return
+	    }
+	    if (uplat > (90 - bzero))
+		uplat = 90 - bzero
+	}
+
+	# Calculate the minimum and maximum values of y in the input image that
+	# we will need to map this output row of pixels and check these
+	# values against the value of the current bottom of the subraster.
+
+	if (bzero > 0) {
+
+	    # Calculate y position in image.
+	    lminusl0 = 90.
+	    latitude = uplat
+	    call getxy (latitude, lminusl0, bzero, el, x, y, skip)
+	    ymax = int(y + .5)
+	    lminusl0 = -90.
+	    latitude = uplat
+	    call getxy (latitude, lminusl0, bzero, el, x, y, skip)
+	    if (int(y + .5) > ymax)
+		ymax = int(y + .5)
+
+	    # Calculate min y position.
+	    lminusl0 = 0.
+	    latitude = downlat
+	    call getxy (latitude, lminusl0, bzero, el, x, y, skip)
+	    ymin = int(y + .5)
+
+	} else {
+
+	    # Calculate y position in image.
+	    lminusl0 = 90.
+	    latitude = downlat
+	    call getxy (latitude, lminusl0, bzero, el, x, y, skip)
+	    ymin = int(y + .5)
+	    lminusl0 = -90.
+	    latitude = downlat
+	    call getxy (latitude, lminusl0, bzero, el, x, y, skip)
+	    if (int(y + .5) < ymin) ymin = int(y + .5)
+
+	    # Calculate max y position.
+	    lminusl0 = 0.
+	    latitude = uplat
+	    call getxy (latitude, lminusl0, bzero, el, x, y, skip)
+	    ymax = int(y + .5)
+	}
+
+	# If ymin or ymax is outside the current subraster, then map in
+	# an appropriate subraster.
+
+	if ((ymin < (inim_subras_bottom + 5)) ||
+	    (ymax > (inim_subras_bottom + 140))) {
+	    if ((ymax - ymin) > 150) {
+		call printf ("Subraster too small(ymax-ymin > 150), bye")
+	    }
+	    if ((ymin +  144) > 2048) {
+		ymin = 2048 - 144
+	    }
+	    if ((ymin - 5) < 1) {
+		skip = true
+		return
+	    }
+	    inim_subras_ptr = imgs2s (inputim, 1, DIM_VTFD, (ymin - 5),
+		(ymin + 144))
+	    inim_subras_bottom = ymin - 5
+	}
+end
+
+
+# ROWMAP -- Map this output row pixel by pixel.
+
+procedure rowmap (inim_subras_bottom, in_subraster, bzero, outputrow,
+	out_subraster, outw_subraster, outa_subraster, el, muzero, meanf,
+	meanaf, numpix helium)
+
+real	bzero						# lat of sub-earth
+real	el[LEN_ELSTRUCT]		# ellipse parameters data structure
+int	inim_subras_bottom				# bottom of current
+int	outputrow					# output row
+short	in_subraster[DIM_VTFD, DIM_IN_RAS]		# subraster
+int	out_subraster[DIM_SQUAREIM, DIM_SQUAREIM]	# output image
+short	outw_subraster[DIM_SQUAREIM, DIM_SQUAREIM]	# output weights
+int	outa_subraster[DIM_SQUAREIM, DIM_SQUAREIM]	# output abs. value
+double	muzero						# mode estimate
+double	meanf						# mean field
+double	meanaf						# mean absolute field
+int	numpix						# number of pixels
+bool	helium						# 10830 flag
+
+int	pixel
+
+errchk	pixelmap
+
+begin
+	# Do all 180 pixels in this output row.
+	do pixel = 1,180 {
+	    call pixelmap (pixel, in_subraster, inim_subras_bottom,
+		bzero, outputrow, out_subraster, outw_subraster, outa_subraster,
+		el, muzero, meanf, meanaf, numpix, helium)
+	}
+end
+
+
+# PIXELMAP -- Sum up and count the input pixels contained inside the
+# given output pixel.  The sum is carried out in the following way:
+#
+# Calculate, on the input image, the position of the center of the
+#    output pixel to be mapped.
+# Calculate the values of the partial derivitives of latitude and
+#    longitude with respect to x and y.
+# Calculate the boundaries of the pixel in the input image and
+#    sum and count all the pixels inside, assign the value
+#    to the output pixel = sum/count.
+
+procedure pixelmap (pixel, in_subraster, inim_subras_bottom,
+	bzero, outputrow, out_subraster, outw_subraster, outa_subraster,
+	el, muzero, meanf, meanaf, numpix, helium)
+
+int	pixel						# which pixel
+short	in_subraster[DIM_VTFD, DIM_IN_RAS]		# subraster
+int	inim_subras_bottom				# bottom of current
+real	bzero						# lat of sub-earth
+int	outputrow					# output row
+int	out_subraster[DIM_SQUAREIM, DIM_SQUAREIM]	# output image
+short	outw_subraster[DIM_SQUAREIM, DIM_SQUAREIM]	# output weights
+int	outa_subraster[DIM_SQUAREIM, DIM_SQUAREIM]	# output abs. value
+real	el[LEN_ELSTRUCT]		# ellipse parameters data structure
+double	muzero						# mode estimate
+double	meanf						# first moment accum.
+double	meanaf						# mean absolute field
+int	numpix						# number of pixels
+bool	helium						# helium flag
+
+real	lat_mid, long_mid, lat_bot, lat_top
+real	long_rite, long_left
+double	sum, sumabs
+int	count
+real	xpixcenter, ypixcenter
+real	dlongdx, dlatdy
+int	xleft,xright,ybottom,ytop,x,y
+int	num_pix_vert, num_pix_horz
+pointer	sp
+pointer	num					# numeric structure pointer
+real	dat
+
+begin
+	call smark (sp)
+	call salloc (num, VT_LENNUMSTRUCT, TY_STRUCT)
+
+	# First obtain the parameters necessary from numeric.
+	call numeric (bzero, el, outputrow, pixel, xpixcenter, ypixcenter, num)
+
+	dlongdx = VT_DLODX(num)
+	dlatdy = VT_DLATDY(num)
+	lat_top = VT_LATTOP(num)
+	lat_bot = VT_LATBOT(num)
+	long_left = VT_LOLEFT(num)
+	long_rite = VT_LORITE(num)
+	lat_mid = VT_LATMID(num)
+	long_mid = VT_LOMID(num)
+
+	if (lat_top == 10000.) {
+	    out_subraster[pixel,outputrow] = 0
+	    outw_subraster[pixel,outputrow] = 0
+	    outa_subraster[pixel,outputrow] = 0
+	    call sfree (sp)
+	    return
+	}
+
+	# Calculate the box of pixels we want.
+	num_pix_horz = int((1.0 / dlongdx) + .5)
+	xleft = xpixcenter - int((.5 / dlongdx) + .5)
+	xright = xleft + num_pix_horz - 1
+	num_pix_vert = int((abs(abs(lat_top) - abs(lat_bot)) / dlatdy) + .5)
+	ybottom = ypixcenter - int(((abs(abs(lat_mid) - abs(lat_bot))) /
+	    dlatdy) + .5) - (inim_subras_bottom - 1)
+	ytop = ybottom + num_pix_vert - 1
+	
+	# Sum up the pixels inside this box.
+	count = 0
+	sum = 0.0
+	sumabs = 0.0
+
+	do x = xleft, xright {
+	    do y = ybottom, ytop {
+		if (and(int(in_subraster[x,y]),17B) >= THRESHOLD+1) {
+		    count = count + 1
+
+		    # Divide by 16 to remove squibby brightness
+		    # Accumulate the various moment data.
+		    dat = real(in_subraster[x,y]/16)
+		    sum = sum + double(dat)
+		    sumabs = sumabs + double(abs(dat - muzero))
+		}
+	    }
+	}
+
+	outw_subraster[pixel,outputrow] = short(count)
+	out_subraster[pixel,outputrow] = int(sum - double(count*muzero) + .5)
+	if (!helium)
+	    outa_subraster[pixel,outputrow] = int(sumabs + .5)
+	meanf = meanf + sum
+	meanaf = meanaf + sumabs
+	numpix = numpix + count
+
+	call sfree (sp)
+end
+
+
+# EMPTYROW -- Set this row in the output image to zero.
+
+procedure emptyrow (outputrow, out_subraster, outw_subraster, outa_subraster)
+
+int	outputrow
+int	out_subraster[DIM_SQUAREIM, DIM_SQUAREIM]
+short	outw_subraster[DIM_SQUAREIM, DIM_SQUAREIM]
+int	outa_subraster[DIM_SQUAREIM, DIM_SQUAREIM]
+
+int	pixel
+
+begin
+	# Do all 180 pixels in this output row.
+	do pixel = 1,180 {
+	    out_subraster[pixel, outputrow] = 0
+	    outw_subraster[pixel, outputrow] = 0
+	    outa_subraster[pixel, outputrow] = 0
+	}
+end
+
+
+double procedure rmap_mode (inputim, histoname, helium)
+
+pointer	inputim			# Input image
+char	histoname[SZ_FNAME]	# Histogram name
+bool	helium
+
+int	count, i, j
+int	dati, hist_middle
+pointer	imline, histim, hiptr
+int	histo[LEN_HISTO]
+
+# Stuff for mrqmin.
+real    a[3], x[LEN_HISTO], y[LEN_HISTO], sig[LEN_HISTO]
+int     lista[3]
+real    alambda, chisq, covar[3,3], alpha[3,3]
+short	k
+
+pointer	imgl2s(), impl1i(), immap()
+short	shifts()
+
+extern	gauss
+
+begin
+	# Initialize.
+	count = 0
+	k = -4
+	do i = 1, LEN_HISTO
+	    histo[i] = 0
+
+	do i = 1, DIM_VTFD, SPACING{
+	    imline = imgl2s (inputim, i)
+	    do j = 1, DIM_VTFD, SPACING {
+		if (and(int(Mems[imline+j-1]),17B) >= THRESHOLD+1) {
+		    count = count + 1
+		    dati = shifts(Mems[imline+j-1], k)
+
+		    # Put the data into a histogram.
+		    hist_middle = (LEN_HISTO-1)/2 + 1
+		    if (abs(dati) <= hist_middle-1)
+			histo[dati+hist_middle] = histo[dati+hist_middle] + 1
+		}
+	    }
+	}
+
+	# Write this histogram out to an image.
+	histim = immap (histoname, NEW_COPY, inputim)
+	IM_NDIM(histim) = 1
+	IM_LEN(histim, 1) = LEN_HISTO
+	IM_PIXTYPE(histim) = TY_INT
+	hiptr = impl1i (histim)
+
+	# Put the histogram into this image.
+	do i = 1, LEN_HISTO
+	    Memi[hiptr+i-1] = histo[i]
+
+	if (!helium) {
+	    # Set up arrays, etc. for gaussian fit.
+	    a[2] = 1.0
+	    a[1] = real(histo[1])
+	    do i = 1, LEN_HISTO {
+	        x[i] = real(i)
+	        y[i] = real(histo[i])
+	        sig[i] = 1.0
+	        if (histo[i] > a[1]) {
+		    a[1] = real(histo[i])
+		    a[2] = real(i)
+	        }
+	    }
+	    a[3] = 15.0
+
+	    do i = 1, 3
+	        lista[i] = i
+
+	    # Fit the gaussian.
+	    alambda = -1.0
+	    call mrqmin (x, y, sig, LEN_HISTO, a, 3, lista, 3, covar, alpha, 3,
+	        chisq, gauss, alambda)
+	    do i = 1, 10 {
+	        call mrqmin (x, y, sig, LEN_HISTO, a, 3, lista, 3, covar,
+		    alpha, 3, chisq, gauss, alambda)
+	    }
+
+	    call imaddr (histim, "GSS_AMPL", a[1])
+	    call imaddr (histim, "GSS_CNTR", a[2])
+	    call imaddr (histim, "GSS_WDTH", a[3])
+
+	    # Put the mode estimate in the header.
+	    call imaddr (histim, "MUZERO", (a[2] - real(hist_middle)))
+	}
+
+	call imunmap (histim)
+
+	if (helium)
+	    return (0.0)
+	else
+	    return (double(a[2] - real(hist_middle)))
+end