/*
 *
 * This is a program to study microstructural evolution using
 * Cahn-Hilliard equations (i.e., for conserved order parameters).
 *
 * Author: M. P. Gururajan
 *
 * Copyright (c) 2003 Computational Materials Science Laboratory,
 * Department of Metallurgy, Indian Institute of Science, Bangalore
 * 560 012. INDIA.
 *
 */

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

#include <gsl/gsl_rng.h>
#include <fftw3.h>

/*
 * local function declarations
 */

#include "functions.h"

static void microstructure_init(int n_x, int n_y, fftw_complex *comp, 
                                gsl_rng *r, double c_0, double noise_str);
static int save_settings(const char *readme, int n_x, int n_y, 
                  double delta_x, double delta_y, double delta_t, 
                  double c_zero, double noise_str, int time_steps);

#define EXIT_SUCCESS 0
#define EXIT_FAILURE 1

int main(void)
{

    int n_x, n_y;               /* The number of nodes in x and y directions */
    double delta_x, delta_y;    /* The grid spacing in x and y directions */
    int time_steps;             /* The total number of time steps */
    double delta_t;             /* The time step */
    fftw_complex *comp;         /* The composition field */
    double c_zero;              /* The nominal alloy composition */
    double noise_str;           /* The strength of initial noise */
    gsl_rng *r;                 /* The random number generator */
    unsigned long seed = 2UL;   /* rng seed */
    const char readme[] = "output/README";

    /*
     * Grid size and cell dimensions
     */

    n_x = 512;
    n_y = 512;
    delta_x = 1.0;
    delta_y = 1.0;

    /*
     * Time step size and the number of steps
     */

    delta_t = 1.0;
    time_steps = 4000;

    /*
     * Average composition and the strength of initial noise
     */

    c_zero = 0.5;
    noise_str = 0.05e-2;

    /*
     * The values used for a particular run are written into the file
     * "README" in the directory named "output/".
     */

    if( save_settings(readme, n_x, n_y, delta_x, delta_y, delta_t, c_zero, noise_str, time_steps) == EXIT_FAILURE)
    {
        fprintf(stderr, "Unable to open %s. Exiting", readme);
        exit (EXIT_FAILURE);
    }

    /*
     * Memory allocation for the microstructrure
     */

    comp = (fftw_complex *) fftw_malloc((size_t) (n_x * n_y) * sizeof(fftw_complex));

    /*
     * Set up the random number generator.
     */

    r = gsl_rng_alloc (gsl_rng_taus);
    gsl_rng_set (r, seed);

    /*
     * Generate the initial random microstructure
     */

    microstructure_init (n_x, n_y, comp, r, c_zero, noise_str);

    /*
     * evolve the microstructure
     */

    evolve(n_x, n_y, delta_x, delta_y, delta_t, time_steps, comp);

    /*
     * Release the allocated memory
     */

    fftw_free(comp);
    gsl_rng_free(r);

    return(EXIT_SUCCESS);
}


/*
 * Generate the initial random microstructure
 */
static void microstructure_init (int n_x, int n_y, fftw_complex *comp, gsl_rng *r, double c_0, double noise_str)
{

    int i, j;
    double av = 0.0;
    fftw_complex del;

    for (i = 0; i < n_x; i++)
    {
        for (j = 0; j < n_y; j++)
        {
            /*
             * creal(comp[j+n_y*i]) = c_0 + noise_str*(0.5-gsl_rng_uniform(r));
             * cimag(comp[j+n_y*i]) = 0.0;
             */
            comp[j + n_y * i] = (fftw_complex) (c_0 + noise_str * (0.5 - gsl_rng_uniform(r)));
            av += creal(comp[j + n_y * i]);
        }
    }
    av /=  (n_x * n_y);
    del = (fftw_complex) (c_0 - av);

    /*
     * correct the deviation of av from c_0
     */
    for (i = 0; i < n_x; i++)
    {
        for (j = 0; j < n_y; j++)
        {
            comp[j + n_y * i] += del;
        }
    }
}

/*
 * The values used for a particular run are written into a file
 */
static int save_settings(const char *readme, int n_x, int n_y, 
                  double delta_x, double delta_y, double delta_t, 
                  double c_zero, double noise_str, int time_steps)
{
    FILE *fpw;

    if ((fpw = fopen(readme, "w")) == NULL)
    {
        return EXIT_FAILURE;
    }

    fprintf(fpw, "n_x = %d\nn_y = %d\n", n_x, n_y);
    fprintf(fpw, "delta_x = %le\ndelta_y = %le\n", delta_x, delta_y);
    fprintf(fpw, "delta_t = %le\n", delta_t);
    fprintf(fpw, "time_steps = %d\n", time_steps);
    fprintf(fpw, "c_zero = %le\n", c_zero);
    fprintf(fpw, "noise_str = %le\n", noise_str);

    (void) fclose(fpw);
    return EXIT_SUCCESS;
}