reaction.h

用Monte Carlo方法

//头文件中写入各种变量、系数、还有反应通道、平衡方程

//参数定义
#define Ve 2.54e-9
#define Ku 1.0e6*Ve
#define Ks 1.0e7*Ve
#define Kf 1.0e6*Ve
#define a0 0.3
#define as 4.0
#define ad 0.3
#define Kd 20
#define Kfold 15.0e3
#define enta(T) ((T>40)?200:40)
#define K(T) ((T>40)?150:75)
#define a 0.2

#define n_channel 5
#define n_parameter 10

//变量定义
double St=40.0,Sf;
double Dt=10000.0,Df;
double Ft=10000.0,Ff;
double Pt=2.0e6;
double SD;
double UD;
double SDF;

int T=37;

//数组
void (*p_channel[n_channel+1])();
double *parameter[n_parameter];
double (*poss[n_channel+1])();
char *id_parameter;

//通道函数
void cf1()
{
	Sf++;
	St++;

}
void cf2()
{
	Sf--;
	St--;
}
void cf3()
{
	SDF--;
	Df++;
	St--;
}
void cf4()
{
	Df++;
	Dt++;
}
void cf5()
{
	Dt--;
	Df--;
}

//反应概率
double pf1()
{
	return enta(T);
}
double pf2()
{
	return St*a0;
}
double pf3()
{
	return SDF*as;
}
double pf4()
{
	return Sf*Kd;
}
double pf5()
{
	return Dt*ad;
}

//初始化
void init()
{
	id_parameter="Time\tSt\tSf\tDt\tDf\tFt\tFf\tPt\tSD\tUD\tSDF\n";

	parameter[0]=&St;
	parameter[1]=&Sf;
	parameter[2]=&Dt;
	parameter[3]=&Df;
	parameter[4]=&Ft;
	parameter[5]=&Ff;
	parameter[6]=&Pt;
	parameter[7]=&SD;
	parameter[8]=&UD;
	parameter[9]=&SDF;

	p_channel[1]=cf1;
	p_channel[2]=cf2;
	p_channel[3]=cf3;
	p_channel[4]=cf4;
	p_channel[5]=cf5;

	poss[1]=pf1;
	poss[2]=pf2;
	poss[3]=pf3;
	poss[4]=pf4;
	poss[5]=pf5;
}
//快反应平衡控制函数
void equilibrium()
{
	Df=K(T)*K(T)+2*Ku*K(T)*K(T)+2*Kfold*Ku*Dt*K(T)+Ku*Ku*K(T)*K(T)*Pt*Pt-2*Ku*Ku*K(T)*Pt*Dt*Kfold+Dt*Dt*Kfold*Kfold*Ku*Ku;
	Df=pow(Df,0.5);
	Df+=Dt*Kfold*Ku-K(T)-Ku*K(T)*Pt;
	Df/=2*Kfold*Ku;
	Sf=St/(1+Ks*Df*(1+Kf*a*Dt));
	SD=Ks*Sf*Df;
	SDF=Kf*a*Dt*SD;
	UD=Ku*Df*K(T)*Pt/(Kfold*Ku*Df+K(T));
}