stdafx.cpp

光脉冲在光纤传输中的自陡峭问题Ｃ语言源程序。

 
#include "stdafx.h"

Nlse1Dim::Nlse1Dim()
{
	NN=1024;
	int dbnn=NN<<1;
	Dt=5/NN;
	data= new double[dbnn];
	Dz=0.01;
	Nz=100;
	OUTF_Time_Win=5;
	OUTF_Freq_Win=5;
	OUTF_Numb=NN;

	Cout_constant();
}
Nlse1Dim::Nlse1Dim(int nn, double ld, double dz, int nz,double Otw, 
				   double Ofw, int On)
{
	NN=nn;
	int dbnn=NN<<1;
	Dt=ld/NN;
	data= new double[dbnn];
	Dz=dz;
	Nz=nz;
	OUTF_Time_Win=Otw;
	OUTF_Freq_Win=Ofw;
	OUTF_Numb=On;

	Cout_constant();
}
Nlse1Dim::~Nlse1Dim()
{
	delete []data;
}
void Nlse1Dim::Cout_constant()
{
	cout<<"\n GVD Sign is "<< GVD_Sign <<";"<<endl;
	cout<<"\n Nonlinear Order is "<<NLS_Coef_NN<<";"<<endl;
	if(GVD_Coef_Thd!=0)
		cout<<"\n Third GVD coeffience is "<<GVD_Coef_Thd <<";"<<endl;
	if(NLS_Coef_SS!=0)
		cout<< "\n Self-Steeping Coeffience is "<<NLS_Coef_SS<<endl;
	if(NLS_Coef_RT!=0)
		cout<<"\n The simplified Retarded Nonlinear Response (SRS) has been considered.\n"<< 
		" Its coeffience is "<<NLS_Coef_RT <<";"<<endl;
}

void Nlse1Dim::FFT1D(double *data1d, int NN, int isign)
//one dimension Fourier Transform; isign=1, fft; isign=-1, inverse fft
{                   
	int i,istep,j,m,mmax,N;                                              
	double tempi,tempr,*Data;
	double theta,wi,wpi,wpr,wr,wtemp;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          

	Data=data1d-1;
	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++){
		data1d[i+i]/=wtemp;
		data1d[i+i+1]/=wtemp;
	}
	
}

///////////////////////////////////////////////////
void  Nlse1Dim::Phase(double *datatm,double *temp)
//obtain the phase after propagation; two erro phase must is less Pi
{
	int  i,n2=NN/2,odd;
	double swap;

	for(i=0;i<=NN-1;i++)
		temp[i]=atan2(datatm[i+i+1],datatm[i+i]);
	odd=0;swap=0;
	for(i=n2;i<NN-1;i++){
		if(temp[i+1]*swap<-Pi) odd+=1;
		swap=temp[i+1];
        temp[i+1] -= (2*Pi*odd);
	}
	temp[0]-=(2*Pi*odd);
  	for(i=0;i<n2-1;i++){
		if(temp[i+1]*swap<-Pi) odd-=1;
		swap=temp[i+1];		 
		temp[i+1] -= (2*Pi*odd);
	}
} 

void Nlse1Dim::Save_ReIm_InPh(char *fn)
{
	int Size=(int)((Dt*NN)/OUTF_Time_Win);
	if (Size<1) Size=1;
	int np=NN/(2*Size);
	int nndot=(int)(np*2/OUTF_Numb);
	if(nndot<1) nndot=1;
	ofstream outfile(fn);

	outfile<<"Time	"<<"Real	"<<"Imag	"<<"Intensity	"<<"Phase	"<<endl;//For Origin
	double intensity, *phase;
	phase=new double[NN];
	Phase(data,phase);

	for(int it=-np;it<np;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]);
			outfile << it*Dt <<"    " <<data[i+i]<<"    " <<data[i+i+1]
				<<"    " << intensity << "    " <<phase[i] << endl;	
		}
	delete []phase;
	outfile.close();
}

void Nlse1Dim::Copydata(double *data1d)
{
	for(int i=0 ;i<2*NN; i++)
		data1d[i]=data[i];
}

