fresnel.c

Fresnel带的验证程序

#include "stdafx.h"
#include <stdio.h>
#include <math.h>
#define max 10

double m=5,l[max],pi=3.1416926,r=1000,R=1000;
static double rd[max] = { 1000.0 };
static double Rd[max] = { 1000.0 };
double K;

double f1(double x)/*f1(x)是被积函数的实部*/
{ 
  K=2*pi/m;

  return x*R*K*sin(K*(sqrt(r*r+x*x)+sqrt(R*R+x*x)))/(sqrt(r*r+x*x)*(R*R+x*x));
}
double f2(double x)/*f2(x)是被积函数的虚部*/
{ 
  K=2*pi/m;

  return x*R*K*cos(K*(sqrt(r*r+x*x)+sqrt(R*R+x*x)))/(sqrt(r*r+x*x)*(R*R+x*x));
}

double f3(double x)/*f3(x)是求带Fresnel半径的函数*/
{   
  return sqrt(r*r+x*x)+sqrt(R*R+x*x)-r-R;
}

void  fresnel_zones_radii ( double l[max] ) /*求Fresnel带的半径*/
{  int i,j;
   double x1,x2,midpoint;
	for ( i=1;i<max;i++) {
		x1=0;
		x2 =sqrt((0.5*i*m+r)*(0.5*i*m+r)-r*r);
        midpoint=0.5*(x1+x2);

		for(j=1;j<100000;j++,midpoint=0.5*(x1+x2)) {
		   if ( f3(midpoint)-0.5*i*m > 0 ) 
             x2=midpoint;
		   else x1=midpoint;
		   if ( fabs(f3(midpoint)-0.5*i*m) <= 0.0000001 ) { l[i] = midpoint; break;}

		}
	rd[i]=sqrt(rd[0]*rd[0]+l[i]*l[i]);
    Rd[i]=sqrt(Rd[0]*Rd[0]+l[i]*l[i]);

	printf("rd[%5d]=%12.6f",i,rd[i]);
    printf("   Rd[%5d]=%12.6f",i,Rd[i]);
    printf("   l[%5d]=%12.6f\n",i,l[i]);
	}
}

double Simp(double a,double b,double (*fun)())
/*形参b,a是积分的上下限，函数指针fun指向被积函数*/
{
  #define NUM 10000  /*等分积分区间*/
  double Simp=0,dx=(b-a)/NUM;
  int i;
  double x1=a,x2=a+dx;
  for(i=0;i<NUM;i++,x1=x2,x2+=dx)
  Simp += ((*fun)(x1)+4*(*fun)((x1+x2)/2)+(*fun)(x2))/6*dx;
  return Simp;
}

void write(double *array,int n,char name[100])
{
  FILE *fp1;
  double *fr=array;
  int i;
  char path[100];
  printf("\n");
  printf("Please input the file name to save the data:\n");
  printf("%s\n",name);
  scanf("%s",path);
  if((fp1=fopen(path,"rb+"))==NULL)
  {
  printf("can't open the file\n");
  exit(0);
  }
  for(i=0;i<n;i++)
	{
	   fprintf(fp1,"%f",fr[i]);
	   fprintf(fp1,"\r\n");
	}
	fclose(fp1);
	printf("\n");
for(i=0;i<n;i++)
{
  printf("%f\t",fr[i]);
}
printf("\n");
}

main()
{ int i,j;
  FILE *fp1,*fp2,*fp3;
  double a1, b1, s1,a[max],b[max],s[max],accuma[max],accumb[max],accums[max];
  
  if((fp1=fopen("Fresnel_zones_each.txt","w+"))==NULL)
  {
  printf("can't open the file\n");
  exit(0);
  }

  if((fp2=fopen("Fresnel_zones_sum.txt","w+"))==NULL)
  {
  printf("can't open the file\n");
  exit(0);
  }
  
  if((fp3=fopen("Zones_radii.txt","w+"))==NULL)
  {
  printf("can't open the file\n");
  exit(0);
  }
  fresnel_zones_radii ( l );

  for (j=1;j<max;j++) {
  a[j]=Simp(l[j-1],l[j],f1);
  b[j]=Simp(l[j-1],l[j],f2);
  s[j]=sqrt(a[j]*a[j]+b[j]*b[j]);
  
  accuma[0]=0;
  accumb[0]=0;
  accums[0]=0;
  accuma[j]=accuma[j-1] + a[j];
  accumb[j]=accumb[j-1] + b[j];
  accums[j]=sqrt(accuma[j]*accuma[j]+accumb[j]*accumb[j]);
  
  }
  
 for (j=1;j<max;j++) {
    if ( b[j] >= 0 )
	{printf("I[%2d]≈%12.8f +%12.8fi  ABS(I[%2d])≈%12.8f\n",j,a[j],b[j],j,s[j]);
	 fprintf(fp1,"%12.8f +%12.8fi  %12.8f\n",a[j],b[j],s[j]);
	}
	else
    {printf("I[%2d]≈%12.8f -%12.8fi  ABS(I[%2d])≈%12.8f\n",j,a[j],b[j]*(-1),j,s[j]);
	 fprintf(fp1,"%12.8f -%12.8fi  %12.8f\n",a[j],b[j]*(-1),s[j]);
	}
}
fclose(fp1);

for (j=1;j<max;j++) {
    if ( accumb[j] >= 0 )
	{ printf("ΣI(1-%d)≈%12.8f +%12.8fi  ABS(ΣI(1-%d)])≈%12.8f\n",j,accuma[j],accumb[j],j,accums[j]);
     fprintf(fp2,"%12.8f +%12.8fi  %12.8f\n",accuma[j],accumb[j],accums[j]);
	}	
    else
    {printf("ΣI(1-%d)≈%12.8f -%12.8fi  ABS(ΣI(1-%d)])≈%12.8f\n",j,accuma[j],accumb[j]*(-1),j,accums[j]);
    fprintf(fp2,"%12.8f -%12.8fi  %12.8f\n",accuma[j],accumb[j]*(-1),accums[j]);
	}
  
  //fprintf(fp1,"\r\n");
}
  fclose(fp2);

for (i=1;i<max;i++) {
    fprintf(fp3,"%12.6f %12.6f %12.6f\n",rd[i],Rd[i],l[i]);
    //fprintf(fp2,"\r\n");
}
  fclose(fp3);

// *计算第一菲涅耳带面积之半的积分*//

  a1=Simp(0.0,0.70710657*l[1],f1);
  b1=Simp(0.0,0.70710678*l[1],f2);
  s1=sqrt(a1*a1+b1*b1);
  printf("%12.8f +%12.8fi  %12.8f\n",a1,b1,s1);
 //write(accums,max,"1.txt");
 //write(l,max,"2.txt");
}