plasma2.cpp

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

	outfile<<"x   " <<" amp  "<<" tham  "<<" dr   "<<" di   "<<endl; //By Origin;
	for(int i=nn/2;i<=nn-1;i++){
		if (i%nndot==0){
			x=((i+n2)%nn-n2)*h;
			mod=sqrt(data[i+i]*data[i+i]+data[i+i+1]*data[i+i+1]);
			dr=(data[i+i+2]-data[i+i])/h;
			di=(data[i+i+3]-data[i+i+1])/h;
			tham=exp(intam(x)/4/nc);
			outfile << x << "    " << mod <<"    " <<tham<<"   "<< dr<<"    " <<di<< endl;
		}
	}	
	for(i=0;i<=nn/2-1;i++){
		if (i%nndot==0){
			x=((i+n2)%nn-n2)*h;
			mod=sqrt(data[i+i]*data[i+i]+data[i+i+1]*data[i+i+1]);
			dr=(data[i+i+2]-data[i+i])/h;
			di=(data[i+i+3]-data[i+i+1])/h;
			tham=exp(intam(x)/4/nc);
			outfile << x << "    " << mod <<"    " <<tham<<"   "<< dr<<"    " <<di<< endl;
		}
	}

	outfile.close();
}

double plasma::intam(double x)
// get the corrective term value by propagation eq. or theory solution. 
{
	double temp2=0,temp1=0, temp3=0, h=z/n,ing=0,hx=ld/nn;
	for(int i=-n/2;i<n/2;i++){
//		temp2=temp1=0;
//		for(int j=-n/2;j<=i;j++){	temp1 += selfne(x,j*h)*h;	temp2 += selfne(x+hx,j*h)*h;}
		temp1 += selfne(x,i*h)*h;
		temp2 += selfne(x+hx,i*h)*h;
		temp3 += selfne(x+2*hx,i*h)*h;
		ing += ((temp3-temp2)/hx-(temp2-temp1)/hx)/hx*h;
	}
	return ing;
}


double plasma::integ(double x) 
// inte the electron density indirect from the assuming ne for theory phase.
{
	double temp=0,h=z/n;
	for(int i=1;i<=n;i++)
		temp += selfne(x,(i-n/2)*h)*h;		
	return temp;
}


void plasma::printphase(char *fn,int nndot, bool correct)
// put out the phase;*fn is file name; nndot is the interval dot.
// out dots number is nn/nndot; If correct== True, it will calculate the corrective term and will
//need much time. 
{
	int  i,odd;
	double *temp,swap;
	int n2=nn/2;

	temp=new double[nn];
	for(i=0;i<=nn-1;i++)
		if (i<n2){
		     temp[i+n2]=atan2(data[i+i+1],data[i+i]);                     					
		}
		else{
		     temp[i-n2]=atan2(data[i+i+1],data[i+i]);
		}
		odd=0;swap=0;

	for(i=0;i<nn/2-1;i++){
        if(temp[i+1]*swap<-Pi) odd+=1;
		swap=temp[i+1];
        temp[i+1]-=(2*Pi*odd);
	}
	temp[nn/2]-=(2*Pi*odd);
  	for(i=nn/2;i<nn-1;i++){
	    	if(temp[i+1]*swap<-Pi) odd-=1;
            swap=temp[i+1];		 
            temp[i+1]-=(2*Pi*odd);
	}
	double norm	= -Pi/wavelength/nc;

	ofstream outfile(fn);
/*    outfile<<"VARIABLES = \"x\"  "<<"\"sim\"  "<<"\"the\"  "<<" \"error \"  ";
	if (correct)
		outfile<<"\" correction\" ";
	outfile<< endl;
	outfile  << "zone T=\"Phase\" "<<endl;
*/// By Tecplot

    outfile<<"x  "<<" sim   "<<" the   "<<"  error    ";
	if (correct)	outfile<<"  correction  ";
	outfile<< endl; //By Origin

	double dx=ld/nn,x,phase, correction,maxph=0,maxcr=0,max=0,max2=0;
	for(i=0;i<=nn-1;i++)	 
		if (i%nndot==0){
			x=(i-nn/2)*dx;
			phase=integ(x)*norm;
			outfile << (i-nn/2)*ld/nn << "    " << temp[i] <<"   "<<phase <<"    " <<phase-temp[i];
		
			if(correct){
				if (x==0) 
					correction=0;
				else correction= -intcr(x)*norm/4/nc;
				outfile << "   " << correction ;
			}
			outfile << endl;
	}
	outfile.close();

	delete temp;	
}

