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

/*
  Every program which uses the CFITSIO interface must include the
  the fitsio.h header file.  This contains the prototypes for all
  the routines and defines the error status values and other symbolic
  constants used in the interface.  
*/
#include "fitsio.h"

int main( int argc, char *argv[] );

void read_write_image( char filename[], int naxis1, int naxis2, char value[], int hdunum );

void printerror( int status);

int main(int argc, char *argv[])
{
/*************************************************************************
   This is a simple main program that calls the following routines:

   
    readimage     - read a FITS image and compute the min and max value
    
**************************************************************************/
   int naxis1, naxis2, hdunum ;
   if(argc != 6)
    {
      printf("\n This Program requires 5 command line arguments..\n") ;
      exit(0) ;
    }
   naxis1 = (int)atoi(argv[2]) ;
   naxis2 = (int)atoi(argv[3] ) ;
   hdunum = (int)atoi(argv[5]) ;
    read_write_image(argv[1], naxis1, naxis2, argv[4], hdunum);
    printf("\n\n");
    return(0);
}

/*--------------------------------------------------------------------------*/
void read_write_image(char filename[], int naxis1, int naxis2, char value[], int hdunum )

    /************************************************************************/
    /* Read a FITS image and determine the minimum and maximum pixel values */
    /************************************************************************/
{
    fitsfile *fptr,*fptr1;       /* pointer to the FITS file, defined in fitsio.h */
    #define buffsize 1000
    int status,  nfound, anynull, buffer2[buffsize], bits_per_pixel, new_val2, hdutype ;
    long naxes[2], fpixel, nbuffer, npixels, ii, new_loc, buffer3[buffsize], new_val3;
    float  nullval, buffer[buffsize], new_val;
    double buffer1[buffsize], new_val1 ;
    char *comment ;
    status = 0;

    if ( fits_open_file(&fptr, filename, READWRITE, &status) )
         printerror( status );

    if ( fits_movabs_hdu(fptr, hdunum, &hdutype, &status) ) ;
           printerror( status );

    if ( fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status) )
         printerror( status );
    
    fits_read_key(fptr, TINT, "BITPIX",&bits_per_pixel,comment, &status ) ;

    npixels  = naxes[0] * naxes[1];         /* number of pixels in the image */
    fpixel   = 1;
    nullval  = 0;                /* don't check for null values in the image */
  
    new_loc = (naxis2 * naxes[0]) + naxis1 ;

    if(bits_per_pixel == -32)
      new_val = (float)atof(value) ;
    
    
    if(bits_per_pixel == -64)
      new_val1 = (double)atof(value) ;

    
    if(bits_per_pixel == 16)
      new_val2 = (int)atoi(value) ;

    
    if(bits_per_pixel == 32)
      new_val3 = (long)atoi(value) ;
    

 while (npixels > 0)
   {
      nbuffer = npixels;
      if (npixels > buffsize)
        nbuffer = buffsize;     /* read as many pixels as will fit in buffer */

      
      /* If the image has got floating point values.*/                      
     if(bits_per_pixel == -32)
     {
      if ( fits_read_img(fptr, TFLOAT, fpixel, nbuffer, &nullval,
                  buffer, &anynull, &status) )
           printerror( status );

      if ( new_loc <= (fpixel + nbuffer)) /* Change the value of the proper location. */
	buffer[new_loc] = new_val ; 

       if ( fits_write_img(fptr, TFLOAT, fpixel, nbuffer,buffer, &status) )
         printerror( status );
      
     }
 
     if(bits_per_pixel == -64)
     {
      if ( fits_read_img(fptr, TDOUBLE, fpixel, nbuffer, &nullval,
                  buffer1, &anynull, &status) )
           printerror( status );

      if (  new_loc <= (fpixel + nbuffer)) /* Change the value of the proper location. */
	buffer1[new_loc] = new_val1 ; 

       if ( fits_write_img(fptr, TDOUBLE, fpixel, nbuffer,buffer1, &status) )
         printerror( status );
     
     }
 
     if(bits_per_pixel == 16)
     {
      if ( fits_read_img(fptr, TINT, fpixel, nbuffer, &nullval,
                  buffer2, &anynull, &status) )
           printerror( status );

      if (new_loc <= (fpixel + nbuffer)) /* Change the value of the proper location. */
	buffer2[new_loc] = new_val2 ; 

       if ( fits_write_img(fptr, TINT, fpixel, nbuffer,buffer2, &status) )
         printerror( status );
    
     
     }

    if(bits_per_pixel == 32)
     {
      if ( fits_read_img(fptr, TLONG, fpixel, nbuffer, &nullval,
                  buffer3, &anynull, &status) )
           printerror( status );

      if (  new_loc <= (fpixel + nbuffer)) /* Change the value of the proper location. */
	buffer3[new_loc] = new_val3 ; 

       if ( fits_write_img(fptr, TLONG, fpixel, nbuffer,buffer3, &status) )
         printerror( status );
     
    }
      npixels -= nbuffer;    /* increment remaining number of pixels */
      fpixel  += nbuffer;    /* next pixel to be read in image */
 }

  

    if ( fits_close_file(fptr, &status) )
         printerror( status );
   
    
    return;
}

/*--------------------------------------------------------------------------*/
void printerror( int status)
{
    /*****************************************************/
    /* Print out cfitsio error messages and exit program */
    /*****************************************************/


    if (status)
    {
       fits_report_error(stderr, status); /* print error report */
       exit( status );    /* terminate the program, returning error status */
    }
    return;
}