resonator.cpp

球形谐振腔的模拟仿真程序

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

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
#define Pi 3.14159265338979
#define w0 1
#define wavelength 1.053e-3
#define index2real(i,NN)  ((i+NN/2)%NN-NN/2)



void simpson(double a, double b, double *data, double *d, int n, int k)
{
	int i;
	double h=(b-a)/(n);
	double s0even=data[0]+data[2*n-2];
	double s0odd=data[1]+data[2*n-1];
	double sr=0,si=0;
	double sor=0,soi=0;           ////奇数项的和
	double ser=0,sei=0;           ////偶数项的和
	for (i=1;i<=n-2;i++)
	{
		
		if (i%2==0)
		{
			ser+=data[i+i];
			sei+=data[i+i+1];
		}
		else
		{
			sor+=data[i+i];
			soi+=data[i+i+1];
		}
	}
	sr=h*(s0even+2*ser+4*sor)/3;
	si=h*(s0odd+2*sei+4*soi)/3;
	d[k+k]=sr;
	d[k+k+1]=si;	
}

void four1(double *Data,unsigned long nn,int isign)
{
   unsigned long n,mmax,m,j,istep,i;
   double wtemp,wr,wpr,wpi,wi,theta;
   double tempr,tempi;
   double *data;
   data=Data-1;
   n=nn<<1;
   j=1;
   for (i=1;i<n;i+=2){
	   if(j>i){
		   SWAP(data[j],data[i]);
           SWAP(data[j+1],data[i+1]);
	   }
	   m=n>>1;
	   while (m>=2&&j>m){
		   j-=m;
		   m>>=1;
	   }
	   j+=m;
   }
   mmax=2;
   while (n>mmax){
	   istep=mmax<<1;
	   theta=isign*2*Pi/mmax;
	   wtemp=sin(0.5*theta);
	   wpr=-2.0*wtemp*wtemp;
	   wpi=sin(theta);
	   wr=1.0;
	   wi=0.0;
	   for (m=1;m<mmax;m+=2){
		   for (i=m;i<=n;i+=istep){
			   j=i+mmax;
			   tempr=wr*data[j]-wi*data[j+1];
			   tempi=wr*data[j+1]+wi*data[j];
               data[j]=data[i]-tempr;
               data[j+1]=data[i+1]-tempi;
			   data[i]+=tempr;
			   data[i+1]+=tempi;
		   }
		   wr=(wtemp=wr)*wpr-wi*wpi+wr;
		   wi=wi*wpr+wtemp*wpi+wi;
	   }
	   mmax=istep;
   }
   double tem=sqrt(nn);
   for(i=1;i<=2*nn;i++)
   {
	   data[i]=data[i]/tem;    
   }
}
   
void Gaussian(double *data,double ld, int n)
{
	double x;
    int i;
	
	double h=2*ld/(n);
	for (i=0;i<n;i++){
		x=h*index2real(i, n)/w0;
		data[i+i]=exp(-x*x);
		data[i+i+1]=0;
	}
}

void theory_gauss(double *data, int n, double ld, double z)
{
	int i;
	double x;
	double h=2*ld/(n);
	double k=2*Pi/wavelength;
	double z0=Pi*w0*w0/wavelength;
	double w=w0*sqrt(1+z*z/z0/z0);
	double r=z+z0*z0/z;
	double alfa=atan2(z,z0)/2;
	double amplitude=sqrt(w0/w);
	for (i=0;i<n;i++){	
		x=h*index2real(i,n);
		double phi=-k*(x*x/2/r)+alfa;
		data[i+i]=amplitude*exp(-x*x/w/w)*cos(phi);
		data[i+i+1]=amplitude*exp(-x*x/w/w)*sin(phi);
	}
}

void print(char *fn,double *data,int n,double Ld)
{
	FILE *cfPtr;
	int  i,nn;
	double h,x;
	h=2*Ld/(n);
	nn=n/2;
	if((cfPtr=fopen(fn,"w"))==NULL)
		printf("File could not be opened\n");
	else
	{
		for(i=0;i<=n-1;i++)
		{
			x=h*index2real(i, n);
			double t=sqrt(data[i+i]*data[i+i]+data[i+i+1]*data[i+i+1]);
			double phi=acos(data[i+i]);
			fprintf(cfPtr,"%-f    %-f   %-f ",
				x/w0, phi,t);	
			fprintf(cfPtr,"\n");
		}
	}		
	fclose(cfPtr);
}