void Nlse1Dim::Save_Freq_ReIm_InPh(char *fn)
{
	int Size=(int)((1/Dt)/OUTF_Freq_Win);
	if (Size<1) Size=1;

	int np=NN/(2*Size);
	int nndot=(int)(np*2/OUTF_Numb);
	if(nndot<1) nndot=1;
	ofstream outfile(fn);

	outfile<<"Frequency	"<<"Real	"<<"Imag	"<<"Intensity	"<<"Phase	"<<endl;//For Origin
	double intensity, *phase, *dtfreq, hf=1/(NN*Dt);
	dtfreq=new double[2*NN];
	phase=new double[NN];
	Copydata(dtfreq);
	FFT1D(dtfreq,NN,1);
	Phase(dtfreq,phase);
	for(int it=-np;it<np;it++)
		if (it%nndot==0){
			int i= (it+NN)%NN;
			intensity=sqrt(dtfreq[i+i]*dtfreq[i+i]+dtfreq[i+i+1]*dtfreq[i+i+1]);
			outfile << it*hf <<"    " <<dtfreq[i+i]<<"    " <<dtfreq[i+i+1]
				<<"    " << intensity << "    " <<phase[i] << endl;	
		}
	delete []dtfreq;
	delete []phase;
	outfile.close();
}
//////////////////////////////////////////////////////////////////////

void Nlse1Dim::Sech(double Tau, double Chirp)
{
	cout<<"\n\n The incident field is soliton (hyperbolic secant),  with the Chirp "<<Chirp<<endl;
	double x;
	double coef=Dt/Tau, ch, tm;
	int i,n2=NN>>1;
	for(i=0;i<NN;i++){
		x= coef*((n2+i)%NN-n2);
		ch= -Chirp*x*x;
		tm=1/cosh(x);
		data[i+i]= tm*sin(ch);
		data[i+i+1]= tm*cos(ch);
	}	
}

void Nlse1Dim::triangle(double Tau)
{
	cout<<"\n\n The incident field is soliton (triangle pulse),  with the Chirp "<<0<<endl;
	double x;
	double coef=Dt/Tau;
	int i,n2=NN>>1;
	for(i=0;i<NN;i++){
		x= coef*((n2+i)%NN-n2);
		if (x>=0){
			data[i+i]=sqrt(3)*(1-x);
		}
		else{
			data[i+i]=sqrt(3)*(1+x);
		}		
		data[i+i+1]=0;
	}	
}


void Nlse1Dim::Bsoliton(double Tau)
{
	cout<<"\n\n The incident field is black soliton, i.e., Tanh."<<endl;
	double x;
	double coef=Dt/Tau;
	int i,n2=NN>>1;
	for(i=0;i<NN;i++){
		x= coef*((n2+i)%NN-n2);
		data[i+i]= tanh(x);
		data[i+i+1]= 0;
	}	
}


void Nlse1Dim::Gauss(int m, double Chirp, double Tau)
{
	cout<<"\n\n The incident field is "<<m <<" order SuperGauss  with the Chirp "<<Chirp<<endl;

	double x;
	double coef= Dt/sqrt(2)/Tau,temp,chtm;
	int i,n2=NN/2;
	for(i=0;i<NN;i++){
		x= coef*((n2+i)%NN-n2);
		x *= x;
		temp= pow(x,m);
		chtm=temp*Chirp;
		temp=exp(-temp);
		data[i+i]= temp*cos(-chtm);
		data[i+i+1]= temp*sin(-chtm);
	}	
}

void Nlse1Dim::Line_Prop(double prop_rang)
{
	int  i,n2=NN>>1;
	double fl=PI2*Pi*prop_rang*GVD_Sign, gvd=GVD_Coef_Thd*Pi*PI2*PI2*prop_rang/3;
	double it,ff,temp,wr,wi;
	double ld=1/(NN*Dt);

	FFT1D(data,NN,1);
	for(i=0;i<=NN-1;i++){
		it=((i+n2)%NN-n2)*ld;
		ff= it*it;
		temp=ff*(fl+gvd*it);
	    wr=cos(temp);
		wi=sin(temp);
		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;
	}
	FFT1D(data,NN, -1);
}

