// PROGRAM analyze
// determine the Fourier coefficients a_k and b_k

#include "TrueBASIC.h"
int main();
void parameters(int *N, int *nmax, double *delta, double *period);
void screen(int nmax, double period, int IC1_, int IC2_);
void plotaxis(int nmax, double ticksize);
void coefficents(int N, int nmax, double delta, double period, int IC1_, 
                 int IC2_);
double f(double t);

int main()
{
    int N, nmax, IC1_ = 1, IC2_ = 2;
    double delta, period;

    GWopen(0);

    parameters(&N, &nmax, &delta, &period);
    screen(nmax, period, IC1_, IC2_);
    coefficents(N, nmax, delta, period, IC1_, IC2_);

    GWquit();
    return 0;
}

void parameters(int *N, int *nmax, double *delta, double *period)
{
    char Stmp1_[_LBUFF_];

    printf("number of data points N (even) = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%d", N);
    printf("sampling time dt = ");
    fgets(Stmp1_, _LBUFF_, stdin);
    sscanf(Stmp1_, "%lg", delta);
    (*period) = (*N)*(*delta);              // assumed period
    // maximum value of mode corresponding to Nyquist frequency
    (*nmax) = (*N)/2;
}

void screen(int nmax, double period, int IC1_, int IC2_)
{
    double ticksize, ymax;
    char Stmp_[_LBUFF_];

    ymax = 2;
    ticksize = ymax/50;
    GWvport(0.0f, 0.5f, GWaspect(1), 1.0f);
    sprintf(Stmp_, "     a_k          frequency interval = %7.5f", 2*pi/period);
    GWindow(-1.0f, (float)(-ymax), (float)(nmax+1), (float)(ymax));
    GWputtxt(0.0f, (float)(ymax*0.7), Stmp_);
    plotaxis(nmax, ticksize);
    GWsetpen(RED, -1, -1, -1);
    GWsavevp(IC1_);
    GWreset();
    GWvport(0.0f, 0.0f, GWaspect(1), 0.5f);
    sprintf(Stmp_, "     b_k");
    GWindow(-1.0f, (float)(-ymax), (float)(nmax+1), (float)(ymax));
    GWputtxt(0.0f, (float)(ymax*0.7), Stmp_);
    plotaxis(nmax, ticksize);
    GWsetpen(RED, -1, -1, -1);
    GWsavevp(IC2_);
}

void plotaxis(int nmax, double ticksize)
{
    double k;

    GWline(0.0f, 0.0f, (float)(nmax), 0.0f);
    for(k = 1; k <= nmax; ++k)
        GWline((float)(k), (float)(-ticksize), (float)(k), (float)(ticksize));
}

void coefficents(int N, int nmax, double delta, double period, int IC1_, 
                 int IC2_)
{
    int i, k;
    double ak, bk, t, wk;

    for(k = 0; k <= nmax; ++k)
    {
        ak = 0;
        bk = 0;
        wk = 2*pi*k/period;
        // rectangular approximation
        for(i = 0; i <= N - 1; ++i)
        {
            t = i*delta;
            ak += f(t)*cos(wk*t);
            bk += f(t)*sin(wk*t);
        }
        ak *= 2.0/N;
        bk *= 2.0/N;
        GWselvp(IC1_);
        GWline((float)(k), 0.0f, (float)(k), (float)(ak));
        GWselvp(IC2_);
        GWline((float)(k), 0.0f, (float)(k), (float)(bk));
    }
}

double f(double t)
{
    double f_, w0;

    w0 = 0.1*pi;
    f_ = sin(w0*t);                         // simple example
    return f_;
}