/*==========================================================================*/
/*                                                                          */
/*       Diffusion Limited Aggregation on the simple square lattice         */
/*                                                                          */
/*       PROGRAM NAME:      dla                                             */
/*                                                                          */
/*       CALLING SEQUENCE:                                                  */
/*                                                                          */
/*         Main program without command line arguments.                     */
/*                                                                          */
/*       RETURN:                                                            */
/*                                                                          */
/*         Exit code      -1  The requested lattice does not fit into the   */
/*                            memory. Resize your request                   */
/*                        -2  The set-up circles are larger than the        */
/*                            lattice.                                      */
/*                                                                          */
/*   Version 1.00                                                           */
/*                                                                          */
/*   (c) Satoru Ozawa and Dieter W. Heermann                                */
/*   Hitachi, Heidelberg July 1994                                          */
/*                                                                          */
/*==========================================================================*/

/* ==== Include files for this program ==== */

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

# include "sgl.h"               /* include for the graphics library */
# include "r250.c"               /* random number generator          */

/* ==== Constants ==== */

# define NMax                   /* largest saw which can be generated */
# define WindowSize  500        /* Size in pixle of the graphics      */
                                /* window                             */
# define dx          0.005       /* height and width, in terms of the  */
                                /* normalize window, to draw          */
# define HIT         1          /* Display the fractal every time the */
                                /* cluster has grown by HIT number of */
                                /* particles                          */

/*--------------------------------------------------------------------------*/
/*                         Main program                                     */
/*--------------------------------------------------------------------------*/
#ifdef  __STDC__
void    main(int argc, char *argv[])
#else
void    main(argc,argv)
        int     argc;
        char    *argv[];
#endif
 {
    /*---------------------------------------------------------- */
    /*                     Declarations                          */
    /*---------------------------------------------------------- */

    /* ====  Area Declarations ==== */

    int      mf[251];
    float    *random;
    int      *grid;

    /* ==== General declarations ==== */

    int      i,j;
    int      index;
    int      L;
    int      sweep, sample_size;
    int      seed;
    int      rad_inc;
    int      return_code;
    int      midx,midy;
    int      x,y;
    int      xn,yn;
    int      xq,yq;
    int      int_radius = 0;
    int      new_radius = 1;
    int      rem_radius;
    int      check_radius;
    int      nearx,neary;
    int      inside = 1;
    int      hit_count = 0;
    int      size;
    int      mass = 1;
    int      icmx, icmy;
    int      xsq, ysq;
    int      irad;

    float    theta;
    float    cmx, cmy;
    float    gyration_rad;
    float    box_radius = 0.0;
    float    PI=3.141592654;

    /*---------------------------------------------------------- */
    /*                End of declares                            */
    /*---------------------------------------------------------- */
 
    /*---------------------------------------------------------- */
    /*      Start with setting up procedures                     */
    /*---------------------------------------------------------- */

    /* ====  Print the header for this test  ==== */

    printf("================================================\n");
    printf("======>  Diffusion Limited Aggregation  <=======\n");
    printf("================================================\n");
    printf("\n\n");
    printf("(c) Satoru Ozawa and Dieter W. Heermann \n");
    printf("Version 1.00 (1994) \n");
    printf("\n\n");

    printf("Linear dimension of the lattice  = ");
    scanf("%d",&L);
    printf("Sample size                      = ");
    scanf("%d",&sample_size);
    printf("Increment for the set-up circle  = ");
    scanf("%d",&rad_inc);
    printf("Seed for Random Number Generator = ");
    scanf("%d",&seed);

    /* ==== Prepare by initializing all necessary quantities ==== */

    size      = (L+2)*(L+2);
    grid = (int *) calloc(size, sizeof(int));
    if ( grid == NULL ) {
        printf("Not enough memory to hold the lattice");
        exit(-1);
    }
 
    /* ==== Grep the storage for the random numbers.        ==== */
    random = (float *) calloc(L, sizeof(float));
    if ( random == NULL ) {
        printf("Not enough memory to hold the lattice");
        free(grid);
        exit(-1);
    }
 
    /* ==== set up the random number generator              ==== */
    seed       = (seed << 2) + 1;
    return_code = init_r250(seed,mf);

    /* ==== Prepare the lattice                             ==== */
    for(i=0; i < L+2; i++) {
       for(j=0; j < L+2; j++) {
        grid[i*(L+2) + j] = 0;
       }
    }

    /* ==== Set the seed for the walker                     ==== */
    midx = (L+1)/2;
    midy = (L+1)/2;
    grid[midx * (L+2) + midy+1] = 1;
    cmx = midx;
    cmy = midy;
    xsq = midx * midx;
    ysq = midy * midy;

    /* ==== initialize the graphics ==== */

    sgl_winsize(WindowSize,WindowSize);
    sgl_init();

    /*---------------------------------------------------------- */
    /*                M o n t e  C a r l o                       */
    /*---------------------------------------------------------- */

    for(sweep = 0; sweep < sample_size; sweep++) {

       /* ==== Display the current Walker number ==== */
       if ( sweep % 100 == 0 ) {
           printf("Walker Nummer %d ",sweep);
       }

       if ( new_radius) {
          new_radius   = 0;
          check_radius = int_radius * int_radius;
          int_radius  += rad_inc;
          rem_radius   = (int_radius + rad_inc) * (int_radius + rad_inc);

          if (int_radius >= (L/2) ) {
             printf("Lattice boundary is reached");
             printf("To close the graphics window type in any key");
             i = getchar();
             i = getchar();
             free(random);
             free(grid);
             exit(-2);
          }

       }  /* --- new radius --- */

       return_code = r250(2,random,mf);
       theta = random[1] * 2.0 * PI;
       x     = int_radius * cos(theta) + midx ;
       y     = int_radius * sin(theta) + midy ;

       return_code = r250(L,random,mf);
       index  = 0;
       inside = 1;
       while( inside ) {

          /* ---- Choose a new nn site ---- */
          i = floor(random[index++]* 4.0);
          switch (i) {
             case 0: xn = x-1;
                     yn = y;
                     break;
             case 1: yn = y-1;
                     xn = x;
                     break;
             case 2: yn = y+1;
                     xn = x;
                     break;
             case 3: xn = x+1;
                     yn = y;
                     break;
          } /* ---- switch i ---- */

          if ( index == L ) {
              return_code = r250(L,random,mf);
              index = 0;
          }  /* ---- index == L ---- */

          /* ---- Check whether out of bounds ---- */
          if ( (xn < 1) || (xn > L) ||
               (yn < 1) || (yn > L) ||
               ((xn-midx)*(xn-midx) + (yn-midy)*(yn-midy) > rem_radius) ) {
            inside = 0;
            continue;
          } /* ---- bounds ---- */

          /* ---- Check whether a nearest neighbour is occupied ---- */
          nearx = xn+1;
          neary = yn;
          if ( grid[nearx*(L+2) + neary +1] == 1 ) {
            grid[xn*(L+2)+yn+1] = 1;
            xq = (xn - midx) * (xn - midx);
            yq = (yn - midy) * (yn - midy);
            if ((xq + yq) >= check_radius) new_radius = 1;
            if ((xq + yq) >  box_radius  ) box_radius = xq + yq;
            mass++;
            cmx += xn;
            cmy += yn;
            xsq += xn * xn;
            ysq += yn * yn;

            inside = 0;
            hit_count++;
            continue;
          }
          nearx = xn;
          neary = yn+1;
          if ( grid[nearx*(L+2) + neary +1] == 1 ) {
            grid[xn*(L+2)+yn+1] = 1;
            xq = (xn - midx) * (xn - midx);
            yq = (yn - midy) * (yn - midy);
            if ((xq + yq) >= check_radius) new_radius = 1;
            if ((xq + yq) >  box_radius  ) box_radius = xq + yq;
            mass++;
            cmx += xn;
            cmy += yn;
            xsq += xn * xn;
            ysq += yn * yn;
            inside = 0;
            hit_count++;
            continue;
          }
          nearx = xn-1;
          neary = yn;
          if ( grid[nearx*(L+2) + neary+1] == 1 ) {
            grid[xn*(L+2)+yn+1] = 1;
            xq = (xn - midx) * (xn - midx);
            yq = (yn - midy) * (yn - midy);
            if ((xq + yq) >= check_radius) new_radius = 1;
            if ((xq + yq) >  box_radius  ) box_radius = xq + yq;
            mass++;
            cmx += xn;
            cmy += yn;
            xsq += xn * xn;
            ysq += yn * yn;
            inside = 0;
            hit_count++;
            continue;
          }
          nearx = xn;
          neary = yn-1;
          if ( grid[nearx*(L+2) + neary+1] == 1 ) {
            grid[xn*(L+2)+yn+1] = 1;
            xq = (xn - midx) * (xn - midx);
            yq = (yn - midy) * (yn - midy);
            xsq += xn * xn;
            ysq += yn * yn;
            if ((xq + yq) >= check_radius) new_radius = 1;
            if ((xq + yq) >  box_radius  ) box_radius = xq + yq;
            mass++;
            cmx += xn;
            cmy += yn;
            inside = 0;
            hit_count++;
            continue;
          }

          x = xn;
          y = yn;

       } /* ---- While ---- */

       if ( hit_count == HIT) {

        /* ---- Graphics part ---- */

        sgl_color_background(SGL_BLACK);
        sgl_color(SGL_WHITE);
        for(i=0;i < L;i++) {
           for(j=0;j < L;j++) {
               if ( grid[i*(L+2)+j+1] == 1) {
                  sgl_rectarea((float)i/(float) L,(float)j / (float) L,dx,dx);
               }
           }
         }
         sgl_display();

           hit_count = 0;
       }
    } /* end walker */

   free(random);
   free(grid);

 exit(0);
} /* ---- end of main ---- */