sec-bj.cpp

光脉冲在光纤中传输需要的分布傅立叶变换的源程序。

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

#define Pi 3.14159265358979

//傅立叶变换，isign=1 正变换 isign=-1 逆变换, nn=2^N
void four1(double *dat, int  nn,int isign)
{                   
	int i,istep,j,m,mmax,n;                                              
    double tempi,tempr,*data;
	data=dat-1;
    double theta,wi,wpi,wpr,wr,wtemp;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          

	n=nn<<1;
    j=1;
    for(i=1;i<n;i+=2){
		if(j>i){
			tempr=data[j];    data[j]=data[i];  data[i]=tempr;
            tempi=data[j+1];  data[j+1]=data[i+1];  data[i+1]=tempi;
		}
        m=n>>1;
		while ((m>=2)&&(j>m)){
			j-=m;
            m/=2;
        }
        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;
	}
	wtemp=sqrt(nn);
	for (i=0;i<nn;i++){
		dat[i+i]/=wtemp;
		dat[i+i+1]/=wtemp;
	}
}

//线性传输部分
void propagation(double *data,int nn,double ld,double zstep,int step)
{
	int  i,n2=nn/2;	
	double fl=-2*Pi*Pi*zstep,ff,temp,wr,wi;
	four1(data,nn,1);
	for(i=0;i<=nn-1;i++){
	    ff=((i+n2)%nn-n2)*((i+n2)%nn-n2)/ld/ld;
		wr=cos(fl*ff);
		wi=sin(fl*ff);
		temp=data[i+i]*wr-data[i+i+1]*wi;				
		data[i+i+1]=data[i+i+1]*wr+data[i+i]*wi;
		data[i+i]=temp;
	}
	four1(data,nn,-1);
}

//非线性传输部分
void nonlinear(double *data,int nn,double ld,double zstep)
{
	double temp,wr,wi;
	int i;
	for(i=0;i<nn;i++){
		temp=data[i+i]*data[i+i]+data[i+i+1]*data[i+i+1];
		temp*=zstep;
		wr=cos(temp);
		wi=sin(temp);
		temp=data[i+i]*wr-data[i+i+1]*wi;
		data[i+i+1]=data[i+i]*wi+data[i+i+1]*wr;
		data[i+i]=temp;
	}
}

//初值
void sech(double *data,int nn,double ld,int N)
{
	double h,x;
	int i,n2=nn/2;
	h=ld/nn;
	for(i=0;i<nn;i++){
		x=h*((n2+i)%nn-n2);
		data[i+i]=N*1/cosh(x);
		data[i+i+1]=0;
	}	
}

//数据输出
void print(char *fn,double *data,int nn,double interval, int nndot)
{	
	ofstream outfile(fn);
	//outfile<<"VARIABLES = \"time\"  "<<"\"Intensity\"  "<<"\" Phase\"  "<< endl;
	//outfile  << "zone T=\"Density\" "<<endl;

	double intensity, phase;
	for(int it=-nn/2;it<nn/2;it++)
		if (it%nndot==0){
			int i= (it+nn)%nn;
			intensity=sqrt(data[i+i]*data[i+i]+data[i+i+1]*data[i+i+1]);
			phase=atan2(data[i+i+1], data[i+i]);
			outfile << it*interval << "    " << intensity << "    " <<phase << endl;	
		}
	outfile.close();
}

//整个NLSE
void NLS(double *data,int nn,double ld,int n,double z,int N)
{	
	double zstep=z/n;
	int step;	
	sech(data,nn,ld,N);
 	for(step=1;step<=n;step++){
		propagation(data,nn,ld,zstep,step);			
		nonlinear(data,nn,ld,zstep);			
	}
}


main()
{
	double ld,z, *data;
	int nn,n;
	char *profile="sec00.dat";
    int N=2;
 
	ld=20;
	z=4.0;	
	n=500;
    nn=1024;
	data=(double *) calloc(nn*2,sizeof(double)); 

	NLS(data,nn,ld,n,z,N);
	print(profile, data, nn, ld/nn, 1);
   
	return 0;
}