c - 我有一个多变量函数，我想将其集成到 C 中(使用梯形规则)

Question

我必须针对 E 进行积分，其下限和上限作为我为每个源打印的 log10(E) 列的第一个值，上限作为同一列的最后一个值。我的函数是:fo E^-spectralindexexp(-tau1)，这只是我在这里定义的 CorrectedDifflux 。 [此函数针对 E]。为了简单起见，这里我将“v”和“w”作为下限和上限，但是我在这里执行的积分并没有为这些提供正确的值。我的代码:

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
float ebl(float zz, float ee) // Function ebl assigned for returning ta3 
value to be used later
{



float e1, e2, z1, z2, tau_a, tau_b, ta1, ta2, ta3, tau[500], z[500], e[500];
  int i, j, ii, jj;
float tauM[501][500];
FILE *fp;

 fp=fopen("tau_dominguez11_cta.txt", "r"); // Reading a file for creating a 
 matrix which has redshift and energy values

 for(i=0; i<501; i++)
   {
    for(j=0; j<400; j++)
   {
 fscanf(fp, "%f", &tauM[i][j]);

     }

   }
 fclose(fp);

 for(j=0; j<400; j++)
  {
  if(tauM[0][j]<=zz)
   {jj=j;}
 }

for(i=0; i<501; i++)
{

  if(tauM[i][0]<=ee)
  {ii=i;}
 }


 ta1=(zz-tauM[0][jj])*(tauM[ii][jj+1] -tauM[ii][jj])/(tauM[0][jj+1]-tauM[0]
 [jj]) + tauM[ii][jj]; //across zz

  ta2=(zz-tauM[0][jj])*(tauM[ii+1][jj+1] - tauM[ii+1][jj]) / (tauM[0][jj+1] 
 - tauM[0][jj]) + tauM[ii+1][jj]; //across zz

 ta3= (ta2 - ta1)*(ee-tauM[ii][0]) /(tauM[ii+1][0] - tauM[ii][0]) + ta1; 
 //across ee

 return ta3;
 }



float t(float E) // Defining to carry out the integration (help needed in 
this part) 
{
float fo, spectralindex, tau1;
return log10(fo * pow(E, -1*spectralindex) * exp(-1*tau1));
}

int main()
{
 FILE *fp1, *fp2,*files;
 int j,i,k;
 float fo;
 char source;
 float a,spectralindex,flux,E1,E2,E3,E4,x,logE,E,z,red,tau1;
 double diffflux,correctedDiffflux;
 float v=0.1, w=10, N=100,sum=0, h, p;

 E1=0.1, E2=100;
 E3=0.01, E4=10;
 float minLogE=log10(E3), maxlogE=log10(E4);

  int n=100;


fp1=fopen("actualvalues.txt", "r"); // reading flux and spectral index values for 3 sources




 for(j=0; j<3; j++)
{

 char filename[300];
 sprintf(filename,"source%d.txt",j);


  files=fopen(filename,"w");


   fscanf(fp1, "%s %e %f %f\n", &source,  &flux, &spectralindex,&red);

  a=-1*(spectralindex)+1;

  fo=(flux*a)/(pow(E2, a) - pow(E1, a));


   x=(log10(E4/E3))/n;
   logE=minLogE;
  for(i=0; i<n; i++)

   {


    logE=logE+x;

   E=pow(10,logE);



   diffflux=fo*(pow(E,-1*(spectralindex)));

  tau1=ebl(red,E);
   correctedDiffflux=diffflux*exp(-1*tau1);

  h = (w-v) / N;
   p = t(v);

   for(k=1; k<N; k++)
  {
     p=p+2*t(v+k*h);
   }
   p = p+t(w);
   p = p*(h/2);  // p is the integration variable

     fprintf(files,"%f    %f    %f    
      %f\n",log10(E),log10(diffflux),log10(correctedDiffflux),p);

      }
     }
       fclose(fp1);
      fclose(fp2);
    fclose(files);
      }

Answer 1

在此功能中，您评论说您需要一些帮助，

float t(float E) // Defining to carry out the integration (help needed in this part) 
{
    float fo, spectralindex, tau1;
    return log10(fo * pow(E, -1*spectralindex) * exp(-1*tau1));
}

这三个局部变量未初始化。它们与其他函数中的类似变量没有任何联系。

如果您启用编译器警告，这应该会通知您。

此外，我建议您停止使用 float，转而使用 double，除非有充分的理由不这样做。二十世纪教科书？