Initial commit

1 files changed, 198 insertions, 0 deletions

diff --git a/src/libcf/cf_astigmatism.c b/src/libcf/cf_astigmatism.c
new file mode 100644
index 0000000..e02f361
--- /dev/null
+++ b/src/libcf/cf_astigmatism.c
@@ -0,0 +1,198 @@
+/*****************************************************************************
+ *              Johns Hopkins University
+ *              Center For Astrophysical Sciences
+ *              FUSE
+ *****************************************************************************
+ *
+ * Description: Modifies X array to correct for instrumental astigmatism.
+ *
+ *		Since we have no astigmatism correction for extended 
+ *		sources, this program does not attempt to apply one.
+ *
+ *              In keeping with the CalFUSE standard, the astigmatism 
+ *		correction is added to the measured X position of each photon.
+ *		Versions of the ASTG_CAL file prior to v009 have the opposite 
+ *		sign and are incompatible with this version of the program.
+ *                    
+ * Returns:     0 upon successful completion.
+ *
+ * History:     12/02/02   1.1   jch	Adapt from cf_astig.c
+ *		12/20/02   1.2   wvd	Install
+ *		02/11/03   1.3   wvd	Apply astig correction only if 
+ *					YCENT# > O for target aperture.
+ *					Add ASTG_KEY keyword.  Comment field:
+ *					Performed? 1=LiF, 2=SiC, 3=Both
+ *		02/13/03   1.4   wvd	Change all indices to type long.
+ *              02/24/03   1.5   peb    Changed include file to calfusettag.h
+ *              03/04/03   1.6   peb    Removed unused variables and math.h.
+ *                                      Corrected sprintf:185 argument error.
+ *		03/11/03   1.7   wvd	Change channel to type char
+ *		03/12/03   1.8   wvd	If wavelength undefined, set
+ *					channel[k] = 0
+ *		04/04/03   1.9   wvd	Test for overflow of astig array.
+ *		04/07/03   1.10  wvd	Delete test for quality of YCENT.
+ *					Delete ASTG_KEY keyword.
+ *					Don't round yoff to nearest integer.
+ *		04/09/03   1.11  wvd	Write name of ASTG_CAL file to trailer.
+ *              05/20/03   1.12  rdr    Added call to cf_proc_check
+ *              09/16/03   1.13  wvd    Change sign of ASTG_CAL file so that
+ *					correction is additive.
+ *              09/22/03   1.14  wvd    Initialize overflow counter.
+ *					Print out number of overflows.
+ *              10/17/03   1.15  wvd    Discard photons that fall outside of 
+ *					the astigmatism-correction window.
+ *              10/31/03   1.16  wvd    Change channel to unsigned char.
+ *              11/11/03   1.17  wvd    Interpolate between tabulated
+ *					wavelength values.  Use cf_read_col
+ *					to read wavecal files.
+ *              02/17/04   1.18  wvd    Replace cf_nint() with cf_nlong()
+ *              03/02/04   1.19  wvd    Change name of assign_wavelengths
+ *					to cf_x2lambda and make it
+ *					generally accessible.
+ *              07/16/04   1.20  wvd    Must change ASTG_COR to COMPLETE or
+ *					SKIPPED by hand.
+ *              02/18/05   1.21  wvd    Add some diagnostic output.
+ *              11/30/05   1.22  wvd    Don't complain when photon events
+ *					fall outside of astig images.
+ *              05/15/06   1.23  wvd    Divide cf_astigmatism_and_dispersion
+ *					into two separate routines.
+ *
+ ****************************************************************************/
+
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "calfuse.h"
+
+static char CF_PRGM_ID[] = "cf_astigmatism";
+static char CF_VER_NUM[] = "1.23";
+
+
+/* 
+ * Loop over the apertures, read appropriate astigmatism correction,
+ * and add to the input X shifts of pixels within extraction window.
+ */
+static int
+add_astig_shifts(fitsfile *infits, long nevents, float *x, float *y,
+	unsigned char *channel)
+{
+    char  	astig_file[FLEN_VALUE], keyword[FLEN_VALUE];
+    int   	status=0;
+    int   	ap, active_ap[2], ast_cen, extno, i, overflow;
+    long  	k, ii, jj, kk, npix, nx, ny;
+    float 	centroid, yoff;
+    float 	*astig=NULL;
+    fitsfile 	*astigfits;
+
+    /* Read target apertures from the file header */
+    (void) cf_source_aper (infits, active_ap);
+
+    /* Open the astigmatism file */
+    FITS_read_key(infits, TSTRING, "ASTG_CAL", astig_file, NULL, &status);
+    FITS_open_file(&astigfits, cf_cal_file(astig_file), READONLY, &status);
+
+    /*   Since we presently have no astig correction for extended
+     *   sources, we apply correction only to target aperture.
+     */
+    for (i = 0; i < 2; i++) {
+	ap = active_ap[i];
+
+	/* Read the centroid of the aperture from the input file. */
+	sprintf(keyword, "YCENT%1d", ap); 
+    	FITS_read_key(infits, TFLOAT, keyword, &centroid, NULL, &status);
+        cf_verbose(3, "Target aperture: %d\t Y centroid: %f", ap, centroid);
+
+	/* Read astigmatism correction for this aperture */
+	extno = ap+1;
+	cf_verbose(3, "Reading extension %d of %s", extno, astig_file);
+	FITS_movabs_hdu(astigfits, extno, NULL, &status);
+	FITS_read_key(astigfits, TINT, "SLIT_CEN", &ast_cen, NULL, &status);
+	FITS_read_key(astigfits, TLONG, "NAXIS1", &nx, NULL, &status);
+	FITS_read_key(astigfits, TLONG, "NAXIS2", &ny, NULL, &status);
+
+	npix = nx * ny;
+	astig = (float *) cf_malloc(sizeof(float) * npix);
+	FITS_read_img(astigfits, TFLOAT, 1L, npix, 0, astig, NULL, &status);
+
+	/* Determine offset between spectrum and astig correction file */
+	yoff = centroid - ast_cen;
+	cf_verbose(3, "Offset between spectrum and astig correction: %f", yoff);
+
+	/*  Go through the event list and find the astigmatism correction
+	 *  appropriate for the photon's (x, y) position.
+	 */
+	overflow = 0;
+	for (k = 0; k < nevents; k++) {
+	    if (channel[k] == ap) {
+		if (((ii = cf_nlong(x[k])) < 0 || ii >= nx) ||
+		    ((jj = cf_nlong(y[k] - yoff)) < 0 || jj >= ny) ||
+		    ((kk = jj * nx + ii) >= npix)) {
+			channel[k] = 0;
+			overflow++;
+			cf_verbose(3, "Cannot correct pixel %05ld (%f, %f)",
+				k, x[k], y[k]);
+		}
+		else x[k] += astig[kk]; 
+	    } 
+	}
+	/* If overflow flag is set, there is a problem either with
+		the ASTG_CAL file or the input photon list.
+	if (overflow)
+    	    cf_if_warning("%d events fell outside of aperture %d ASTG_CAL window",
+		overflow, ap);
+	Comment this out, as it is always triggered by cf_bad_pixels. - wvd (11/30/05) */
+
+        /* Space for astig array is allocated in each loop. */
+	free(astig);
+    	cf_verbose(3, "End of loop for aperture %d", ap);
+    }
+
+    FITS_close_file(astigfits, &status);
+    return status;
+}
+
+
+/*
+ *  If APERTURE = RFPT or target is not a point source, exit without
+ *  applying an astigmatism correction.
+ */
+int cf_astigmatism(fitsfile *infits, long nevents, float *x,
+	float *y, unsigned char *channel)
+{
+    char 	aperture[FLEN_VALUE];
+    char	source_type[FLEN_VALUE];
+    int 	errflg=0, fileok=TRUE, status=0;
+
+    /* Enter a timestamp into the log. */
+    cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+
+    if ((errflg = cf_proc_check(infits, CF_PRGM_ID))) return errflg;
+
+    FITS_read_key(infits, TSTRING, "APERTURE", aperture, NULL, &status);
+    if (!strncmp(aperture, "RFPT", 4)) {
+        cf_verbose(1, "Aperture is RFPT. No astig correction applied.");
+	fileok = FALSE;
+    }
+
+    /* Check source type */
+    FITS_read_key(infits, TSTRING, "SRC_TYPE", source_type, NULL, &status);
+    if (strncmp(source_type, "P", 1)) {
+	cf_verbose(1, "Source type is %s. No astig correction applied.",
+	    source_type);
+	fileok = FALSE;
+    }
+
+    if (fileok) {
+        add_astig_shifts(infits, nevents, x, y, channel);
+	cf_proc_update(infits, CF_PRGM_ID, "COMPLETE");
+    }
+    else
+	cf_proc_update(infits, CF_PRGM_ID, "SKIPPED");
+
+    /* Update processing flags. */
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done Processing");
+
+    return fileok;
+}