Repatch (from linux) of OSX IRAF

1 files changed, 762 insertions, 0 deletions

diff --git a/noao/imred/vtel/merge.x b/noao/imred/vtel/merge.x
new file mode 100644
index 00000000..79aa23eb
--- /dev/null
+++ b/noao/imred/vtel/merge.x
@@ -0,0 +1,762 @@
+include <mach.h>
+include	<imhdr.h>
+include	"vt.h"
+
+# MERGE -- Put together all appropriate daily grams to produce a full
+# carrington rotation map.  This is done both for the average input images
+# and for the absolute value input images.  The output of the program is
+# 4 images, average image, absolute value image, weight image, ratio of
+# first image to second image.
+
+procedure t_merge()
+
+char	mergelist[SZ_FNAME]		# list of images to be merged
+
+int	wavelength, listfd
+char	inputimage[SZ_FNAME]
+pointer	inputim
+
+pointer	immap()
+int	imgeti(), open(), fscan()
+errchk	immap, open
+
+begin
+	# Get the image name file from the cl and open it.
+	call clgstr ("mergelist", mergelist, SZ_FNAME)
+	listfd = open (mergelist, READ_ONLY, TEXT_FILE)
+
+	# Get the wavelength from the first image in the mergelist.
+	if (fscan (listfd) != EOF) {
+	    call gargwrd (inputimage, SZ_FNAME)
+	    inputim = immap (inputimage, READ_ONLY, 0)
+	    wavelength = imgeti (inputim, "WV_LNGTH")
+	    call close (listfd)
+	} else {
+	    call error (0, "No images in 'mergelist'")
+	    call close (listfd)
+	    return
+	}
+
+	if (wavelength == 8688)
+	    call mergem (mergelist, wavelength)
+	else
+	    call mergeh (mergelist, wavelength)
+end
+
+
+# MERGEM -- MERGE Magnetograms.
+
+procedure mergem (mergelist, wavelength)
+
+char	mergelist[SZ_FNAME]		# list of images to be merged
+int	wavelength			# wavelength of images
+
+pointer	outputim, outw, outa, outr
+pointer outptr, outwptr, outaptr, outrptr
+char	outputimage[SZ_FNAME], outweight[SZ_FNAME]
+char	outabs[SZ_FNAME], outratio[SZ_FNAME]
+real	longout, weight_tbl[SZ_WTBL], bzeroave
+int	i, mapmonth, mapday, mapyear
+
+real	clgetr()
+int	clgeti()
+pointer	immap(), imps2r()
+errchk	immap, imps2r
+
+begin
+	# Get parameters from the cl.
+
+	# Output images.
+	call clgstr ("outputimage", outputimage, SZ_FNAME)
+	call clgstr ("outweight", outweight, SZ_FNAME)
+	call clgstr ("outabs", outabs, SZ_FNAME)
+	call clgstr ("outratio", outratio, SZ_FNAME)
+
+	# Longitude of center of output Carrington rotation map.
+	longout = clgetr ("longout")
+
+	# Month, day, and year of the center of the output map.
+	mapmonth = clgeti ("mapmonth")
+	mapday = clgeti ("mapday")
+	mapyear = clgeti ("mapyear")
+
+	# Open output image.
+	outputim = immap (outputimage, NEW_IMAGE, 0)
+
+	# Define some parameters for the output images.
+	IM_NDIM(outputim) = 2
+	IM_LEN(outputim, 1) = DIM_XCARMAP
+	IM_LEN(outputim, 2) = DIM_SQUAREIM
+	IM_PIXTYPE(outputim) = TY_REAL
+
+	# Open the rest of the output images.
+	outw = immap (outweight, NEW_COPY, outputim)
+	outa = immap (outabs, NEW_COPY, outputim)
+	outr = immap (outratio, NEW_COPY, outputim)
+
+	# Map the outputimages into memory.
+	outptr = imps2r (outputim, 1, DIM_XCARMAP, 1, DIM_SQUAREIM)
+	outwptr = imps2r (outw, 1, DIM_XCARMAP, 1, DIM_SQUAREIM)
+	outaptr = imps2r (outa, 1, DIM_XCARMAP, 1, DIM_SQUAREIM)
+	outrptr = imps2r (outr, 1, DIM_XCARMAP, 1, DIM_SQUAREIM)
+
+	# Create weight table.
+	do i = 1,SZ_WTBL
+	    weight_tbl[i] = (cos((real(i-91)+.5)*3.1415926/180.))**4
+
+	call mmall (mergelist, Memr[outptr], Memr[outwptr], Memr[outaptr],
+	    outputim, outw, outa, outr, wavelength, weight_tbl, longout,
+	    mapmonth, mapday, mapyear, bzeroave)
+
+	# Fill the ratio image.
+	call imratio (Memr[outptr],Memr[outaptr],Memr[outrptr],DIM_XCARMAP,
+	    DIM_SQUAREIM)
+
+	# Write some information out to the image headers.
+	call imaddr (outputim, "AV_BZERO", bzeroave)
+	call imaddi (outputim, "WV_LNGTH", wavelength)
+	call imaddr (outw, "AV_BZERO", bzeroave)
+	call imaddr (outw, "WV_LNGTH", wavelength)
+	call imaddb (outw, "WEIGHTS", TRUE)
+	call imaddr (outa, "AV_BZERO", bzeroave)
+	call imaddr (outr, "AV_BZERO", bzeroave)
+	call imaddr (outa, "WV_LNGTH", wavelength)
+	call imaddr (outr, "WV_LNGTH", wavelength)
+	call imaddb (outa, "ABS_VALU", TRUE)
+	call imaddb (outr, "POLARITY", TRUE)
+	
+	# Weight the data image and the abs image.
+	call imratio (Memr[outptr],Memr[outwptr],Memr[outptr],DIM_XCARMAP,
+	    DIM_SQUAREIM)
+	call imratio (Memr[outaptr],Memr[outwptr],Memr[outaptr],DIM_XCARMAP,
+	    DIM_SQUAREIM)
+
+	# Close images
+	call imunmap (outputim)
+	call imunmap (outw)
+	call imunmap (outa)
+	call imunmap (outr)
+end
+
+
+# MERGEH -- MERGE Helium 10830 grams.
+
+procedure mergeh (mergelist, wavelength)
+
+char	mergelist[SZ_FNAME]		# list of images to merge
+int	wavelength			# wavelength of observation
+
+pointer	outputim, outw
+pointer outptr, outwptr
+char	outputimage[SZ_FNAME], outweight[SZ_FNAME]
+real	longout, weight_tbl[SZ_WTBL], bzeroave
+int	i, mapmonth, mapday, mapyear
+
+real	clgetr()
+int	clgeti()
+pointer	immap(), imps2r()
+errchk	immap, imps2r
+
+begin
+	# Get parameters from the cl.
+
+	# Output images.
+	call clgstr ("outputimage", outputimage, SZ_FNAME)
+	call clgstr ("outweight", outweight, SZ_FNAME)
+
+	# Longitude of center of output Carrington rotation map.
+	longout = clgetr ("longout")
+
+	# Month, day, and year of the center of the output map.
+	mapmonth = clgeti ("mapmonth")
+	mapday = clgeti ("mapday")
+	mapyear = clgeti ("mapyear")
+
+	# Open output image.
+	outputim = immap (outputimage, NEW_IMAGE, 0)
+
+	# Define some parameters for the output images.
+	IM_NDIM(outputim) = 2
+	IM_LEN(outputim, 1) = DIM_XCARMAP
+	IM_LEN(outputim, 2) = DIM_SQUAREIM
+	IM_PIXTYPE(outputim) = TY_REAL
+
+	# Open the other output image.
+	outw = immap (outweight, NEW_COPY, outputim)
+
+	# Map the outputimages into memory.
+	outptr = imps2r (outputim, 1, DIM_XCARMAP, 1, DIM_SQUAREIM)
+	outwptr = imps2r (outw, 1, DIM_XCARMAP, 1, DIM_SQUAREIM)
+
+	# Create weight table.
+	do i = 1,SZ_WTBL
+	    weight_tbl[i] = (cos((real(i-91)+.5)*3.1415926/180.))**4
+
+	call mhall (mergelist, Memr[outptr], Memr[outwptr], outputim,
+	    outw, wavelength, weight_tbl, longout, mapmonth,
+	    mapday, mapyear, bzeroave)
+
+	# Write some information out to the image headers.
+	call imaddr (outputim, "AV_BZERO", bzeroave)
+	call imaddi (outputim, "WV_LNGTH", wavelength)
+	call imaddr (outw, "AV_BZERO", bzeroave)
+	call imaddr (outw, "WV_LNGTH", wavelength)
+	call imaddb (outw, "WEIGHTS", TRUE)
+	
+	# Weight the data image.
+	call imratio (Memr[outptr],Memr[outwptr],Memr[outptr],DIM_XCARMAP,
+	    DIM_SQUAREIM)
+
+	# Close images.
+	call imunmap (outputim)
+	call imunmap (outw)
+end
+
+
+# MMALL -- Merge Magnetograms ALL.
+# Map in each input image, weight it, figure out where it goes
+# and add it to the output image.
+
+procedure mmall (mergelist, outarray, outarrayw, outarraya, outputim,
+	outw, outa, outr, wavelength, weight_tbl, longout, mapmonth, mapday,
+	mapyear, bzeroave)
+
+char	mergelist[SZ_FNAME]			# list of images to be merged
+int	wavelength				# wavelength of observations
+real	outarray[DIM_XCARMAP, DIM_SQUAREIM]	# output data array
+real	outarrayw[DIM_XCARMAP, DIM_SQUAREIM]	# output weights array
+real	outarraya[DIM_XCARMAP, DIM_SQUAREIM]	# output absolute value array
+pointer	inputim					# pointer to input image
+pointer	outputim				# pointer to output image
+pointer	outw					# pointer to weight image
+pointer	outa					# pointer to abs value image
+pointer	outr					# pointer to ratio image
+int	mapmonth, mapday, mapyear		# date of output map
+real	weight_tbl[SZ_WTBL]			# weight table
+real	longout					# longitude of map center
+real	bzeroave				# average b-zero for map
+
+char	inputimage[SZ_FNAME], inweight[SZ_FNAME], inabs[SZ_FNAME]
+pointer	inw, ina, inptr, inwptr, inaptr
+int	listfd, month, day, year, count
+real	longin, bzero, bzerosum
+int	obsdate, temp, i, j
+char	ltext[SZ_LINE]
+
+int	open(), fscan(), imgeti()
+real	imgetr()
+pointer	immap(), imgs2i(), imgs2s()
+errchk	open, immap, imgs2i, imgs2s
+
+begin
+	count = 0
+	bzerosum = 0.0
+	listfd = open (mergelist, READ_ONLY, TEXT_FILE)
+
+	# Zero the output images.
+	do i = 1, DIM_XCARMAP {
+	    do j = 1, DIM_SQUAREIM {
+		outarray[i,j] = 0.0
+		outarrayw[i,j] = 0.0
+		outarraya[i,j] = 0.0
+	    }
+	}
+
+	# Get inputimages from the mergelist until they are all used up.
+	while (fscan (listfd) != EOF) {
+	    call gargwrd (inputimage, SZ_FNAME)
+
+	    # Get absolute value image.
+	    if(fscan (listfd) != EOF)
+	        call gargwrd (inabs, SZ_FNAME)
+	    else
+	        call error (0, "wrong number of file names in mergelist")
+
+	    # Get weight image.
+	    if(fscan (listfd) != EOF)
+	        call gargwrd (inweight, SZ_FNAME)
+	    else
+	        call error (0, "wrong number of file names in mergelist")
+
+	    # Open input image, its corresponding weight map, and its
+	    # corresponding absolute value map.
+
+	    inputim = immap (inputimage, READ_ONLY, 0)
+	    inw = immap (inweight, READ_ONLY, 0)
+	    ina = immap (inabs, READ_ONLY, 0)
+
+	    bzero = imgetr (inputim, "B_ZERO")
+	    bzerosum = bzerosum + bzero
+	    longin = imgetr (inputim, "L_ZERO")
+	    obsdate = imgeti (inputim, "OBS_DATE")
+
+	    # Check to see that the date is same on the three input images.
+	    temp = imgeti (inw, "OBS_DATE")
+ 	    if (temp != obsdate) {
+		call eprintf ("ERROR: date on weight image differs from that ")
+		call eprintf ("on data image!\n")
+		break
+	    }
+
+	    temp = imgeti (ina, "OBS_DATE")
+	    if (temp != obsdate) {
+		call eprintf ("ERROR: date on abs image differs from that ")
+	        call eprintf ("on data image!\n")
+	        break
+	    }
+
+	    # Decode month, day, year.
+	    month = obsdate/10000
+	    day = obsdate/100 - 100 * (obsdate/10000)
+	    year = obsdate - 100 * (obsdate/100)
+
+	    # Pack a name for this date and longitude and then put them out
+	    # into the outputimages' headers.
+
+	    count = count + 1
+	    call sprintf (ltext, SZ_LINE, "DATE%04d")
+		call pargi (count)
+	    call imaddi (outputim, ltext, obsdate)
+	    call imaddi (outw, ltext, obsdate)
+	    call imaddi (outa, ltext, obsdate)
+	    call imaddi (outr, ltext, obsdate)
+
+	    call sprintf (ltext, SZ_LINE, "LONG%04d")
+		call pargi (count)
+	    call imaddr (outputim, ltext, longin)
+	    call imaddr (outw, ltext, longin)
+	    call imaddr (outa, ltext, longin)
+	    call imaddr (outr, ltext, longin)
+	    
+	    # Map the inputimage, the weight map, and abs_image into memory.
+	    inptr = imgs2i (inputim, 1, DIM_SQUAREIM, 1, DIM_SQUAREIM)
+	    inwptr = imgs2s (inw, 1, DIM_SQUAREIM, 1, DIM_SQUAREIM)
+	    inaptr = imgs2i (ina, 1, DIM_SQUAREIM, 1, DIM_SQUAREIM)
+
+	    # Weight this image and add it to the output image.
+	    call addmweight (Memi[inptr],Mems[inwptr],Memi[inaptr],outarray,
+	        outarrayw, outarraya, weight_tbl, longin, longout,
+		month, day, year, mapmonth, mapday, mapyear)
+
+	    # Close this input image.
+	    call imunmap (inputim)
+	    call imunmap (inw)
+	    call imunmap (ina)
+
+	} # end of do loop on input images
+
+	bzeroave = bzerosum/real(count)
+	call close (listfd)
+end
+
+
+# MHALL -- Merge Heliumgrams ALL.
+# Map in each input image, weight it, figure out where it goes
+# and add it to the output image.
+
+procedure mhall (mergelist, outarray, outarrayw, outputim,
+	outw, wavelength, weight_tbl, longout, mapmonth, mapday,
+	mapyear, bzeroave)
+
+char	mergelist[SZ_FNAME]			# list of images to be merged
+int	wavelength				# wavelength of observations
+real	outarray[DIM_XCARMAP, DIM_SQUAREIM]	# output data array
+real	outarrayw[DIM_XCARMAP, DIM_SQUAREIM]	# output weights array
+pointer	inputim					# pointer to input image
+pointer	outputim				# pointer to output image
+pointer	outw					# pointer to weight image
+int	mapmonth, mapday, mapyear		# date of output map
+real	weight_tbl[SZ_WTBL]			# weight table
+real	longout					# longitude of map center
+real	bzeroave				# average b-zero for map
+
+char	inputimage[SZ_FNAME], inweight[SZ_FNAME]
+pointer	inw, inptr, inwptr
+int	listfd, month, day, year, count
+real	longin, bzero, bzerosum
+int	obsdate, temp, i, j
+char	ltext[SZ_LINE]
+
+real	imgetr()
+int	open(), fscan(), imgeti()
+pointer	immap(), imgs2i(), imgs2s()
+errchk	open, immap, imgs2i, imgs2s
+
+begin
+	count = 0
+	bzerosum = 0.0
+	listfd = open (mergelist, READ_ONLY, TEXT_FILE)
+
+	# Zero the output images.
+	do i = 1, DIM_XCARMAP {
+	    do j = 1, DIM_SQUAREIM {
+		outarray[i,j] = 0.0
+		outarrayw[i,j] = 0.0
+	    }
+	}
+
+	# Get inputimages from the mergelist until they are all used up.
+	while (fscan (listfd) != EOF) {
+	    call gargwrd (inputimage, SZ_FNAME)
+
+	    # Get weight image.
+	    if (fscan (listfd) != EOF)
+	        call gargwrd (inweight, SZ_FNAME)
+	    else
+	        call error (0, "wrong number of file names in mergelist")
+
+	    # Open input image, its corresponding weight map, and its
+	    # corresponding absolute value map.
+
+	    inputim = immap (inputimage, READ_ONLY, 0)
+	    inw = immap (inweight, READ_ONLY, 0)
+
+	    bzero = imgetr (inputim, "B_ZERO")
+	    bzerosum = bzerosum + bzero
+	    longin = imgetr (inputim, "L_ZERO")
+	    obsdate = imgeti (inputim, "OBS_DATE")
+
+	    # Check to see that the date is same on the three input images.
+	    temp = imgeti (inw, "OBS_DATE")
+ 	    if (temp != obsdate) {
+		call eprintf ("ERROR: date on weight image differs from that ")
+		call eprintf ("on data image!\n")
+		break
+	    }
+
+	    # Decode month, day, year.
+	    month = obsdate/10000
+	    day = obsdate/100 - 100 * (obsdate/10000)
+	    year = obsdate - 100 * (obsdate/100)
+
+	    # Pack a name for this date and longitude and then put them out
+	    # into the outputimages' headers.
+
+	    count = count + 1
+	    call sprintf (ltext, SZ_LINE, "DATE%04d")
+		call pargi (count)
+	    call imaddi (outputim, ltext, obsdate)
+	    call imaddi (outw, ltext, obsdate)
+
+	    call sprintf (ltext, SZ_LINE, "LONG%04d")
+		call pargi (count)
+	    call imaddr (outputim, ltext, longin)
+	    call imaddr (outw, ltext, longin)
+	    
+	    # Map the inputimage, the weight map, and abs_image into memory.
+	    inptr = imgs2i (inputim, 1, DIM_SQUAREIM, 1, DIM_SQUAREIM)
+	    inwptr = imgs2s (inw, 1, DIM_SQUAREIM, 1, DIM_SQUAREIM)
+
+	    # Weight this image and add it to the output image.
+	    call addhweight (Memi[inptr], Mems[inwptr], outarray, outarrayw,
+		weight_tbl, longin, longout, month, day, year, mapmonth,
+		mapday, mapyear)
+
+	    # Close this input image.
+	    call imunmap (inputim)
+	    call imunmap (inw)
+
+	} # end of do loop on input images
+
+	bzeroave = bzerosum/real(count)
+	call close (listfd)
+end
+
+
+# ADDMWEIGHT -- Weight input image by cos(longitude - (L-L0))**4, and add
+# it to the output image in the proper place.
+
+procedure addmweight (inim, inwim, inaim, outim, outwim, outaim,
+	weight_tbl, longin, longout, month, day, year, mapmonth, mapday,
+	mapyear)
+
+int	inim[DIM_SQUAREIM, DIM_SQUAREIM]	# input image
+short	inwim[DIM_SQUAREIM, DIM_SQUAREIM]	# input image weights
+int	inaim[DIM_SQUAREIM, DIM_SQUAREIM]	# input absolute image
+real	outim[DIM_XCARMAP, DIM_SQUAREIM]	# outputimage
+real	outwim[DIM_XCARMAP, DIM_SQUAREIM]	# output image weights
+real	outaim[DIM_XCARMAP, DIM_SQUAREIM]	# output absolute image
+int	month, day, year			# date of input image
+int	mapmonth, mapday, mapyear		# date of output image
+real	weight_tbl[DIM_SQUAREIM]		# weight table
+real	longin, longout				# longitudes of images
+
+int	p1offset, p2offset, firstpix, lastpix, column, row
+int	offset, datein, dateout, temp, temp2
+int	d1900()
+
+begin
+	# Translate the two dates into julian day numbers to make comparisons
+	# simpler.
+
+	datein = d1900 (month, day, year)
+	dateout = d1900 (mapmonth, mapday, mapyear)
+
+	# Figure out the pixel offset between the first pixel of the input
+	# image and the first pixel of ther output image.
+	# Actually, there may be two pixel offsets for a particular image
+	# corresponding to the correct position of the image and the 360
+	# degree offset position.
+
+	p1offset = mod(abs(int(longin - longout + .5)), 360)    # This is one.
+	p2offset = 360 - p1offset                         # This is the other.
+
+	# Determine which side of the output image center is each of these
+	# offsets.
+
+	if (datein > dateout) {
+	    if (longout > 180) {
+	        if (((longin >= longout) && (longin <= 360))  || 
+		    (longin <= mod((longout + 180.),360.))) {
+		    if (p1offset < 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+	        } else {
+		    if (p1offset >= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+		}
+	    } else {
+	        if ((longin >= longout) && (longin <= (longout + 180))) {
+		    if (p1offset <= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+	        } else {
+		    if (p1offset >= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+		}
+	    }
+        } else {
+	    if (longout < 180) {
+	        if (((longin >= (180 + longout)) && (longin <= 360))  || 
+		    (longin <= longout)) {
+		    if (p1offset < 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+	        } else {
+		    if (p1offset >= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+		}
+	    } else {
+	        if ((longin < longout) && (longin > (longout - 180))) {
+		    if (p1offset < 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+	        } else {
+		    if (p1offset >= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+		}
+	    }
+	}
+
+        # Make sure the sign is right
+        if (datein > dateout)
+	    offset = -offset
+
+        # Check for the case that the two longitudes are equal.
+        if (longin == longout) {
+	    if (abs(datein - dateout) <= 1) {
+	        offset = 0
+	    } else {
+	        call eprintf ("input day too far from center of output map\n")
+	        return
+	    }
+        }
+	    
+        # Check for the case that the two dates are equal.
+        if (datein == dateout)
+	    offset = longin - longout
+
+        # If the offset is too large then do not use this image.
+        if (abs(offset) > 240) {
+	    call eprintf ("input day too far from center of output map\n")
+	    return
+        }
+
+        # Determine what part, if not all, of the input image will lie on the
+        # output image.
+
+        firstpix = 1
+        if (offset < -90)
+	    firstpix = abs(offset+90)
+        lastpix = DIM_SQUAREIM
+        if (offset > 90) 
+	    lastpix = 180 - (offset - 90)
+
+
+        # Do all 180 columns in the image.
+	if (offset <= 0)
+	    temp = 91
+	else
+	    temp = 90
+
+	do column = firstpix,lastpix {
+	    do row = 1, DIM_SQUAREIM {
+	        temp2 = column + temp + offset
+	        outim[temp2,row] = outim[temp2, row] +
+		    inim[column, row] * weight_tbl[column]
+	        outwim[temp2,row] = outwim[temp2, row] +
+		    inwim[column, row] * weight_tbl[column]
+	        outaim[temp2,row] = outaim[temp2, row] +
+		    inaim[column, row] * weight_tbl[column]
+	    }
+	}
+end
+
+
+# ADDHWEIGHT -- Weight input image by cos(longitude - (L-L0))**4, and add
+# it to the output image in the proper place. (For 10830 grams)
+
+procedure addhweight (inim, inwim, outim, outwim, weight_tbl, longin, longout,
+	month, day, year, mapmonth, mapday, mapyear)
+
+int	inim[DIM_SQUAREIM, DIM_SQUAREIM]	# input image
+short	inwim[DIM_SQUAREIM, DIM_SQUAREIM]	# input image weights
+real	outim[DIM_XCARMAP, DIM_SQUAREIM]	# outputimage
+real	outwim[DIM_XCARMAP, DIM_SQUAREIM]	# output image weights
+int	month, day, year			# date of input image
+int	mapmonth, mapday, mapyear		# date of output image
+real	weight_tbl[DIM_SQUAREIM]		# weight table
+real	longin, longout				# longitudes of images
+
+int	p1offset, p2offset, firstpix, lastpix, column, row
+int	offset, datein, dateout, temp, temp2
+int	d1900()
+
+begin
+	# Translate the two dates into julian day numbers to make comparisons
+	# simpler.
+
+	datein = d1900 (month, day, year)
+	dateout = d1900 (mapmonth, mapday, mapyear)
+
+	# Figure out the pixel offset between the first pixel of the input
+	# image and the first pixel of ther output image.
+	# Actually, there may be two pixel offsets for a particular image
+	# corresponding to the correct position of the image and the 360
+	# degree offset position.
+
+	p1offset = mod(abs(int(longin - longout + .5)), 360)    # this is one.
+	p2offset = 360 - p1offset                         # this is the other.
+
+	# Determine which side of the output image center is each of these
+	# offsets.
+
+	if (datein > dateout) {
+	    if (longout > 180) {
+	        if (((longin >= longout) && (longin <= 360))  || 
+		    (longin <= mod((longout + 180.),360.))) {
+		    if (p1offset < 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+	        } else {
+		    if (p1offset >= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+		}
+	    } else {
+	        if ((longin >= longout) && (longin <= (longout + 180))) {
+		    if (p1offset <= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+	        } else {
+		    if (p1offset >= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+		}
+	    }
+        } else {
+	    if (longout < 180) {
+	        if (((longin >= (180 + longout)) && (longin <= 360))  || 
+		    (longin <= longout)) {
+		    if (p1offset < 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+	        } else {
+		    if (p1offset >= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+		}
+	    } else {
+	        if ((longin < longout) && (longin > (longout - 180))) {
+		    if (p1offset < 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+	        } else {
+		    if (p1offset >= 180)
+		        offset = p2offset
+		    else
+		        offset = p1offset
+		}
+	    }
+	}
+
+        # Make sure the sign is right.
+        if (datein > dateout)
+	    offset = -offset
+
+        # Check for the case that the two longitudes are equal.
+        if (longin == longout) {
+	    if (abs(datein - dateout) <= 1) {
+	        offset = 0
+	    } else {
+	        call eprintf ("Input day too far from center of output map.\n")
+	        return
+	    }
+        }
+	    
+        # Check for the case that the two dates are equal.
+        if (datein == dateout)
+	    offset = longin - longout
+
+        # If the offset is too large then do not use this image.
+        if (abs(offset) > 240) {
+	    call eprintf ("input day too far from center of output map\n")
+	    return
+        }
+
+        # Determine what part, if not all, of the input image will lie on the
+        # output image.
+
+        firstpix = 1
+        if (offset < -90)
+	    firstpix = abs(offset+90)
+        lastpix = DIM_SQUAREIM
+        if (offset > 90) 
+	    lastpix = 180 - (offset - 90)
+
+
+        # Do all 180 columns in the image.
+	if (offset <= 0)
+	    temp = 91
+	else
+	    temp = 90
+
+	do column = firstpix, lastpix {
+	    do row = 1, DIM_SQUAREIM {
+	        temp2 = column + temp + offset
+	        outim[temp2,row] = outim[temp2, row] +
+		    inim[column, row] * weight_tbl[column]
+	        outwim[temp2,row] = outwim[temp2, row] +
+		    inwim[column, row] * weight_tbl[column]
+	    }
+	}
+end