/* Compile using `gcc -lm thing.c'. In addition to displaying the coupled map
lattice being studied on the screen, this program also outputs the second
projection to a graphics file `thing.bmp'. If you ensure that the bare bones
bitmap file of that name is present in the working directory on a given
iteration, and you remove it before iterating again, you will have a
viewable RGB bitmap. */

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

main()
  {
    FILE *thing;
    double rnd = atan(1);
    float PI, angle, d, e, r, z, xold[151][151], yold[151][151], xnew[151][151], ynew[151][151], ang[151][151], temp1, temp2, temp3, temp4;
    char ch;
    int c, gap, m, n, mmax, nmax, time = 0;
    PI = 4 * atan(1);

    printf("\nThis program generates three different symbolic representations of a \n");
    printf("coupled map lattice. The first shows the locations of clockwise and \n");
    printf("anticlockwise spiral centres, labelled a and b respectively. The \n");
    printf("second shows in what quadrant of local state space (the unit disc) \n");
    printf("each site variable lies, by using the symbols ., \\, * and /. The \n");
    printf("third and final projection shows the modulus of each site variable, \n");
    printf("with the symbols ., 1, 2, 3, 4 and * representing complex numbers of \n");
    printf("increasing magnitude (approaching the unit circle). \n");
    printf("\nEpsilon (coupling parameter) (0 <= mx <= 1): ");
    scanf("%f", &e);
    printf("\nDelta (local map parameter): ");
    scanf("%f", &d);
    printf("\nBest to maximise this window and enter lattice dimensions that will \n");
    printf("fit onto your screen three times horizontally and once vertically. \n");
    mmax = 100;
    nmax = 100;

    printf("\nAfter how many iterations should I start displaying the lattices and \n");
    printf("waiting for you to press enter at each Nth iteration (much slower than \n");
    printf("doing it invisibly and automatically)? ");
    scanf("%d", &c);
    printf("\nEnter step size N: ");
    scanf("%d", &gap);
    printf("\nEnter a random number between 0 and 1 to generate initial conditions: ");
    scanf("%f", &rnd);
    ch = getchar();
    printf("\nPlease wait... \n \n");

    for (n = 0; n < nmax + 2; n++)
      {
        for (m = 0; m < mmax + 2; m++)
          {
            rnd = rnd - floor(rnd);
            xold[m][n] = (rnd - 0.5) / 5;
            rnd = rnd * 12345;
            rnd = rnd - floor(rnd);
            yold[m][n] = (rnd - 0.5) / 5;
            rnd = rnd * 12345;
            rnd = rnd - floor(rnd);
            xnew[m][n] = (rnd - 0.5) / 5;
            rnd = rnd * 12345;
            rnd = rnd - floor(rnd);
            ynew[m][n] = (rnd - 0.5) / 5;
            rnd = rnd * 12345;
          }
      };

    do
      {
        for (n = 1; n < nmax + 1; n++)
          {
            for (m = 1; m < mmax + 1; m++)
              {
                r = xnew[m][n] * xnew[m][n] + ynew[m][n] * ynew[m][n];
                angle = d * r * (1 - r);
                xold[m][n] = (xnew[m][n] * cos(angle) - ynew[m][n] * sin(angle)) * (11 - r) / 10;
                yold[m][n] = (xnew[m][n] * sin(angle) + ynew[m][n] * cos(angle)) * (11 - r) / 10;
              }
          }

        for (n = 1; n < nmax + 1; n++)
          {
            for (m = 1; m < mmax + 1; m++)
              {
                xold[0][n] = xold[mmax][n];
                xold[mmax + 1][n] = xold[1][n];
                xold[m][0] = xold[m][nmax];
                xold[m][nmax + 1] = xold[m][1];
                yold[0][n] = yold[mmax][n];
                yold[mmax + 1][n] = yold[1][n];
                yold[m][0] = yold[m][nmax];
                yold[m][nmax + 1] = yold[m][1];
              }
          }

        xold[0][0] = xold[mmax][nmax];
        xold[mmax + 1][nmax + 1] = xold[1][1];
        xold[0][nmax + 1] = xold[mmax][1];
        xold[mmax + 1][0] = xold[1][nmax];
        yold[0][0] = yold[mmax][nmax];
        yold[mmax + 1][nmax + 1] = yold[1][1];
        yold[0][nmax + 1] = yold[mmax][1];
        yold[mmax + 1][0] = yold[1][nmax];

        for (n = 0; n < nmax + 2; n++)
          {
            for (m = 0; m < mmax + 2; m++)
              {
                xnew[m][n] = (1 - e) * xold[m][n] + e * (xold[m+1][n+1]/8 + xold[m+1][n]/8 + xold[m][n+1]/8 + xold[m-1][n-1]/8 + xold[m-1][n]/8 + xold[m][n-1]/8 + xold[m+1][n-1]/8 + xold[m-1][n+1]/8);
                ynew[m][n] = (1 - e) * yold[m][n] + e * (yold[m+1][n+1]/8 + yold[m+1][n]/8 + yold[m][n+1]/8 + yold[m-1][n-1]/8 + yold[m-1][n]/8 + yold[m][n-1]/8 + yold[m+1][n-1]/8 + yold[m-1][n+1]/8);
              }
          }

        /* angle = 0;
        if (ynew[1][1] > 0)
          angle = PI/2;
        if (ynew[1][1] < 0)
          angle = 3*PI/2;
        if (xnew[1][1] > 0)
          angle = atan(ynew[1][1]/xnew[1][1]);
        if (xnew[1][1] < 0)
          angle = atan(ynew[1][1]/xnew[1][1]) + PI; */

        for (n = 0; n < nmax + 2; n++)
          {
            for (m = 0; m < mmax + 2; m++)
              {
                xold[m][n] = xnew[m][n];
                yold[m][n] = ynew[m][n];
                /* xnew[m][n] = xold[m][n] * cos(angle) + yold[m][n] * sin(angle);
                ynew[m][n] = -xold[m][n] * sin(angle) + yold[m][n] * cos(angle); */
                ang[m][n] = 0;
                if (ynew[m][n] == 0)
                  ang[m][n] = PI;
                if (ynew[m][n] > 0)
                  ang[m][n] = PI/2;
                if (ynew[m][n] < 0)
                  ang[m][n] = -PI/2;
                if (xnew[m][n] > 0)
                  ang[m][n] = atan(ynew[m][n]/xnew[m][n]);
                if (xnew[m][n] < 0 && ynew[m][n] > 0)
                  ang[m][n] = atan(ynew[m][n]/xnew[m][n]) + PI;
                if (xnew[m][n] < 0 && ynew[m][n] < 0)
                  ang[m][n] = atan(ynew[m][n]/xnew[m][n]) - PI;
	      }
	  }

        if (time > c && time == (time / gap) * gap)
          {
            for (n = 1; n < nmax + 1; n++)
              {
                for (m = 1; m < mmax + 1; m++)
                  {
                    temp1 = ang[m+1][n] - ang[m][n];
                    if (temp1 > PI)
                      temp1 = ang[m+1][n] - ang[m][n] - 2*PI;
                    if (temp1 < -PI)
                      temp1 = ang[m+1][n] - ang[m][n] + 2*PI;
                    temp2 = ang[m+1][n+1] - ang[m+1][n];
                    if (temp2 > PI)
                      temp2 = ang[m+1][n+1] - ang[m+1][n] - 2*PI;
                    if (temp2 < -PI)
                      temp2 = ang[m+1][n+1] - ang[m+1][n] + 2*PI;
                    temp3 = ang[m][n+1] - ang[m+1][n+1];
                    if (temp3 > PI)
                      temp3 = ang[m][n+1] - ang[m+1][n+1] - 2*PI;
                    if (temp3 < -PI)
                      temp3 = ang[m][n+1] - ang[m+1][n+1] + 2*PI;
                    temp4 = ang[m][n] - ang[m][n+1];
                    if (temp4 > PI)
                      temp4 = ang[m][n] - ang[m][n+1] - 2*PI;
                    if (temp4 < -PI)
                      temp4 = ang[m][n] - ang[m][n+1] + 2*PI;
                    yold[m][n] = temp1 + temp2 + temp3 + temp4;
                    if (yold[m][n]/PI < 0.0001 && yold[m][n]/PI > -0.0001)
                      printf(".");
                    if (yold[m][n]/PI < 2.0001 && yold[m][n]/PI > 1.9999)
                      printf("a");
                    if (yold[m][n]/PI < -1.9999 && yold[m][n]/PI > -2.0001)
                      printf("b");
                  }
                printf("|");
                thing = fopen("thing.bmp", "a");
                for (m = 1; m < mmax + 1; m++)
                  {
                    if (ang[m][n] == 0)
                      {
                        printf("/");
                        fprintf(thing, "1");
		      }
                    if (ang[m][n] > -PI && ang[m][n] < -PI/2)
                      {
                        printf(".");
                        fprintf(thing, "n");
		      }
                    if (ang[m][n] > -PI/2 && ang[m][n] < 0)
                      {
                        printf("\\");
                        fprintf(thing, "g");
		      }
                    if (ang[m][n] > 0 && ang[m][n] < PI/2)
                      {
                        printf("*");
                        fprintf(thing, "Z");
		      }
                    if (ang[m][n] > PI/2 && ang[m][n] < PI)
                      {
                        printf("/");
                        fprintf(thing, "1");
		      }
                  }
                fclose(thing);
                printf("|");
                for (m = 1; m < mmax + 1; m++)
                  {
                    r = sqrt(xnew[m][n] * xnew[m][n] + ynew[m][n] * ynew[m][n]);
                    if (r < 0.9)
                      printf(".");
                    if (r > 0.9 && r < 0.92)
                      printf("1");
                    if (r > 0.92 && r < 0.94)
                      printf("2");
                    if (r > 0.94 && r < 0.96)
                      printf("3");
                    if (r > 0.96 && r < 0.98)
                      printf("4");
                    if (r > 0.98)
                      printf("*");
                  }
                printf("\n");
              }

            printf("\nz (1, 1) = %f + %fi \n", xnew[1][1], ynew[1][1]);
            printf("\nz (2, 1) = %f + %fi \n", xnew[2][1], ynew[2][1]);
            printf("\nz (1, 2) = %f + %fi \n", xnew[1][2], ynew[1][2]);
            printf("\nz (2, 2) = %f + %fi \n", xnew[2][2], ynew[2][2]);
            printf("\nTime = %d. Press enter to iterate. \n \n", time);
            ch = getchar();
	  }

        time++;

        for (n = 1; n < nmax + 1; n++)
          {
            for (m = 1; m < mmax + 1; m++)
              {
                xold[m][n] = xnew[m][n];
                yold[m][n] = ynew[m][n];
              }
          }
      }
    while (1 != 2);
  }