📄 soanoise.cpp
字号:
/*******************************************************************
This program is to investigate the noise characteristics of
semiconductor optical amplifiers.
*******************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <malloc.h>
#include <math.h>
#include <complex.h>
#include <iostream.h>
int const N=4096;
double const PI=3.14159265;
double const ph=1.05e-34; //Plank constant
double const c=3.0e8; //light velocity
double an,achc,achv,ashb,ep,epchc,epchv,epshb,ts,tchc,tchv,tshb;
double kappa,Psat,Esat,A,L,gama,vg,gn,N0,G0;
//G0: unsaturated power gain
//an,achc,achv,ashb: the linewidth enhancement factors
//ts: Carrier lifetime
//tchc,tchv,tshb: characteristic time constant with CH and SHB
//epchc,epchv,epshb: gain compression factors
//ep: total Gain compression factor
//kappa: conversion factor from photon density to power, kappa=h*w*Aeff*vg
//Psat,Esat: the saturation power and energy. Psat=Isat*Aeff, Esat=Psat*ts
//A,L: the active area and the length of SOA
//gama: mode confinement factor
//vg: light velocity in SOA
//gn: differential gain coefficient
//N0: the carrier density at transparency
void soagain1(double, double, double *, double *, double *);
double gfunc(double, double, double, double);
void soagain2(double, double, double *, double *, double *);
double soagain3(double);
/*****************************************************************
The main program
*****************************************************************/
void main()
{
int i,j,ch;
double Tp,Tpb,tdy,cc,En1,En2,P01,P02,T1,T2,Pav1,Pav2,t,dt,omiga,re,im,w0;
//Tp: the pump pulse width
//Tpb: the probe pulse width
//tdy: the time dalay between pump and probe, tdy<0 corresponding to probe pulse precede to pump pulse
//cc: phenomenological parameter compensate nonplane waveguide
//En1,En2: the input pulse energy energy of the pump and probe signal
//Pav1,Pav2: the input average power of the pump and probe signal
//P01,P02: the input pulse peak power of the pump and probe signal
//T1,T2: the period of the pump and probe pulses
//y0: the initial gain coefficient
//Tw: the time window width
//omiga: the beat frequency
double pue[N],pre[N],pin[N],dpin[N],g[N],g1[N],g3,nx[N],Wase,Pase;
double s1,s2,t1,B,g00,G,ef,SNR;
m_complex x1,x2,x3,pout1[N],f[N],k[N];
//pue[],pre[]: store pump pulse and probe pulse respecively
//pin[]: pin[]=pue[]+pre[]
//dpin[N]: store the differential of pin
//pout1[]: store conjugatd signal
//g1[],g2[]: store the SOA gain
//nx[]: store the noise parameter for SOA
//Wase: the output noise power spectral density of the SOA
//Pase: the output noise power of the SOA, Pase=Wase*B
//s1,s2,s3,ft,t1: temporary variables
FILE *fp;
//Set the parameters of SOA
ts=200.0; //in ps
tchc=0.7;
tchv=0.0;
tshb=0.07;
an=4.0;
achc=3.0;
achv=0.0;
ashb=0.0;
epchc=1.0e-23;
epchv=0.0;
epshb=0.5e-23;
ep=epchc+epchv+epshb;
// epchc=epchv=epshb=0.0;
// ep=0.0;
A=0.4e-12; //in m*m, that is 0.4 um*um
L=2.0e-4; //200um
gama=0.3;
gn=3.0e-20; //in m*m
N0=1.0e24; //in 1/(m*m*m)
Psat=28.4*1.0e-3;
Esat=Psat*ts; //in pJ
cc=1.0;
w0=2.0*PI*c/1.55e-6;
vg=c/3.56;
kappa=ph*w0*A/gama*vg;
B=c*(1.0/1.55e-6-1.0/1.551e-6); //the bandwidth of the optical bandpass filter is 3nm
/*
do
{
cout<<"The input signal is continuous wave or pulses:"<<endl;
cout<<"1.continuous wave; 2.pulses; Please choose:";
cin>>ch;
}while(ch!=1&&ch!=2);
*/
ch=2;
fp=fopen("c:\\temp\\x1.dat","wt");
if(fp==NULL)
{
printf("File eta.dat open error!!\n");
exit(1);
}
Tpb=10.0/1.665;
Tp=10.0/1.665;
T1=100.0/Tp; //T1=100ps
T2=T1;
dt=T1/N;
j=0;
tdy=0.0;
omiga=1.5*2.0*PI; //in THz
// Pav1=10.0;
Pav1=10.0;
for(Pav2=-20.0;Pav2<=20.0;Pav2=Pav2+1.0)
// for(G=2.0;G<=40.0;G=G+1.0)
{
G=20.0;
G0=pow(10.0,G/10.0);
g00=log(G0)/L;
//calculate ASE noise
if(ch==1) //pump is the continuous wave
{
g3=soagain3(pow(10.0,(Pav1-30.0)/10.0));
Wase=(g00+gama*gn*N0)/g00*(g3-1)+gama*gn*N0/g00*g3*log(g3)*pow(10.0,(Pav1-30.0)/10.0)/Psat;
Wase=(g00+gama*gn*N0)/g00*(g3-1); //constant nsp
Wase=ph*w0*Wase;
P01=pow(10.0,(Pav1-30.0)/10.0); //in W
En2=pow(10.0,(Pav2-30.0)/10.0)*T2*Tp; //in pJ
P02=En2/(Tpb*sqrt(PI)); //in W
for(i=-N/2;i<N/2;i++)
{
t=dt*(double)i;
pue[i+N/2]=P01;
pre[i+N/2]=P02*exp(-t*t);
pin[i+N/2]=pue[i+N/2]+pre[i+N/2];
}
/*
for(i=0;i<N;i++)
{
if(i==0)
s1=pin[N-1];
else
s1=pin[i-1];
if(i==N-1)
s2=pin[0];
else
s2=pin[i+1];
dpin[i]=(s2-s1)/(2.0*dt);
}
soagain1(Tp,dt,pin,dpin,g1);
Wase=0.0;
for(i=0;i<N;i++)
Wase=Wase+g1[i]-1.0+gama*gn*N0*L/log(g1[i])*(ep/kappa*pin[i]*g1[i]*(log((1+ep/kappa*pin[i]*g1[i])/(1+ep/kappa*pin[i]))-log(g1[i]))+g1[i]-1.0);
Wase=ph*w0*Wase/N;
*/ }
else //pump is a train of optical pulses
{
En1=pow(10.0,(Pav1-30.0)/10.0)*T1*Tp; //in pJ
P01=En1/(Tp*sqrt(PI));
En2=pow(10.0,(Pav2-30.0)/10.0)*T2*Tp; //in pJ
P02=En2/(Tpb*sqrt(PI)); //in W
for(i=-N/2;i<N/2;i++)
{
t=dt*(double)i;
t1=t*Tp/Tpb-tdy;
pue[i+N/2]=P01*exp(-t*t);
pre[i+N/2]=P02*exp(-(t1*t1));
pin[i+N/2]=pue[i+N/2]+pre[i+N/2];
}
for(i=0;i<N;i++)
{
if(i==0)
s1=pin[N-1];
else
s1=pin[i-1];
if(i==N-1)
s2=pin[0];
else
s2=pin[i+1];
dpin[i]=(s2-s1)/(2.0*dt);
}
soagain1(Tp,dt,pin,dpin,g1);
Wase=0.0;
for(i=0;i<N;i++)
Wase=Wase+g1[i]-1.0+gama*gn*N0*L/log(g1[i])*(ep/kappa*pin[i]*g1[i]*(log((1+ep/kappa*pin[i]*g1[i])/(1+ep/kappa*pin[i]))-log(g1[i]))+g1[i]-1.0);
// Wase=Wase+(g1[i]-1.0)*(g00+gama*gn*N0)/g00; //nsp=constant
Wase=ph*w0*Wase/N;
// soagain2(Tp,dt,pin,g1,nx); //numerically calculate the noise property
// Wase=0.0;
// for(i=0;i<N;i++)
// Wase=Wase+g1[i]-1+nx[i]*g1[i];
// Wase=ph*w0*Wase/N;
}
Pase=Wase*B;
//calculate FWM signal
if(ch==1)
{
for(i=0;i<N;i++)
// g[i]=log(g1[i]);
g[i]=log(g3);
for(i=0;i<N;i++)
{
x1=(exp(g[i])-1.0)/m_complex(1.0+exp(g[i])*ts/Esat*pue[i],-omiga*ts)*ts/Esat;
x3=-an*x1;
x2=m_complex(1.0,-achc)*epchc/kappa*Psat/m_complex(1.0,-omiga*tchc);
x2=x2+m_complex(1.0,-achv)*epchv/kappa*Psat/m_complex(1.0,-omiga*tchv);
x2=x2+m_complex(1.0,-ashb)*epshb/kappa*Psat/m_complex(1.0,-omiga*tshb);
x2=(exp(g[i])-1.0)*x2*ts/Esat;
f[i]=-0.5*(x1+x2+m_complex(0.0,1.0)*x3);
}
for(i=0;i<N;i++)
{
x1=m_complex(1.0,-an)*g[i]/2.0;
x1=exp(x1);
x1=x1*cc*f[i];
s1=pue[i]*sqrt(pre[i]);
pout1[i]=s1*x1;
}
s1=0.0;
s2=0.0;
for(i=0;i<N;i++)
{
s1=s1+pre[i]*dt*Tp;
s2=s2+norm(pout1[i])*dt*Tp;
}
}
else
{
for(i=0;i<N;i++)
g[i]=log(g1[i]);
for(i=0;i<N;i++)
{
x1=(exp(g[i])-1.0)/m_complex(1.0+exp(g[i])*ts/Esat*pue[i],-omiga*ts)*ts/Esat;
x3=-an*x1;
x2=m_complex(1.0,-achc)*epchc/kappa*Psat/m_complex(1.0,-omiga*tchc);
x2=x2+m_complex(1.0,-achv)*epchv/kappa*Psat/m_complex(1.0,-omiga*tchv);
x2=x2+m_complex(1.0,-ashb)*epshb/kappa*Psat/m_complex(1.0,-omiga*tshb);
x2=(exp(g[i])-1.0)*x2*ts/Esat;
f[i]=-0.5*(x1+x2+m_complex(0.0,1.0)*x3);
}
for(i=0;i<N;i++)
{
x1=m_complex(1.0,-an)*g[i]/2.0;
x1=exp(x1);
x1=x1*cc*f[i];
s1=pue[i]*sqrt(pre[i]);
pout1[i]=s1*x1;
}
s1=0.0;
s2=0.0;
for(i=0;i<N;i++)
{
s1=s1+pre[i]*dt*Tp;
s2=s2+norm(pout1[i])*dt*Tp;
}
}
Wase=Pase*T1*Tp;
ef=10.0*log10(s2/s1);
SNR=10.0*log10(s2/Wase);
// for(i=0;i<N;i++)
// fprintf(fp,"%lf\t%lf\n",dt*(i-N/2)*Tp,exp(g[i]));
fprintf(fp,"%lf\t%lf\t%lf\n",Pav2,ef,SNR);
// printf("Pave=%lf has been calculated!\n",Pav1);
// fprintf(fp,"%lf\t%lf\n",Pav1,10.0*log10(Pase*1.0e3));
// printf("G0=%lf has been calculated!\n",G);
}
// for(i=0;i<N;i++)
// fprintf(fp,"%lf\t%lf\n",(double)(i-N/2)*dt*Tp,10.0*log10(g2[i]));
fclose(fp);
printf("OK!\n");
}
/*******************************************************************
The end of main program
*******************************************************************/
/*******************************************************************
This procedure is to calculate the gain of SOA according to the
input power in pulsed input case.
The boundary condition is a[0]=a[n-1]
*******************************************************************/
void soagain1(double T0, double dt, double *pin, double *dpin, double *a)
//g0: input parameter, the initial gain coefficient
//T0: input parameter, the pulse width
//dt: input parameter, the step
//pin[]: input parameter, store the input signal power profile pin(t)
//dpin[]: input parameter, store the differential of the input signal power profile pin(t)
//a[]: output parameter, store the gain coefficient of SOA
{
int i, i0;
double k1, k2, k3, k4, g0, dg0, g, dg;
g=log(G0);
dg0=g/2.0;
g0=g;
a[0]=g0;
i0=0;
for(i=2;i<N;i=i+2)
{
k1=2.0*dt*gfunc(g0,T0,pin[i0],dpin[i0]);
k2=2.0*dt*gfunc(g0+0.5*k1,T0,pin[i0+1],dpin[i0+1]);
k3=2.0*dt*gfunc(g0+0.5*k2,T0,pin[i0+1],dpin[i0+1]);
k4=2.0*dt*gfunc(g0+k3,T0,pin[i0+2],dpin[i0+2]);
g0=g0+(k1+2.0*k2+2.0*k3+k4)/6.0;
i0=i0+2;
a[i]=g0;
}
a[1]=(a[0]+a[2])/2.0;
i0=1;
g0=a[1];
for(i=3;i<N;i=i+2)
{
k1=2.0*dt*gfunc(g0,T0,pin[i0],dpin[i0]);
k2=2.0*dt*gfunc(g0+0.5*k1,T0,pin[i0+1],dpin[i0+1]);
k3=2.0*dt*gfunc(g0+0.5*k2,T0,pin[i0+1],dpin[i0+1]);
k4=2.0*dt*gfunc(g0+k3,T0,pin[i0+2],dpin[i0+2]);
g0=g0+(k1+2.0*k2+2.0*k3+k4)/6.0;
i0=i0+2;
a[i]=g0;
}
dg=fabs(a[0]-a[N-1]);
while(dg>0.001)
{
if(a[0]>a[N-1])
g=g-dg0;
else
g=g+dg0;
dg0=dg0/2.0;
g0=g;
i0=0;
a[0]=g0;
for(i=2;i<N;i=i+2)
{
k1=2.0*dt*gfunc(g0,T0,pin[i0],dpin[i0]);
k2=2.0*dt*gfunc(g0+0.5*k1,T0,pin[i0+1],dpin[i0+1]);
k3=2.0*dt*gfunc(g0+0.5*k2,T0,pin[i0+1],dpin[i0+1]);
k4=2.0*dt*gfunc(g0+k3,T0,pin[i0+2],dpin[i0+2]);
g0=g0+(k1+2.0*k2+2.0*k3+k4)/6.0;
i0=i0+2;
a[i]=g0;
}
a[1]=(a[0]+a[2])/2.0;
i0=1;
g0=a[1];
for(i=3;i<N;i=i+2)
{
k1=2.0*dt*gfunc(g0,T0,pin[i0],dpin[i0]);
k2=2.0*dt*gfunc(g0+0.5*k1,T0,pin[i0+1],dpin[i0+1]);
k3=2.0*dt*gfunc(g0+0.5*k2,T0,pin[i0+1],dpin[i0+1]);
k4=2.0*dt*gfunc(g0+k3,T0,pin[i0+2],dpin[i0+2]);
g0=g0+(k1+2.0*k2+2.0*k3+k4)/6.0;
i0=i0+2;
a[i]=g0;
}
dg=fabs(a[0]-a[N-1]);
}
for(i=0;i<N;i++)
a[i]=exp(a[i]);
}
/**************************************************************************/
/*************************************************************************
the SOA gain differential equation when t0<=tc
dy/dt=f(t)=gfunc. The parameters are quoted from ref.2
*************************************************************************/
double gfunc(double y0, double T0, double pin, double dpin)
//g0: the initial gain coefficient of SOA
//T0: the width of the input pulse, that is, the normalized time
//pin: the input pulse power
//dpin: the differential of the input pulse power
{
double y,g00;
g00=log(G0);
y=T0/ts*(g00-y0)-ep/kappa*(exp(y0)-1.0)*dpin-(exp(y0)-1.0)*pin*(ep/kappa*T0/ts+T0/Esat);
y=y/(1.0+ep/kappa*exp(y0)*pin);
return(y);
}
/**************************************************************************/
/*****************************************************************************
This procedure is to calculate the gain of SOA by solving the following
equations:
dg/dt=(g0-g)/ts-g/(1+ep*P)*P/Esat
dP/dz=g/(1+ep*P)*P
Here P and g are the functions of z and t
the gain G(t)=P(L,t)/P(0,t), where L is the lengtg of SOA
*****************************************************************************/
void soagain2(double Tp,double dt, double *pin, double *G, double *nx)
//Tp,dt: the input parameters
//pin[]: the input parameter, the optical input power profile
//G[]: the output parameter, the gain of SOA
//nx[]: the output parameter for calaulation of noise of SOA
{
int i;
double g00,g0[N],g1[N],h[N],P0[N],P1[N],z,dL,x0,flag;
g00=log(G0)/L;
dL=L/4000.0;
for(i=0;i<N;i++)
{
g0[i]=g00;
P0[i]=pin[i];
}
for(i=0;i<N;i++)
{
nx[i]=0.0;
h[i]=0.0;
}
for(z=0;z<=L;z=z+dL)
{
for(i=0;i<N;i++)
P1[i]=P0[i]+g0[i]*dL*P0[i]/(1.0+ep/kappa*P0[i]);
/*
g1[0]=g0[0];
for(i=1;i<N;i++)
g1[i]=g1[i-1]+((g00-g1[i-1])/ts-g1[i-1]/(1.0+ep/kappa*P1[i-1])*P1[i-1]/Esat)*dt*Tp;
*/
x0=g0[N-1];
g1[0]=x0;
for(i=1;i<N;i++)
g1[i]=g1[i-1]+((g00-g1[i-1])/ts-g1[i-1]/(1.0+ep/kappa*P1[i-1])*P1[i-1]/Esat)*dt*Tp;
x0=g1[N-1]+((g00-g1[N-1])/ts-g1[N-1]/(1.0+ep/kappa*P1[N-1])*P1[N-1]/Esat)*dt*Tp;
flag=(g1[0]-x0)/g1[0];
while(flag>1.0e-4)
{
g1[0]=x0;
for(i=1;i<N;i++)
g1[i]=g1[i-1]+((g00-g1[i-1])/ts-g1[i-1]/(1.0+ep/kappa*P1[i-1])*P1[i-1]/Esat)*dt*Tp;
x0=g1[N-1]+((g00-g1[N-1])/ts-g1[N-1]/(1.0+ep/kappa*P1[N-1])*P1[N-1]/Esat)*dt*Tp;
flag=(g1[0]-x0)/g1[0];
}
for(i=0;i<N;i++)
{
g0[i]=g1[i];
P0[i]=P1[i];
}
for(i=0;i<N;i++)
h[i]=h[i]+g0[i]/(1.0+ep/kappa*P0[i])*dL;
for(i=0;i<N;i++)
{
nx[i]=nx[i]+dL*exp(-h[i])/(1.0+ep/kappa*P0[i]);
}
}
for(i=0;i<N;i++)
G[i]=exp(h[i]);
for(i=0;i<N;i++)
nx[i]=gama*gn*N0*nx[i];
}
/*****************************************************************************************
The end of procedure soagain2
*****************************************************************************************/
/*****************************************************************************************
This procedure is to calculate the gain of SOA in the injection of the continuous wave
*****************************************************************************************/
double soagain3(double p)
//p is the input parameter, the power of the input continuous wave
//This procedure retuen the gain of SOA
{
double z,dL,g00,p0,G,g1;
/*
FILE *fp;
fp=fopen("e:\\temp\\n1.dat","wt");
if(fp==NULL)
{
printf("File eta.dat open error!!\n");
exit(1);
}
*/
g00=log(G0)/L;
dL=L/4000.0;
p0=p;
for(z=0;z<L;z=z+dL)
{
g1=g00/(1.0+p0/Psat);
p0=p0+g1*dL*p0;
// fprintf(fp,"%lf\t%lf\n",z/L,(g1+gama*gn*N0)/g1);
}
// fclose(fp);
G=p0/p;
return(G);
}
/*****************************************************************************************
The end of procedure soagain3
*****************************************************************************************/⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -