/* Driver for routine fourfs */#include <stdio.h>#include <math.h>#define NRANSI#include "nr.h"#include "nrutil.h"#define SWAP(a,b) {iswap=(a);(a)=(b);(b)=iswap;}#define FSWAP(a,b) {flswap=(a);(a)=(b);(b)=flswap;}#define NX 8#define NY 32#define NZ 4#define NDAT (2*NX*NY*NZ)#define TEMPFILE1 "frfstmp1"#define TEMPFILE2 "frfstmp2"#define TEMPFILE3 "frfstmp3"#define TEMPFILE4 "frfstmp4"#if defined(MSDOS) || defined(_MSDOS) || defined(_MSDOS_) || defined(__MSDOS__)#define BINREADWRITE "wb+"#else#define BINREADWRITE "w+"#endifchar *fnames[4]={TEMPFILE1,TEMPFILE2,TEMPFILE3,TEMPFILE4};int main(void){	int cc,nnv=3,nwrite;	long idum=(-23);	unsigned long dim[4],i,j,k,l,ll,iswap;	float diff,smax,sum,sum1=0.0,sum2=0.0,tot,*data1,*data2,*data1p,*data2p;	FILE *flswap,*file[5];	data1=vector(1,NDAT);	data2=vector(1,NDAT);	dim[1]=NX;	dim[2]=NY;	dim[3]=NZ;	tot=(float)NX*(float)NY*(float)NZ;	for (j=1;j<=4;j++)		if ((file[j]=fopen(fnames[j-1],BINREADWRITE)) == NULL)			nrerror("Couldn't open temporary file");	for (i=1;i<=dim[3];i++)		for (j=1;j<=dim[2];j++)			for (k=1;k<=dim[1];k++) {				l=k+(j-1)*dim[1]+(i-1)*dim[2]*dim[1];				l=(l<<1)-1;				data2[l]=data1[l]=2*ran1(&idum)-1;				l++;				data2[l]=data1[l]=2*ran1(&idum)-1;			}	nwrite=NDAT >> 1;	cc=fwrite(&data1[1],sizeof(float),nwrite,file[1]);	if (cc != nwrite) nrerror("write error in xfourfs");	cc=fwrite(&data1[nwrite+1],sizeof(float),nwrite,file[2]);	if (cc != nwrite) nrerror("write error in xfourfs");	rewind(file[1]);	rewind(file[2]);	printf("**************** now doing fourfs *********\n");	fourfs(file,dim,nnv,1);	for (j=1;j<=4;j++) rewind(file[j]);	cc=fread(&data1[1],sizeof(float),nwrite,file[3]);	if (cc != nwrite) nrerror("read error in xfourfs");	cc=fread(&data1[nwrite+1],sizeof(float),nwrite,file[4]);	if (cc != nwrite) nrerror("read error in xfourfs");	printf("**************** now doing fourn *********\n");	fourn(data2,dim,nnv,1);	sum=smax=0.0;	for (i=1;i<=dim[3];i++)		for (j=1;j<=dim[2];j++)			for (k=1;k<=dim[1];k++) {				l=k+(j-1)*dim[1]+(i-1)*dim[2]*dim[1];				l=(l<<1)-1;				ll=i+(j-1)*dim[3]+(k-1)*dim[3]*dim[2];				ll=(ll<<1)-1;				diff=sqrt(SQR(data2[ll]-data1[l])+SQR(data2[ll+1]-data1[l+1]));				sum2 += SQR(data1[l])+SQR(data1[l+1]);				sum += diff;				if (diff > smax) smax=diff;			}	sum2=sqrt(sum2/tot);	sum=sum/tot;	printf("(r.m.s.) value, (max,ave) discrepancy= %12.7f %12.7f %12.7f\n",		sum2,smax,sum);	/* now check the inverse transforms */	SWAP(dim[1],dim[3]);	/* This step swap step is conceptually a reversal, but for	three dimensions a swap accomplishes that. */	FSWAP(file[1],file[3])	FSWAP(file[4],file[2])	for (j=1;j<=4;j++) rewind(file[j]);	printf("**************** now doing fourfs *********\n");	fourfs(file,dim,nnv,-1);	for (j=1;j<=4;j++) rewind(file[j]);	cc=fread(&data1[1],sizeof(float),nwrite,file[3]);	if (cc != nwrite) nrerror("read error in xfourfs");	cc=fread(&data1[nwrite+1],sizeof(float),nwrite,file[4]);	if (cc != nwrite) nrerror("read error in xfourfs");	SWAP(dim[1],dim[3]);	printf("**************** now doing fourn *********\n");	fourn(data2,dim,nnv,-1);	sum=smax=0.0;	data1p=data1;	data2p=data2;	for (j=1;j<=NDAT;j+=2) {		sum1 += SQR(data2p[1])+SQR(data2p[2]);		diff=sqrt(SQR(data2p[1]-data1p[1])+SQR(data2p[2]-data1p[2]));		sum += diff;		if (diff > smax) smax=diff;		data1p += 2;		data2p += 2;	}	sum=sum/tot;	sum1=sqrt(sum1/tot);	printf("(r.m.s.) value, (max,ave) discrepancy= %12.7f %12.7f %12.7f\n",		sum1,smax,sum);	printf("ratio of r.m.s. values, expected ratio= %12.6f %12.6f\n",		sum1/sum2,sqrt(tot));	for (j=1;j<=4;j++)		if (fclose(file[j]) == EOF) nrerror("Couldn't close temporary file");	free_vector(data2,1,NDAT);	free_vector(data1,1,NDAT);	for (j=1;j<=4;j++)		if (remove(fnames[j-1])) nrerror("Couldn't delete temporary file");	return 0;}#undef NRANSI