/* This is a template program for reading the 2dFgg database object fits files   */
/* For the filename provided on the command line, the DSS image is read into     */
/* an array, and APM and DSS data are read from primary header keywords          */
/* The number of spectral extensions ( observations ) are determined by          */
/* stepping though the FITS file. The user is prompted to select one of the      */
/* spectral extensions to read the data from. The object spectra, variance and   */
/* sky are loaded into arrays and simple operations performed on the data.       */
/* Data related to the spectra are read from keywords in the header of the       */
/* chosen spectral extension ( redshift, quality, observed ra & dec. The file    */
/* is closed and allocated memory released. Written by Ian Price, June 1999.     */ 
/*                                                                               */
/* To compile, the cfitsio include files "fitsio.h" and "longnam.h", and the     */
/* library libcfitsio.a are required. Typical compile line would be :            */
/*                                                                               */
/* cc -o 2dfgrs_example 2dfgrs_example.c -I/path_to_include_files -lcfitsio -lm  */
/*                                                                               */
/* Minor changes: more long names,main stuff in subroutine, mean stats, etc.     */
/*   Karl Glazebrook, June 1999.                                                 */
/* Minor updates: new user def keywords, Carole Jackson, Apr 2001.               */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "fitsio.h"  /* normally just "fitsio.h" */

#define BUFFSIZE 1024    /* buffer size when reading from fits file */

/* Declare C routines used in main() */

double calc_mean (float*, int);
void do_stuff_to_data (float*, float*, float*, int, double, double, double );

