1 files changed, 106 insertions, 0 deletions
diff --git a/vendor/voclient/voapps/zzwcs.c b/vendor/voclient/voapps/zzwcs.c
new file mode 100644
index 00000000..e681db51
--- /dev/null
+++ b/vendor/voclient/voapps/zzwcs.c
@@ -0,0 +1,106 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "fitsio.h"
+
+
+int voc_imgQParams (char *name, double *ra, double *dec, double *rad);
+
+
+int main (int argc, char *argv[])
+{
+    double  ra, dec, rad;
+    int     status = 0;
+
+    status = voc_imgQParams (argv[1], &ra, &dec, &rad);
+    printf ("RA=%f  DEC=%f  SZ=%f\n", ra, dec, rad);
+}
+
+
+int voc_imgQParams (char *name, double *ra, double *dec, double *rad)
+{
+    fitsfile *fptr;
+
+    double  xpos=0.0, ypos=0.0, xpix=0.0, ypix=0.0;
+    double  xrval=0.0, yrval=0.0, xrpix=0.0, yrpix=0.0;
+    double  xinc=0.0, yinc=0.0, rot=0.0;
+    double  ll_x=0.0, ll_y=0.0, ur_x=0.0, ur_y=0.0, c_x=0.0, c_y=0.0;
+    long    naxes[3], pcount, gcount;
+    int     status=0, extend, simple, naxis, bitpix;
+    char    ctype[5];
+
+
+    /*  Open the FITS file.
+     */
+    if (ffopen (&fptr, name, READWRITE, &status) > 0) {
+	fprintf (stderr, "Error: open status = %d\n", status);
+	exit (0);
+    }
+
+    /*  Get the primary header keywords.
+     */
+    ffghpr (fptr, 99, &simple, &bitpix, &naxis, naxes, &pcount,
+           &gcount, &extend, &status);
+
+fprintf (stderr, "gcount = %d  pcount = %d  extend = %d \n", gcount, pcount, extend);
+
+    /*  Get the header WCS keywords.
+     */
+    ffgics (fptr, &xrval, &yrval, &xrpix,
+               &yrpix, &xinc, &yinc, &rot, ctype, &status);
+    if (status != 506)
+        fprintf (stderr, "Read WCS keywords with ffgics status = %d\n",status);
+
+
+    xpix = 0.5;
+    ypix = 0.5;
+    status = 0;
+    if (ffwldp (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, ctype,
+               &ll_x, &ll_y, &status) > 0) {
+        fprintf (stderr, "Calculated sky coordinate with ffwldp status = %d\n",
+	    status);
+
+    } else {
+        status = 0;
+        xpix = (double) naxes[0];
+        ypix = (double) naxes[1];
+        ffwldp (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, ctype,
+           &ur_x, &ur_y, &status);
+
+        status = 0;
+        xpix = (double) naxes[0] / 2.;
+        ypix = (double) naxes[1] / 2.;
+        ffwldp (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, ctype,
+           &c_x, &c_y, &status);
+    }
+
+    *ra   = c_x;
+    *dec  = c_y;
+    *rad  = sqrt ((c_x - ll_x)*(c_x - ll_x) + (c_y - ll_y)*(c_y - ll_y));
+
+#ifdef DBG_WCS
+    printf ("  CRVAL1, CRVAL2 = %16.12f, %16.12f\n", xrval, yrval);
+    printf ("  CRPIX1, CRPIX2 = %16.12f, %16.12f\n", xrpix, yrpix);
+    printf ("  CDELT1, CDELT2 = %16.12f, %16.12f\n", xinc, yinc);
+    printf ("  Rotation = %10.3f, CTYPE = %s\n", rot, ctype);
+    printf ("  Pixels (%8.4f,%8.4f) --> (%11.6f, %11.6f) Sky\n",
+        xpix, ypix, ll_x, ll_y);
+
+    fprintf (stderr, "RA=%f  DEC=%f  SZ=%f\n", c_x, c_y,
+	sqrt ((c_x - ll_x)*(c_x - ll_x) + (c_y - ll_y)*(c_y - ll_y)));
+
+    ffxypx (xpos, ypos, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, ctype,
+           &xpix, &ypix, &status);
+    printf ("Calculated pixel coordinate with ffxypx status = %d\n", status);
+    printf ("  Sky (%11.6f, %11.6f) --> (%8.4f,%8.4f) Pixels\n",
+            xpos, ypos, xpix, ypix);
+#endif
+
+
+    if (ffclos (fptr, &status) > 0) {
+	fprintf (stderr, "Error: close status = %d\n", status);
+	return (-1);
+    }
+
+    return (0);
+}