#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <fitsio.h>
#include <longnam.h>
#include <wcshdr.h>
#include <wcs.h>
#include <math.h>

main(int argc, char **argv){
  FILE *fp;
  char infile[256], inputlist[256], PRJ[4];
  char outfile[256], RAproj[9], DECproj[9];
  long outimagesizeX, outimagesizeY;
  fitsfile *fptr_in, *fptr_out;
  int status=0, ncards, relax, nreject, nwcs, ctrl, i, j, anynul;
  long naxes[2];
  float *inimage, *outimage, X_in, Y_in, X_out, Y_out;
  char *header;
  struct wcsprm *wcs_in,*wcs_out;
  double pixcrd_in[1][2], imgcrd[1][2], phi[1], theta[1], world[1][2], alpha, delta;
  double crpix1, crpix2, pixcrd_out[1][2];
  int i_out, j_out, nfile;
  long imagesize;
  double pixelsize;

  if(argc !=9){
    printf("usage: convertWCS inputfilelist PRJ outimagesizeX outimagesizeY crpix1 crpix2 pixelsize outfitsname\n");
    printf("where PRJ must be one of the following; CAR, AIT, ZEA, CSC, PCO.\n");
    exit(1);
  }else
    {
      strcpy(inputlist, argv[1]);
      strcpy(PRJ, argv[2]);
      sscanf(argv[3], "%ld", &outimagesizeX);
      sscanf(argv[4], "%ld", &outimagesizeY);
      sscanf(argv[5], "%lf",  &crpix1);
      sscanf(argv[6], "%lf",  &crpix2);
      sscanf(argv[7], "%lf",  &pixelsize);
      strcpy(outfile, argv[8]);
  }
  printf("input file list = %s\n", inputlist);
  if(strcmp(PRJ, "CAR")!=0 &&strcmp(PRJ, "AIT")!=0&&strcmp(PRJ, "ZEA")!=0&&strcmp(PRJ, "CSC")!=0&&strcmp(PRJ, "PCO")!=0){
    printf("PRJ must be one of the following: CAR, AIT, ZEA, CSC, PCO\n");
    exit(1);
  }

  /* Read the input image ascii list */
  if(NULL==(fp=fopen(inputlist,"r"))){
    printf("cannot open %s\n", inputlist);
    exit(1);
  }

  inimage=malloc(512*512*4);
  outimage=malloc(outimagesizeX*outimagesizeY*4);
  for(i=0;i<outimagesizeX;i++){
    for(j=0;j<outimagesizeY;j++){
      outimage[i*outimagesizeY+j]=0.0;
    }
  }

  /* put the output wcs keywords in the structure */
  relax=1; /* 0 means recognize only FITS keywords defined by the 
	      published WCS standard. 1 is more benign.*/
  ctrl =0; /* If 0, not report any rejected header cards.
	     ctrl=1,2,3 gives detailed error messages in this order */
 
  /* Create  the output FITS file */
  if (fits_create_file(&fptr_out, outfile, &status)) {
    fits_report_error(stderr, status);
    return 1;
  }
  naxes[0]=outimagesizeX;  naxes[1]=outimagesizeY;
  fits_create_img(fptr_out, FLOAT_IMG, 2L, naxes, &status);
  fits_write_key_str(fptr_out, "RADESYS", "FK5", "Coordinate system", &status);
  fits_write_key_dbl(fptr_out, "EQUINOX", 2000.0, 3, "Equinox",&status);
  sprintf(RAproj, "RA---%s", PRJ);
  fits_write_key_str(fptr_out, "CTYPE1", RAproj, "DEC projection",&status);
  sprintf(DECproj, "DEC--%s", PRJ);
  fits_write_key_str(fptr_out, "CTYPE2", DECproj, "DEC projection",&status);
  fits_write_key_dbl(fptr_out, "CRPIX1", crpix1, 13, "X reference pixel",&status);
  fits_write_key_dbl(fptr_out, "CRPIX2", crpix2, 13, "Y reference pixel",&status);
  fits_write_key_dbl(fptr_out, "CRVAL1", 0.0, 13, "RA of reference pixel",&status);
  fits_write_key_dbl(fptr_out, "CRVAL2", 0.0, 13, "DEC of reference pixel",&status);
  fits_write_key_dbl(fptr_out, "CDELT1", -pixelsize, 13, "X pixel increment (degree)",&status);
  fits_write_key_dbl(fptr_out, "CDELT2",  pixelsize, 13, "Y pixel increment (degree)",&status);

  /* read the output FITS header */
  if (fits_hdr2str(fptr_out, 1, NULL, 0, &header, &ncards, &status)) {
    fits_report_error(stderr, status);
    return 1;
  }

  /* put the output wcs keywords out the structure */
  if (status = wcspih(header, ncards, relax, ctrl, &nreject, &nwcs, &wcs_out)) {
    fprintf(stderr, "wcspih ERROR %d: %s.\n", status, wcs_errmsg[status]);
  }
  if (status = wcsset(wcs_out)) {
    fprintf(stderr, "wcsset ERROR %d: %s.\n", status,
	    wcs_errmsg[status]);
  }
  wcsprt(wcs_out); 
  free(header);

  nfile=0;
  while(fscanf(fp, "%s", infile)==1){
    nfile++;
    /* open the input FITS file */
    if (fits_open_file(&fptr_in, infile, READONLY, &status)) {
      fits_report_error(stderr, status);
      return 1;
    }
    printf("####### file %d: %s\n", nfile, infile);
    
    /* read the input FITS header */
    if (fits_hdr2str(fptr_in, 1, NULL, 0, &header, &ncards, &status)) {
      fits_report_error(stderr, status);
      return 1;
    }
    /* Put the input image into the image array */
    fits_read_2d_flt(fptr_in, 0, 0.0, 512, 512, 512, inimage, &anynul, &status);
    
    /* put the input wcs keywords in the structure */
    if (status = wcspih(header, ncards, relax, ctrl, &nreject, &nwcs, &wcs_in)) {
      fprintf(stderr, "wcspih ERROR %d: %s.\n", status, wcs_errmsg[status]);
    }
    free(header);
    
    if (status = wcsset(wcs_in)) {
      fprintf(stderr, "wcsset ERROR %d: %s.\n", status,
	      wcs_errmsg[status]);
    }
    /*    wcsprt(wcs_in); */
    /* Read the input image, 
       convert the original pixel coordinates to world coordinates*/
    for(i=0;i<512;i++){
      for(j=0;j<512;j++){
	/* Define the image coordinages (X,Y), where (1,1) 
	   is the lower left and X increases to left, Y increases 
	   to up, then X=j+1, Y=i+1. */
	X_in = (float)j+1.0; Y_in=(float)i+1.0;
	/* calculate the World Coordinates (right ascention and declination)
	   using wcslib. */
	pixcrd_in[0][0]=(double) X_in; pixcrd_in[0][1]=(double) Y_in;
	wcsp2s(wcs_in,1,2,pixcrd_in[0],imgcrd[0],phi,theta,world[0],&status);
	alpha = world[0][0]; delta=world[0][1];
	/* calculate the X,Y pixel values in the output file using wcslib. */
	wcss2p(wcs_out,1,2,world[0],phi,theta,imgcrd[0],pixcrd_out[0],&status);
	X_out=(float) pixcrd_out[0][0]; Y_out=(float) pixcrd_out[0][1];
	/*      printf("X=%f Y=%f alpha=%f dec=%f flux=%f status=%d Xout=%f Yout=%f\n", 
		X_in, Y_in, alpha, delta,inimage[i*512+j], status,X_out, Y_out);*/
	i_out = floorf(X_out)-1;j_out=floor(Y_out)-1;
	if(i_out>=0&&i_out<outimagesizeX&&j_out>=0&&j_out<outimagesizeY) {
	  outimage[i_out+j_out*outimagesizeX]=outimage[i_out+j_out*outimagesizeX] + inimage[i*512+j];
	  /*	printf("outimage %d %d %f\n",i_out, j_out, outimage[i_out*outimagesizeX+j_out]);*/
	}
      }
    }
    /* close the input file */
    fits_close_file(fptr_in, &status);
    wcsfree(wcs_in);
  }
  imagesize=outimagesizeX*outimagesizeY;
  fits_write_2d_flt(fptr_out, 0L, outimagesizeX, outimagesizeX, outimagesizeY, outimage, &status);

  /* close the output file*/
  fits_close_file(fptr_out, &status);            
  wcsfree(wcs_out);

}