傅立葉級數係數​​gnuplot

Question

這個程序應該找到傅里葉級數的係數，然後繪製從計算值建立的函數。

分段爲噸 H（T）= 1-（T/PI）對於t> 0

#include <stdio.h> 
#include <math.h> 
#include "comphys.c" 
#include "comphys.h" 
#define pi 3.141592653589793 

double trap1 (double l, double u, int m, int N); 
double trap(double l,double u,int m,int N); // Use trapezoidal method 
double inte(double x,int m); 
double integrand(double x, int m); // This is the integrand of the bms 
double w=1.0; //omega (here - w-2pi/Period where Period=2pi) 
int main() 
{ 
    FILE *out,*rsp; 
    int m; // Index of the b coefficients 
    int M=5; // Number of coeeficients to compute 
    int T=201; // Number of time steps 
    int t; // This is a counter to keep track of the h(t) function to plot 
    double *b,*h,j,dj,*a,a0; 
    double ll,ul; 
    a=dvector(1,M); 
    b=dvector(1,M); // coefficients 
    h=dvector(1,T); // function built from Fourier terms 
    ll=-pi; // lower time limit 
    ul=pi; // upper time limit 

    // Prepare to write for plotting 
    printf("\nBrute force complex Fourier transform algorithm"); 
    if ((out=fopen("fs.out","w"))==NULL) { 
    printf("\nCannot open file for output\n"); 
    getchar(); 
    return(0);   
    } 

    // Loop through the wave numbers  
    for (m=1;m<=M;m++) 
    { 
     b[m]=(1.0/pi)*trap(ll,ul,m,T); 
     a[m]=(1.0/pi)*trap1(ll,ul,m,T); 
     printf("\nb[%d]=%f a[%d]=%f",m,b[m],m,a[m]); 

     // Build the h(t) function as we go 
     t=0; 
     dj=(ul-ll)/(double)T; 
     for (j=ll;j<=ul;j+=dj) 
    { 
     t++; 
     h[t]+=(b[m]*sin((double)m*w*j))+(a[m]*cos((double)m*w*j)); 
    } 
    } 

    // Write out for plotting 
    t=0; 
    for (j=ll;j<=ul;j+=dj) 
    { 
     t++; 
     fprintf(out,"%f %f\n",j,h[t]); 
    } 
    fclose(out); 
    printf("\nProgram complete without known error.\n"); 

    if((rsp = fopen("gnuplot.rsp","w")) == NULL) 
    { 
     printf("\nCannot open gnuplot.rsp for writing\n"); 
     exit(1); 
    } 
    fprintf(rsp,"plot '%s' using 1:2 with lines\n",out); 
    fprintf(rsp,"pause mouse\n"); 
    fprintf(rsp,"replot\n"); 
    fclose(rsp); 
    if(system("gnuplot gnuplot.rsp") == -1) 
    { 
     printf("\nCommand could not be executed\n"); 
     exit(1); 
    } 
    return; 
} 

double trap (double l, double u, int m, int N) 

{ 
    double x,integral,dx; 
    int i; 

    if (u==l) return (0.0); 
    dx=(u-l)/(double)(N-1); 

    integral=0.0; 
    integral+=0.5*(integrand(u,m)-integrand(l,m)); 
    x=l+dx; 
    for (i=1;i<N;i++) 
    { 
     integral += integrand(x,m); 
     x+=dx; 
    } 
    integral*=dx; 
    printf("\n%f\n",integral); 
    return(integral); 
} 

double trap1 (double l, double u, int m, int N) 
{ 
    double x,integral,dx; 
    int i; 

    if (u==l) return (0.0); 
    dx=(u-l)/(double)(N-1); 

    integral=0.0; 
    integral+=0.5*(inte(u,m)-inte(l,m)); 
    x=l+dx; 
    for (i=1;i<N;i++) 
    { 
     integral += inte(x,m); 
     x+=dx; 
    } 
    integral*=dx; 
    return(integral); 
} 


double integrand(double x, int m) 
{ 
    if (x<0) return (sin((double)m*w*x)+((x/pi)*sin((double)m*w*x))); 
    return(sin((double)m*w*x)-((x/pi)*sin((double)m*w*x))); 
} 

double inte(double x,int m) 
{ 
    if (x<0) return (cos((double)m*w*x)+((x/pi)*cos((double)m*w*x))); 
    return(cos((double)m*w*x)-((x/pi)*cos((double)m*w*x))); 
} 

H（噸）= 1 +（T/PI）

我有的問題是，它編譯罰款和列出係數（不知道，如果他們是正確的或任何東西），但最後當試圖繪製它給我

「不能讀取數據文件」？， 「 」gnuplot.rsp「，行：1 util.c：沒有這樣的文件或目錄」

而我不知道如何解決這個問題。我用gnuplot繪製了大量的程序，所以它不應該是一個丟失的文件或任何東西。

Answer 1

這裏：

fprintf(rsp,"plot '%s' using 1:2 with lines\n",out); 

你應該直接使用字符串「fs.out」，而不是「出」，這是文件指針（也已關閉）。