关于奇偶震荡项的一个代码，不知道它哪里出了问题

Question

```#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:[e-mail]       *\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")

Answer 1

下面是修改编译错误后的代码和执行结果。
以只读方式打开文件，若文件不存在，则不会新建文件，因此此处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:[e-mail]       *\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);
}