📄 求功率.txt
字号:
#include <iostream.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
double gr(double);
double maxD();
double D[8],X[8];
double S[8][8];
double L=0.1,Len=40.0,Keff=2.0,Aeff=50*1e-12,ak=0.3;
double pw[9]; /*w*/
double v[9]; //HZ
double I[9];
double wl[9]; //nm
double Istar[9]; //w*m
double maxd;
int N;
void Ifun(); //求功率积分
void Dfun(); //打靶结果偏差向量
void Sfun(); //雅可比行列式
int Xfun(); //打靶调整因子
double maxD(); //判断打靶的结束
typedef struct node
{
double data; //存储元素的值
int x; //存储元素的横坐标
int y; //存储元素的纵坐标
} array;
double sum; //全局变量，存储行列式的值
void Create(double H[][8], double X[]); //构造一个行列式
void Solve(const double H[][8], array S[], int i, int NiXu);//采用递归方式求行列式的值
bool Judge(const array S[], int line, int len); //判断行列式的元素的纵坐标是否重复
void main()
{
int i;
wl[1]=1409.5;wl[2]=1416.2;wl[3]=1431.4;wl[4]=1449.8;wl[5]=1485.9;wl[6]=1487.8;wl[7]=1489.8;wl[8]=1491;//nm
Istar[1]=1206.66;Istar[2]=748.308;Istar[3]=576.977;Istar[4]=475.042;Istar[5]=236.958;Istar[6]=285.721;Istar[7]=224.716;Istar[8]=384.069;
N=int(Len/L);
for(i=1;i<9;i++)
{
v[i]=3e5/wl[i];
}
for(i=1;i<9;i++)
pw[i]=(Istar[i]*ak/4.343)/((1-exp(-ak*Len/4.343))*1000); //只考虑衰减的输入泵浦功率
/*pw[1]=0.084;pw[2]=0.050;pw[3]=0.045;pw[4]=0.030;pw[5]=0.020;pw[6]=0.023;pw[7]=0.018;pw[8]=0.025;*/
for(i=1;i<9;i++)
cout<<pw[i]<<" pw "<<endl;
// *****************************
Ifun(); //第一次打靶的功率积分
Dfun(); //打靶结果偏差向量
maxD(); //打靶结果最大偏差
cout<<maxd<<"********maxd"<<endl;
while(maxd>0.01)
{
Sfun(); //雅可比行列式
/*for(i=0;i<8;i++)
{for(int j=0;j<8;j++)
{cout<<S[i][j]<<' ';}
cout<<endl;}*/
Xfun(); //打靶调整因子
for(i=0;i<8;i++)
cout<<X[i]<<"////X[i]"<<endl;
for(i=1;i<9;i++)
pw[i]=pw[i]+X[i-1]; //下一次打靶输入功率
Ifun();
Dfun();
maxD();
cout<<maxd<<"********maxd"<<endl;
};
for(i=1;i<9;i++)
cout<<pw[i]<<" pw"<<endl;
// ******************
}
void Ifun() //第一次打靶功率积分
{
double p[9][4000],id[9][4000];
for(int m=1;m<9;m++)
I[m]=0.0;
for(int i=1;i<9;i++)
{
p[i][0]=pw[i];
}
for(int j=0;j<N-1;j++)
{
p[1][j+1]=p[1][j]*exp(-ak*L/4.343)*exp(-gr(v[1]-v[2])*v[1]*p[2][j]*L/(v[2]*Keff*Aeff)-gr(v[1]-v[3])*v[1]*p[3][j]*L/(v[3]*Keff*Aeff)\
-gr(v[1]-v[4])*v[1]*p[4][j]*L/(v[4]*Keff*Aeff)-gr(v[1]-v[5])*v[1]*p[5][j]*L/(v[5]*Keff*Aeff)\
-gr(v[1]-v[6])*v[1]*p[6][j]*L/(v[6]*Keff*Aeff)-gr(v[1]-v[7])*v[1]*p[7][j]*L/(v[7]*Keff*Aeff)\
-gr(v[1]-v[8])*v[1]*p[8][j]*L/(v[8]*Keff*Aeff));
p[2][j+1]=p[2][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[2])*p[1][j]*L/(Keff*Aeff))*exp(-gr(v[2]-v[3])*v[2]*p[3][j]*L/(v[3]*Keff*Aeff)\
-gr(v[2]-v[4])*v[2]*p[4][j]*L/(v[4]*Keff*Aeff)-gr(v[2]-v[5])*v[2]*p[5][j]*L/(v[5]*Keff*Aeff)\
-gr(v[2]-v[6])*v[2]*p[6][j]*L/(v[6]*Keff*Aeff)-gr(v[2]-v[7])*v[2]*p[7][j]*L/(v[7]*Keff*Aeff)\
-gr(v[2]-v[8])*v[2]*p[8][j]*L/(v[8]*Keff*Aeff));
p[3][j+1]=p[3][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[3])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[3])*p[2][j]*L/(Keff*Aeff))*\
exp(-gr(v[3]-v[4])*v[3]*p[4][j]*L/(v[4]*Keff*Aeff)-gr(v[2]-v[5])*v[2]*p[5][j]*L/(v[5]*Keff*Aeff)\
-gr(v[2]-v[6])*v[2]*p[6][j]*L/(v[6]*Keff*Aeff)-gr(v[2]-v[7])*v[2]*p[7][j]*L/(v[7]*Keff*Aeff)\
-gr(v[2]-v[8])*v[2]*p[8][j]*L/(v[8]*Keff*Aeff));
p[4][j+1]=p[4][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[4])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[4])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[4])*p[3][j]*L/(Keff*Aeff))*exp(-gr(v[4]-v[5])*v[4]*p[5][j]*L/(v[5]*Keff*Aeff)-gr(v[4]-v[6])*v[4]*p[6][j]*L/(v[6]*Keff*Aeff)\
-gr(v[4]-v[6])*v[4]*p[6][j]*L/(v[6]*Keff*Aeff)-gr(v[4]-v[7])*v[4]*p[7][j]*L/(v[7]*Keff*Aeff)\
-gr(v[4]-v[8])*v[4]*p[8][j]*L/(v[8]*Keff*Aeff));
p[5][j+1]=p[5][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[5])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[5])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[5])*p[3][j]*L/(Keff*Aeff)+gr(v[4]-v[5])*p[4][j]*L/(Keff*Aeff))*exp(-gr(v[5]-v[6])*v[5]*p[6][j]*L/(v[6]*Keff*Aeff)\
-gr(v[5]-v[7])*v[5]*p[7][j]*L/(v[7]*Keff*Aeff)-gr(v[5]-v[8])*v[5]*p[8][j]*L/(v[8]*Keff*Aeff));
p[6][j+1]=p[6][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[6])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[6])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[6])*p[3][j]*L/(Keff*Aeff)+gr(v[4]-v[6])*p[4][j]*L/(Keff*Aeff)+gr(v[5]-v[6])*p[5][j]*L/(Keff*Aeff))\
*exp(-gr(v[6]-v[7])*v[6]*p[7][j]*L/(v[7]*Keff*Aeff)-gr(v[6]-v[8])*v[6]*p[8][j]*L/(v[8]*Keff*Aeff));
p[7][j+1]=p[7][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[7])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[7])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[7])*p[3][j]*L/(Keff*Aeff)+gr(v[4]-v[7])*p[4][j]*L/(Keff*Aeff)+gr(v[5]-v[7])*p[5][j]*L/(Keff*Aeff)\
+gr(v[6]-v[7])*p[6][j]*L/(Keff*Aeff))*exp(-gr(v[7]-v[8])*v[7]*p[8][j]*L/(v[8]*Keff*Aeff));
p[8][j+1]=p[8][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[8])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[8])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[8])*p[3][j]*L/(Keff*Aeff)+gr(v[4]-v[8])*p[4][j]*L/(Keff*Aeff)+gr(v[5]-v[8])*p[5][j]*L/(Keff*Aeff)\
+gr(v[6]-v[8])*p[6][j]*L/(Keff*Aeff)+gr(v[7]-v[8])*p[7][j]*L/(Keff*Aeff));
}
for(i=1;i<9;i++)
for(j=0;j<N;j++)
{
id[i][j]=p[i][j]*L*1000;
}
for(i=1;i<9;i++)
for(j=0;j<N;j++)
I[i]+=id[i][j];
}
double gr(double y) //拉曼增益系数
{
double x=fabs(y);
double gr;
if((x>0)&&(x<=0.5))
gr=0+0.34*x;
else if((x>0.5)&&(x<=11.7))
gr=-0.56708+1.60247*x-0.39691*x*x+0.04341*x*x*x-0.00147*x*x*x*x;
else if((x>11.7)&&(x<=13.9))
gr=5.53278+2.3999/(1.79291*sqrt(3.14/2))*exp(-2*(x-13.14066)*(x-13.14066)/(1.79291*1.79291));
else if((x>13.9)&&(x<=14.4))
gr=194.37371-26.77718*x+0.95297*x*x;
else if((x>14.4)&&(x<=16.0))
gr=2.1840+6.5874*exp(-2*(x-14.491)*(x-14.491)/(1.2349*1.2349))/(1.2349*sqrt(3.14/2));
else if((x>16.0)&&(x<=17.6))
gr=175.6676-20.52019*x+0.6056*x*x;
else if((x>17.6)&&(x<=20.0))
gr=0.87774+1.9890*exp(-2*(x-17.874)*(x-17.874)/(1.1425*1.1425))/(1.1425*sqrt(3.14/2));
else if((x>20.0)&&(x<24.0))
gr=0.02375*x*x-1.0475*x+12.31;
else cout<<"Out of region";
return gr*1e-14;
}
void Dfun() //打靶结果偏差向量
{
for(int i=1;i<9;i++)
cout<<(D[i-1]=I[i]-Istar[i])<<" ***D"<<endl;
}
void Sfun() //雅可比行列式
{
double p[9][4000],id[9][4000];
for(int m=1;m<9;m++)
{
pw[m]+=0.001; /*第一个泵浦功率增加1mw*/
for(int i=1;i<9;i++)
{
p[i][0]=pw[i];
}
for(int j=0;j<N-1;j++)
{
p[1][j+1]=p[1][j]*exp(-ak*L/4.343)*exp(-gr(v[1]-v[2])*v[1]*p[2][j]*L/(v[2]*Keff*Aeff)-gr(v[1]-v[3])*v[1]*p[3][j]*L/(v[3]*Keff*Aeff)\
-gr(v[1]-v[4])*v[1]*p[4][j]*L/(v[4]*Keff*Aeff)-gr(v[1]-v[5])*v[1]*p[5][j]*L/(v[5]*Keff*Aeff)\
-gr(v[1]-v[6])*v[1]*p[6][j]*L/(v[6]*Keff*Aeff)-gr(v[1]-v[7])*v[1]*p[7][j]*L/(v[7]*Keff*Aeff)\
-gr(v[1]-v[8])*v[1]*p[8][j]*L/(v[8]*Keff*Aeff));
p[2][j+1]=p[2][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[2])*p[1][j]*L/(Keff*Aeff))*exp(-gr(v[2]-v[3])*v[2]*p[3][j]*L/(v[3]*Keff*Aeff)\
-gr(v[2]-v[4])*v[2]*p[4][j]*L/(v[4]*Keff*Aeff)-gr(v[2]-v[5])*v[2]*p[5][j]*L/(v[5]*Keff*Aeff)\
-gr(v[2]-v[6])*v[2]*p[6][j]*L/(v[6]*Keff*Aeff)-gr(v[2]-v[7])*v[2]*p[7][j]*L/(v[7]*Keff*Aeff)\
-gr(v[2]-v[8])*v[2]*p[8][j]*L/(v[8]*Keff*Aeff));
p[3][j+1]=p[3][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[3])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[3])*p[2][j]*L/(Keff*Aeff))*\
exp(-gr(v[3]-v[4])*v[3]*p[4][j]*L/(v[4]*Keff*Aeff)-gr(v[2]-v[5])*v[2]*p[5][j]*L/(v[5]*Keff*Aeff)\
-gr(v[2]-v[6])*v[2]*p[6][j]*L/(v[6]*Keff*Aeff)-gr(v[2]-v[7])*v[2]*p[7][j]*L/(v[7]*Keff*Aeff)\
-gr(v[2]-v[8])*v[2]*p[8][j]*L/(v[8]*Keff*Aeff));
p[4][j+1]=p[4][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[4])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[4])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[4])*p[3][j]*L/(Keff*Aeff))*exp(-gr(v[4]-v[5])*v[4]*p[5][j]*L/(v[5]*Keff*Aeff)-gr(v[4]-v[6])*v[4]*p[6][j]*L/(v[6]*Keff*Aeff)\
-gr(v[4]-v[6])*v[4]*p[6][j]*L/(v[6]*Keff*Aeff)-gr(v[4]-v[7])*v[4]*p[7][j]*L/(v[7]*Keff*Aeff)\
-gr(v[4]-v[8])*v[4]*p[8][j]*L/(v[8]*Keff*Aeff));
p[5][j+1]=p[5][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[5])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[5])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[5])*p[3][j]*L/(Keff*Aeff)+gr(v[4]-v[5])*p[4][j]*L/(Keff*Aeff))*exp(-gr(v[5]-v[6])*v[5]*p[6][j]*L/(v[6]*Keff*Aeff)\
-gr(v[5]-v[7])*v[5]*p[7][j]*L/(v[7]*Keff*Aeff)-gr(v[5]-v[8])*v[5]*p[8][j]*L/(v[8]*Keff*Aeff));
p[6][j+1]=p[6][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[6])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[6])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[6])*p[3][j]*L/(Keff*Aeff)+gr(v[4]-v[6])*p[4][j]*L/(Keff*Aeff)+gr(v[5]-v[6])*p[5][j]*L/(Keff*Aeff))\
*exp(-gr(v[6]-v[7])*v[6]*p[7][j]*L/(v[7]*Keff*Aeff)-gr(v[6]-v[8])*v[6]*p[8][j]*L/(v[8]*Keff*Aeff));
p[7][j+1]=p[7][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[7])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[7])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[7])*p[3][j]*L/(Keff*Aeff)+gr(v[4]-v[7])*p[4][j]*L/(Keff*Aeff)+gr(v[5]-v[7])*p[5][j]*L/(Keff*Aeff)\
+gr(v[6]-v[7])*p[6][j]*L/(Keff*Aeff))*exp(-gr(v[7]-v[8])*v[7]*p[8][j]*L/(v[8]*Keff*Aeff));
p[8][j+1]=p[8][j]*exp(-ak*L/4.343)*exp(gr(v[1]-v[8])*p[1][j]*L/(Keff*Aeff)+gr(v[2]-v[8])*p[2][j]*L/(Keff*Aeff)\
+gr(v[3]-v[8])*p[3][j]*L/(Keff*Aeff)+gr(v[4]-v[8])*p[4][j]*L/(Keff*Aeff)+gr(v[5]-v[8])*p[5][j]*L/(Keff*Aeff)\
+gr(v[6]-v[8])*p[6][j]*L/(Keff*Aeff)+gr(v[7]-v[8])*p[7][j]*L/(Keff*Aeff));
}
double ist[9]={0,0,0,0,0,0,0,0,0};
for(i=1;i<9;i++)
for(j=0;j<N;j++)
{
id[i][j]=p[i][j]*L*1000;
}
for(i=1;i<9;i++)
for(j=0;j<N;j++)
ist[i]+=id[i][j];
/*for(i=1;i<9;i++)
cout<<ist[i]<<" I"<<m<<"+0.001"<<endl;*/
for(i=0;i<8;i++)
S[i][m-1]=(ist[i+1]-Istar[i+1]-D[i])/0.001; //矩阵S
/*S[i][m-1]=ist[i+1]-Istar[i+1]-D[i]; */
pw[m]-=0.001;
}
}
int Xfun(void) //打靶调整因子
{
array SL[8]; //栈，存储行列式的每一个乘积项的元素（因子）
double H[8][8], CH[8][8];//存储系数行列式
double B[8]; //存储常数项
double d; //分别存储系数行列式的值 和 方程组的解
int i, j, k;
Create(H,B); //构造一个行列式
/*for(i=0;i<8;i++)
{
for(j=0;j<8;j++)
{cout<<H[i][j]<<"///";}
cout<<B[i]<<endl;
}*/
sum = 0;
Solve(H, SL, 0, 0); //采用递归方式求行列式的值
d = sum;/*cout<<sum<<"*****";*/
//输出该行列式的值
for(k=0; k<8; k++)
{
for(i=0; i<8; i++) //复制行列式
{
for(j=0; j<8; j++)
CH[i][j] = H[i][j];
CH[i][k] = B[i]; //把系数行列式D中第k列的元素用方程组右端的常数项代替
}
sum = 0;
Solve(CH, SL, 0, 0); //采用递归方式求行列式的值
/*cout<<sum<<"++++++++";*/
X[k] = sum / d;
}
/*for(i=0; i<8; i++)
printf("X%d = %f\n", i+1, X[i]); //输出该行列式的值 */
return 0;
}
void Create(double H[][8], double x[])
{
int i, j;
for(i=0; i<8; i++)
{
for(j=0; j<8; j++)
{H[i][j]=S[i][j]; //输入系数项
/*cout<<H[i][j]<<"*";*/}//
x[i]=-D[i]; //输入常数项
/*cout<<x[i]<<endl;//*/
fflush(stdin);
}
}
void Solve(const double H[][8], array s[], int i, int NiXu)//采用递归方式求行列式的值
{
array CS[8]; //栈，存储S[]的拷贝
int j, k, top = i;
double mul; //存储每一个乘积项的值
int CNiXu; //累积每一个乘积项的逆序数
for(j=0; j<8; j++)
{
if(Judge(s, j, top))//如果当前元素的纵坐标不与栈中存储的元素重复，将其入栈
{
s[top].x = i;
s[top].y = j;
s[top].data = H[i][j];
CNiXu = NiXu; //把逆序数复制到CNiXu
for(k=0; k<top; k++)
{
if(j < s[k].y) //累积逆序数
CNiXu++;
}
for(k=0; k<=top; k++) //复制栈
CS[k] =s[k];
if(i<7) //如果未分析到该乘积项的最后一个元素，递归继续分析
Solve(H, CS, i+1, CNiXu);
else //否则计算该乘积项的值，并存储到栈中
{
for(mul=1, k=0; k<=top; k++)
mul *= s[k].data;
if(CNiXu%2==0) //如果逆序数为偶数，该乘积项为正
sum += mul;
else //否则为负
sum -= mul;
}
}
}
}
bool Judge(const array s[], int line, int len)
{
int i;
for(i=0; i<len; i++)
if(line == s[i].y)
return 0;
return 1;
}
double maxD()
{
maxd=fabs(D[0]);
for(int i=0;i<8;i++)
{
if(maxd<fabs(D[i]))
maxd=fabs(D[i]);
}
return maxd;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -