#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <conio.h>

#define MMAX 300

double x[1025] , y[1025] , fpe[MMAX] , r[MMAX] , rr[MMAX] , rfpe[MMAX];
double rm , wnmin , wnmax , sum1 , sum2 , sum , sumn , sumd ,wnint;
int    isw , nd , mmax , isw , ns , minm;
double dt , pmm , pm , fpemin , av;
double ci = 2*3.14159;


void gs1880(void)
{
	int i;
	
	minm = mmax;
	pmm = pm;
	for(i = 1 ; i <= mmax ; i++)
	{
		rfpe[i] = r[i];
	}
}

void gs1800(int m)
{
	int i;
	
	fpe[m] = (double)(nd+m+1) / (double)(nd-m-1) * pm;
	if(fpe[m] > fpemin) return;
	fpemin = fpe[m];
	minm = m;
	pmm = pm;
	for(i = 1 ; i <= m ; i++)
	{
		rfpe[i] = r[i];
	}
}

void gs1340(void)
{
	char nm[256];
	FILE *fp;
	int i;
	double ccc , bbb;
	
	printf(" ***   ÃÞ°À̧²ÙÒ²       =");
	scanf("%s" , nm);
	if((fp = fopen(nm , "r")) == NULL)
	{
		perror(nm);
		exit(0);
	}
	for(i = 1 ; i <= nd ; i++)
	{
		fscanf(fp , "%lf %lf %lf" , &(x[i]) , &bbb , &ccc);
		printf("I=%d  X=%g\n" , i , x[i]);
	}
	fclose(fp);
}

void gs1420(void)
{
	double  z;
	int i , m;
	
	sum = 0.0;
	for(i = 0 ; i <= nd ; i++) sum += x[i];
	av = sum / (double)nd;
	sum = 0.0;
	for(i = 0 ; i <= nd ; i++)
	{
		z = x[i] - av;
		x[i] = z;
		y[i-1] = z;
		sum += z*z;
	}
	pm = sum / (double)nd;
	fpemin = (double)(nd + 1) / (double)(nd - 1) * pm;
	fpe[0] = fpemin;
	for(m = 1 ; m <= mmax ; m++)
	{
		sumn = 0.0;
		sumd = 0.0;
		for(i = 1 ; i <= nd-m ; i++)
		{
			sumn += x[i] * y[i];
			sumd += x[i] * x[i] + y[i] * y[i];
		}
		rm = -2.0 * sumn / sumd;
		r[m] = rm;
		pm *= 1.0 - rm * rm;
		if(m == 1) goto label1680;
		for(i = 1 ; i <= m-1 ; i++)
		{
			r[i] = rr[i] + rm * rr[m-i];
		}
label1680:;
		for(i = 1 ; i <= m ; i++)
		{
			rr[i] = r[i];
		}
		for(i = 1 ; i <= nd-m-1 ; i++)
		{
			x[i] = x[i] + rm * y[i];
			y[i] = y[i+1] + rm * x[i+1];
		}
		if(isw == 0) gs1800(m);
	}
	if(isw == 1) gs1880();
}

void gs1960(void)
{
	double f;
	int i , j;
	
	ci *= dt;
	i = 0;
	for(f = wnmin ; f <= wnmax ; f += wnint)
	{
		sum1 = 1.0;
		sum2 = 0.0;
		for(j = 1 ; j <= minm ; j++)
		{
			sum1 += rfpe[j] * cos(ci * f * j);
			sum2 += rfpe[j] * sin(ci * f * j);
		}
		x[i] = pmm * dt / (sum1 * sum1 + sum2 * sum2);
		i = i + 1;
	}
}


void gs2110(void)
{
	char nm[256];
	FILE *fp;
	int i;
	
	printf("  Out Data File Name ");
	scanf("%s" , nm);
	if((fp = fopen(nm , "w")) == NULL)
	{
		perror(nm);
		exit(0);
	}
	for(i = 1 ; i <= ns ; i++)
	{
		fprintf(fp , "%g\n" , x[i]);
	}
	fclose(fp);
}

int mem(void)
{
	char an , p[256] , n[256];
	
	
	printf("      ***** MEM *****\n\n");
	printf(" ***  »ÝÌßÙ Ãݽ³       =");
	scanf("%d" , &nd);
	printf("nd:%d\n" ,nd);
	printf(" ***   »ÝÌßÙ ¶Ý¶¸       =");
	scanf("%lf" , &dt);
	printf("dt:%g\n" , dt);
label1130:;
	printf(" ÓÃÞÙ ¼Þ½³ ¦ ±À´Ï½¶ (Y/N)\n");
	an = getch();
	if(an == 'Y' || an =='y')
	{
		isw = 1;
		strcpy(p , " ***   ÓÃÞÙ ¼Þ½³           =");
	}
	else if(an == 'N' || an == 'n')
	{
		isw = 0;
		strcpy(p , " ***   »²ÀÞ² ÓÃÞÙ ¼Þ½³     =");
	}
	else goto label1130;
	printf("%s" , p);
	scanf("%d" , &mmax);
	printf("mmax:%d\n" ,mmax);
	gs1340();
	ns = nd;
	wnmin = 0.0;
	wnmax = nd;
	strcpy(n , " Data ");
/*	gs2110();*/
	gs1420();
label1220:;
	printf("   »²¼®³ ¼­³Ê½³ =");
	scanf("%lf" , &wnmin);
	printf("   »²ÀÞ² ¼­³Ê½³ =");
	scanf("%lf" , &wnmax);
	printf("wnmin,max:%g  ,  %g\n" , wnmin , wnmax);
	printf("   ½Í߸ÄÙ Ãݽ³  =");
	scanf("%d" , &ns);
	wnint = (wnmax - wnmin) / (double)ns;
	printf("ns:%d   wnint:%g\n" , ns , wnint);
	gs1960();
	strcpy(n , " MEM Spectrum ");
	gs2110();
label1270:;
	printf("  Ó³²ÁÄÞ ¼Ï½¶ (Y/N)");
	an = getch();
	if(an == 'Y' || an == 'y') goto label1220;
	else if(an != 'N' && an != 'n') goto label1270;
	return(1);
}


void main(void)
{
	mem();
}