/* etude.c     (c) Kenrou Adachi 2022 */
/* etude system for TEK4014 terminal emulator */
/* using 4.3bsd plot(1G), plot(3X), plot(3F), plot(5G) */
/* gcc -o etude -lm -lplot -l4014 */
/* ./etude > etude.out */
/* plot -Ttek etude.out */

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

#define MAXX 1024
#define MAXY 780
#define RAND_MAX 2147483647; /* 2**31-1 */

int xs, ys;

void drawing(double a, double b, double eps, double h, int l);
void axes(void);
void repeat(double *x, double *y,
            double a, double b, double eps, double h, double e);
double fx(double x, double y, double r);
double gy(double x, double y, double r);


int main(void)
{
    double a, b, eps, h;
    double x[40], y[40];
    int    m, l;
    int    n;
 
    a = 8.0;
    b = a;
    eps = 0.1;
    m = 5;
    l = 150;

    h = 2.0 * sqrt((a / 3120.0) * (a / 3120.0) +(b / 3120.0) * (b / 3120.0));
    openpl();
    axes();
    for (n = 0; n < 3; n++)
    {
        drawing(a, b, eps, h, l);
    }
    closepl();
    return 0;
}


void drawing(double a, double b, double eps, double h, int l)
{
    int    j, count;
    double rnd, e1, e2, e, x[40], y[40], xx[40], yy[40];

    j = 0;
    for (e1 = -1.0; e1 <= 0.9; e1 += 0.25)
        {
        for (e2 = -1.0; e2 <= 0.9; e2 += 0.4)
        {
            rnd = (double)(rand()) / (double)RAND_MAX;
            *(xx + j) = a * (e1 + 0.25 * rnd);
            rnd = (double)(rand()) / (double)RAND_MAX;
            *(yy + j) = b * (e2 + 0.4 * rnd);
            j++;
        }
    }
    e = 1.0;
    do
    {
        for (j = 0; j <= 39; j++)
        {
            *(x + j) = *(xx + j);
            *(y + j) = *(yy + j);
        }
        for (count = 0; count < 300; count++)
        {
            repeat(x, y, a, b, eps, h, e);
        }
        e = -e;
    } while (e < 0);
}


void repeat(double *x, double *y,
            double a, double b, double eps, double h, double e)
{
    register int j;
    double x0, y0, xr0, x1, y1, xr1, wx, wy, wx1, wy1;
    double u, v, u1, v1, q, q_h;
    int px, py, px1, py1;

    for (j = 0; j <= 39; j++)
    {
        x0 = *(x + j);
        y0 = *(y + j);

        if ((fabs(x0) <= 2.0 * a) && (fabs(y0) <= 2.0 * b))
        {
            xr0 = sqrt(x0 * x0 + y0 * y0);
            if (xr0 >= h)
            {
                u = fx(x0, y0, xr0);
                v = gy(x0, y0, xr0);
                q = sqrt(u * u + v * v);
                if (q >= eps)
                {
                    q_h = h / q;
                    x1 = x0 + 0.5 * e * u * q_h;
                    y1 = y0 + 0.5 * e * v * q_h;
                    xr1 = sqrt(x1 * x1 + y1 * y1);
                    u1 = fx(x1, y1, xr1);
                    v1 = gy(x1, y1, xr1);
                    *(x + j) = x0 + e * u1 * q_h;
                    *(y + j) = y0 + e * v1 * q_h;
                    wx = 1560.0 * x0 / a;
                    wy = 1560.0 * y0 / b;
                    wx1 = 1560.0 * (*(x + j)) / a;
                    wy1 = 1560.0 * (*(y + j)) / b;
                    px = (int)(wx + 1560.0);
                    py = (int)(wy + 1560.0);
                    px1 = (int)(wx1 + 1560.0);
                    py1 = (int)(wy1 + 1560.0);
                    if ( (px > 0) && (px1 > 0) && (px < 3120) && (px1 < 3120)
                      && (py > 0) && (py1 > 0) && (py < 3120) && (py1 < 3120))
                        line((int)(wx + 1560.0), (int)(wy + 1560.0),
                             (int)(wx1 + 1560.0), (int)(wy1 + 1560.0));
                }
            }
        }
    }
}


void axes(void)
{
  line(   0, 1560, 3120, 1560);
  line(1560, 3120, 1560, 3120);
}


double fx(double x, double y, double r)
{
    return cos(x) - x * sin(y);
}


double gy(double x, double y, double r)
{
    return x * sin(y) + y * sin(x);
}