plasma2.cpp

光脉冲在等离子体中传输的模拟仿真程序。用于电子密度诊断。

#include "plasma2.h"

plasma::plasma():Pi(3.14159265358979)
{
	z=5;
	nn=1024;
	n=1000;
	w0=1.0;
	data=new double[2*nn];
	databak=new double[nn];
	nebak=new double[nn/2];
	ld=5;
	rd=2;
	nc=1000;
	wavelength=15.5e-6;
	double x=ld/nn;
	for(int i=0;i<nn/2;i++)
		nebak[i]=selfne(x*i,0);
}

plasma::plasma(double Z, double LD, int NN, int N,double NC, 
			   double W0, double RD,double WL):Pi(3.14159265358979)
{
	z=Z;
	nn=NN;
	n=N;
	data = new double[2*NN];
	databak =new double[NN];
	nebak=new double[NN/2];
	ld=LD;
	rd=RD;
	w0=W0;
	nc=NC;
	wavelength=WL;
	double x=ld/nn;
	for(int i=0;i<nn/2;i++)
		nebak[i]=selfne(x*i,0);
}

plasma::~plasma()
{
	delete []data;
	delete []databak;
	delete []nebak;
}

void plasma::four1(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=data-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++){
		data[i+i]/=wtemp;
		data[i+i+1]/=wtemp;
	}
}


void  plasma::initvalue() // initial the data, real =1, imag=0;
{
	int i;

	for(i=0;i<=nn-1;i++){
		data[i+i]=1;
		data[i+i+1]=0;
	}	
}

void  plasma::linpropagate()// linear propagation in linear media; in frequency domain
{
	int  i,n2=nn/2;	
	double fl,ff,temp,wr,wi;
	double zstep=z/n;
	double tmd=1.0/ld/ld;
	
	fl=-Pi*wavelength*zstep;
	for(i=0;i<=nn-1;i++){
	    ff=((i+n2)%nn-n2)*((i+n2)%nn-n2)*tmd;
		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;
	}	
}

void  plasma::nonlinear(int step)// nonlinear propagation; in space domain
{
	double temp,wr,wi,h,x,z1;
	int i,n2=nn/2;
	double ne,tt;

 	double zstep=z/n;
       
	z1=step*zstep;
	tt=-Pi*zstep/wavelength/nc;
    h=ld/nn;
	for(i=0; i<=nn-1; i++){
		x=((n2+i)%nn-n2)*h;
		ne=tt*selfne(x,z1);
		wr=cos(ne);
		wi=sin(ne);
		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  plasma::nls() //solve the nls schrodinger eq. by linear and nonlinear propagation
{
	int step,n2=n/2;
	for(step=-n2;step<=n2;step++){
		four1(1);
		linpropagate();
		four1(-1);		
	    nonlinear(step);		
	}
}

void  plasma::phase()//obtain the phase after propagation; two erro phase must is less Pi
{
	int  i,n2=nn/2,odd;
	double *temp,swap;
	double norm	= -Pi/wavelength/nc;

	temp=new double[nn];

	for(i=0;i<=nn-1;i++)
		temp[i]=atan2(data[i+i+1],data[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);
	}

	for(i=0;i<nn;i++){
		data[i+i]=temp[i]/norm;
		databak[i]=data[i+i];// bak the phase data for corrective term.
		data[i+i+1]=0;
	}	   
	delete temp;	
};

double  plasma::bessj1(double x)//computer 1st-Bessel function for hankel transform;
{
	double ax,z;
	double xx,y,ans,ans1,ans2;
	
	if ((ax=fabs(x)) < 8.0) {
	y=x*x; 
	ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
		+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
	 ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
		+y*(99447.43394+y*(376.9991397+y*1.0))));
	 ans=ans1/ans2;
	} else {
	z=8.0/ax; 
	y=z*z; 
	xx=ax-2.356194491; 
	ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4 
		+y*(0.2457520174e-5+y*(-0.240337019e-6)))); 
	ans2=0.04687499995+y*(-0.2002690873e-3 
		+y*(0.8449199096e-5+y*(-0.88228987e-6 
		+y*0.105787412e-6))); 
	ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2); 
	if (x < 0.0) ans = -ans;
	 }
	 return ans;
 };

void  plasma::hankel() //slowly hankel transform; directly 
{
	double rt,k1,k2,x,x1,temp,tm;
	int i,j,in,jn,nn2=nn/2;
	double r,h,*ne,*bssj;

    ne=new double[nn];
	bssj=new double[nn2*nn2];

	x=2*Pi;
	rt=rd/nn2;
	h=ld/nn2;
	for(j=0;j<nn;j++) 
		ne[j]=0;
	for(j=0;j<nn2-1;j++){			
		k1=j*h;
		k2=k1+k1+h;
		data[j+j]=data[j+j]+data[j+j+2];
		data[j+j+1]=(data[j+j+1]+data[j+j+3]);
		temp=data[j+j]*k2;  
		ne[0]+=temp;
		temp=data[j+1+j]*k2;  
		ne[1]+=temp; 
	}
	ne[0]*=x*h/4;
	ne[1]*=x*h/4;
	
	for(i=0;i<nn2;i++){
		r=x*rt*i;
		k2=i*h;
		in=i*nn2;
		for(j=0;j<=i;j++){
			k1=j*h;
			x1=k1*r;
			jn=j*nn2;
			temp=bessj1(x1);
			bssj[in+j]=k1*temp;
			bssj[i+jn]=k2*temp;
		}
	}	
	for(i=1;i<nn2;i++){
		r=rt*i;
		in=i*nn2;
		for(j=0;j<nn2-1;j++){
			tm=(bssj[in+j+1]-bssj[in+j]);	
			temp=data[j+j]*tm;	
			ne[i+i]+=temp;
			temp=data[j+j+1]*tm;
			ne[i+i+1]+=temp;	
		}
		tm=2*r;
		ne[i+i]/=tm;
		ne[i+i+1]/=tm;               	
	}
	
	for(i=0;i<nn2;i++){
		data[i+i]=ne[i+i];
		data[i+i+1]=ne[i+i+1];
	}
	delete  bssj;
	delete  ne;
}


void  plasma::abel() // abel transform, need hankel transform
{
	int i,n2=nn/2;
	double max=0;
	for( i=0;i<nn;i++)
		if (fabs(data[i+i])>max) max=data[i+i];	
	for(i=0;i<nn;i++)
		data[i+i]/= max;
    four1(1);
	 for(i=0;i<nn;i++){
		if(i<n2){
			data[i+i]*=ld /sqrt(nn) ;
			data[i+i+1]*=ld /sqrt(nn) ;		
		}else{
	        data[i+i]=0 ;
			data[i+i+1]=0 ;		
		}
	}
	ld=nn/2/ld;
	hankel();
	ld=nn/2/ld;
	for(i=0;i<nn/2;i++)
		data[i+i]*= max;

}

void  plasma::printamp(char *fn,int nndot)
//put out the amplitude; *fn is file name; nndot is the interval dot. out dots number is nn/nndot;
{
	int n2=nn/2;
	double h=ld/nn,mod,dr,di,x;
	double tham=0;
	ofstream outfile(fn);
/*   outfile<<"VARIABLES =\"x \" " <<"\"amp\" "<<"\"tham\" "<<"\"dr\"  "<<"\"di\"  "<<endl;
	outfile  << "zone T=\" "<<"amplitude \" "<<endl;
*///by Tecplot