www.pudn.com > hilert.rar > hilert.cpp

/*   main.c  
/******************************************/  
   
#include   
#include   
   
#define N 500  
#define PI 3.14159265  
   
#define sq(X) ((X)*(X))  
 
void hilbert(int, double[], double[]);    
 
int     main()  
{ 
	FILE *fp1,*fp2,*fp3,*fp4; 
	fp1=fopen("input.txt","w+"); 
	fp2=fopen("output.txt","w+"); 
	fp3=fopen("phase.txt","w+"); 
	fp4=fopen("angel.txt","w+"); 
 
        int             i;  
        double          x;  
        double          delta[N];  
        double          kappa[N];  
        double          y;  
        double          xmin = -10.;  
        double          xmax = 10.;  
        double          cd = -1.;  
        double          w = 2.;  
        double          h = (xmax - xmin) / N;  
   
        for (i = 0; i < N; i++)  
        {  
                x = 2. * (xmin + i * h - cd) / w;  
                if (x < -1.)  
                        delta[i] = 0.;  
                else if (x < 1.)  
                        delta[i] = sqrt(1. - sq(x));  
                else  
                        delta[i] = 0.;  
//				fprintf(fp1,"%.3f  ",delta[i]); 
        }  
		for (i = 0; i < N; i++) 
		{ 
			delta[i]=cos(4.*PI*i/N); 
			fprintf(fp1,"%.3f  ",delta[i]); 
		} 
   
        hilbert(N, delta, kappa);  
   
        for (i = 0; i < N; i++)  
        { 
                x = 2. * (xmin + i * h - cd) / w;  
                if (x < -1.)  
                        y = x + sqrt(sq(x) - 1.);  
                else if (x < 1.)  
                        y = x;  
                else  
                        y = x - sqrt(sq(x) - 1.);  
 
				fprintf(fp2,"%.3f  ",kappa[i]); 
				fprintf(fp3,"%.3f  ",xmin + i * h); 
				fprintf(fp4,"%.3f  ",y); 
 
                (void) printf("%.3f %.3f %.3f %.3f\n", xmin + i * h, delta[i], kappa[i], y);  
        }  
   
        return 0;  
 
		fclose(fp1); 
		fclose(fp2); 
		fclose(fp3); 
		fclose(fp4); 
}  
   
/******************************************/  
/*   hilbert.c  
/******************************************/  
   
 
   
/******************************************/  
/*   hilbert.c  
/******************************************/  
 
void hilbert(int n, double delta[], double kappa[])  
{  
        double d1, d2, d3, d4;  
        int i1, i2;  
   
        for (i1 = 0; i1 < n; i1++)  
        {  
                kappa[i1] = 0.;  
                for (i2 = 1; i2 < n; i2++)  
                {  
                        d1 = (i1+i2=0)? delta[i1-i2]: 0.;  
                        d3 = (i1+i2+1=0)? delta[i1-i2-1]: 0.;  
   
                        kappa[i1] -= 0.5 * (d1-d2) / i2 + 0.5 * (d3 - d4) / (i2+1);  
                }  
                kappa[i1] /= PI;  
        }  
}