/*
 * This function evolves a given microstructure using the Cahn-Hilliard
 * equation. It uses FFTW for carrying out the Fourier transforms.
 *
 * 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_math.h>
#include <fftw3.h>

#include "functions.h"

#define NAMELEN 50
#define LABELLEN 7

static void local_contrib(int n_x, int n_y, fftw_complex *comp, fftw_complex *g)
{
    int i, j;
    double compos;
    for (i = 0; i < n_x; i++)
    {
        for (j = 0; j < n_y; j++)
        {
            compos = creal(comp[j + n_y * i]);
            g[j + n_y * i] = 2 * compos * (1. - compos) * (1. - 2 * compos);
        }
    }
}

void evolve(int n_x, int n_y, double delta_x, double delta_y,
            double delta_t, int time_steps, fftw_complex *comp)
{

    int step;
    int stone;
    int i, j;
    int half_nx, half_ny;
    double kx, ky, delta_kx, delta_ky;
    double kx2, ky2, k2, k4;
    double denom;
    double *c;
    double norm;

    const char nameformat[] = "output/time%06d%s";
    const char labelformat[] = "%06d";
    const char ext[] = ".gif";

    char label[LABELLEN];
    char name[NAMELEN];

    fftw_complex *g;
    fftw_plan plan1, plan2, plan3;

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

    plan1 = fftw_plan_dft_2d (n_x, n_y, comp, comp, FFTW_FORWARD, FFTW_ESTIMATE);
    plan2 = fftw_plan_dft_2d (n_x, n_y, g, g, FFTW_FORWARD, FFTW_ESTIMATE);
    plan3 = fftw_plan_dft_2d (n_x, n_y, comp, comp, FFTW_BACKWARD, FFTW_ESTIMATE);


    c = (double *) malloc((size_t) (n_x * n_y) * sizeof(double));

    half_nx = (int) n_x / 2;
    half_ny = (int) n_y / 2;

    delta_kx = (2. * M_PI) / (n_x * delta_x);
    delta_ky = (2. * M_PI) / (n_y * delta_y);

    norm = 1./(n_x*n_y);

    /*
     * Make a snapshot of the initial microstructure
     */

    for (i = 0; i < n_x; i++)
    {
        for (j = 0; j < n_y; j++)
        {
            c[j + n_y * i] = creal(comp[j + n_y * i]);
        }
    }

    step = 0;
    snprintf (name, NAMELEN, nameformat, step, ext);
    snprintf (label, LABELLEN, labelformat, step);

    generate_snapshot(name, label, n_x, n_y, c);

    /*
     * Evolve the microstructure in time for the given number of time steps.
     * n_x/2 and n_y/2 are needed for imposing the periodic boundary conditions.
     */

    for (step = 1; step <= time_steps; step++)
    {

        /* calculate g */

        local_contrib(n_x, n_y, comp, g);
        /*
        for (i = 0; i < n_x; i++) {
            for (j = 0; j < n_y; j++) {
                compos = creal(comp[j + n_y * i]);
                g[j + n_y * i] = 2 * compos * (1. - compos) * (1. - 2 * compos);
            }
        }
        */

        /* Fourie transforms of g and comp */

        fftw_execute(plan1);
        fftw_execute(plan2);

        /* evolve the composition profile */

        for (i = 0; i < n_x; i++)
        {
            kx = (i < half_nx) ? i * delta_kx : (i - n_x) * delta_kx;
            kx2 = kx * kx;

            for (j = 0; j < n_y; j++)
            {
                ky = (j < half_ny) ? j * delta_ky : (j - n_y) * delta_ky;
                ky2 = ky * ky;

                k2 = kx2 + ky2;
                k4 = k2 * k2;

                denom = 1. + 2 * k4 * delta_t;
                comp[j + n_y * i] = (comp[j + n_y * i] - k2 * delta_t * g[j + n_y * i]) / denom;

            }
        }

        /* Transform comp back to real space */

        fftw_execute(plan3);

        for (i = 0; i < n_x; i++)
        {
            for (j = 0; j < n_y; j++)
            {
                comp[j + n_y * i] *= norm;
            }
        }

        /*
         * Store the microstructure
         */

        if (step < 150)
        {
            stone = 1;
        }
        else if (step < 300)
        {
            stone = 2;
        }
        else if (step < 600)
        {
            stone = 3;
        }
        else if (step < 1000)
        {
            stone = 5;
        }
        else if (step < 2000)
        {
            stone = 10;
        }
        else
        {
            stone = step/100;
        }

        if (step % stone == 0)
        {

            for (i = 0; i < n_x; i++)
            {
                for (j = 0; j < n_y; j++)
                {
                    c[j + n_y * i] = creal(comp[j + n_y * i]);
                }
            }

            snprintf(name, NAMELEN, nameformat, step, ext);
            snprintf(label, LABELLEN, labelformat, step);

            generate_snapshot(name, label, n_x, n_y, c);

            printf("finished step %7d out of %7d\r", step, time_steps);
            fflush(stdout);
        }

    }
    printf("Done                         \n");

    free(c);

    fftw_free(g);

    fftw_destroy_plan(plan1);
    fftw_destroy_plan(plan2);
    fftw_destroy_plan(plan3);
}