void main(int argc, char **argv)
{

  int status,anynull,header_unit;
  float *spec_data,*vari_data,*sky_data,sp_nullval,bjsel;
  float *dss_image, mean;
  fitsfile *fptr;
  char comment[FLEN_COMMENT];
  short dss_buff[BUFFSIZE],im_nullval;
  long im_naxis1,im_naxis2,sp_naxis1,sp_naxis2,nbuffer,fpixel;
  long npixel,serial,igal,quality;
  int i,num_spec;
  double ra,dec,cdelt1,crpix1,crval1,obsra,obsdec;
  float plate_scale,z,z_helio,snr,eta;

  if (argc==1)
    {
      printf("\n Usage : template 2dF_database_fits_file\n\n");
      exit(1);
    }

  /* FITS status flag is 0 until an error occurs. If not 0, the fits */
  /* routine will  be skipped over.                                  */

  status = 0;

  /* open the database object FITS file for reading */

  fits_open_file(&fptr,argv[1],READONLY,&status);

  if (status)
    {
      printf(" File NOT Found - %s\n",argv[1]);
      exit(1);
    }

/* Read APM and DSS image related keywords from primary FITS header        */

/* Available keywords are :-                                               */
/*                                                                         */
/* FITS descriptive keywords                                               */
/* SIMPLE   BITPIX   NAXIS    NAXIS1   NAXIS2   EXTEND   BSCALE   BZERO    */
/*                                                                         */
/* 2dFGRS/APM descriptive keywords				           */
/* SEQNUM   NAME	IMAGE	 RA	  DEC	   EQUINOX  BJSEL    PROB  */
/* PARK     PARMU	IGAL	 JON	  ORIENT   ECCENT   AREA     X_BJ  */ 
/* Y_BJ     DX          DY	 BJG	  RMAG     PMAG     FMAG     SMAG  */
/* IFIELD   IGFIELD     REGION   OBJEQNX  OBJRA    OBJDEC   PLTSCALE XPIXELSZ*/
/* YPIXELSZ OBJPLTX     OBJPLTY  DATAMAX  DATAMIN  BJSELOLD BJG_OLD        */
 
  /* get the DSS image size keywords */

  fits_read_key_lng(fptr,"NAXIS1",&im_naxis1,comment,&status);
  fits_read_key_lng(fptr,"NAXIS2",&im_naxis2,comment,&status);

  /* allocate memory to store the image */

  dss_image = (float *)calloc(im_naxis1*im_naxis2,sizeof(float));

  /* read the DSS image via a buffer - probably overkill for the small 2dFgg files */

  /* do not perform checking for null data values */
  im_nullval = 0.0;

  npixel = im_naxis1*im_naxis2;
  fpixel = 1;

  while (npixel>0)
    {
      nbuffer = npixel;
      if (nbuffer>BUFFSIZE) nbuffer = BUFFSIZE;

      /* read in a buffer full of data from the image array */
      
      fits_read_img(fptr,TSHORT,fpixel,nbuffer,&im_nullval,dss_buff,&anynull,&status);

      /* copy the buffer - useful for type casts ( ie. read float image, and store */
      /* as doubles                                                                */

      for (i=0;i<nbuffer;++i)
	{
	  *(dss_image+fpixel+i-1) = (float) dss_buff[i];
	}
      npixel -= nbuffer;
      fpixel += nbuffer;
    }
    
  mean = calc_mean(dss_image,   im_naxis1*im_naxis2);
  printf("Mean of DSS image = %f counts\n",mean);

  /* read APM keywords - BJSEL RA DEC SEQNUM IGAL */

  fits_read_key_dbl(fptr,"RA",&ra,comment,&status);
  fits_read_key_dbl(fptr,"DEC",&dec,comment,&status);
  fits_read_key_flt(fptr,"BJSEL",&bjsel,comment,&status);
  fits_read_key_lng(fptr,"SEQNUM",&serial,comment,&status);
  fits_read_key_lng(fptr,"IGAL",&igal,comment,&status);
  fits_read_key_flt(fptr,"PLTSCALE",&plate_scale,comment,&status);

  if (!status) 
	{
	  printf("\nAPM parameters:\n\n");
	  printf("BJSEL = %f\n",bjsel);
	  printf("RA = %f\n",ra);
	  printf("DEC = %f\n",dec);
	  printf("SEQNUM = %ld\n",serial);
	  printf("IGAL = %ld\n",igal);
	  printf("PLTSCL = %f\n\n",plate_scale);
	  
	}
      else
    {
      printf("Failed to read data and keywords from the primary image\n");
    }

  /* Finished with the primary header and DSS image - move on to the spectra */

  
  num_spec = 0;

  /* find out how many spectral extensions there are.   */
  /* step through the headers until status error occurs */

  while (!status)
    {
      fits_movrel_hdu(fptr,1,&header_unit,&status);
      if (!status) ++num_spec;
    }

  printf("%d spectral extensions found\n",num_spec);

  
  /* reset status to zero, and select spectral extension interegate */

  if (num_spec)
    {
      i = 0;
      while (i<1 || i>num_spec)
	{
	  printf("Which spectrum do you want to read from file (1..%d)\n",num_spec);
	  status = scanf("%d",&i);
	  if (status != 1) i = 0;
	}

      /* resetting the status flag */

      status = 0;

      /* move to the selected extension */

      fits_movabs_hdu(fptr,i+1,&header_unit,&status);
      
      if (status)
	printf("Failed to move to the chosen extension\n");


      /* read the spectra/variance/sky image dimension keywords */

      fits_read_key_lng(fptr,"NAXIS1",&sp_naxis1,comment,&status);
      fits_read_key_lng(fptr,"NAXIS2",&sp_naxis2,comment,&status);

      /* check that all components exist - could also check sp_naxis1==1024 */

      if (sp_naxis2 != 3)
	{
	  status = 1;
	  printf("fits file has invalid structure\n");
	}
	
      printf("Dimensions of spectrum array = %ld\n", sp_naxis1);

      /* allocate memory for the spectra, variance and sky */

      spec_data = (float *)calloc(sp_naxis1,sizeof(float));
      vari_data = (float *)calloc(sp_naxis1,sizeof(float));
      sky_data  = (float *)calloc(sp_naxis1,sizeof(float));

      /* read the wavelength scale keywords */

      /* wavelength defined as CRVAL1 + CDELT1*(pixel - CRPIX1), where pixel */
      /* starts at 1, and finishes at NAXIS1                                 */

      fits_read_key_dbl(fptr,"CDELT1",&cdelt1,comment,&status);
      fits_read_key_dbl(fptr,"CRVAL1",&crval1,comment,&status);
      fits_read_key_dbl(fptr,"CRPIX1",&crpix1,comment,&status);

      /* read the spectral keywords Z Z_HELIO OBSRA OBSDEC QUALITY */

      fits_read_key_flt(fptr,"Z",&z,comment,&status);
      fits_read_key_flt(fptr,"Z_HELIO",&z_helio,comment,&status);
      fits_read_key_dbl(fptr,"OBSRA",&obsra,comment,&status);
      fits_read_key_dbl(fptr,"OBSDEC",&obsdec,comment,&status);
      fits_read_key_lng(fptr,"QUALITY",&quality,comment,&status);

      if (!status)
	{
	  printf("\n\nSpectrum parameters:\n\n");
	  printf("Z = %f\n",z);
	  printf("Z_HELIO = %f\n",z_helio);
	  printf("OBSRA = %f\n",obsra);
	  printf("OBSDEC = %f\n",obsdec);
	  printf("QUALITY = %ld\n",quality);
	}
      else
	printf("Failed to read spectral keywords\n");

      /* read optional spectra keywords SNR and ETA */

      fits_read_key_flt(fptr,"SNR",&snr,comment,&status);
      if (!status)
	{
	  printf("SNR = %f\n",snr);
        }

      fits_read_key_flt(fptr,"ETA",&eta,comment,&status);
      if (!status)
	{
	  printf("ETA = %f\n",eta);
        }

      /* load the data without buffering, and replace null data values with */
      /* -1.0e20                                                            */

      sp_nullval = -1.0e20;

      fits_read_img(fptr,TFLOAT,1,sp_naxis1,&sp_nullval,spec_data,&anynull,&status);
      fits_read_img(fptr,TFLOAT,sp_naxis1+1,sp_naxis1,&sp_nullval,vari_data,&anynull,&status);
      fits_read_img(fptr,TFLOAT,2*sp_naxis1+1,sp_naxis1,&sp_nullval,sky_data,&anynull,&status);
      
      /* Do some stuff with the data */
      
      if (!status) {
         do_stuff_to_data (spec_data, vari_data, sky_data, sp_naxis1, crval1, crpix1, cdelt1);
      } 
 
      /* reclaim the allocated memory from the spectra, variance and sky */

      free(spec_data);
      free(vari_data);
      free(sky_data);
    }

  /* reset the status flag before closing to ensure it closes */

  status = 0;
  
  fits_close_file(fptr,&status);

  /* reclaim the memory allocated for the DSS image */

  free(dss_image);

}
 
