/* imgfftio.c -- routines fimgr and fimgw called by filter_img.f
	W. H. F. Smith, 24 June 1996
	filename string must be fixed with '\0' before call.
	fimgw opens for append, so main program should remove old
		file before starting a new run.

	Modified 27 June 1996 by W H F Smith to add remove and restore
		of mean values features.
*/

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

void fimgr_(savex, nximg, nyfft, jrowseek, ikillm, ifilnam, len)
float	*savex;
int	*nximg, *nyfft, *jrowseek, *ikillm, *len;
char	*ifilnam;
{
	int	ntot, jseek, k;
	FILE	*fp = NULL;
	short int	dummy;
	double 	sum;

	sum += 0.0;

	if ((fp = fopen(ifilnam, "r")) == NULL) {
		fprintf(stderr,"fimgr ERROR cannot open r %s\n", ifilnam);
		exit(-1);
	}

	jseek = *jrowseek;
	ntot = (*nximg) * (*nyfft);

	if (jseek && fseek(fp, jseek*(*nximg)*2, 0) ) {
		fprintf(stderr,"fimgr ERROR seeking %d rows in %s\n", jseek,
					ifilnam);
		exit(-1);
	}

	for (k = 0; k < ntot; k++) {
		if ( (fread((char *)&dummy, 2, 1, fp)) != 1) {
			fprintf(stderr,"fimgr ERROR reading from %s\n",
				ifilnam);
			exit(-1);
		}
		savex[k] = (float)dummy;
		sum += (double)dummy;
	}
	
	fclose(fp);

	if (*ikillm) {
		sum /= ntot;
		for (k = 0; k < ntot; k++) savex[k] -= sum;
	}

	return;
}

void fimgw_(blend, nximg, nyfft2, iaddx, x, ofilnam, len)
float	*blend;
int	*nximg, *nyfft2, *iaddx, *len;
double	*x;
char	*ofilnam;
{
	int	ntot, k, idummy, nint();
	FILE	*fp = NULL;
	short int	dummy;

	if ((fp = fopen(ofilnam, "a")) == NULL) {
		fprintf(stderr,"fimgw ERROR cannot open a %s\n", ofilnam);
		exit(-1);
	}

	ntot = (*nximg) * (*nyfft2);

	if (*iaddx) {
		for (k = 0; k < ntot; k++) {
			idummy = nint((*x) + (double)blend[k]);
			if (abs(idummy) > 32767) {
				fprintf(stderr,"fimgw ERROR data exceeds 32767\n");
			}
			dummy = (short)idummy;
			if ( (fwrite((char *)&dummy, 2, 1, fp)) != 1) {
				fprintf(stderr,"fimgw ERROR writing to %s\n",
					ofilnam);
				exit(-1);
			}
		}
	}
	else {
		for (k = 0; k < ntot; k++) {
			idummy = nint((double)blend[k]);
			if (abs(idummy) > 32767) {
				fprintf(stderr,"fimgw ERROR data exceeds 32767\n");
				dummy = (idummy < 0) ? -32767 : 32767;
			}
			else {
				dummy = (short)idummy;
			}
			if ( (fwrite((char *)&dummy, 2, 1, fp)) != 1) {
				fprintf(stderr,"fimgw ERROR writing to %s\n",
					ofilnam);
				exit(-1);
			}
		}
	}
	
	fclose(fp);
	return;
}


void findgm_(gmean, nximg, nyimg, ifilnam, len)
double	*gmean;
int	*nximg, *nyimg, *len;
char	*ifilnam;
{
	int	ntot, k;
	FILE	*fp = NULL;
	short int	dummy;
	double	sum = 0.0;

	if ((fp = fopen(ifilnam, "r")) == NULL) {
		fprintf(stderr,"findgm ERROR cannot open r %s\n", ifilnam);
		exit(-1);
	}

	ntot = (*nximg) * (*nyimg);

	for (k = 0; k < ntot; k++) {
		if ( (fread((char *)&dummy, 2, 1, fp)) != 1) {
			fprintf(stderr,"findgm ERROR reading from %s\n",
				ifilnam);
			exit(-1);
		}
		sum += (double)dummy;
	}
	
	fclose(fp);

	*gmean = sum/ntot;
	return;
}