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

/***************************************************/
/***          Divide one spectrum file by another **/
/*      FFT format data files only   *****/
/***************************************************/

void main(int argc, char *argv[]) {
    FILE *files_in[16],*file_out;
    float *data[16], *output;
    float temp, mag, phase;
    float PI = 4.*atan( 1. ) ;
    long i, j, length, num_files;
    
    if (argc < 3)        {
        printf("useage: average infile.fft <infile2.fft ...> outfile.fft\n");
        exit(0);
    }
    files_in[0] = fopen(argv[1],"rb");
    if (files_in[0])   {
        j = 0;
        while (fread(&temp,4,1,files_in[0])) j++;
    }
    length = j;
    fseek(files_in[0], 0, 0);
    data[0] = (float *) malloc(4 * length);
    fread(data[0],4,i,files_in[0]);
    for (i=2;i<argc-1;i++)    {
        files_in[i-1] = fopen(argv[i],"rb");
	if (files_in[i-1])  {
    	    j = 0;
    	    while (fread(&temp,4,1,files_in[i-1])) j++;
	    if (j != length) {
		printf("Files Not Same Length!!\n");
		exit(0);
	    }
	    fseek(files_in[i-1], 0, 0);
	    data[i-1] = (float *) malloc(4 * length);
	    fread(data[i-1],4,length,files_in[i-1]);
	}
	else	{
	    printf("File Problems!!\n");
	    exit(0);
	}
    }
    num_files = argc-2;
    printf("%i\n", num_files);
    file_out = fopen(argv[argc-1],"wb");
    output = (float *) malloc(4 * length);
    output[0] = 0.0;
    output[1] = 0.0;
    for (i=0;i<num_files;i++)	{
	output[0] += data[i][0];
	output[1] += data[i][1];
    }
    output[0] /= num_files;
    output[1] /= num_files;
    for (j=2;j<length;j+=2) {
        output[j] = 0.0;
        output[j+1] = 0.0;
	for (i=0;i<num_files;i++)   {
	    mag = data[i][j]*data[i][j] 
		    + data[i][j+1]*data[i][j+1];
	    mag = sqrt(mag);
	    phase = atan2(data[i][j+1],data[i][j]);
	    output[j] += mag;
	    output[j+i] += phase;
	}
        mag = output[j] / num_files;
	phase = output[j+1];
	while (phase > PI) phase -= PI;
	while (phase < -PI) phase += PI;	
	output[i] = mag * cos(phase);
	output[i+1] = mag * sin(phase);
    }
    fwrite(output,4,length,file_out); 
    
    for (i=0;i<num_files;i++)    {
	free(data);
        fclose(files_in[i]);
    }
    free(output);
    fclose(file_out);
}