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

/***************************************************/
/***     Invert a spectrum file  **/
/*      FFT format data files only   *****/
/***************************************************/

void main(int argc, char *argv[]) {
    FILE *file_in,*file_in2,*file_out;
    float *num, *den, *output;
    float temp, peaknum, peakout, magnum, magden, phase, norm = 1.0;
    long i,j;
    if (argc < 3)        {
        printf("useage: spectinv infile.fft outfile <norm> \n");
        printf("        where optional <norm> is a normalization factor \n");
        exit(0);
    }
    if (argc==4)    {
        norm = atof(argv[3]);
	printf("Norm = %f\n", norm);
    }
    file_in = fopen(argv[1],"rb");
    if (file_in)        {
	j = 0;
    	while (fread(&temp,4,1,file_in)) j++;
	num = (float *) malloc(4 * j);
	den = (float *) malloc(4 * j);
	file_out = fopen(argv[2],"wb");
	fseek(file_in, 0, 0);
	fread(den,4,j,file_in);
	output = (float *) malloc(4 * j);
	output[0] = 1.0 / den[0];
	output[1] = 1.0 / den[1];
	peaknum = 1.0;
        peakout = fabs(output[0]);
	for (i=2;i<j;i+=2)	{
	    magnum = 1;
	    if (magnum > peaknum) peaknum = magnum;
/*	    printf("%f ",magnum);	*/
	    magden = den[i]*den[i] + den[i+1]*den[i+1];
	    magden = sqrt(magden);
	    magnum /= magden;
	    if (magnum > peakout) peakout = magnum;
	    phase = atan2(num[i+1],num[i]);
	    phase -= atan2(den[i+1],den[i]);
/*	    printf("%f ",phase);	*/
	    output[i] = magnum * cos(phase);
	    output[i+1] = magnum * sin(phase);

/*	    printf("%f %f\n",magden,magnum);
*/
        }
	printf ("%f %f\n",peaknum, peakout);
	norm = norm * peaknum / peakout;
	printf("Norm = %f\n", norm);
	for (i=0;i<j;i+=2)	{
	    output[i] *= norm;
	    output[i+1] *= norm;
	    magnum = output[i]*output[i] + output[i+1]*output[i+1];
	    magnum = sqrt(magnum);
/*	    printf("%f %f %f \n",output[i],output[i+1],magnum);
*/
	}	   

	fwrite(output,4,j,file_out); 

	free(num);
	free(den);
	free(output);
        fclose(file_in);
        fclose(file_out);
    }
}