/* This program may be time-consuming. Smaller values of
transient, tMax, or paraNum result in shorter execution
time. */ //使用较小的值可以提高该程序运行时间
#include<stdio.h> //头文件
#include<math.h>
#define PI 3.141592653589793238
#define C 2.99792458e8
#define M 6 // number of equations 方程数量
// variable parameter values 可变参量值
double r3Dri = 0.008;
// reflectivity of external mirror of drive 驱动激光器外部反射镜的反射率；K与r3之间存在关系
double JRatio = 1.3;
// normalized injection current 标准注入电流
double LDri = 0.6;
// external cavity length (one-way) for drive [m] 驱动激光器外腔长度（one-way）
double LInj = 1.2;
// distance between drive and response lasers [m]	DL与RL之间的距离
double detun = 0.0e9;
// optical frequency detuning [Hz]  光失谐频率
double kapInjRatio = 1.0;
// injection strength ratio to drive feedback strength 注入强度与驱动反馈强度之比
// fixed parameter values  固定参数值
double r2 = 0.556;
// reflectivity of internal mirror 内镜反射率
double tauIn = 8.0e-12;
// internal cavity round-trip time [s] 内腔往返时间
double lambda = 1.537e-6;
// optical wavelength [m] 光学波长 在选择传输波长时，主要综合考虑光纤损耗和散射。目的是通过向最远的距离、以最小的光纤损耗来传输最多的数据。在传输中信号强度的损耗就是衰减。衰减度与波形的长度有关，波形越长，衰减越小。光纤中使用的光在850、1310、1550nm处的波长较长，故此光纤的衰减较小，这也导致较少的光纤损耗。并且这三个波长几乎具有零吸收，最为适合作为可用光源在光纤中传输。
double Gn = 8.4e-13; // gain coefficient  增益系数
double N0 = 1.4e24;
// carrier density at transparency 透明载流子浓度
double tauP = 1.927e-12; // photon lifetime [s] 光子寿命
double tauS = 2.04e-9; // carrier lifetime [s] 载流子寿命
double alpha = 3.0; // alpha parameter aipha参数 线宽增强因子
double Nth, Jth;
double J, kapDri, kapInj, tauDri, tauInj;
double omega0, phaseShiftDri, phaseShiftInj;

#define DELAY_MAX 100000
// maximum array size for delay 延迟的最大数组大小
double eDelayDri[DELAY_MAX], phiDelayDri[DELAY_MAX];
double eDelayInj[DELAY_MAX], phiDelayInj[DELAY_MAX];
int delayDriNum, delayInjNum;
int delayDriIndex, delayInjIndex;
#define WAVE_MAX 200000
// maximum array size to save waveforms 最大阵列大小保存波形
double wave[3][WAVE_MAX];
// calculation of parameter values 参数值的计算
void calcParameter(double h)  //计算参数
{
	Nth = N0 + 1.0 / tauP / Gn;
	// carrier density at threshold 阈值载流子密度
	Jth = Nth / tauS;
	// injection current at threshold 阈值注入电流
	J = JRatio * Jth; // injection current 注入电流
	kapDri = (1 - r2 * r2) * r3Dri / r2 / tauIn;
	// feedback strength of drive DL反馈强度
		kapInj = kapInjRatio * kapDri;
	// injection strength from drive to response D到L的注入强度
	tauDri = 2.0 * LDri / C;
	// external-cavity round-trip time for drive  D 外腔往返时间
	tauInj = LInj / C;
	// coupling time from drive to response  从D-RL的耦合时间
	omega0 = 2.0 * PI * C / lambda;
	// optical angular frequency  光学角频率
	phaseShiftDri = fmod(omega0 * tauDri, 2.0 * PI);
	// initial phase shift for drive  DL初始相移
	phaseShiftInj = fmod(omega0 * tauInj, 2.0 * PI);
	// initial phase shift for coupling  耦合的初始相移
	delayDriNum = (int)(tauDri / h);
	delayInjNum = (int)(tauInj / h);
}
void initializeDelay(double a[])
{
	int i;
	delayDriIndex = delayInjIndex = 0;//延时驱动指数；延时注入指数
	for (i = 0; i < DELAY_MAX; i++) {
		eDelayDri[i] = a[0];
		phiDelayDri[i] = a[1];
	}
	for (i = 0; i < DELAY_MAX; i++) {
		eDelayInj[i] = 0.0;
		phiDelayInj[i] = 0.0;
	}
}
void laser(double x[], double b[], double theta[]) //unidirectionally激光器L-K方程
{
	b[0] = 1.0 / 2.0 * (Gn * (x[2] - N0)
		- 1.0 / tauP) * x[0] + kapDri
		* eDelayDri[delayDriIndex]
		* cos(theta[0]);    //Eqs.6.1  X[2]为N，X[0]为E
	b[1] = alpha / 2.0 * (Gn * (x[2] - N0)
		- 1.0 / tauP) - kapDri
		* eDelayDri[delayDriIndex]
		/ x[0] * sin(theta[0]);  //Eqs.6.3 振幅微分eqs
	b[2] = J - x[2] / tauS - Gn * (x[2] - N0)
		* x[0] * x[0];//载流子密度方程

		b[3] = 1.0 / 2.0 * (Gn * (x[5] - N0)
			- 1.0 / tauP) * x[3] + kapInj
		* eDelayInj[delayInjIndex]
		* cos(theta[1]);     //响应激光器电场微分方程
	b[4] = alpha / 2.0 * (Gn * (x[5] - N0)
		- 1.0 / tauP) - kapInj
		* eDelayInj[delayInjIndex]
		/ x[3] * sin(theta[1]);  //振幅
	b[5] = J - x[5] / tauS - Gn * (x[5] - N0)
		* x[3] * x[3];  //载流子密度
}
void rungeKutta(double a[], double h, double t)  //龙格库塔常微分方程积分
{
	int i, j;
	double x[M], b[4][M];
	double theta[2];
	theta[0] = fmod(phaseShiftDri + a[1]
		- phiDelayDri[delayDriIndex], 2.0 * PI);
	theta[1] = fmod(phaseShiftInj + a[4]
		- phiDelayInj[delayInjIndex]
		- 2.0 * PI * detun * t, 2.0 * PI);
	for (i = 0; i < 4; i++) {
		for (j = 0; j < M; j++) {
			if (i == 0) x[j] = a[j];
			if (i == 1) x[j] = a[j] + h * b[0][j] / 2.0;
			if (i == 2) x[j] = a[j] + h * b[1][j] / 2.0;
			if (i == 3) x[j] = a[j] + h * b[2][j];
		}
		laser(x, b[i], theta);
	}
	for (i = 0; i < M; i++) {
		a[i] += h * (b[0][i] + 2.0 * b[1][i]
			+ 2.0 * b[2][i] + b[3][i]) / 6.0;
	}
	// Renew arrays for delay  为延迟更新数组
	eDelayDri[delayDriIndex] = a[0];
	eDelayInj[delayInjIndex] = a[0];
	phiDelayDri[delayDriIndex] = a[1];
	phiDelayInj[delayInjIndex] = a[1];
		delayDriIndex = (delayDriIndex + 1)
		% delayDriNum;
	delayInjIndex = (delayInjIndex + 1)
		% delayInjNum;
}