double plasma::intcr(double x)
// get the corrective term value by propagation eq. or theory solution. 
{
	double temp2,temp1,h=z/n,ing=0,hx=ld/nn;
	double tm;
	
	temp2=temp1=tm=0.0;
	for(int i=-n/2;i<n/2;i++){
//		temp2=temp1=0;
//		for(int j=-n/2;j<=i;j++){
//		temp1 += selfne(x,j*h)*h;
//		temp2 += selfne(x+hx,j*h)*h;			
//		}
		temp1 = selfne(x,i*h);
		temp2 = selfne(x+hx,i*h);
		tm += ((temp2-temp1)/hx)*h;
		ing += tm*tm*h;
	}
	return ing;
}


void  plasma::printne(char *fn,int nndot)
//put out the electron density of simulation and assume;
//*fn is file name; nndot is the interval dot. out dots number is nn/nndot;
{
	int  i,rn=nn/2;
	double h=rd/rn;
 
	ofstream outfile(fn);
/*	outfile<<"VARIABLES = \"r\"  "<<"\"abel\"  "<<"\" imag\"  "<<" \"assum \"  "<< endl;
	outfile  << "zone T=\"Density\" "<<endl;
*/ //By Tecplot

	outfile<<"r   "<<" abel   "<<"  imag   "<<"  assum    "<< endl; //By Origin
	for(i=0;i<=rn-1;i++)
		if (i%nndot==0)
			outfile << i*h << "    " << data[i+i] << "    " 
			   << data[i+i+1] << "    " << selfne(h*i,0) << endl;		

	outfile.close();
}


void plasma::derivative()
// derivate the electron density in corrective term by limitted difference.
{
	double hr=2*rd/nn;
	for(int i=0;i<nn/2;i++){
		nebak[i]=data[i+i];
	}
//	printreim("ff0.dat",4);
	data[0]=(nebak[1]-nebak[0])/hr;
	for(i=1;i<nn/2-1;i++){
		data[i+i] = (nebak[i+1]-nebak[i-1])/hr/2;
	}
	data[nn-2]=(nebak[nn/2-2]-nebak[nn/2-1])/hr;
	printreim("ff2.dat",4);
	double max=0, min=0;
	int m_max=0, m_min=0;
	for(i=1;i<nn/2-1;i++){// smoothing
		if(data[i+i]>data[i+i-2] && data[i+i]>data[i+i+2]){
			if(fabs(data[i+i]-data[i+i-2])>0.2 && fabs(data[i+i]-data[i+i+2])>0.2) 
				max=(data[i+i-2]+data[i+i]+data[i+i+2])/3;
			else 
				max=data[i+i];
			m_max=i;
		}
		if(data[i+i]<data[i+i-2] && data[i+i]<data[i+i+2]){
			if(fabs(data[i+i]-data[i+i-2])>0.2 && fabs(data[i+i]-data[i+i+2])>0.2) 
				max= (data[i+i-2]+data[i+i]+data[i+i+2])/3;
			else 
				min=data[i+i];
			m_min =i;
		}
		if(max !=0.0 && min !=0.0 && abs(m_max-m_min)<10 && data[i]!=0)
			data[i+i]= (max+min)/2;
	}
	printreim("fftsm.dat",4);
}

