home *** CD-ROM | disk | FTP | other *** search
/ QBasic & Borland Pascal & C / Delphi5.iso / C / Samples / C-SSP.ARJ / FORIE.C < prev    next >
Encoding:
Text File  |  1984-08-30  |  2.4 KB  |  97 lines
  1.    forie(n,nk,r,k,e,a1,b1)
  2.  
  3.      int k[],n,nk;
  4.      float r[],e[],a1[],b1[];
  5.  
  6.    {
  7.  
  8.      int h,i,ix,j,l,kmax,m;
  9.      float a,a0,b,d,sum,suma,sumb;
  10.      float t[200],s[200],c[200],f[200];
  11.      extern double cos(),sin();
  12.  
  13.       ix = 1;
  14.       kmax = k[nk-1];
  15.       if(n <= 0 || n > 200) goto e2;
  16.       if(kmax > (n/2)) goto e3;
  17.       if(k[0] != 0) goto e1;
  18.       for(i = 0; i <= n-1; i++)
  19.       {
  20.        a1[i] = 0.;
  21.        b1[i] = 0.;
  22.       }
  23.       l = 1;
  24.       i = n - 1;
  25.       d = i;
  26. reit: m = k[l-1];
  27.       sumb = 0.;
  28.       if(m <  0) goto e4;
  29.       if (m <= 0)            /* compute average level */
  30.       {
  31.        sum = (r[0]+r[n-1])/2.;
  32.        for(j = 1; j <= i-1; j++)
  33.         sum = sum + r[j];
  34.        a = sum/d;
  35.        a0 = a;
  36.        b = 0.;
  37.       }
  38.        else
  39.        {
  40.         h = m;          /* compute sine-cosine */
  41.         s[1] = sin((6.2831854/d) * h);
  42.         c[1] = cos((6.2831854/d) * h);
  43.  
  44.         for(j = 2; j <= i-1; j++)
  45.         {
  46.          s[j] = s[j-1] * c[1] + c[j-1] * s[1];
  47.          c[j] = c[j-1] * c[1] - s[j-1] * s[1];
  48.         }
  49.         suma = (r[0] + r[n-1] -2. * a0)/2.;
  50.         for(j = 1; j <= i-1; j++)
  51.         {
  52.          suma = suma + (r[j]-a0) * c[j];
  53.          sumb = sumb + (r[j]-a0) * s[j];
  54.         }
  55.         a = (2.0 * suma)/d;
  56.         b = (2.0 * sumb)/d;
  57.                                         /*compute t vector*/
  58.         for(j = 1; j <= i-1; j++)
  59.          t[j] = a * c[j] + b * s[j];
  60.  
  61.         t[0] = a;
  62.         t[n-1] = a;
  63.        }
  64.       a1[ix-1] = a;
  65.       b1[ix-1] = b;
  66.       ix = ix + 1;
  67.  
  68.       if(m <= 0)                         /*compute f*/
  69.       {
  70.        for(j = 0; j <= n-1; j++)
  71.          f[j] = a0;
  72.       }
  73.        else
  74.        {
  75.         for(j = 0; j <= n-1; j++)
  76.          f[j] = f[j] + t[j];
  77.        }
  78.                                        /*increment harmonic index*/
  79.       if (m != kmax)
  80.       {
  81.        l = l + 1;
  82.        goto reit;
  83.       }
  84.                                        /*compute error*/
  85.       for(j = 0; j <= n-1; j++)
  86.        e[j] = r[j] - f[j];
  87.       return;
  88.  
  89.  e1:  printf("error in -Forie...first harmonic index not 0\n");
  90.       return;
  91.  e2:  printf("error in -Forie...illegal data point total\n");
  92.       return;
  93.  e3:  printf("error in -Forie...high harmonic index too large\n");
  94.       return;
  95.  e4:  printf("error in -Forie...illegal harmonic index\n");
  96.     }
  97.