/*************************************************************************/
/* etude_tek.c     (c) Kenrou Adachi 2022                                */
/* etude for TEK4014 terminal emulator with 4.3BSD on simh vax simulator */
/* gcc -O2 -o etude etude_tek.c -lm                                      */
/*************************************************************************/

#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);
void   startDraw(int x1, int y1);
void   draw(int x2, int y2);
void   endDraw(void);


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

    a = 10.24;
    b = 7.8;
    eps = 0.1;
    m = 5;
    l = 150;

    h = 2.0 * sqrt((a / 512.0) * (a / 512.0) +(b / 390.0) * (b / 390.0));
    axes();
    for (n = 0; n < 5; n++)
    {
        drawing(a, b, eps, h, l);
    }
    printf("The Program has been terminated!\n");
    printf("Put Enter Key!\n");
    getchar();
    startDraw(1,1);
    endDraw();
    return 0;
}


/* for Tektronics 4041 Graphics Terminal */ 

void startDraw(int x1, int y1)
{
    if (x1 < 0)      x1 = 0;
    if (x1 >= MAXX)  x1 = MAXX - 1;
    if (y1 < 0)      y1 = 0;
    if (y1 >= MAXY)  y1 = MAXY - 1;
    putchar(29);
    putchar((y1 >> 5) + 32);
    putchar((y1 & 31) + 96);
    putchar((x1 >> 5) + 32);
    putchar((x1 & 31) + 64);
    xs = x1;
    ys = y1;
}


void draw(int x2,int y2)
{
    int hxchange, lychange;

    if (x2 < 0) x2 = 0;
    if (x2 >= MAXX) x2 = MAXX - 1;
    if (y2 < 0) y2 = 0;
    if (y2 >= MAXY) y2 = MAXY - 1;

    if ((y2 >> 5) != (ys >> 5))        /* if high y has changed */
        putchar((y2 >> 5) + 32);
    hxchange = (x2 >> 5) != (xs >> 5);
    lychange = (y2 & 31) != (ys & 31);
    if (hxchange || lychange) putchar((y2 & 31) + 96);
    if (hxchange)
         putchar((x2 >> 5) + 32);
    putchar((x2 & 31) + 64);
    xs = x2;
    ys = y2;
}


void endDraw()
{
    putchar(31);
    fflush(stdout);
}


void line(int x1, int y1, int x2, int y2)
{
     startDraw(x1, y1);
     draw(x2,y2);
     endDraw();
}


/* End of Tektronics Routines */


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 < 150; 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;

    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 = 512.0 * x0 / a;
                    wy = 390.0 * y0 / b;
                    wx1 = 512.0 * (*(x + j)) / a;
                    wy1 = 390.0 * (*(y + j)) / b;
                    line((int)(wx + 512.0), (int)(wy + 390.0),
                         (int)(wx1 + 512.0), (int)(wy1 + 390.0));
                }
            }
        }
    }
}


void axes(void)
{
  line(0, 390, 1023, 390);
  line(512, 0, 512,  779);
}


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);
}