void plasma::correct()
//inte the derivative value of the electron density.
// obtain the value of the corrective term for the abel transform value.
{
	double *ne= new double[nn/2];
	double zstep=z/(n);
	double drd=2*rd/nn;
	double tempr, r,r2,x,x2,zn,z2,z1;

	derivative();
	for(int i=0; i<nn/2; i++){
		x=i*ld/nn;
		x2=x*x;
		ne[i]=0;
		for(int j=-n/2; j<n/2; j++){
			zn=zstep*j;
			tempr=0;
			for(int ir=0; ir<nn/2-1; ir ++){
				r = ir*drd;
				r2 = r+drd;
				if (x==0 && r==0){
					z2=sqrt(r2);
					if (zn>0){
						tempr += data[ir+ir+2]/r2*z2;
						if(z2<zn)
							tempr += data[ir+ir+2]/r2*z2;
						else 
							tempr += data[ir+ir+2]/(r2)*zn;
					}else{
						if(z2>zn)
							tempr += data[ir+ir+2]/r2*(z2-zn);
					}
				}
				if(r==0 || r<x || r2<x) continue;
				z1=sqrt(r*r-x2);
				z2=sqrt(r2*r2-x2);
				if (zn>=0){
					if(z2<z/2)
						tempr += (data[ir+ir]/r+data[ir+ir+2]/r2)*(z2-z1)/2;//cal -L/2~ 0;
					else if (z2>z/2 && z1<z/2)
						tempr += data[ir+ir]/r*(z/2-z1);
					if (z2<=zn){
						tempr += (data[ir+ir]/r+data[ir+ir+2]/r2)*(z2-z1)/2;//cal 0~ z;
					}else if(z2>zn && z1<zn)
						tempr += data[ir+ir]/r*(zn-z1);// cal z1<z;
				}else {
					if(z2<=z/2 && z1>=zn){
						tempr += (data[ir+ir]/r+data[ir+ir+2]/r2)*(z2-z1)/2;//cal -L/2~ z;
					}else{
						if(z2>zn && z1<zn)
							tempr += data[ir+ir+2]/(r2)*(z2-zn);
						if(z2>z/2 && z1<z/2)
							tempr += data[ir+ir]/r*(z/2-z1);
					}			
				}
			}
			ne[i] += tempr*tempr;
		}
		
		ne[i] *= (x2*zstep);
		ne[i] /= (4.0*nc);
	}

	double max=0, min=0;
	for(i=1;i<nn/2-1;i++){//SMOOTH
		x=i*ld/nn;
		if (ne[i]>ne[i-1] && ne[i]>ne[i+1]){
			max=ne[i];
//			m_max=i;
		}
		if (ne[i]<ne[i-1] && ne[i]<ne[i+1]){
			min=ne[i];
//			m_min=i;
		}
		if (min!=0 && max!=0 && ne[i]!=0) 
			ne[i]=(max+min)/2;
	}
	
	for(i=0;i<nn/2;i++)
		databak[i]=databak[i]+ne[i];
	for(i=nn-1;i>=nn/2;i--)
		databak[i] = databak[i]+ne[nn-1-i];

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

bool plasma::judgene(double val)
// By val to judge whether need loop for correct.
{
	double max=0;
	double temp;
	for(int i=0;i<nn/8;i++){
		temp=fabs(nebak[i]-data[i+i]);
		if(temp>max) max=temp;
	}	cout<<max<<endl;
	if (max>val) return true;
	else return false;

}

void plasma::printreim(char* fn ,int nndot)
// put out the real and image part of the electron density.
{
	int  i;
	double h=rd/nn;
 
	ofstream outfile(fn);
/*	outfile<<"VARIABLES = \"x\"  "<<"\"real\"  "<<"\" imag\"  "<< endl;
	outfile  << "zone T=\"Density\" "<<endl;
*/ //by Tecplot
	
	outfile<< "x  "   <<"real  "<<"imag  "<< endl; // by origin
	for(i=0;i<=nn-1;i++)
		if (i%nndot==0)
			outfile << i*h << "    " << data[i+i] << "    " << data[i+i+1] << endl;	

	outfile.close();
}