Initial commit

1 files changed, 736 insertions, 0 deletions

diff --git a/vo/votools/t_vodata.x b/vo/votools/t_vodata.x
new file mode 100644
index 00000000..e6b3de58
--- /dev/null
+++ b/vo/votools/t_vodata.x
@@ -0,0 +1,736 @@
+#
+# VODATA -- Query and Access VO Data Services
+
+include <ctype.h>
+include	<error.h>
+include <imhdr.h>
+include <mwset.h>
+include <tbset.h>
+
+
+define 	DEBUG		FALSE
+define 	DEF_BASE   	"res"
+
+define  BP_NAMES    "|any|radio|millimeter|infrared|optical|uv|x-ray|gamma-ray|"
+define	BP_ANY		1
+define	BP_RADIO	2
+define	BP_MILLIMETER	3
+define	BP_INFRARED	4
+define	BP_OPTICAL	5
+define	BP_UV		6
+define	BP_XRAY		7
+define	BP_GAMMARAY	8
+
+define  TYP_NAMES	"|any|catalog|image|spectral|"
+define	TYP_ANY		1
+define  TYP_CATALOG	2
+define  TYP_IMAGE	3
+define  TYP_SPECTRAL	4
+
+define	OUT_FMTS	"|ascii|csv|tsv|html|raw|fits|"
+define	FMT_ASCII	1
+define	FMT_CSV		2
+define	FMT_TSV		3
+define	FMT_HTML	4
+define	FMT_RAW		5
+define	FMT_FITS	6
+
+
+#    resources	list of Resources to query (or "" or ServiceURL)
+#    objects	list of Object names/images to query (or VOTable of positions)
+#    ra		RA of query position
+#    dec	Dec of query position
+#    sr		Search radius (degrees, or min or sec)
+#
+#    bandpass	bandpass constraint (radio|infrared|optical|uv|xray|gamma)
+#    type	service type to query (catalog|image)
+#    count	print only a count of the results
+#    all	query all data sources
+#    extract	extract positions and URLs
+#
+#    format	output format (|fits|ascii|csv|tsv|html|raw)
+#    output	output filename (or 'samp' to broadcast table)
+#    sequential	use sequential file numbers (or object name)
+#
+#    verbose	verbose output?
+#
+
+
+procedure t_vodata ()
+
+char	resources[SZ_LINE], objects[SZ_FNAME]
+char	format[SZ_FNAME], verb[SZ_FNAME], output[SZ_FNAME]
+char	bandpass[SZ_FNAME], svctype[SZ_FNAME], dbgflag[SZ_FNAME]
+char	fmt[SZ_FNAME], sr[SZ_FNAME], allflag[SZ_FNAME]
+char	resdb[SZ_FNAME], resfile[SZ_FNAME], objfile[SZ_FNAME]
+char	pos[SZ_FNAME], size[SZ_FNAME], ivorn[SZ_FNAME], listflag[SZ_FNAME]
+char	sname[SZ_FNAME], url[SZ_FNAME], type[SZ_FNAME]
+
+pointer	sp, name, tout, old, new, extn
+pointer	cmd, bpass, stype
+int	fd, verbose, ndat
+int	objlist, reslist, nres, nobjs, n
+bool	quiet, do_samp, count, all
+double	ra, dec, srad
+
+int	clgeti(), strdic(), open(), imaccess(), vod_rename()
+int	access(), fntopnb(), fntgfnb(), fntlenb()
+real	clgetr()
+bool	clgetb(), streq(), rdb_lookup()
+
+begin
+	call smark (sp)
+	call salloc (name, SZ_FNAME, TY_CHAR)
+	call salloc (tout, SZ_PATHNAME, TY_CHAR)
+	call salloc (old,  SZ_PATHNAME, TY_CHAR)
+	call salloc (new,  SZ_PATHNAME, TY_CHAR)
+	call salloc (cmd,  SZ_PATHNAME, TY_CHAR)
+	call salloc (extn, SZ_PATHNAME, TY_CHAR)
+
+	call salloc (bpass, SZ_FNAME, TY_CHAR)
+	call salloc (stype, SZ_FNAME, TY_CHAR)
+
+
+	# Get the task parameters.
+	call clgstr ("resources", resources, SZ_LINE)	# query params
+	all        = clgetb ("all")
+	if (!all)
+	    call clgstr ("objects", objects, SZ_LINE)
+	call clgstr ("size", sr, SZ_FNAME)
+
+	count      = clgetb ("count")
+	verbose    = clgeti ("verbose")
+	quiet      = clgetb ("quiet")
+	do_samp    = clgetb ("samp")
+
+	call clgstr ("output", output, SZ_FNAME)	# output params
+	call clgstr ("format", format, SZ_FNAME)
+	call clgstr ("bandpass", bandpass, SZ_FNAME)	# constraints
+	call clgstr ("type", svctype, SZ_FNAME)
+	call clgstr ("resdb", resdb, SZ_FNAME)
+
+
+	################################
+	# Build up arguments
+	################################
+
+	# Get the resources to query.
+	call mktemp ("/tmp/res", resfile, SZ_FNAME)
+	fd = open (resfile, NEW_FILE, TEXT_FILE)
+
+	if (resources[1] == EOS || streq (resources, "any")) {
+	    call fprintf (fd, "any\n")
+	    nres = INDEFI
+	    call strcpy ("+l", listflag, SZ_FNAME)
+	} else {
+	    reslist = fntopnb (resources, NO)
+	    nres = fntlenb (reslist)
+            while (fntgfnb (reslist, Memc[name], SZ_FNAME) != EOF) {
+		if (rdb_lookup (resdb, Memc[name], type, sname, ivorn, url)) {
+		    # Found an alias, print the unambiguous IVORN.
+	            call fprintf (fd, "%s\n")
+		        call pargstr (ivorn)
+		} else {
+		    # Not an alias, let the task resolve it.
+	            call fprintf (fd, "%s\n")
+		        call pargstr (Memc[name])
+		}
+	    }
+	    call fntclsb (reslist)
+	    call strcpy ("+n", listflag, SZ_FNAME)
+	}
+	call close (fd)
+
+
+	# Get the object names to query.
+	call mktemp ("/tmp/obj", objfile, SZ_FNAME)
+	fd = open (objfile, NEW_FILE, TEXT_FILE)
+
+	if (objects[1] != EOS) {
+	    objlist = fntopnb (objects, NO)
+	    nobjs = fntlenb (objlist)
+            while (fntgfnb (objlist, Memc[name], SZ_FNAME) != EOF) {
+
+		if (imaccess (Memc[name], READ_ONLY) == YES) {
+		    call vod_img_query (Memc[name], ra, dec, srad, pos, size)
+
+	    	    call fprintf (fd, "%g %g\n")
+			call pargd (ra)
+			call pargd (dec)
+	    	    call sprintf (sr, SZ_FNAME, "%.6g")
+			call pargd (srad)
+
+		} else {
+		    # Object name, let the task resolve it.
+	            call fprintf (fd, "%s\n")
+			call pargstr (Memc[name])
+		}
+            }
+	    call fntclsb (objlist)
+
+	} else {
+	    # No names given, use the explicit position in the parameters.
+	    nobjs = 1
+	    call fprintf (fd, "%g %g\n")
+		call pargr (clgetr ("ra"))
+		call pargr (clgetr ("dec"))
+	}
+	call close (fd)
+
+
+	# FIXME -- In this release we're limited to image services for blind
+	#	   queries due to the catalog services querying non-std 
+	#	   interfaces in VODATA>
+	# until we figure out the threading error problem.
+	if (IS_INDEFI(nres)) {
+	    if (svctype[1] == EOS || !streq (svctype, "image")) {
+	      call eprintf (
+		"Error: Blind resource queries limited to 'image' services\n")
+	      return
+	    }
+	}
+
+
+	# Check constraints.
+	call aclrc (Memc[bpass], SZ_FNAME)
+	if (bandpass[1] == EOS || bandpass[1] == NULL) {
+	    call strcpy ("any", Memc[bpass], SZ_FNAME)
+	} else {
+	    n = strdic (bandpass, Memc[bpass], SZ_FNAME, BP_NAMES)
+	    if (!quiet && n == 0) {
+	        call eprintf ("Invalid bandpass constraint '%s', ignoring.\n")
+		    call pargstr (bandpass)
+	    }
+	}
+
+	call aclrc (Memc[stype], SZ_FNAME)
+	if (svctype[1] == EOS || svctype[1] == NULL) {
+	    call strcpy ("any", Memc[stype], SZ_FNAME)
+	} else {
+	    n = strdic (svctype, Memc[stype], SZ_FNAME, TYP_NAMES)
+	    if (!quiet && n == 0) {
+	        call eprintf ("Invalid service constraint '%s', ignoring.\n")
+		    call pargstr (svctype)
+	    }
+	}
+
+	# Get output processing options.
+	call vod_get_format (format, fmt, Memc[extn])		
+	if (IS_INDEFI(nres) || do_samp || all) {
+	    # Unknown multiple resources.
+	    call strcpy (output, Memc[tout], SZ_FNAME)
+	} else
+	    call strcpy ("/tmp/vod", Memc[tout], SZ_FNAME)
+	if (output[1] == EOS)
+	    call strcpy ("vod", output, SZ_FNAME)	# set default name
+
+
+	call strcpy ("+n", verb, SZ_FNAME)		# query verbosity
+	if (verbose > 1)
+	    call strcpy ("-v", verb, SZ_FNAME)
+	if (verbose > 2)
+	    call strcpy ("-vv", verb, SZ_FNAME)
+
+	if (all)					# query all data svcs?
+	    call strcpy ("-a", allflag, SZ_FNAME)
+	else
+	    call strcpy ("+n", allflag, SZ_FNAME)
+
+	if (access ("/tmp/VOC_VDEBUG", 0, 0) == YES)	# debug flag
+	    call strcpy ("+d", dbgflag, SZ_FNAME)
+	else
+	    call strcpy ("+n", dbgflag, SZ_FNAME)
+
+
+	######################################
+	# Call the VODATA task interface.
+	######################################
+	call vx_vodata (16,    				# argc
+			"-t", Memc[stype],		# type constraint
+	        	"-b", Memc[bpass], 		# bandpass constraint
+	        	"-O", Memc[tout], 		# output filename (temp)
+			verb,				# query verbosity
+			allflag,			# query all data?
+			dbgflag,			# debug flag
+			listflag,			# svc_list flag
+			"-q",				# suppress output
+			"-N",				# numeric output name
+			"-R",				# output format
+			resfile,			# resources to query
+			objfile, sr)			# objects to query
+
+
+	if (!quiet && (nres > 1 || IS_INDEFI(nres)) && streq(output,"STDOUT")) {
+	    call printf("# Rows  Cols  File         \tResource Title\n")
+	    call printf("#\n")
+	}
+
+	# Clean up the output filenames.  For a single resource/object we
+	# still end up with a name file 'foo_000_000.xml', use a more 
+	# logical naming scheme.
+	if (!all) {
+	    ndat = vod_rename (output, Memc[tout], format, Memc[extn], 
+	        nres, nobjs, quiet, count, do_samp)
+	}
+
+	if (!quiet && (nres > 1 || IS_INDEFI(nres)) && streq(output,"STDOUT")) {
+	    call printf ("#\n# Data found for %d of %d resources queried.\n")
+	        call pargi (ndat)
+	        call pargi (nres)
+	}
+
+	# Clean up the temporary resource and object/position files.
+	if (access (resfile, 0, 0) == YES)
+	    call delete (resfile)
+	if (access (objfile, 0, 0) == YES)
+	    call delete (objfile)
+end
+
+
+#  VOD_RENAME -- Clean up the output filenames.  For a single resource/object
+#  we still end up with a name file 'foo_000_000.xml', use a more logical 
+#  naming scheme.
+
+int procedure vod_rename (output, tout, format, extn, nres, nobjs, 
+				quiet, count, samp)
+
+char	output[ARB]				#i output base name
+char	tout[ARB]				#i temp output name
+char	format[ARB]				#i output format
+char	extn[ARB]				#i filename extension
+int	nres					#i number of resources
+int	nobjs					#i number of objects
+bool	quiet					#i print summary?
+bool	count					#i print result count only?
+bool	samp					#i broadcast as SAMP message?
+
+pointer	sp, line, svc, obj, dir, new, old
+int	i, j, ip, sfd, nch, nfound, nrows, ncols, stat
+
+int	access(), stridx(), strldx(), getline(), open(), ctoi()
+int	vot_convert()
+bool	streq()
+
+begin
+	call smark (sp)
+	call salloc (svc, SZ_FNAME, TY_CHAR)
+	call salloc (obj, SZ_FNAME, TY_CHAR)
+	call salloc (new, SZ_FNAME, TY_CHAR)
+	call salloc (old, SZ_FNAME, TY_CHAR)
+	call salloc (line, SZ_LINE, TY_CHAR)
+	call salloc (dir, SZ_PATHNAME, TY_CHAR)
+
+	nrows = 0
+	ncols = 0
+	nfound = 0
+
+	# Get the cwd.
+	call fpathname (".", Memc[dir], SZ_PATHNAME)
+	ip = stridx ('/', Memc[dir])
+	call strcpy (Memc[dir+ip-1], Memc[dir], SZ_PATHNAME)
+
+	call sprintf (Memc[obj], SZ_FNAME, "%s.objects")
+	    call pargstr (output)
+
+	sfd = NULL
+	if (IS_INDEFI(nres)) {
+	    # Get the number of services queried.
+	    call sprintf (Memc[svc], SZ_FNAME, "%s.services")
+		call pargstr (output)
+	    sfd = open (Memc[svc], READ_ONLY, TEXT_FILE)
+	    nch = getline (sfd, Memc[line])
+
+	    ip = 20
+	    nch = ctoi (Memc[line], ip, nres)
+
+	    nch = getline (sfd, Memc[line]) 		# Skip ahead
+	    nch = getline (sfd, Memc[line])
+	}
+
+	# Process the output files as needed.
+
+	for (i=0; i < nres; i=i+1) {
+	    if (IS_INDEFI(nres))
+	        nch = getline (sfd, Memc[line])
+
+	    for (j=0; j < nobjs; j=j+1) {
+		# Construct the vodata filename.
+	    	call sprintf (Memc[old], SZ_FNAME, "%s_%03d_%03d.xml")
+		    call pargstr (tout)
+		    call pargi (i)
+		    call pargi (j)
+
+	    	if (streq ("STDOUT", output) || streq ("STDERR", output)) {
+		    call strcpy (output, Memc[new], SZ_FNAME)
+		    if (access (Memc[old], 0, 0) == YES) {
+	    	        call fcopy (Memc[old], Memc[new])
+	        	call delete (Memc[old])
+			nfound = nfound + 1
+		    }
+
+	    	} else {
+
+		    if (nres == 1 && nobjs == 1) {
+	    	        call sprintf (Memc[new], SZ_FNAME, "%s%s%s")
+		    	    call pargstr (Memc[dir])
+			    call pargstr (output)
+			    call pargstr (extn)
+		    } else if (nres == 1) {
+	    	        call sprintf (Memc[new], SZ_FNAME, "%s%s_%03d%s")
+		    	    call pargstr (Memc[dir])
+			    call pargstr (output)
+			    call pargi (j)
+			    call pargstr (extn)
+		    } else if (nobjs == 1) {
+	    	        call sprintf (Memc[new], SZ_FNAME, "%s%s_%03d%s")
+		    	    call pargstr (Memc[dir])
+			    call pargstr (output)
+			    call pargi (i)
+			    call pargstr (extn)
+		    } else {
+	    	        call sprintf (Memc[new], SZ_FNAME, "%s%s_%03d_%03d%s")
+		    	    call pargstr (Memc[dir])
+			    call pargstr (output)
+			    call pargi (i)
+			    call pargi (j)
+			    call pargstr (extn)
+		    }
+
+		    if (access (Memc[old], 0, 0) == YES) {
+			if (!quiet) {
+			    call vod_tinfo (Memc[old], nrows, ncols)
+			    ip = strldx ('/', Memc[new])
+			    call eprintf (" %5d  %4d  %s")
+			        call pargi (nrows)
+			        call pargi (ncols)
+			        call pargstr (Memc[new+ip])
+			    if (nres > 1) {
+			        call eprintf ("  %s")
+			            call pargstr (Memc[new+22])
+			    } else
+			        call eprintf ("\n")
+			}
+
+			if (count) {
+			    # If ony printing a count, delete saved results.
+		            if (access (Memc[old], 0, 0) == YES)
+	    	                call delete (Memc[old])
+			} else {
+
+		            if (access (Memc[new], 0, 0) == NO) {
+			        if (!streq (extn, "xml")) {
+            			    stat = vot_convert (Memc[old], Memc[new],
+					format)
+	    	                    call delete (Memc[old])
+			        } else 
+	    	                    call frename (Memc[old], Memc[new])
+			    }
+
+	    	    	    if (samp) 			# broadcast the table
+			        call vod_bcast_table (Memc[new])
+			}
+
+			nfound = nfound + 1
+		    }
+	    	}
+		    
+		if (DEBUG) {
+		    call eprintf ("%d/%d: %s -> %s\n")
+			call pargi (nres) ; call pargi (nobjs)
+			call pargstr (Memc[old])
+			call pargstr (Memc[new])
+		}
+	    }
+	}
+
+			
+	if (count || streq (output, "STDOUT")) {
+	    # If ony printing a count, delete saved results.
+	    if (access (Memc[svc], 0, 0) == YES)
+		call delete (Memc[svc])
+	    if (access (Memc[obj], 0, 0) == YES)
+		call delete (Memc[obj])
+	}
+	if (sfd != NULL)
+	    call close (sfd) 
+	call sfree (sp)
+
+	return (nfound)
+end
+
+
+# VOD_TINFO -- Get table dimensions.
+
+procedure vod_tinfo (tblname, nrows, ncols)
+
+char	tblname[ARB]				#i table name
+int	nrows, ncols				#o table dimensions
+
+pointer	tp, tbtopn()
+int	tbpsta()
+
+begin
+        iferr {
+            tp = tbtopn (tblname, READ_ONLY, 0)
+        } then {
+            call eprintf ("Error: can't open '%s'\n")
+                call pargstr (tblname)
+	    return
+	}
+ 
+        nrows   = tbpsta (tp, TBL_NROWS)
+        ncols   = tbpsta (tp, TBL_NCOLS)
+
+	call tbtclo (tp)
+end
+
+
+# VOD_BCAST_TABLE -- Broadcast a table name for loading.
+
+procedure vod_bcast_table (tblname)
+
+char	tblname[ARB]				#i table name
+
+pointer	sp, cwd, cmd, name
+int	status
+
+begin
+	call smark (sp)
+	call salloc (cwd, SZ_PATHNAME, TY_CHAR)
+	call salloc (cmd, SZ_PATHNAME, TY_CHAR)
+	call salloc (name, SZ_PATHNAME, TY_CHAR)
+
+	# Get the current working dir
+	call zfgcwd (Memc[cwd], SZ_PATHNAME, status)
+       	call strupk (Memc[cwd], Memc[cwd], SZ_PATHNAME)
+
+	# Broadcast the table
+	call sprintf (Memc[cmd], SZ_PATHNAME, "samp loadVOTAble file://%s\n")
+	    call pargstr (tblname)
+
+	#  FIXME -- Need to be using the real SAMP interface here.
+	call clcmd (Memc[cmd])
+
+	call sfree (sp)
+end
+
+
+# VOD_GET_FORMAT -- Convert the format parameter to the flag and file extn.
+
+procedure vod_get_format (format, fmt, extn)
+
+char	format[ARB]				#i format parameter
+char	fmt[ARB]				#o format flag
+char	extn[ARB]				#o output extension
+
+int	ofmt
+int	strdic()
+
+begin
+	ofmt = strdic (format, fmt, SZ_FNAME, OUT_FMTS)
+	switch (ofmt) {
+	case FMT_ASCII:   
+	    call strcpy ("-A", fmt, SZ_FNAME)
+	    call strcpy (".asv", extn, SZ_FNAME)
+	case FMT_CSV: 	  
+	    call strcpy ("-C", fmt, SZ_FNAME)
+	    call strcpy (".csv", extn, SZ_FNAME)
+	case FMT_TSV: 	  
+	    call strcpy ("-T", fmt, SZ_FNAME)
+	    call strcpy (".tsv", extn, SZ_FNAME)
+	case FMT_HTML: 	  
+	    call strcpy ("-H", fmt, SZ_FNAME)
+	    call strcpy (".html", extn, SZ_FNAME)
+	case FMT_RAW: 	  
+	    call strcpy ("-R", fmt, SZ_FNAME)
+	    call strcpy (".xml", extn, SZ_FNAME)
+	case FMT_FITS: 	  
+	    call strcpy ("-F", fmt, SZ_FNAME)
+	    call strcpy (".fits", extn, SZ_FNAME)
+	default: 	  
+	    call strcpy ("-R", fmt, SZ_FNAME)
+	    call strcpy (".xml", extn, SZ_FNAME)
+	}
+end
+
+
+# VOD_TBL_COPY -- Copy a single table to a new file.  If the file exists, append
+# as a new extension.
+
+procedure vod_tbl_copy (oldfile, newfile)
+
+char	oldfile[ARB]	# i: current file name
+char	newfile[ARB]	# i: new file name
+
+int	phu_copied	# set by tbfpri and ignored
+pointer	sp, oldname, newname
+
+bool	use_fcopy	# true if we should copy the file with fcopy
+
+bool	streq()
+int	tbtacc()
+errchk	tbfpri, tbtcpy
+
+begin
+	call smark (sp)
+	call salloc (oldname, SZ_FNAME, TY_CHAR)
+	call salloc (newname, SZ_FNAME, TY_CHAR)
+
+	# Check to make sure the copy is legal
+
+	use_fcopy = false
+	if (streq (oldfile, newfile)) {
+	    call eprintf ("Cannot copy table to itself:  %s\n")
+	    call pargstr (oldfile)
+
+	    if (tbtacc (oldfile) == YES)
+		use_fcopy = true
+
+	    if (use_fcopy) {
+		call tbtext (oldfile, Memc[oldname], SZ_FNAME)
+		call tbtext (newfile, Memc[newname], SZ_FNAME)
+
+		iferr (call fcopy (Memc[oldname], Memc[newname])) {
+		    call erract (EA_WARN)
+		}
+	    }
+
+	} else {
+	    # Table extensions are copied by the table
+	    # library function tbtcpy
+	    iferr {
+		call tbfpri (oldfile, newfile, phu_copied)
+		call tbtcpy (oldfile, newfile)
+	    } then {
+		call erract (EA_WARN)
+	    }
+	}
+
+	call sfree (sp)
+	return
+end
+
+
+#  VOD_IMG_QUERY -- Get the query params for the named image.
+
+procedure vod_img_query (imname, ra, dec, sr, pos, size)
+
+char	imname[ARB]				#i image name
+double	ra, dec, sr				#o RA/DEC/SR params
+char	pos[ARB]				#o POS string
+char	size[ARB]				#o SIZE string
+
+double  r[IM_MAXDIM], w[IM_MAXDIM], cd[2,2]
+double	xrot, yrot, rot, scale
+
+pointer	im, mw, ctw, co
+int	stat
+double	cx, cy, ra0, dec0
+
+pointer	immap()
+int	sk_decim(), mw_stati()
+pointer	mw_sctran()
+
+begin
+	#  Open the image.
+  	iferr {
+            im = immap (imname, READ_ONLY, 0)
+	} then {
+            # Unable to decode image WCS
+	    call eprintf ("Error: Cannot open image '%s'\n")
+		call pargstr (imname)
+	    return
+	}
+
+
+	#  Get the WCS.
+  	iferr {
+            stat = sk_decim (im, "world", mw, co)
+            if (stat == ERR || mw == NULL) {
+            	# Unable to decode image WCS
+		call eprintf ("Error: No image WCS present.\n")
+		return
+
+            } else if (mw != NULL) {
+                ctw = mw_sctran (mw, "logical", "world", 03B)
+		cx  = IM_LEN(im,1) / 2
+		cy  = IM_LEN(im,2) / 2
+        	call mw_c2trand (ctw, cx,  cy,  ra0, dec0)
+
+		# Get the CD matrix for scale and rotation
+           	call vod_gfterm (mw, r, w, cd, mw_stati (mw, MW_NPHYSDIM))
+        	scale = 3600. * sqrt ((cd[1,1]**2 + cd[2,1]**2 +
+				       cd[1,2]**2 + cd[2,2]**2) / 2.)
+
+        	xrot  = abs (atan2 ( cd[2,1], cd[1,1]))
+        	yrot  = abs (atan2 (-cd[1,2], cd[2,2]))
+        	rot   = (xrot + yrot) / 2.0		# NOT USED
+
+
+		ra = ra0
+		dec = dec0
+		sr =  max (((scale * IM_LEN(im,1)) / 3600.),
+	      		    (scale * IM_LEN(im,2)) / 3600.)
+
+		call sprintf (pos, SZ_FNAME, "%.6g,%.6g")
+		    call pargd (ra0)
+		    call pargd (dec0)
+		call sprintf (size, SZ_FNAME, "%.4g,%.4g")
+		    call pargd ((scale * IM_LEN(im,1)) / 3600.)
+		    call pargd ((scale * IM_LEN(im,2)) / 3600.)
+	    }
+
+	} then {
+            # Unable to decode image WCS
+	    call eprintf ("Error: Unable to decode WCS.\n")
+	}
+
+	call imunmap (im)
+end
+
+
+# VOD_GFTERM -- Compute the output FITS CRPIX, CRVAL, and CD arrays from the
+# MWCS LTERM and WTERM. Note that the CD matrix terms are still transposed
+# from the usual Fortran order.
+
+procedure vod_gfterm (mw, crpix, crval, cd, ndim)
+
+pointer mw              			#i the input mwcs pointer
+double  crpix[ndim]     			#o the output FITS CRPIX array
+double  crval[ndim]     			#o the output FITS CRVAL array
+double  cd[ndim,ndim]   			#o the output FITS CD matrix
+int     ndim            			#i the dimensionality of the wcs
+
+pointer sp, r, wcd, ltv, ltm, iltm
+int	i
+
+errchk	mw_gwtermd, mw_gltermd
+
+begin
+	call smark (sp)
+	call salloc (r, ndim, TY_DOUBLE)
+	call salloc (wcd, ndim * ndim, TY_DOUBLE)
+	call salloc (ltv, ndim, TY_DOUBLE)
+	call salloc (ltm, ndim * ndim, TY_DOUBLE)
+	call salloc (iltm, ndim * ndim, TY_DOUBLE)
+
+	iferr {
+	    call mw_gwtermd (mw, Memd[r], crval, Memd[wcd], ndim)
+	    call mw_gltermd (mw, Memd[ltm], Memd[ltv], ndim)
+	    call mwvmuld (Memd[ltm], Memd[r], crpix, ndim)
+	    call aaddd (crpix, Memd[ltv], crpix, ndim)
+	    call mwinvertd (Memd[ltm], Memd[iltm], ndim)
+	    call mwmmuld (Memd[wcd], Memd[iltm], cd, ndim)
+
+	} then {
+	    # Set up a default value.
+	    call aclrd (cd, ndim*ndim)
+	    for (i=1; i <= ndim; i=i+1) {
+	        crpix[i] = 1.0d0
+	        crval[i] = 1.0d0
+	        cd[i,i] = 1.0d0
+	    }
+	}
+
+	call sfree (sp)
+end