1 files changed, 177 insertions, 0 deletions
diff --git a/noao/imred/vtel/numeric.x b/noao/imred/vtel/numeric.x
new file mode 100644
index 00000000..640778c8
--- /dev/null
+++ b/noao/imred/vtel/numeric.x
@@ -0,0 +1,177 @@
+include <mach.h>
+include "vt.h"
+include "numeric.h"
+
+# NUMERIC -- calculate some of the necessary information, including
+# the partial derivitives of latitude and longitude with respect
+# to x and y, that we need to do the projection.
+
+procedure numeric (bzero, el, outputrow, pixel, xpixcenter, ypixcenter, num)
+
+real	bzero				# latitude of subearth point
+real	el[LEN_ELSTRUCT]		# ellipse parameters data structure
+int	outputrow			# which output row are we working on
+int	pixel				# which pixel in that output row
+real	xpixcenter, ypixcenter		# coordinates of center of pixel
+pointer	num				# numeric structure pointer
+
+real	dlatdy, dlongdx			# partial derivitives
+real	lat_top, lat_bot		# latitude of top and bottom of pix
+real	long_left, long_rite		# longitude of left and right of pix
+real	lat_mid, long_mid		# latitude and longitude of middle
+real	lat1, long1, lat3, long3, lat4, long4, lat5, long5
+real	x1, y1, x3, y3, x4, y4, x5, y5
+bool	skip
+
+begin
+	skip = false
+
+	# First calculate lats, longs for this pixel.
+	lat_top   = 180./3.1415926*asin(real(outputrow - 90)/90.)
+	lat_bot   = 180./3.1415926*asin(real(outputrow - 91)/90.)
+	long_left = real(pixel - 1) - 90.
+	long_rite = real(pixel) - 90.
+	lat_mid   = .5 * (lat_top + lat_bot)
+	long_mid  = .5 * (long_left + long_rite)
+
+	# Check the proximity of this pixel to the image boundary, if its
+	# too close, set that output pixel to zero.
+
+	if (abs(abs(lat_mid) - 90.0) < abs(bzero)) {
+	    if (abs(abs(lat_top) - 90.0) >= abs(bzero)) {
+		lat_bot = -90.0 + bzero
+		lat_mid = .5 * (lat_top + lat_bot)
+	    } else {
+		if (abs(abs(lat_bot) - 90.0) >= abs(bzero)) {
+		    lat_top = 90.0 + bzero
+		    lat_mid = .5 * (lat_top + lat_bot)
+		} else {
+	    	    # Nothing to map!
+		    # Flag to pixelmap marking zero pixel.
+		    VT_LATTOP(num) = 10000.
+	    	    return
+		}
+	    }
+	} else {
+	    if (abs(abs(lat_top) - 90.0) < abs(bzero))
+		lat_top = 90.0 + bzero
+	    else
+		if (abs(abs(lat_bot) - 90.0) < abs(bzero))
+		    lat_bot = -90.0 + bzero
+	}
+
+	# Now that we have the pixel we want defined, calculate the partial
+	# derivitives we need numerically.  First calculate the latitude and
+	# longitude of the centers of the 4 adjacent pixels.
+
+	lat1 = lat_mid
+	if (pixel == 1)
+	    long1 = long_mid
+	else
+	    long1 = long_mid - 1.0
+
+	lat3 = lat_mid
+	if (pixel == 180)
+	    long3 = long_mid
+	else
+	    long3 = long_mid + 1.0
+
+	long5 = long_mid
+	if (outputrow == 1)
+	    lat5 = lat_mid
+	else
+	    lat5 = 180./3.1415926*((asin(real(outputrow - 92)/90.) +
+		    asin(real(outputrow - 91)/90.))/2.)
+
+	long4 = long_mid
+	if (outputrow == 180)
+	    lat4 = lat_mid
+	else
+	    lat4 = 180./3.1415926*((asin(real(outputrow - 89)/90.) +
+	            asin(real(outputrow - 90)/90.))/2.)
+
+	# Given these latitudes and longitudes, find out where in xy coords
+	# they are. Get xpixcenter and ypixcenter then the x#s and y#s.
+
+	call getxy (lat_mid, long_mid, bzero, el, xpixcenter, ypixcenter, skip)
+	if (skip) {
+
+	    # Off the limb or behind the sun.
+	    # Flag to pixelmap marking zero pixel.
+	    VT_LATTOP(num) = 10000.
+	    return
+	}
+
+	call getxy (lat1, long1, bzero, el, x1, y1, skip)
+	call getxy (lat3, long3, bzero, el, x3, y3, skip)
+	call getxy (lat4, long4, bzero, el, x4, y4, skip)
+	call getxy (lat5, long5, bzero, el, x5, y5, skip)
+
+	# Calculate the partials.
+	if (x3 == x1)
+	    dlongdx = 9999.
+	else
+	    dlongdx = (long3 - long1) / (x3 - x1)
+
+	if (y4 == y5)
+	    dlatdy = 9999.
+	else
+	    dlatdy = (lat4 - lat5) / (y4 - y5)
+
+	VT_DLODX(num) = dlongdx 
+	VT_DLATDY(num) = dlatdy 
+	VT_LATTOP(num) = lat_top 
+	VT_LATBOT(num) = lat_bot 
+	VT_LOLEFT(num) = long_left 
+	VT_LORITE(num) = long_rite 
+	VT_LATMID(num) = lat_mid 
+	VT_LOMID(num) = long_mid 
+end
+
+
+# GETXY -- Given the latitude and longitude of a point and the image
+# parameters, return the x and y position of that point.
+
+procedure getxy (lat, lml0, b0, el, x, y, skip)
+
+real	lat			# latitude of point on image
+real	lml0			# distance in longitude from disk center
+real	b0			# latitude of sub earth point
+real	el[LEN_ELSTRUCT]	# ellipse parameters data structure
+real	x, y			# returned position
+bool	skip			# skip flag
+
+real	sinlat, coslat, sinbzero, cosbzero, sinlminusl0, coslminusl0
+real	cosrho, sinrho, sinpminustheta, cospminustheta 
+real	latitude, lminusl0, bzero
+
+begin
+	skip = false
+	lminusl0 = lml0*3.1415926/180.
+	bzero = b0*3.1415926/180.
+	latitude = lat*3.1415926/180.
+	sinlat = sin(latitude)
+	coslat = cos(latitude)
+	sinbzero = sin(bzero)
+	cosbzero = cos(bzero)
+	sinlminusl0 = sin(lminusl0)
+	coslminusl0 = cos(lminusl0)
+	cosrho = sinbzero * sinlat + cosbzero * coslat * coslminusl0
+
+	# If we are behind limb return skip = true.
+	if (cosrho <= 0.00) skip = true
+	sinrho = (1. - cosrho**2)**.5
+	if (sinrho >= EPSILONR) {
+	    sinpminustheta = (coslat/sinrho) * sinlminusl0
+	    cospminustheta = (coslat/sinrho) * (cosbzero * tan(latitude) -
+	        sinbzero * coslminusl0)
+	} else {
+	    sinpminustheta = 0.000001
+	    cospminustheta = 0.000001
+	}
+
+	x = E_XSEMIDIAMETER(el) * sinrho * sinpminustheta
+	y = E_YSEMIDIAMETER(el) * sinrho * cospminustheta
+	x = x + real(E_XCENTER(el))
+	y = y + real(E_YCENTER(el))
+end