/* This is a solver for the equation 
 * x''(t) + sin(x(t)) = amp sin(omega t)
 * or alternatively
 * x'(t) = y
 * y'(t) = -sin(x(t)) + amp sin(omega t)
 * with ICs x(0) = x'(0) = 0, and any value for omega.
 *
 * Compilation:
 *      gcc -o diffsolve diffsolve.c
 *
 * Usage:
 *      ./diffsolve amp omega tmax dt x0 y0
 * 
 * (We say that we 'pass' the arguments amp, omega, etc. to the executable
 * diffsolve. On a Unix-like system, you type the above into the directory where
 * diffsolve has been compiled. I don't know how you do it on Windows or in an
 * IDE, sorry.)
 */

#include<stdio.h>
#include<math.h>
#include<stdlib.h>

typedef struct {
    double amp;
    double omega;
} problem_t;

double forcing(double t, problem_t params) 
{
    return params.amp * sin( params.omega * t );
}

double rhs_x(double t, double x, double y, problem_t params)
{
    return y;
}

double rhs_y(double t, double x, double y, problem_t params)
{
    return -sin(x) + forcing(t, params);
}

int main(int argc, char* argv[]) 
{
    double x;
    double y;
    double tmax;
    double dt;
    problem_t params;

    if (argc > 6) 
    {
        params.amp   = atof(argv[1]);
        params.omega = atof(argv[2]);
        tmax  = atof(argv[3]);
        dt    = atof(argv[4]);
        x     = atof(argv[5]);
        y     = atof(argv[6]);
    }
    else
    {
        fprintf(stderr, "Usage: %s amp omega tmax dt x0 y0\n", argv[0]);
        exit(-1);
    }

    double t = 0;
    while (t < tmax) 
    {
        // If the time from the current time to tmax is greater than dt, then
        // advance things by dt. Otherwise, advance just by tmax-t.
        dt = fmin(dt, tmax-t); 

        // Print the current result, to six decimal places of accuracy.
        printf("%.6f %.6f %.6f\n", t, x, y); 
        
        // Advance the solution using the Runge-Kutta method. Note that 'the'
        // RK4 method (as given on Wikipedia) is designed for a single
        // first-order equation; here we have a pair of couple equations so we
        // need to do a little bit more work.

        double k1, k2, k3, k4;
        double l1, l2, l3, l4;
        k1 = rhs_x(t,          x,             y,             params);
        l1 = rhs_y(t,          x,             y,             params);
        k2 = rhs_x(t + dt*0.5, x + k1*dt*0.5, y + l1*dt*0.5, params);
        l2 = rhs_y(t + dt*0.5, x + k1*dt*0.5, y + l1*dt*0.5, params);
        k3 = rhs_x(t + dt*0.5, x + k2*dt*0.5, y + l2*dt*0.5, params);
        l3 = rhs_y(t + dt*0.5, x + k2*dt*0.5, y + l2*dt*0.5, params);
        k4 = rhs_x(t + dt,     x + k3*dt,     y + l3*dt,     params);
        l4 = rhs_y(t + dt,     x + k3*dt,     y + l3*dt,     params);

        x += (dt / (double)6) * (k1 + 2*k2 + 2*k3 + k4);
        y += (dt / (double)6) * (l1 + 2*l2 + 2*l3 + l4);

        // Advance the time. 
        t += dt;
    }

    // Print the final state of the system. 
    printf("%.6f %.6f %.6f\n", t, x, y); 

    return 0;
}