/*Subtract SKYVAL background from 2MASS fits image
R. Stiening 09/06/2006
Usage: ./sub_backgd 2mass.fits
Output file will value of SKYVAL set to zero.
Output file name is input file name with leading s_ added.
Compile command: gcc -Wall sub_backgd.c -lcfitsio -lm -osub_backgd 
*/
#include <string.h>
#include <stdio.h>
#include "fitsio.h"
int main(int argc, char *argv[])
{
    fitsfile *afptr, *outfptr;  /* FITS file pointers */
    int status = 0;  /* CFITSIO status value MUST be initialized to zero! */
    int anaxis, keytype, check = 1, ii;
    long npixels = 1, firstpix[3] = {1,1,1},anaxes[3] = {1,1,1};
    double *apix,bkgd;
    char card[FLEN_CARD],newcard[FLEN_CARD],oldvalue[FLEN_VALUE], comment[FLEN_COMMENT];
    char outfilename[20];
    if (argc != 2) { 
      printf("USAGE: sub_backgd inimage \n");
      return(0);
    }
    fits_open_file(&afptr, argv[1], READONLY, &status); /* open input image */
    fits_get_img_dim(afptr, &anaxis, &status);  /* read dimensions */
    fits_get_img_size(afptr, 3, anaxes, &status);
    if (status) {
       fits_report_error(stderr, status); /* print error message */
       return(status);
    }
    if (anaxis > 3) {
       printf("Error: images with > 3 dimensions are not supported\n");
       check = 0;
    }
	strcpy(outfilename,"s_");
        strcat(outfilename,argv[1]);
    remove(outfilename);
    if (check && !fits_create_file(&outfptr, outfilename, &status) )
    {
      /* copy all the header keywords from first image to output file */
      fits_copy_header(afptr, outfptr, &status);
      npixels = anaxes[0];  /* no. of pixels to read in each row */
      apix = (double *) malloc(npixels * sizeof(double)); /* mem for 1 row */
      if (apix == NULL) {
        printf("Memory allocation error\n");
        return(1);
      }
      fits_read_card(afptr,"SKYVAL", card, &status);
      fits_parse_value(card, oldvalue, comment, &status);
      strcpy(newcard, "SKYVAL");     /* copy keyword name */
      strcat(newcard, " = ");       /* '=' value delimiter */
      strcat(newcard, "0.0");     /* new value */
      if (*comment) {
         strcat(newcard, " / ");  /* comment delimiter */
         strcat(newcard, comment);     /* append the comment */
         }
      fits_parse_template(newcard, card, &keytype, &status);
      fits_update_card(outfptr, "SKYVAL", card, &status);
      /*printf("Old Value= %s Comment= %s \n", oldvalue,comment); */
      bkgd=atof(oldvalue);
      //*printf("%g\n",bkgd); */
      /* loop over all planes of the cube (2D images have 1 plane) */
      for (firstpix[2] = 1; firstpix[2] <= anaxes[2]; firstpix[2]++)
      {
        /* loop over all rows of the plane */
        for (firstpix[1] = 1; firstpix[1] <= anaxes[1]; firstpix[1]++)
        {
          if (fits_read_pix(afptr, TDOUBLE, firstpix, npixels, NULL, apix, NULL, &status)) break;
          for(ii=0; ii< npixels; ii++)  
                {
		if(!isnan(apix[ii])) apix[ii]=apix[ii]-bkgd;
		}
        fits_write_pix(outfptr, TDOUBLE, firstpix, npixels,apix, &status);
        }
      }
      fits_close_file(outfptr, &status);
      free(apix);
    }
    fits_close_file(afptr, &status);
    if (status) fits_report_error(stderr, status); /* print any error message */
    return(status);
}