/**************************************************************/
/*  motion in central force of f(r) = a * r^n                 */
/*        r = sqrt(x * x + y * y)                             */
/*  Original Program in N88BASIC by Katsumi Sakurai, 1992     */
/*  Ported to 4.3BSD plot version by Kenrou Adachi,           */
/*  1993,2004,2022                                            */
/**************************************************************/

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


float  x, y, vx, vy, a, n, t, h;
int    xr, yr;

void  axes(void);
void  runge(void);
float fx(float x0);
float gx(float x0, float y0, float n0, float a0);
float fy(float vy0);
float gy(float x0, float y0, float n0, float a0);


int main(void)
{
    register int  j;
    float y1;

    a = 2.0;
    n = -2.0;
    y = -18.0;
    y1 = y;
    openpl();
    erase();
    axes();
    for (j = 1; j <= 72; j++)
    {
        vx = 1.0;
        vy = 0.0;
        x  = -20.0;
        y  = y1;
        t  = 0.0;
        h  = 0.5;
        xr = (int)(9.0 * x) + 200;
        yr = -((int)(9.0 * y)) + 200;
        runge();
        y1 += 0.5;
    }
    closepl();
    return 0;
}


void axes()
{
    circle(1000, 1000, 25);
    line( 100,  100, 1900,  100);
    line( 100,  100,  100, 1900);
    line( 100, 1900, 1900, 1900);
    line(1900,  100, 1900, 1900);
} 


void runge(void)
{
    register int    i;
    int    xt, yt;
    float kx1, kx2, kx3, kx4, lx1, lx2, lx3, lx4;
    float ky1, ky2, ky3, ky4, ly1, ly2, ly3, ly4;

    for (i = 0; i <= 2000; i++)
    {
        xt = (int)(9.0 * x) + 200;
        yt = -((int)(9.0 * y)) + 200;
        line(xr * 5, yr * 5, xt * 5, yt * 5);
        xr = xt;
        yr = yt;

        kx1 = h * fx(vx);
        lx1 = h * gx(x, y, n, a);
        kx2 = h * fx(vx + lx1 / 2.0);
        lx2 = h * gx(x + kx1 /2.0, y, n, a);
        kx3 = h * fx(vx + lx2 / 2.0);
        lx3 = h * gx(x + kx2 / 2.0, y, n, a);
        kx4 = h * fx(vx + lx3);
        lx4 = h * gx(x + kx3, y, n, a);
        x  += (kx1 + 2.0 * kx2 + 2.0 * kx3 + kx4) / 6.0;
        vx += (lx1 + 2.0 * lx2 + 2.0 * lx3 + lx4) / 6.0;

        ky1 = h * fy(vy);
        ly1 = h * gy(x, y, n, a);
        ky2 = h * fy(vy + ly1 / 2.0);
        ly2 = h * gy(x, y + ky1 /2.0, n, a);
        ky3 = h * fy(vy + ly2 / 2.0);
        ly3 = h * gy(x, y + ky2 / 2.0, n, a);
        ky4 = h * fy(vy + ly3);
        ly4 = h * gy(x, y + ky3, n, a);
        y  += (ky1 + 2.0 * ky2 + 2.0 * ky3 + ky4) / 6.0;
        vy += (ly1 + 2.0 * ly2 + 2.0 * ly3 + ly4) / 6.0;

        t  += h;

        if (fabs(x) > 20.0 || fabs(y) > 20.0)
            break;
    }
}


float fx(float vx0)
{
    return vx0;
}


float gx(float x0, float y0, float n0, float a0)
{
    return a0 * x0 * pow(x0 * x0 + y0 * y0, n0 / 2.0 - 0.5);
}


float fy(float vy0)
{
    return vy0;
}


float gy(float x0, float y0, float n0, float a0)
{
    return a0 * y0 * pow(x0 * x0 + y0 * y0, n0 / 2.0 - 0.5);
}