// PROGRAM synthesize

#include "TrueBASIC.h"
int main();
void plotf(float xmin, float xmax, float ymin, float ymax);
void fourier(int N);
void coefficients(int N, double a[], double b[]);

int main()
{
    GWopen(0);

    plotf(0.0f, 0.5f, 0.5f, 1.0f);
    plotf(0.0f, 0.5f, 0.0f, 0.5f);
    plotf(0.5f, 1.0f, 0.5f, 1.0f);
    plotf(0.5f, 1.0f, 0.0f, 0.5f);

    GWquit();
    return 0;
}

void plotf(float xmin, float xmax, float ymin, float ymax)
{
    int N;
    char Stmp1_[_LBUFF_];

    GWvport(xmin*GWaspect(1), ymin, xmax*GWaspect(1), ymax);
    GWindow(-4.0f, -2.0f, 4.0f, 2.0f);
    GWsetpen(BLACK, -1, -1, -1);
    GWline((float)(-pi), 0.0f, (float)(pi), 0.0f);
    GWline(0.0f, -1.5f, 0.0f, 1.5f);
    GWinput("number of modes = ", Stmp1_, _LBUFF_);
    sscanf(Stmp1_, "%d", &N);
    GWsetpen(-1, -1, -1, -1);
    GWsetpen(RED, -1, -1, -1);
    sprintf(Stmp1_, "number of modes = %d", N);
    GWputtxt((float)(-pi), 1.5f, Stmp1_);
    fourier(N);
}

void fourier(int N)
// フーリエ級数の計算と関数の表示
{
    int k, nplot;
    double a[1001], b[1001], dt, f, t;

    coefficients(N, a, b);
    nplot = 100;
    t = -pi;
    dt = pi/nplot;
    while(t <= pi)
    {
        f = a[0]/2;
        for(k = 1; k <= N; ++k)
        {
            if(a[k])
                f += a[k]*cos(k*t);
            if(b[k])
                f += b[k]*sin(k*t);
        }
        GWline2((float)(t), (float)(f));
        t += dt;
    }
}

void coefficients(int N, double a[], double b[])
// 特別な場合のフーリエ係数の生成
{
    int k;

    a[0] = 0;
    for(k = 1; k <= N; ++k)
    {
        a[k] = 0;
        if(k % 2)
            b[k] = 2/(k*pi);
        else
            b[k] = 0;
    }
}