dfblaser.cpp

dfb laaser c++ novel

/*
写出函数Ref_Trans.dat中，A:波长，B:透射谱％，C:反射谱％，D:透射谱DB，E:反射谱DB，F:透射时延，G:反射时延。
  用传输矩阵计算光纤的透射谱，反射谱，及实验特性。
*/

#include <iostream>
#include <complex>
#include <fstream>
#include <cstdlib>
#include <iomanip>
using namespace std;
#define PI 3.1415926535898
#define PLANCK 6.626196e-27
#define OPC 299792458

void DFBLaser(long double,long double,int,int);
extern void CMatrMuti2222(complex<long double> [2][2],complex<long double> [2][2],complex<long double> [2][2]);

void DFBLaser(long double wave_max,long double wave_min,int n_wave,int n_grating)
{	
	int i,j;
	char ent;//ent用于等待屏幕输入“Enter”以便确认，i，j为计数变量
	long double wave_step,dwave=1e-10;//波长步长
	long double *wave,*z,*dcm,*acm,*k,*phase,*dz;//*inver_x
	long double kappa,betab;//传输矩阵中所用到的参数
	complex<long double> c_i(0,1),c_2(2,0);//c_i表示虚数i
	complex<long double> c_f3[2][2],c_f2[2][2],c_f1[2][2];
	complex<long double> c_dbeta,c_gamma;//复系数deltabeta和r
	complex<long double> c_eb0,c_efz,c_eb0t,c_efzt;//定义传输矩阵中的四个矩阵元,ef为forward方向，eb为back方向，
	long double gain=0;//暂时保留gain的位置，这个值可能是数组或者复数等 多种形式。
	long double ref,trans,refdb,transdb,refps,transps,reftd,transtd,refpst,transpst;//dB表示的,ps:phase

	wave=(long double *)calloc(n_wave, sizeof(long double));  

	z=(long double *)calloc(n_grating, sizeof(long double));//返回一个指针，该指针指向的空闲空间足以容纳由n_grating个long double长的对象组成的数组
	dcm=(long double *)calloc(n_grating, sizeof(long double));
    acm=(long double *)calloc(n_grating, sizeof(long double));
	k=(long double *)calloc(n_grating, sizeof(long double));	
	phase=(long double *)calloc(n_grating, sizeof(long double));
   	dz=(long double *)calloc(n_grating, sizeof(long double));
	//inver_x=(long double *)calloc(n_grating, sizeof(long double));
	//gain=(long double *)calloc(n_grating, sizeof(long double));   这一行留给增益用

	if((z==NULL)||(dcm==NULL)||(acm==NULL)||(k==NULL)||(phase==NULL)||(dz==NULL))
	{       
		printf("不能分配内存！\n回车键退出：");
		scanf("%c",&ent);
		exit(0); 
	} 

	wave_step=(wave_max-wave_min)/n_wave;
	for(i=0;i<n_wave;i++) 
		wave[i]=wave_min+wave_step*i; //波长排布
	
	ifstream in("FBG_struct.dat");
	for(i=0;i<n_grating;i++)
	{	
     	in>>z[i]>>dcm[i]>>acm[i]>>k[i]>>phase[i]>>dz[i];// 光栅的结构参数输入 
		z[i]*=1e6;
		dz[i]*=1e6;
	}
	ofstream out("Ref_Trans.dat");

	for(i=0;i<n_wave;i++)//i的每个值代表一个波长
	{
		c_f1[0][0]=1;
		c_f1[0][1]=0;
		c_f1[1][0]=0;
		c_f1[1][1]=1;
		for(j=0;j<n_grating;j++)//j的每个值代表一个光栅位置，此循环结束以后将得到所有传输矩阵相乘的结果。
		{
			//计算对应于wave[i]的透射反射谱，及其相位

			kappa=PI/wave[i]*acm[j];
			betab=k[j]/2;
			c_dbeta=dcm[j]*2*PI/wave[i]-betab;//+c_i*gain/c_j;//c_j=2;
			c_gamma=sqrt(kappa*kappa-c_dbeta*c_dbeta);

			if((j!=0)&&(phase[j]!=phase[j-1]))//引入相移
			{
				c_f2[0][0]=exp(-c_i*phase[j]/c_2);
				c_f2[0][1]=0;
				c_f2[1][0]=0;
				c_f2[1][1]=exp(c_i*phase[j]/c_2);
				
				CMatrMuti2222(c_f3,c_f2,c_f1);
				
				c_f1[0][0]=c_f3[0][0];//将c_t矩阵转移到c_temp矩阵中，以便下次调用；
				c_f1[0][1]=c_f3[0][1];
				c_f1[1][0]=c_f3[1][0];
				c_f1[1][1]=c_f3[1][1];
			}

			c_f2[0][0]=(cosh(c_gamma*dz[j])-c_i*c_dbeta/c_gamma*sinh(c_gamma*dz[j]));//*exp(-c_i*betab*dz[j]);//分别对四个矩阵元赋值
			c_f2[0][1]=-c_i*kappa/c_gamma*sinh(c_gamma*dz[j]);//*exp(-c_i*phase[j]);//*exp(-c_i*(betab*dz[j]+phase[j]));
			c_f2[1][0]=c_i*kappa/c_gamma*sinh(c_gamma*dz[j]);//*exp(-c_i*phase[j]);//*exp(c_i*(betab*dz[j]+phase[j]));
			c_f2[1][1]=(cosh(c_gamma*dz[j])+c_i*c_dbeta/c_gamma*sinh(c_gamma*dz[j]));//*exp(c_i*betab*dz[j]);

			CMatrMuti2222(c_f3,c_f2,c_f1);

			c_f1[0][0]=c_f3[0][0];//将c_t矩阵转移到c_temp矩阵中，以便下次调用；
			c_f1[0][1]=c_f3[0][1];
			c_f1[1][0]=c_f3[1][0];
			c_f1[1][1]=c_f3[1][1];
		}

		c_eb0=-c_f3[1][0]/c_f3[1][1];//这两个方程是求解二元一次方程组的结果
		c_efz=c_f3[0][0]+c_f3[0][1]*c_eb0;

		trans=abs(c_efz)*abs(c_efz);
		ref=abs(c_eb0)*abs(c_eb0);

		transdb=10*log10(trans);
		refdb=10*log10(ref);

		transps=atan(imag(c_efz)/real(c_efz));
		refps=atan(imag(c_eb0)/real(c_eb0));

		//下面一段计算加上wave[i]微小波长后的反透射相位。		

		wave[i]+=dwave;

		c_f1[0][0]=1;
		c_f1[0][1]=0;
		c_f1[1][0]=0;
		c_f1[1][1]=1;
		for(j=0;j<n_grating;j++)//j的每个值代表一个光栅位置，此循环结束以后将得到所有传输矩阵相乘的结果。
		{
			kappa=PI/wave[i]*acm[j];
			betab=k[j]/2;
			c_dbeta=dcm[j]*2*PI/wave[i]-betab;//+c_i*gain/c_j;//c_j=2;
			c_gamma=sqrt(kappa*kappa-c_dbeta*c_dbeta);

			if((j!=0)&&(phase[j]!=phase[j-1]))//引入相移
			{
				c_f2[0][0]=exp(-c_i*phase[j]/c_2);
				c_f2[0][1]=0;
				c_f2[1][0]=0;
				c_f2[1][1]=exp(c_i*phase[j]/c_2);
				
				CMatrMuti2222(c_f3,c_f2,c_f1);
				
				c_f1[0][0]=c_f3[0][0];//将c_t矩阵转移到c_temp矩阵中，以便下次调用；
				c_f1[0][1]=c_f3[0][1];
				c_f1[1][0]=c_f3[1][0];
				c_f1[1][1]=c_f3[1][1];
			}

			c_f2[0][0]=(cosh(c_gamma*dz[j])-c_i*c_dbeta/c_gamma*sinh(c_gamma*dz[j]));//*exp(-c_i*betab*dz[j]);//分别对四个矩阵元赋值
			c_f2[0][1]=-c_i*kappa/c_gamma*sinh(c_gamma*dz[j]);//*exp(-c_i*phase[j]);//*exp(-c_i*(betab*dz[j]+phase[j]));
			c_f2[1][0]=c_i*kappa/c_gamma*sinh(c_gamma*dz[j]);//*exp(-c_i*phase[j]);//*exp(c_i*(betab*dz[j]+phase[j]));
			c_f2[1][1]=(cosh(c_gamma*dz[j])+c_i*c_dbeta/c_gamma*sinh(c_gamma*dz[j]));//*exp(c_i*betab*dz[j]);

			CMatrMuti2222(c_f3,c_f2,c_f1);

			c_f1[0][0]=c_f3[0][0];//将c_t矩阵转移到c_temp矩阵中，以便下次调用；
			c_f1[0][1]=c_f3[0][1];
			c_f1[1][0]=c_f3[1][0];
			c_f1[1][1]=c_f3[1][1];
		}
		wave[i]-=dwave;

		c_eb0t=-c_f3[1][0]/c_f3[1][1];//这两个方程是求解二元一次方程组的结果
		c_efzt=c_f3[0][0]+c_f3[0][1]*c_eb0t;

		transpst=atan(imag(c_efzt)/real(c_efzt));
		refpst=atan(imag(c_eb0t)/real(c_eb0t));

		//防止相位突变
		if(fabs(transpst-transps)>PI/2)
			(transpst>transps)?(transps+=PI):(transpst+=PI);

		if(fabs(refpst-refps)>PI/2)
			(refpst>refps)?(transps+=PI):(refpst+=PI);
		
		transtd=wave[i]*wave[i]/2/PI/OPC*(transpst-transps)/dwave*1e3;
		reftd=wave[i]*wave[i]/2/PI/OPC*(refpst-refps)/dwave*1e3;

		out<<setprecision(14)<<wave[i]<<" "<<trans<<" "<<ref<<" "<<transdb<<" "<<refdb<<" "<<transtd<<" "<<reftd<<endl;
	}
}