/* Trvial subroutine to calculate the mean of an array */
 
double calc_mean (float *a, int nelem) {
   int i;
   double sum;
   sum = 0.0;
   for (i=0; i<nelem; i++)  {
      sum += a[i];
   }
   return (sum/nelem);
}

void do_stuff_to_data (float *data, float*var, float *sky, int nelem,
                       double crval1, double crpix1, double  cdelt1) {

   
     /* Do some more stuff with the data e.g. */
      /* perform operations on the data - add 42 to the spectra */
      /* take square root of variance                           */
      /* turn sky counts into flux in arbitrary units           */
      /* Work out mean of good data */
      
      int i,n;
      float lambda;
      double sum;

      n=0; sum=0.0;
      for (i=0;i<nelem;++i)
     	{
	  /* Wavelength of each pixel */
	
     	  lambda = crval1 + cdelt1*( (double)(i+1) - crpix1 );
	  if (i==0) 
	     printf("Wavelength min = %7.2f A\n", lambda);
	  if (i==nelem-1) 
	     printf("Wavelength max = %7.2f A\n", lambda);
	     

     	  if ( var[i]>0.0 ) {  /* Good data has positive variance */
	      var[i] = sqrt( var[i] );
	      sum += data[i];
	      n++;
	  }
     	  else {
	      var[i] = -1.0;
          }
	  
     	  data[i] += 42.0; /* Add 42 to data (example of modification) */
     	  sky[i] /= 0.01 * lambda;
     	}
	
	printf("Mean of good data = %f units\n", sum/n);

}