void Nlse1Dim::SelfSteep(int i, double *SS)
{
	int I2=i+i;
	double ssr=0., ssi=0.;
	double dt2=0.0,dt3=0.0,dt1=0.0;
	double temp=0., dtdb= 2*Dt, tmss;

	tmss= (NLS_Coef_SS) / (dtdb);//imag
	if(i!=NN-1) dt3=data[I2+2]*data[I2+2]+data[I2+3]*data[I2+3];	
	else dt3= data[0]*data[0]+data[1]*data[1];
	if(i!=0) dt1=data[I2-2]*data[I2-2]+ data[I2-1]*data[I2-1];
	else dt1= data[NN+NN-2]*data[NN+NN-2]+ data[NN+NN-1]*data[NN+NN-1];

	ssr=(dt3-dt1); ssi=0.;
	
	if(i!=0 && i!=NN-1){
		dt3=data[I2+2]-data[I2-2];
		dt1=data[I2+3]-data[I2-1];
	}else{
		if(i==0){
			dt3= data[2]-data[NN+NN-2];
			dt1= data[3]-data[NN+NN-1];
		}	
		if(i==NN-1){
			dt3= data[0]-data[NN+NN-4];
			dt1= data[1]-data[NN+NN-3];
		}
	}
	ssr += (dt3*data[I2]+ dt1*data[I2+1]);
	ssi += (dt1*data[I2]- dt3*data[I2+1]);

	SS[0]= -(ssi*tmss);
	SS[1]=  (ssr*tmss);
}

void Nlse1Dim::NonLn_Prop(double prop_rang)
{
	bool Self_steep=true, Raman_Diff=true;
	if(NLS_Coef_RT==0) Raman_Diff=false;
	else Raman_Diff=true;
	if(NLS_Coef_SS==0) Self_steep=false;
	else Self_steep=true;

	double temp,wr,wi,tmpr,tmpi;
	int i;
	double retard=0.0, ssr=0., ssi=0.;
	double dt2=0.0,dt3=0.0,dt1=0.0;
	double tmnls=(NLS_Coef_NN)*(NLS_Coef_NN)*prop_rang,dtdb= 2*Dt;

	double *SS;
	SS= new double[2];

	for(i=0;i<NN;i++){
		dt2=data[i+i]*data[i+i]+data[i+i+1]*data[i+i+1];
		if(Self_steep){
			SelfSteep(i, SS);
			ssr=SS[0];ssi=SS[1];
		}else
		{ssr=0.; ssi=0.;}
		if(Raman_Diff){
			if(i<NN-1) dt3=data[i+i+2]*data[i+i+2]+data[i+i+3]*data[i+i+3];	
			else dt3= data[0]*data[0]+data[1]*data[1];
			if(i!=0) dt1=data[i+i-2]*data[i+i-2]+data[i+i-1]*data[i+i-1];
			else dt1= data[NN+NN-2]*data[NN+NN-2]+data[NN+NN-1]*data[NN+NN-1];
			
			retard= (-NLS_Coef_RT)*(dt3-dt1)/(dtdb);
		}else 
			retard=0.0;
		tmpi = (dt2+retard+ ssr)*tmnls;
		tmpr = -ssi*tmnls;
		temp =exp(tmpr);
		wr=cos(tmpi)*temp;
		wi=sin(tmpi)*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;
	}
	delete []SS;
}

void Nlse1Dim::LN_NLS(int nzdot)
{	

	cout<<"\n\n The program is used to the long pulse propagation by simplited SRS."<<endl;
	if(nzdot>26*26||nzdot<1) nzdot=26*26;
	int step, dot= Nz/nzdot;//中间隔dot步输出一次电场
	if(dot<1) dot=1;
	char fn[]="riaa.dat", fq[]="fqaa.dat";
		
	Line_Prop(Dz/2);
	for(step=1;step< Nz;step++){			
		NonLn_Prop(Dz);
		Line_Prop(Dz);		
		if(step%dot==0){
        	Save_ReIm_InPh(fn);
			Save_Freq_ReIm_InPh(fq);
			if(fn[3]=='z'){
				fn[2]+=1;fq[2]++;
				fn[3]='a';fq[3]='a';
			}
			else {
				fn[3]+=1;fq[3]++;
			}
		}
	}
	Line_Prop(Dz/2);
}