void read(char *fn,double *data,int n)
{
	FILE *cfPtr;	
	if((cfPtr=fopen(fn,"r+"))==NULL)
		printf("File could not be opened\n");
	else
	{
		fread(data,n+n,sizeof(double),cfPtr);
	}		
	fclose(cfPtr);
}


void propagation(double *data, int n, double ld, double z)
{
	int i,nn=n/2;
	double t=0;
	double phi=2*Pi*z/wavelength;
	t=-wavelength*Pi*z/(ld*ld*4);
	double u,pr,pi,temp;
	four1(data,n,1);
	for(i=0;i<n;i++){
		u=index2real(i, n);
    	pr=cos(t*u*u);
		pi=sin(t*u*u);
		temp=pr*data[i+i]-pi*data[i+i+1];
		data[i+i+1]=pi*data[i+i]+pr*data[i+i+1];
		data[i+i]=temp;
	}
    four1(data,n,-1);
}


/*void propagation1to2(double *f1, double *f2, int n, 
					 double ld, double z)
{
	int i,nn=n/2;
	double t=0;
	t=-wavelength*Pi*z/(ld*ld*4);
	double u,pr,pi,temp;
	
	for (i=0;i<2*n;i++){
		f2[i]=f1[i];
	}

	four1(f2,n,1);

	for(i=0;i<n;i++){
		u=((i+nn)%n-nn);
    	pr=cos(t*u*u);
		pi=sin(t*u*u+Pi);
		temp=pr*f2[i+i]-pi*f2[i+i+1];
		f2[i+i+1]=(pi*f2[i+i]+pr*f2[i+i+1]);
		f2[i+i]=temp;
	}

    four1(f2,n,-1);
}*/

void propagation1to2(double *fin,double *fout, int n,double a1, 
			   double a2 ,double L)
{
	int i=0,j=0;
	double *temp;
	temp=(double*)calloc(n+n,sizeof(double));
	double t1=-Pi*L/wavelength;
	double t2=-Pi/wavelength/L;	
	double x1=0,x2=0;
	double dx1=2*a1/(n);
	double dx2=2*a2/(n);
	double pr=0,pi=0;
	for (i=0;i<n;i++)
	{
		x2=index2real(i,n)*dx2;
		for (j=0;j<=n-1;j++)
		{
			x1=index2real(j,n)*dx1;
			double t=(x2-x1)*(x2-x1);		
			temp[j+j]=cos(t2*t)*fin[j+j]-sin(t2*t)*fin[j+j+1];
			temp[j+j+1]=cos(t2*t)*fin[j+j+1]+sin(t2*t)*fin[j+j];
		}
		simpson(-a1, a1, temp, fout,  n, i);		
	}
	free(temp);
}








void normalize(double *data, int n)
{
	int i=0;
	double max=0,t;
	for (i=0;i<n;i++){
		t=sqrt(data[i+i]*data[i+i]+data[i+i+1]*data[i+i+1]);
		if (t>max)
			max=t;
	}
	for (i=0;i<2*n;i++){
		data[i]=data[i]/max;
	}
}




void mirror(double *fin, double *phase, int nn)
{
	int i;
	for (i=0;i<nn;i++){
		double t=fin[i+i]*cos(phase[i])-fin[i+i+1]*sin(phase[i]);
		fin[i+i+1]=fin[i+i]*sin(phase[i])+fin[i+i+1]*cos(phase[i]);
		fin[i+i]=t;
	}
}

void propagation2to1(double *fin,double *fout, int n,double a1, 
			   double a2 ,double L)
{
	int i=0,j=0;
	double *temp;
	temp=(double*)calloc(n+n,sizeof(double));
	double t1=-Pi*L/wavelength;
	double t2=-Pi/wavelength/L;	
	double x1=0,x2=0;
	double dx1=2*a1/(n);
	double dx2=2*a2/(n);
	double pr=0,pi=0;
	for (i=0;i<n;i++)
	{
		x1=index2real(i,n)*dx1;
		for (j=0;j<=n-1;j++)
		{
			x2=index2real(j,n)*dx2;
			double t=(x2-x1)*(x2-x1);		
			temp[j+j]=cos(t2*t)*fin[j+j]-sin(t2*t)*fin[j+j+1];
			temp[j+j+1]=cos(t2*t)*fin[j+j+1]+sin(t2*t)*fin[j+j];
		}
		simpson(-a2, a2, temp, fout,  n, i);		
	}
	free(temp);
}
void plane(double *data,int n)
{
	int i;
	//nn=n/2;
	for (i=0;i<n;i++){
		data[i+i]=1;
		data[i+i+1]=0;
	}
}

