/* R. Seaman, NOAO, 12 July 2006
 * based on CFITSIO imcopy.c
 */

#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include "fitsio.h"

#define	SZ_PATH	161

int main(int argc, char *argv[])
{
    fitsfile *infptr, *outfptr;
    int status = 0, ii = 1, iteration = 0, single = 0, multiple = 0, hdupos;
    int hdutype, bitpix, bytepix, naxis = 0, nkeys, datatype = 0, anynul;
    int comptype, narg;
    long naxes[9] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
    long first, totpix = 0, npix;
    double *array, bscale = 1.0, bzero = 0.0, nulval = 0.;
    char card[81], outfile[SZ_PATH], outstring[SZ_PATH];

    if (argc < 2) {
	printf ("Usage:  newimcopy inputImage{s} {'outputImage[compress]'}\n");
	return (0);

    } else if (argc == 3) {
	/* preserve previous command line behavior
	 * this prototype can't compress exactly two FITS files in place
	 */
	strcpy (outstring, argv[2]);

    } else if (argc > 3)
	multiple = 1;

    for (narg = 1; narg < argc; narg++) {
	if (! multiple && narg > 1) break;

	/* arrange a new temporary file name for each
	 * note the devil-may-care usage of strcpy and mktemp
	 */
	strcpy (outfile, "tmp.XXXXXX");
	mktemp (outfile);
	sprintf (outstring, "%s[compress]", outfile);

	/* Open the input file and create output file */
	fits_open_file (&infptr, argv[narg], READONLY, &status);
	fits_create_file (&outfptr, outstring, &status);

	if (status != 0) {    
	    fits_report_error (stderr, status);
	    return (status);
	}

	fits_get_hdu_num (infptr, &hdupos);

	/* Copy only a single HDU if a specific extension was given */ 
	if (hdupos != 1 || strchr (argv[narg], '[')) single = 1;

	for (; !status; hdupos++) {
	    fits_get_hdu_type (infptr, &hdutype, &status);

	    if (hdutype == IMAGE_HDU) {
		for (ii = 0; ii < 9; ii++)
		naxes[ii] = 1;

		fits_get_img_param (infptr, 9, &bitpix, &naxis, naxes, &status);

		totpix = naxes[0] * naxes[1] * naxes[2] * naxes[3] * naxes[4] *
		    naxes[5] * naxes[6] * naxes[7] * naxes[8];
	    }

	    if (hdutype != IMAGE_HDU || naxis == 0 || totpix == 0) { 
		/* just copy tables and null images, update checksums */
		fits_copy_hdu (infptr, outfptr, 0, &status);
		fits_write_chksum (outfptr, &status);

	    } else {
		/* Explicitly create new image, to support compression */
		fits_create_img (outfptr, bitpix, naxis, naxes, &status);

		/* copy all the user keywords (not the structural keywords) */
		fits_get_hdrspace (infptr, &nkeys, NULL, &status); 

		for (ii = 1; ii <= nkeys; ii++) {
		    fits_read_record (infptr, ii, card, &status);
		    if (fits_get_keyclass (card) > TYP_CMPRS_KEY)
			fits_write_record (outfptr, card, &status);
		}

		switch (bitpix) {
		    case BYTE_IMG:
			datatype = TBYTE;
			break;
		    case SHORT_IMG:
			datatype = TSHORT;
			break;
		    case LONG_IMG:
			datatype = TINT;
			break;
		    case FLOAT_IMG:
			datatype = TFLOAT;
			break;
		    case DOUBLE_IMG:
			datatype = TDOUBLE;
			break;
		}

		bytepix = abs (bitpix) / 8;

		npix = totpix;
		iteration = 0;

		/* try to allocate memory for the entire image */
		/* use double type to force memory alignment */
		array = (double *) calloc (npix, bytepix);

		/* if allocation failed, divide size by 2 and try again */
		while (! array && iteration < 10)  {
		    iteration++;
		    npix = npix / 2;
		    array = (double *) calloc (npix, bytepix);
		}

		if (! array)  {
		    printf ("Memory allocation error\n");
		    return (0);
		}

		/* turn off any scaling so that we copy the raw pixel values */
		fits_set_bscale (infptr,  bscale, bzero, &status);
		fits_set_bscale (outfptr, bscale, bzero, &status);

		first = 1;
		while (totpix > 0 && !status) {
		    /* read image in whole or part, write it to the output */
		    fits_read_img (infptr, datatype, first, npix, 
		    &nulval, array, &anynul, &status);

		    fits_write_img (outfptr, datatype, first, npix, array,
			&status);
		    totpix = totpix - npix;
		    first  = first  + npix;
		}

		fits_write_chksum (outfptr, &status);

		free(array);
	    }

	    if (single) break;

	    /* try to move to next HDU */
	    fits_movrel_hdu (infptr, 1, NULL, &status);
	}

	if (status == END_OF_FILE) status = 0;

	fits_close_file (outfptr, &status);
	fits_close_file (infptr, &status);

	if (status)
	   fits_report_error (stderr, status);

	if (argc != 3)
	    rename (outfile, argv[narg]);
    }

    return(status);
}