```#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* Fragmentation reaction CS code: FRACS V1.0 */
/* See details in paper: PRC 95, 034608 (2017) */
double fracs(double Ap, double Zp, double At, double Zt, double A, double Z, double Ep);
main(int argc, char *argv[])
{
double Ap, Zp, Ep, At, Zt, A, Z, result;
printf("\n");
printf("***************************************\n");
printf("* FRACS version 1.0 *\n");
printf("* Empirical parametrization *\n");
printf("* for projectile fragmentation *\n");
printf("* cross sections *\n");
printf("* PRC 95, 034608 (2017) *\n");
printf("* Email:meibo@impcas.ac.cn *\n");
printf("***************************************\n");
printf("\n");
//Add Input: Projectile, Projectile Energy (MeV/nucleon), and Target.
Ap=58;
Zp=28;
Ep=140;
At=181;
Zt=73;
//End of Input.
//Read nuclides from file and write CS output to file
FILE *ifp, *ofp, *ifsp, *ofsp;
char *mode = "r";
char outputFilename[] = "out.list";
int n;
ifp = fopen("in.list", mode);
if (ifp == NULL) {
fprintf(stderr, "Can't open input file in.list!\n");
exit(1);
}
ofp = fopen(outputFilename, "w");
if (ofp == NULL) {
fprintf(stderr, "Can't open output file %s!\n",
outputFilename);
exit(1);
}
while (fscanf(ifp, "%lf %lf", &A, &Z) != EOF) {
result = fracs(Ap, Zp, At, Zt, A, Z, Ep);
fprintf(ofp, " %f %f ++ %f %f -> %f %f %e b \n",Ap,Zp,At,Zt,A,Z,result);
}
}
/********************************************************************************
* Calculate cross section (in barn) for fragment (A,Z) produced by
* projectile fragmentation of projectile (Ap,Zp) on target (At,Zt) at energy Ep
*********************************************************************************/
double fracs(double Ap, double Zp, double At, double Zt, double A, double Z, double Ep)
{
double R,r_0,p,f_pt,f_mod_y,zprob,zbeta,dzbeta_p,delta,dq;
double f_mod,slope,z0,a2,expo_fac,offset,dzprob,u_p,u_n;
double expo,fract,yield_a,zbeta_p;
double cd[9],cp[6],cr[9],cn[3],cq[3],un[2],up[3];
double cs[4],cl[3],cy[3],cb[3],norm[2],result;
double A1, Z1, OES;
FILE *ifsp, *ofsp;
result = -999.0 ;
double At13, Ap13;
Ap13 = pow(Ap,0.33333);
At13 = pow(At,0.33333);
/************************* Constants of FRACS V1.0 *******************************/
cp[1] = 2.0 ; /* mass yield slope P */
cp[2] = 0.6 ;
cp[3] = 15.0 ;
cp[4] = 5.0E-5 ;
cp[5] = -0.01 ;
cs[1] = 0.272 ; /* scaling factor */
cs[2] = -1.60 ;
cs[3] = -8.0;
cs[4] = 4.0E-5;
cy[1] = 0.95 ; /* correction of mass yield */
cy[2] = 10.0 ;
cd[1] = -1.4E+00 ; /* Zprob shift Delta */
cd[2] = +3.447E-02 ;
cd[3] = +2.1353E-04 ;
cd[4] = +7.1350E+01 ;
cd[5] = -24.0 ;
cd[6] = 0.80 ;
cd[7] = +3.5E-05 ;
cd[8] = -1.5E-07 ;
cq[1] = -10.25 ; /* memory effect p-rich projectiles */
cq[2] = +10.25 ;
cn[1] = 0.400 ; /* memory effect n-rich projectiles */
cn[2] = 0.600 ;
cr[1] = 2.78 ; /* width parameter R */
cr[2] = -0.015 ;
cr[3] = 3.20E-05 ;
cr[4] = 0.0412 ;
cr[5] = 0.124 ;
cr[6] = 29.0;
cr[7] = 0.855 ;
cr[8] = -11 ;
cr[9] = -0.12 ;
un[1] = 1.26+(A/Ap)*0.6-Ap*0.0006; /* slope for n-rich side. */
up[1] = 2.1; /* slope parameter for p-rich side */
up[2] = -0.16;
cb[1] = 2.3e-03 ; /* brute force correction for very n-rich */
cb[2] = 2.4 ;
cl[1] = 1.20 ; /* transition point for p-rich side */
cl[2] = 0.647 ;
/*****************************************************************************
/* slope parameterization */
p = cp[1]/Ap+cp[2]/(At13*pow(Ap,0.5))+cp[3]/(Ep*At)+cp[4]*(log10(Ep/100))*At+cp[5];
/* calculate mass yield */
f_pt = cs[1]*p*(Ap13 + At13 + cs[2] + cs[3]*(Ap-2*Zp)*Zt/(Ap*At) + cs[4]*Ap*At*(log10(Ep/100))) ;
yield_a = f_pt * exp(-p * (Ap - A)) ;
/* Mass yield correction near projectile */
f_mod_y=1.0 ;
if (A/Ap > cy[1]) {
f_mod_y=1+cy[2]*Ap*pow(A/Ap-cy[1],2) ;
}
yield_a = yield_a * f_mod_y ;
/* calculate Zprob */
zbeta = A/(1.98+0.0155*pow(A,2./3.)) ;
zbeta_p = Ap/(1.98+0.0155*pow(Ap,2./3.)) ;
dzbeta_p = Zp - zbeta_p ;
if(A > cd[4]) {
delta = cd[1]+cd[2]*A+ cd[7]*A*A + cd[8]*A*A*A;
}
else
{
delta = cd[3]*A*A ;
}
f_mod = 1.0 ;
if(A/Ap > cd[6]) {
f_mod = cd[5] * pow(A/Ap-cd[6],2) +1.0 ;
}
delta = delta * f_mod ;
zprob = zbeta + delta ;
if((Zp-zbeta_p) > 0) { /* Additional correction for p-rich */
dq = exp(cq[1]+cq[2]*(A/Ap)) ;
}
else /* Additional correction for n-rich */
{
dq = A/Ap*(A/Ap*(cn[1]+A/Ap*(A/Ap*cn[2]))) ;
}
dzprob = dq*(Zp-zbeta_p) ;
zprob = zprob +dzprob;
/* calculate width parameter */
if ((Zp-zbeta_p)<0) { /* n-rich */
r_0 = cr[1]* exp(cr[4]*dzbeta_p) ;
}
else
{ /* p-rich */
r_0 = cr[1]* exp(cr[5]*dzbeta_p) ;
}
R = r_0 * exp(cr[2]*A + cr[3]*A*A + (Zp-zbeta_p-1.5)*(A/(4.9*Ap)+cr[9])) ;
f_mod = 1.0 ;
if (A/Ap > cr[7]) {
f_mod = exp((cr[6]+cr[8]*(Zp-zbeta_p-1.5))*sqrt(Ap) * pow(A/Ap-cr[7],3.2)) ;
}
R = R * f_mod ;
/* Exponent parameters for n and p rich */
u_n = un[1];
u_p = up[1]+up[2]*log(R);
if((zprob-Z) > 0) {
/* neutron-rich */
expo = -R * pow(fabs(zprob-Z),u_n) ;
fract = exp(expo) * sqrt(R/3.14159) ;
}
else {
//Transition to exponential slope for p-rich
slope = cl[1] + cl[2] * pow(A/2.,0.3) ;
z0 = zprob+( (cl[1]+cl[2] * pow(A/2.,0.3) ) *log(10.)/(2.*R)) ;
expo = -R * pow(fabs(zprob-Z),u_p) ;
fract = exp(expo) * sqrt(R/3.14159) ;
if( Z > z0 ) {
expo = -R * pow(fabs(zprob-z0),u_p) ;
a2 = exp(expo) * sqrt(R/3.14159) ;
fract = a2/pow( pow(10,slope),(Z-z0) ) ;
}
}
norm[1] = 1.0 ;
/* "brute force factor" for n-rich projectiles */
if((zbeta_p-Zp)>0) {
expo_fac = -cb[1] * pow(fabs(Zp-zbeta_p),1) ;
if ((zbeta-Z)>(cb[2]+Zp-zbeta_p)) {
offset = Zp - zbeta_p + cb[2] ;
norm[1]= pow(10.0,expo_fac * pow(zbeta-Z+offset,3));
}
else
{
norm[1]=1.0 ;
}
}
// Cross section calculations
result = norm[1] * fract * yield_a ;
//Read correction factors for the OES effect
ifsp = fopen("OES.dat", "r");
if (ifsp == NULL) {
fprintf(stderr, "Can't open input file OES.dat!\n");
exit(1);
}
//Correct CS with OES factors
while (fscanf(ifsp, "%lf %lf %lf", &A1, &Z1, &OES) != EOF) {
if(A1==A && Z1==Z) result=exp(log(result)+pow(-1,Z)*OES);
}
fclose(ifsp);
return(result);
}
![img](https://img-mid.csdnimg.cn/release/static/image/mid/ask/683622898056111.jpg "#left")
关于奇偶震荡项的一个代码,不知道它哪里出了问题
- 写回答
- 好问题 0 提建议
- 追加酬金
- 关注问题
- 邀请回答
-
2条回答 默认 最新
- bostonAlen 2022-04-26 00:22关注
下面是修改编译错误后的代码和执行结果。
以只读方式打开文件,若文件不存在,则不会新建文件,因此此处exit了#define _CRT_SECURE_NO_WARNINGS #include <stdio.h> #include <stdlib.h> #include <math.h> /* Fragmentation reaction CS code: FRACS V1.0 */ /* See details in paper: PRC 95, 034608 (2017) */ double fracs(double Ap, double Zp, double At, double Zt, double A, double Z, double Ep); int main(int argc, char *argv[]) { double Ap, Zp, Ep, At, Zt, A, Z, result; printf("\n"); printf("***************************************\n"); printf("* FRACS version 1.0 *\n"); printf("* Empirical parametrization *\n"); printf("* for projectile fragmentation *\n"); printf("* cross sections *\n"); printf("* PRC 95, 034608 (2017) *\n"); printf("* Email:meibo@impcas.ac.cn *\n"); printf("***************************************\n"); printf("\n"); //Add Input: Projectile, Projectile Energy (MeV/nucleon), and Target. Ap = 58; Zp = 28; Ep = 140; At = 181; Zt = 73; //End of Input. //Read nuclides from file and write CS output to file FILE *ifp, *ofp, *ifsp, *ofsp; const char *mode = "r"; char outputFilename[] = "out.list"; int n; ifp = fopen("in.list", mode); if (ifp == NULL) { fprintf(stderr, "Can't open input file in.list!\n"); exit(1); } ofp = fopen(outputFilename, "w"); if (ofp == NULL) { fprintf(stderr, "Can't open output file %s!\n", outputFilename); exit(1); } while (fscanf(ifp, "%lf %lf", &A, &Z) != EOF) { result = fracs(Ap, Zp, At, Zt, A, Z, Ep); fprintf(ofp, " %f %f ++ %f %f -> %f %f %e b \n", Ap, Zp, At, Zt, A, Z, result); } } /******************************************************************************** * Calculate cross section (in barn) for fragment (A,Z) produced by * projectile fragmentation of projectile (Ap,Zp) on target (At,Zt) at energy Ep *********************************************************************************/ double fracs(double Ap, double Zp, double At, double Zt, double A, double Z, double Ep) { double R, r_0, p, f_pt, f_mod_y, zprob, zbeta, dzbeta_p, delta, dq; double f_mod, slope, z0, a2, expo_fac, offset, dzprob, u_p, u_n; double expo, fract, yield_a, zbeta_p; double cd[9], cp[6], cr[9], cn[3], cq[3], un[2], up[3]; double cs[4], cl[3], cy[3], cb[3], norm[2], result; double A1, Z1, OES; FILE *ifsp, *ofsp; result = -999.0; double At13, Ap13; Ap13 = pow(Ap, 0.33333); At13 = pow(At, 0.33333); /************************* Constants of FRACS V1.0 *******************************/ cp[1] = 2.0; /* mass yield slope P */ cp[2] = 0.6; cp[3] = 15.0; cp[4] = 5.0E-5; cp[5] = -0.01; cs[1] = 0.272; /* scaling factor */ cs[2] = -1.60; cs[3] = -8.0; cs[4] = 4.0E-5; cy[1] = 0.95; /* correction of mass yield */ cy[2] = 10.0; cd[1] = -1.4E+00; /* Zprob shift Delta */ cd[2] = +3.447E-02; cd[3] = +2.1353E-04; cd[4] = +7.1350E+01; cd[5] = -24.0; cd[6] = 0.80; cd[7] = +3.5E-05; cd[8] = -1.5E-07; cq[1] = -10.25; /* memory effect p-rich projectiles */ cq[2] = +10.25; cn[1] = 0.400; /* memory effect n-rich projectiles */ cn[2] = 0.600; cr[1] = 2.78; /* width parameter R */ cr[2] = -0.015; cr[3] = 3.20E-05; cr[4] = 0.0412; cr[5] = 0.124; cr[6] = 29.0; cr[7] = 0.855; cr[8] = -11; cr[9] = -0.12; un[1] = 1.26 + (A / Ap)*0.6 - Ap * 0.0006; /* slope for n-rich side. */ up[1] = 2.1; /* slope parameter for p-rich side */ up[2] = -0.16; cb[1] = 2.3e-03; /* brute force correction for very n-rich */ cb[2] = 2.4; cl[1] = 1.20; /* transition point for p-rich side */ cl[2] = 0.647; /***************************************************************************** /* slope parameterization */ p = cp[1] / Ap + cp[2] / (At13*pow(Ap, 0.5)) + cp[3] / (Ep*At) + cp[4] * (log10(Ep / 100))*At + cp[5]; /* calculate mass yield */ f_pt = cs[1] * p*(Ap13 + At13 + cs[2] + cs[3] * (Ap - 2 * Zp)*Zt / (Ap*At) + cs[4] * Ap*At*(log10(Ep / 100))); yield_a = f_pt * exp(-p * (Ap - A)); /* Mass yield correction near projectile */ f_mod_y = 1.0; if (A / Ap > cy[1]) { f_mod_y = 1 + cy[2] * Ap*pow(A / Ap - cy[1], 2); } yield_a = yield_a * f_mod_y; /* calculate Zprob */ zbeta = A / (1.98 + 0.0155*pow(A, 2. / 3.)); zbeta_p = Ap / (1.98 + 0.0155*pow(Ap, 2. / 3.)); dzbeta_p = Zp - zbeta_p; if (A > cd[4]) { delta = cd[1] + cd[2] * A + cd[7] * A*A + cd[8] * A*A*A; } else { delta = cd[3] * A*A; } f_mod = 1.0; if (A / Ap > cd[6]) { f_mod = cd[5] * pow(A / Ap - cd[6], 2) + 1.0; } delta = delta * f_mod; zprob = zbeta + delta; if ((Zp - zbeta_p) > 0) { /* Additional correction for p-rich */ dq = exp(cq[1] + cq[2] * (A / Ap)); } else /* Additional correction for n-rich */ { dq = A / Ap * (A / Ap * (cn[1] + A / Ap * (A / Ap * cn[2]))); } dzprob = dq * (Zp - zbeta_p); zprob = zprob + dzprob; /* calculate width parameter */ if ((Zp - zbeta_p) < 0) { /* n-rich */ r_0 = cr[1] * exp(cr[4] * dzbeta_p); } else { /* p-rich */ r_0 = cr[1] * exp(cr[5] * dzbeta_p); } R = r_0 * exp(cr[2] * A + cr[3] * A*A + (Zp - zbeta_p - 1.5)*(A / (4.9*Ap) + cr[9])); f_mod = 1.0; if (A / Ap > cr[7]) { f_mod = exp((cr[6] + cr[8] * (Zp - zbeta_p - 1.5))*sqrt(Ap) * pow(A / Ap - cr[7], 3.2)); } R = R * f_mod; /* Exponent parameters for n and p rich */ u_n = un[1]; u_p = up[1] + up[2] * log(R); if ((zprob - Z) > 0) { /* neutron-rich */ expo = -R * pow(fabs(zprob - Z), u_n); fract = exp(expo) * sqrt(R / 3.14159); } else { //Transition to exponential slope for p-rich slope = cl[1] + cl[2] * pow(A / 2., 0.3); z0 = zprob + ((cl[1] + cl[2] * pow(A / 2., 0.3)) *log(10.) / (2.*R)); expo = -R * pow(fabs(zprob - Z), u_p); fract = exp(expo) * sqrt(R / 3.14159); if (Z > z0) { expo = -R * pow(fabs(zprob - z0), u_p); a2 = exp(expo) * sqrt(R / 3.14159); fract = a2 / pow(pow(10, slope), (Z - z0)); } } norm[1] = 1.0; /* "brute force factor" for n-rich projectiles */ if ((zbeta_p - Zp) > 0) { expo_fac = -cb[1] * pow(fabs(Zp - zbeta_p), 1); if ((zbeta - Z) > (cb[2] + Zp - zbeta_p)) { offset = Zp - zbeta_p + cb[2]; norm[1] = pow(10.0, expo_fac * pow(zbeta - Z + offset, 3)); } else { norm[1] = 1.0; } } // Cross section calculations result = norm[1] * fract * yield_a; //Read correction factors for the OES effect ifsp = fopen("OES.dat", "r"); if (ifsp == NULL) { fprintf(stderr, "Can't open input file OES.dat!\n"); exit(1); } //Correct CS with OES factors while (fscanf(ifsp, "%lf %lf %lf", &A1, &Z1, &OES) != EOF) { if (A1 == A && Z1 == Z) result = exp(log(result) + pow(-1, Z)*OES); } fclose(ifsp); return(result); }
本回答被题主选为最佳回答 , 对您是否有帮助呢?解决 2无用
悬赏问题
- ¥100 set_link_state
- ¥15 虚幻5 UE美术毛发渲染
- ¥15 CVRP 图论 物流运输优化
- ¥15 Tableau online 嵌入ppt失败
- ¥100 支付宝网页转账系统不识别账号
- ¥15 基于单片机的靶位控制系统
- ¥15 真我手机蓝牙传输进度消息被关闭了,怎么打开?(关键词-消息通知)
- ¥15 装 pytorch 的时候出了好多问题,遇到这种情况怎么处理?
- ¥20 IOS游览器某宝手机网页版自动立即购买JavaScript脚本
- ¥15 手机接入宽带网线,如何释放宽带全部速度