void Gauss(double *data,int n, double ld)
{
	int i;
	double x=0,nn=n/2;
	double h=2*ld/(n);
	for (i=0;i<n;i++){		
		x=index2real(i,n)*h/w0;
		data[i+i]=exp(-x*x);
		data[i+i+1]=0;
	}
}

void SuperGauss(double *data,int n, double ld, int m)
{
	int i;
	double x=0,nn=n/2;
	double h=2*ld/(n);
	for (i=0;i<n;i++){		
		x=index2real(i,n)*h/w0;
		x=pow(x,m);
		data[i+i]=exp(-x*x);
		data[i+i+1]=0;
	}
}

void ChG(double *data,int n, double ld, double omega)
{
	int i;
	double x;
	double b=w0*omega;
	double h=2*ld/(n);
	for (i=0;i<n;i++){		
		x=index2real(i,n)*h/w0;
		data[i+i]=exp(-x*x)*cosh(b*x);
		data[i+i+1]=0;
	}
}




void phase(double *phase, int n)
{
	int i;
	double temp=0;
	for (i=0;i<n;i++){
		temp=sqrt(phase[i+i]*phase[i+i]+phase[i+i+1]*phase[i+i+1]);
		phase[i+i]=phase[i+i]/temp;
		phase[i+i+1]=phase[i+i+1]/temp;	
	}
}

void mirror_ChG(double *data, int nn, double ld,
			     double z, double omega)
{
	int j;
	double b=omega*w0;
	double z0=Pi*w0*w0/wavelength;
	double g=z/z0;
	double delta=2*ld/(nn);
	double temp,wr,wi,x;
	double r=z+z0*z0/z;
	double k=2*Pi/wavelength;	
	for(j=0;j<nn;j++){		
		x=index2real(j,nn)*delta;		
		double alfa=-(x*x+b*b/4)*g/(1+g*g)+atan(tan(g*b*x/(1+g*g))*tanh(g*g*b*x/(1+g*g)));
		wr=cos(alfa);
		wi=sin(alfa);
		temp=data[j+j]*wr-data[j+j+1]*wi;
		data[j+j+1]=data[j+j]*wi+data[j+j+1]*wr;
		data[j+j]=temp;
	}
	
}

void semiconfocal(double *phase, int n, int times,double a1,
				  double a2, double L,double F)
{
	int i;	
	double *f1,*f2;
	f1=(double*)calloc(n+n,sizeof(double));
	f2=(double*)calloc(n+n,sizeof(double));
	//ChG(f1, n, a1,2);
	//Gauss(f1,n,a1);
	plane(f1,n);
    for (i=1;i<=times;i++)
	{
		cout<<"Time is "<<i<<endl;
		if (i%2==1)
		{
			propagation1to2(f1, f2, n, a1, a2, L);
			//propagation1to2(f1, f2, n, a2, L);
			
			mirror(f2,phase,n);
			normalize(f2, n);
		}
		else
		{			
			propagation2to1(f2, f1, n, a1, a2, L);
			normalize(f1,n);
		}
		/*
		if (i==100){
			printfull("mirror2_100.dat",f2, 2*a2/n, n);
	        printfull("mirror1_100.dat",f1, 2*a1/n, n);
		}
		if (i==200){
			printfull("mirror2_00.dat",f2, 2*a2/n, n);
	        printfull("mirror1_200.dat",f1, 2*a1/n, n);
		}
		if (i==300){
			printfull("mirror2_300.dat",f2, 2*a2/n, n);
	        printfull("mirror1_300.dat",f1, 2*a1/n, n);
		}*/
	}

	print("mirror2.dat",f2, n, a2);
	print("mirror1.dat",f1, n, a1);
	free (f1);
	free (f2);

}



void main()

{
	int n=1024;              //n是在反射镜上取点的个数.
	int times=200;              //k是光束在谐振腔里反射的次数.
	double k=2*Pi/wavelength;
	double z0=Pi*w0*w0/wavelength;	
	double L=100*wavelength;         //L是谐振腔两镜面相距的长度		
	double a1=25*wavelength,a2=25*wavelength;
	double F=a1*a2/wavelength/L;
	cout<<"The fresnel number is "<<F<<endl;
	double *data;
	data=(double*)calloc(n+n,sizeof(double));	
	read("ChGPhase(z0).dat",data, n);
	semiconfocal(data, n, times, a1, a2, L,F);
	delete[] data;
	
}