#include <stdio.h>
#include <math.h>

#define Q315
/*#define Q270*/

#ifdef Q315
#define NPIXELS 384
#endif

#ifdef Q270
#define NPIXELS 260
#endif

/* Read a NPIXELSxNPIXELS pixel correction file and an XDS_ASCII.HKL file and
 * then correct the XDS_ASCII.HKL file and write a new one.
 */

int idata[NPIXELS*NPIXELS];

/*#define NX 6144 Number of pixels in CCD image */
/*#define NX 3072 Number of pixels in CCD image */

main(argc, argv)
int argc;
char *argv[];
{
	int NX;
	char str[256];
	FILE *fpin;

	if (argc != 4) {
		fprintf(stderr,"Usage:\n\n %s 3072 correction.raw XDS_ASCII.HKL > XDS_ASCII_cor.HKL\n\n",argv[0]);
		fprintf(stderr,"or for unbinned data:\n\n %s 6144 correction.raw XDS_ASCII.HKL > XDS_ASCII_cor.HKL\n",argv[0]);
		exit(-1);
	}

	NX = atoi(argv[1]);
	/*
	if ((NX != 3072) && (NX != 6144)) {
		fprintf(stderr,"Error: Number of pixels must be 3072 or 6144 (now %d)\n",NX);
		exit(-1);
	}
	*/
	fprintf(stderr,"Correcting data from %dx%d pixel images.\n",NX,NX); fflush(stderr);

	if ((fpin = fopen(argv[2],"r"))==NULL) {
		fprintf(stderr,"error: can't open %s\n",argv[2]);
		exit(-1);
	}

	if (fread(idata, sizeof(int), NPIXELS*NPIXELS, fpin) != NPIXELS*NPIXELS) {
		fprintf(stderr,"error: reading %s\n",argv[2]);
		exit(-1);
	}
	fclose(fpin);

	if ((fpin = fopen(argv[3],"r"))==NULL) {
		fprintf(stderr,"error: can't open %s\n",argv[3]);
		exit(-1);
	}

	while (fgets(str, sizeof(str), fpin) != NULL) {
		if (!strncmp(str,"!END_OF_HEADER",strlen("!END_OF_HEADER"))) {
			fprintf(stdout,"!MODIFIED BY XDSFIX\n"); 
			fprintf(stdout,"%s",str); 
			break;
		}
		fprintf(stdout,"%s",str); 
	}

	while (fgets(str, sizeof(str), fpin) != NULL) {
		int ix, iy, ih, ik, il;
		char junk[256];
		double iobs,isig,xpos,ypos,dcor;
		double interp_int();

		if (!strncmp(str,"!END_OF_DATA",strlen("!END_OF_DATA"))) {
			fprintf(stdout,"%s",str); 
			break;
		}
		sscanf(str,"%d%d%d%lf%lf%lf%lf%[^\n]",&ih,&ik,&il,&iobs,&isig,&xpos,&ypos,junk);

		if (xpos < 0) xpos = 0;
		if (xpos >= NX) xpos = NX-1;

		if (ypos < 0) ypos = 0;
		if (ypos >= NX) ypos = NX-1;

		dcor = interp_int(xpos*NPIXELS/NX,ypos*NPIXELS/NX,NPIXELS,NPIXELS,idata);
		if (dcor == 0.0)
			dcor = 1000.0;

		iobs = iobs * 1000.0/dcor;
		isig = isig * 1000.0/dcor;

		/*
		ix = xpos * (double)NPIXELS/NX + 0.5;
		iy = ypos * (double)NPIXELS/NX + 0.5;

		if (ix < 0) ix = 0;
		if (ix >= NPIXELS) ix = NPIXELS-1;

		if (iy < 0) iy = 0;
		if (iy >= NPIXELS) iy = NPIXELS-1;

		iobs = iobs * 1000.0/idata[iy*NPIXELS + ix];
		isig = isig * 1000.0/idata[iy*NPIXELS + ix];
		*/

		fprintf(stdout,"%6d%6d%6d%11.3e%11.3e%8.1f%8.1f%s\n",ih,ik,il,iobs,isig,xpos,ypos,junk);
	}
}

/* Return interpolated value of f[x][y] using
 * bilinear interpolation. From Numerical Recipes.
 *
 * The data to be interpolated is of type integer and stored in
 * an array of size width x height. The value is interpolated at
 * the point x,y. If x,y is not in the array then 0 is returned.
 */

#define Z1	(f[j*width + i])
#define Z2	(f[(j+1)*width + i])
#define Z3	(f[(j+1)*width + i+1])
#define Z4	(f[j*width + i+1])

double
interp_int(x,y,width,height,f)
register double x, y;
register int width,height;
register int *f;
{
	register int i=x, j=y;
	register double u, t;

	if ((i < 0) || (i >= width) || (j < 0) || (j >= height))
		return(0.0);

	t = (y - j);
	u = (x - i);

	if (i == width-1) {
		if (j == height-1)
			return((double)f[j*width + i]);
		else
			return((double)(1.0-t)*f[j*width + i] + t*f[(j+1)*width + i]);
	}
	else
	if (j == height-1)
		return((double)(1.0-u)*f[j*width + i] + u*f[j*width + i+1]);
	
	return((1.0-t)*(1.0-u)*Z1 + t*(1.0-u)*Z2 + t*u*Z3 + (1.0-t)*u*Z4);
}