From Tauwiki
[edit] Program Operation
[edit] Input Variables
- <Filename>: Contains list of all Swift files to be coadded.
- <gridsize>: Total size of field in degrees.
- <gridspace>: Pixel size in degrees. If left out or 0, uses spacing from first file.
- min_time: Default is 0
[edit] Additional files
- star_list.txt should contain two columns with decimal degrees of ra and dec.
[edit] What does the program do?
- Converts all stellar coordinates from ra and dec to pixel values.
- Read an extension.
- Zero out all the stars. I have assumed a box of 20 pixels on either side of the star.
- Add into the grid. Note that we can select the gridsize.
- Repeat for all extensions.
- Write out a file with two extensions: data in counts/s/pixel and exposure time per pixel.
/*
Will sum all the pixels over extensions in a Swift data file.
We have to work in sky space rather than pixel space.
If a file called stars.txt exists the data around that star will be blanked out before adding.
Libraries needed are CFITSIO, COORD_TRANS, and NRUTIL
Last Modification: 02/09/06
*/
/*Header Definitions*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <fitsio.h>
/*I think these are unused
#include <nrutil.h>
#include <malloc.h>
#include <nr.h>
#include <coord_trans.h>
*/
#include <parse_par.h>
#define BIGERR 1000000.
#define NODATA -1e6
#define SMO 16
#define MIN_EXP_TIME 5
#define r2d (90./atan2(1,0))
#define MAX_STARS 1000 /*Maximum number of point sources*/
#define STAR_SIZE 20 /*Number of pixels to black out*/
#define min(x,y) (((x)<(y)) ? (x) : (y))
#define max(x,y) (((x)>(y)) ? (x) : (y))
struct WCS_HEAD {
double xrefval;
double yrefval;
double xrefpix;
double yrefpix;
double xinc;
double yinc;
double rot;
char coordtype[5];
};
/*Function Definitions*/
/*********************************************************************/
/*Function to optimize centers of grids*/
static void OPTIMIZE_POS(fitsfile *fptr, double *cgridx, double *cgridy,
int hdu_num)
{
double d_tmp;
int status = 0, i_hdr;
char comment[FLEN_COMMENT];
*cgridx = *cgridy = 0;
for (i_hdr = 2; i_hdr <= hdu_num; ++i_hdr){/*Do all headers*/
status = 0;
fits_movabs_hdu(fptr, i_hdr, NULL, &status);
fits_read_key_dbl(fptr, "CRVAL1", &d_tmp, comment, &status);
*cgridx += d_tmp;
fits_read_key_dbl(fptr, "CRVAL2", &d_tmp, comment, &status);
*cgridy += d_tmp;
} /*Endfor*/
*cgridx /= (hdu_num - 1); /*Average positions*/
*cgridy /= (hdu_num - 1);
printf("Center of field is %g %g\n", *cgridx, *cgridy);
}/*End OPTIMIZE_POS*/
/*********************************************************************/
static float ADD_PIXEL(float *grid, float *data, float x, float deltax,
float y, float deltay, long *naxes)
{
int xp = 0, yp = 0, xint, yint;
float xf, yf, data_val = 0, tdata, npoints = 0;
xint = (int) x;
yint = (int) y;
while((yp < deltay) && ((y + yp) < naxes[1]) && ((y - yp) > 0)) {
while((xp < deltax) && ((x + xp) < naxes[0]) && ((x - xp) > 0)) {
xf = min(1, deltax - xp);
yf = min(1, deltay - yp);
tdata = data[(xint - xp) + (yint - yp)*naxes[0]];
if ((isfinite(tdata) != 0) && (tdata > NODATA)){/*If the data are valid*/
data_val += tdata * (xf * yf);
npoints += xf*yf;
}
++xp;
}/*End while*/
xp = 0;
++yp;
}/*End while*/
if (npoints > 0)
*grid += data_val/npoints;
return(npoints);
}/* End ADD_PIXEL*/
/*********************************************************************/
static void SETUP_COORDS(char gal[], char proj[], float gridsize,
long *ngridx, long *ngridy,
struct WCS_HEAD wcs_in, struct WCS_HEAD *wcs_out)
{
double g_lii, g_bii;
if (fabs(wcs_out->xrefval) > 360)
wcs_out->xrefval = wcs_in.xrefval;
if (fabs(wcs_out->yrefval) > 360)
wcs_out->yrefval = wcs_in.yrefval;
if (wcs_out->xinc == 0) {
wcs_out->xinc = -fabs(wcs_in.xinc);
wcs_out->yinc = wcs_in.yinc;
} else {
wcs_out->xinc = -fabs(wcs_out->xinc);/*Reverse coordinates*/
wcs_out->yinc = fabs(wcs_out->xinc);/*Double reversal*/
} /* endif */
strcpy(wcs_out->coordtype, proj);
*ngridx = (long) (gridsize/fabs(wcs_out->xinc));
*ngridy = (long) (gridsize/fabs(wcs_out->yinc));
wcs_out->xrefpix = (long) *ngridx/2.;
wcs_out->yrefpix = (long) *ngridy/2.;
wcs_out->rot = wcs_in.rot;
}/*End SETUP_COORDS*/
/*********************************************************************/
/*---------------------------Begin Main Progam------------------------------*/
int main(int argc, char *argv[], char *envp[])
{
fitsfile *fptr, *out_fptr;
FILE *star_file_ptr=NULL;
int status = 0, anynull, i_arg;
char gal[3];
int num_hdu, i_hdr, nfound;
long ipixel, jpixel, nbad_pts;
long ngrid, ngridx, ngridy, new_gridx, new_gridy;
float npoints = 0;
long naxes[2], npixels, group = 0, out_axes[2];
float *grid = NULL, *tgrid = NULL, *new_grid = NULL, *new_tgrid = NULL;
float nullval = 0, *array = NULL, max_data, max_exp_time;
float f_tmp, f_tmp1, f_tmp2;
double *x_star, *y_star, crpix1, crpix2;
double *ra_star, *dec_star, tmp1, tmp2;
float gridsize, gridspace;
double *ra, *dec, tra, tdec;
double datamin, datamax, exp_time = 0;
double x, y, d_tmp1, d_tmp2;
char str_out[FLEN_KEYWORD];
char out_file_name[120], inp_filename[120];
char status_str[FLEN_STATUS], proj[5];
char comment_str[FLEN_COMMENT];
char date_obs[FLEN_KEYWORD];
int update_coords = 0, istar, nstar;
float min_time;
int min_i, max_i, min_j, max_j;
struct WCS_HEAD wcs_in, wcs_out;
/*Begin Initialization*/
if (argc < 2) {/*Parameter list*/
printf("Parameter order is:\n");
printf("Size of grid in degrees: default 5\n");
printf("Grid spacing: default is that of first file\n");
} /* endif */
/*Read parameters if present*/
i_arg = 2;
if (PARSE_PAR_CHAR(argc, argv, inp_filename, &i_arg)){
printf("Please use filename on command line");
exit(EXIT_FAILURE);
}
if (PARSE_PAR_FLOAT(argc, argv, &gridsize, &i_arg))
gridsize = 1.;
if (PARSE_PAR_FLOAT(argc, argv, &gridspace, &i_arg))
gridspace = 0; /*Default is to use native format*/
if (PARSE_PAR_FLOAT(argc, argv, &min_time, &i_arg))
min_time = 0; /*Default is to use native format*/
printf("Total Extent of integration %g degrees\n", gridsize);
wcs_out.xinc = wcs_out.yinc = gridspace;
wcs_out.xrefval = wcs_out.yrefval = wcs_out.rot = -1000;
strcpy(gal, "eq");/*Use equatorial coordinates*/
strcpy(proj, "-TAN");/*Tangent projection*/
/*Begin processing files*/
status = 0;
fits_open_file(&fptr, inp_filename, READONLY, &status);
fits_get_num_hdus(fptr, &num_hdu, &status);/*Number of HDUS*/
if (status){
printf("Problem opening file with status %i\n", status);
exit(EXIT_FAILURE);
}
/*If ra and dec are not specified then use file headers to get the position*/
if ((wcs_out.xrefval == (-1000)) || (wcs_out.yrefval == (-1000)))
OPTIMIZE_POS(fptr, &wcs_out.xrefval, &wcs_out.yrefval, num_hdu);
status = 0;
x_star = (double *) calloc(sizeof(double), MAX_STARS);
y_star = (double *) calloc(sizeof(double), MAX_STARS);
ra_star = (double *) calloc(sizeof(double), MAX_STARS);
dec_star = (double *) calloc(sizeof(double), MAX_STARS);
/*Read stellar coordinates*/
star_file_ptr = fopen("star_list.txt", "r");
nstar = -1;
if (star_file_ptr != NULL){
while (feof(star_file_ptr) == 0){
fscanf(star_file_ptr, "%f %f\n", &f_tmp1, &f_tmp2);
++nstar;
ra_star[nstar] = f_tmp1;
dec_star[nstar] = f_tmp2;
}
fclose(star_file_ptr);
} else printf("couldn't find file: star_list.txt\n");
for (i_hdr = 2; i_hdr < num_hdu; ++i_hdr){/*Work through each header*/
status = 0;
printf("Processing header: %i\n", i_hdr);
fits_movabs_hdu(fptr, i_hdr, NULL, &status);
/*Get exposure times*/
fits_read_key_dbl(fptr, "EXPOSURE", &exp_time, comment_str, &status);
fits_read_key(fptr, TSTRING, "DATE-OBS", date_obs, comment_str,
&status);
/*WCS keywords*/
ffgics(fptr, &wcs_in.xrefval, &wcs_in.yrefval, &wcs_in.xrefpix,
&wcs_in.yrefpix, &wcs_in.xinc, &wcs_in.yinc, &wcs_in.rot,
&wcs_in.coordtype[0], &status);
if (status > 0) break;
nbad_pts = 0;
for (istar = 0; istar <= nstar; ++istar){/*Convert to pixel coords*/
(fits_world_to_pix(ra_star[istar], dec_star[istar], wcs_in.xrefval,
wcs_in.yrefval, wcs_in.xrefpix, wcs_in.yrefpix,
wcs_in.xinc, wcs_in.yinc, wcs_in.rot,
wcs_in.coordtype, &d_tmp1, &d_tmp2, &status));
status = 0;
x_star[istar] = d_tmp1;
y_star[istar] = d_tmp2;
}/*End convert to pixels*/
/*Set up reference information for output file*/
if (update_coords == 0) { /*We have not yet set up coordinate reference*/
SETUP_COORDS(gal, proj, gridsize, &ngridx, &ngridy, wcs_in, &wcs_out);
/*Inititalize Grid*/
ngrid = ngridx * ngridy;
grid = (float *) malloc(sizeof(float) * ngrid);
memset(grid, 0, sizeof(float) * ngrid);
tgrid = (float *) malloc(sizeof(float) * ngrid);
memset(tgrid, 0, sizeof(float) * ngrid);
printf("Size of array is %i by %i\n",ngridx, ngridy);
if (grid == NULL){
printf("%s\n", "Could not allocate memory");
exit(EXIT_FAILURE);
}/*Endif*/
update_coords = 1;
ra = (double *) malloc(sizeof(double) * ngrid);
dec = (double *) malloc(sizeof(double) * ngrid);
for (jpixel = 0; jpixel < ngridy; ++jpixel) {/*Set up ra and dec array*/
for (ipixel = 0; ipixel < ngridx; ++ipixel) {
status = 0;
fits_pix_to_world(ipixel + 1, jpixel + 1, wcs_out.xrefval,
wcs_out.yrefval, wcs_out.xrefpix, wcs_out.yrefpix,
wcs_out.xinc, wcs_out.yinc, wcs_out.rot,
wcs_out.coordtype, &tra, &tdec, &status);
ra[ipixel + jpixel*ngridx] = tra;
dec[ipixel + jpixel*ngridy] = tdec;
}
}/*setting up ra and dec loop*/
} /* Finished setting up coordinate information for output file*/
/*Read data from file*/
fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status);
npixels = naxes[0] * naxes[1]; /* number of pixels in the image */
array = (float *) malloc(sizeof(float) * npixels);
if (array == NULL) {
printf("%s\n","Could not allocate memory");
exit(EXIT_FAILURE);
} /*Endif*/
if (fits_read_2d_flt(fptr, group, nullval, naxes[0], naxes[0],
naxes[1], array, &anynull, &status))
break;
max_data = -1e6;
max_exp_time = 0;
/*Zero out stars*/
for (istar = 0; istar <= nstar; ++istar){
tmp1 = (y_star[istar] - STAR_SIZE);
tmp2 = (y_star[istar] + STAR_SIZE);
if (tmp1 < 0) tmp1 = 0;
if (tmp2 > naxes[1]) tmp2 = naxes[1];
for (jpixel = max(y_star[istar] - STAR_SIZE, 0);
jpixel < min(y_star[istar] + STAR_SIZE, naxes[1]);
++jpixel){
for (ipixel = max(x_star[istar] - STAR_SIZE, 0);
ipixel < min(x_star[istar] + STAR_SIZE, naxes[0]);
++ipixel)
array[ipixel + jpixel*naxes[0]] = NODATA;
}/*Zero out data*/
}
/*Loop through every pixel in the output image*/
for (jpixel = 0; jpixel < ngridy; ++jpixel) {
for (ipixel = 0; ipixel < ngridx; ++ipixel) {
status = 0;
tra = ra[ipixel + jpixel*ngridx];
tdec = dec[ipixel + jpixel*ngridx];
if (fits_world_to_pix(tra, tdec, wcs_in.xrefval,
wcs_in.yrefval, wcs_in.xrefpix, wcs_in.yrefpix,
wcs_in.xinc, wcs_in.yinc, wcs_in.rot,
wcs_in.coordtype, &x, &y, &status) == 0){
npoints = ADD_PIXEL(grid + ipixel + jpixel*ngridx,
array, x, fabs(wcs_out.xinc/wcs_in.xinc),
y, fabs(wcs_out.yinc/wcs_in.yinc), naxes);
if (npoints > 0)
tgrid[ipixel + jpixel*ngridx] += exp_time;
}/*End add pixel inner loop*/
} /* endfor Ipixel*/
} /* endfor Jpixel*/
free(array);
}/*End hdu read loop*/
/* Write out the file into a bunch of extensions */
/*Get datamin and datamax*/
datamin = 1e6;
datamax = -1e6;
min_i = ngridx;
max_i = 0;
min_j = ngridy;
max_j = 0;
for (jpixel = 0; jpixel < ngridy; ++jpixel) {
for (ipixel = 0; ipixel < ngridx; ++ipixel) {
if (tgrid[ipixel + jpixel*ngridx] > 0)/*Only if there is any exposure*/
grid[ipixel + jpixel*ngridx] /= tgrid[ipixel + jpixel*ngridx];
}/*ipixel*/
}/*End datamin datamax loop*/
for (jpixel = 0; jpixel < ngridy; ++jpixel) {
for (ipixel = 0; ipixel < ngridx; ++ipixel) {
if (tgrid[ipixel + jpixel*ngridx] > min_time){
datamin = min(datamin, grid[ipixel + jpixel*ngridx]);
datamax = max(datamax, grid[ipixel + jpixel*ngridx]);
min_i = min(min_i, ipixel);
max_i = max(max_i, ipixel);
min_j = min(min_j, jpixel);
max_j = max(max_j, jpixel);
}/*Only if the time is greater than the minimum*/
}/*ipixel*/
}/*End datamin datamax loop*/
/*Open output file*/
sprintf(out_file_name,"file1.fit");
remove(out_file_name);
fits_create_file(&out_fptr, out_file_name, &status);
new_gridx = max_i - min_i + 1;
new_gridy = max_j - min_j + 1;
new_grid = (float *) calloc(sizeof(float), new_gridx*new_gridy);
new_tgrid = (float *) calloc(sizeof(float), new_gridx*new_gridy);
for (jpixel = min_j; jpixel < max_j; ++jpixel){
for (ipixel = min_i; ipixel < max_i; ++ipixel){
new_grid[(ipixel - min_i) + (jpixel - min_j)*new_gridx] =
grid[ipixel + jpixel * ngridx];
new_tgrid[(ipixel - min_i) + (jpixel - min_j)*new_gridx] =
tgrid[ipixel + jpixel * ngridx];
}/*ipixel*/
}/*jpixel*/
/*Create 1st HDU Keywords: Remember that the array has been cut short*/
out_axes[0] = new_gridx;
out_axes[1] = new_gridy;
crpix1 = wcs_out.xrefpix - min_i;
crpix2 = wcs_out.yrefpix - min_j;
fits_create_img(out_fptr, FLOAT_IMG, 2, out_axes, &status);
fits_write_key(out_fptr, TDOUBLE, "CRVAL1", &wcs_out.xrefval, "REF POINT VALUE IN DEGREES", &status);
fits_write_key(out_fptr, TDOUBLE, "CRPIX1", &crpix1, "REF POINT PIXEL LOCATION", &status);
fits_write_key(out_fptr, TDOUBLE, "CDELT1", &wcs_out.xinc, "DEGREES PER PIXEL", &status);
fits_write_key(out_fptr, TDOUBLE, "CROTA1", &wcs_out.rot, "ROTATION FROM ACTUAL AXIS", &status);
fits_write_key(out_fptr, TSTRING, "CTYPE1", "RA---TAN", "COORDINATE TYPE", &status);
fits_write_key(out_fptr, TDOUBLE, "CRVAL2", &wcs_out.yrefval, "REF POINT VALUE IN DEGREES", &status);
fits_write_key(out_fptr, TDOUBLE, "CRPIX2", &crpix2, "REF POINT PIXEL LOCATION", &status);
fits_write_key(out_fptr, TDOUBLE, "CDELT2", &wcs_out.yinc, "DEGREES PER PIXEL", &status);
fits_write_key(out_fptr, TDOUBLE, "CROTA2", &wcs_out.rot, "ROTATION FROM ACTUAL AXIS", &status);
fits_write_key(out_fptr, TSTRING, "CTYPE2", "DEC--TAN", "COORDINATE TYPE", &status);
fits_write_key(out_fptr, TSTRING, "EPOCH", "2000", "EQUINOX", &status);
fits_write_key(out_fptr, TDOUBLE, "DATAMIN", &datamin, "Minimum Data Pixel", &status);
fits_write_key(out_fptr, TDOUBLE, "DATAMAX", &datamax, "Minimum Data Pixel", &status);
fits_write_key(out_fptr, TDOUBLE, "EXPOSURE", &exp_time, "Exposure", &status);
fits_write_key(out_fptr, TSTRING, "DATE-OBS", date_obs, "DATE_OBS", &status);
/*Write image data into first HDU */
fits_write_2d_flt(out_fptr, 1, new_gridx, new_gridx, new_gridy, new_grid, &status);
fits_create_img(out_fptr, FLOAT_IMG, 2, out_axes, &status);
fits_write_key(out_fptr, TDOUBLE, "CRVAL1", &wcs_out.xrefval, "REF POINT VALUE IN DEGREES", &status);
fits_write_key(out_fptr, TDOUBLE, "CRPIX1", &crpix1, "REF POINT PIXEL LOCATION", &status);
fits_write_key(out_fptr, TDOUBLE, "CDELT1", &wcs_out.xinc, "DEGREES PER PIXEL", &status);
fits_write_key(out_fptr, TDOUBLE, "CROTA1", &wcs_out.rot, "ROTATION FROM ACTUAL AXIS", &status);
fits_write_key(out_fptr, TSTRING, "CTYPE1", "RA---TAN", "COORDINATE TYPE", &status);
fits_write_key(out_fptr, TDOUBLE, "CRVAL2", &wcs_out.yrefval, "REF POINT VALUE IN DEGREES", &status);
fits_write_key(out_fptr, TDOUBLE, "CRPIX2", &crpix2, "REF POINT PIXEL LOCATION", &status);
fits_write_key(out_fptr, TDOUBLE, "CDELT2", &wcs_out.yinc, "DEGREES PER PIXEL", &status);
fits_write_key(out_fptr, TDOUBLE, "CROTA2", &wcs_out.rot, "ROTATION FROM ACTUAL AXIS", &status);
fits_write_key(out_fptr, TSTRING, "CTYPE2", "DEC--TAN", "COORDINATE TYPE", &status);
fits_write_key(out_fptr, TSTRING, "EPOCH", "2000", "EQUINOX", &status);
fits_write_key(out_fptr, TDOUBLE, "DATAMIN", &datamin, "Minimum Data Pixel", &status);
fits_write_key(out_fptr, TDOUBLE, "DATAMAX", &datamax, "Minimum Data Pixel", &status);
fits_write_key(out_fptr, TDOUBLE, "EXPOSURE", &exp_time, "Exposure", &status);
fits_write_key(out_fptr, TSTRING, "DATE-OBS", date_obs, "DATE_OBS", &status);
/*Write times into second HDU */
fits_write_2d_flt(out_fptr, 1, new_gridx, new_gridx, new_gridy, new_tgrid, &status);
fits_close_file(out_fptr, &status);
printf("%s\n", "done with program");
/*Free variables*/
if (grid != NULL) free(grid);
if (tgrid != NULL) free(tgrid);
free(new_tgrid);
free(new_grid);
free(x_star);
free(y_star);
free(ra_star);
free(dec_star);
return(EXIT_SUCCESS);
}