//混沌同步质量：互相关

double crossCorrelation(double wave1[],
	double wave2[], int n)
{
	int i;
	double sum1, sum2;
	double ave1, ave2;
	double dev1, dev2;
	double var1, var2;
	double cov;
	sum1 = sum2 = 0.0;
	for (i = 0; i < n; i++) {
		sum1 += wave1[i];//对时域波形图各参考点求和，+=求和后的值作为输出值
		sum2 += wave2[i];
	}
	ave1 = sum1 / n; //求平均，取值点个数为n
	ave2 = sum2 / n;
	var1 = var2 = cov = 0.0; 
	for (i = 0; i < n; i++) {
		dev1 = wave1[i] - ave1; 
		dev2 = wave2[i] - ave2;
		var1 += dev1 * dev1; 
		var2 += dev2 * dev2;
		cov += dev1 * dev2;
	}
	var1 /= n;   //标准偏差的平方
	var2 /= n;
	cov /= n;     //互相关值中的尖括号内容，remains time average
	if (var1 == 0.0 || var2 == 0.0) return 0.0;
	return cov / (sqrt(var1) * sqrt(var2));








}
int main()
{
	int i, p;
	double a[M], t;
		double h = 5.0e-12; // calculation step 计算步骤

	double transient = 5000.0e-9; // transient time 瞬态时间
	double tMax = 1000.0e-9; // plot time  
	int trans = (int)(transient / h);
	int n = (int)(tMax / h);
	// range for parameter change 参数更改范围
	double paraMin = 0.0;
	double paraMax = 12.0;
	int paraNum = 100;
	double paraDiv = (paraMax - paraMin) / paraNum;
	// Initial conditions 初始条件
	a[0] = 1.3e10;
	// electric field amplitude for drive DL电场振幅
	a[1] = 0.0; // electric field phase for drive  DL电场相位
	a[2] = 1.90e24; // carrier density for drive  DL载流子密度
	a[3] = 1.4e10;
	// electric field amplitude for response  RL电场振幅				
	a[4] = 0.0; // electric field phase for response  RL电场相位
	a[5] = 1.85e24; // carrier density for response  RL载流子密度
	initializeDelay(a);




	for (p = 0; p < paraNum; p++) {
		kapInjRatio = paraMin + paraDiv * p;
		// parameter change
		calcParameter(h);
		// calculation for transient
		for (i = 0; i < trans; i++) {
			t = h * i;
			rungeKutta(a, h, t);
		}
		for (i = 0.0; i < n; i++) {
			t = h * (trans + i);
			wave[0][i] = a[0] * a[0];
			// Drive intensity
			wave[1][i] = a[3] * a[3];
			// Response intensity  响应强度
			wave[2][i] = eDelayDri[delayDriIndex]
				* eDelayDri[delayDriIndex];
			// Delayed drive intensity
			rungeKutta(a, h, t);
		}
			printf("%e\t", kapInj * 1e-9);//注入强度
		printf("%e\t", crossCorrelation(wave[0],
		wave[1], n));
		// Identical synchronization  相同同步（完全同步）
		printf("%e\n", crossCorrelation(wave[2],
			wave[1], n));
		// Generalized synchronization 广义同步
	}
	return 0;
}