/************************************************************
This function evolves a given microstructure using the Cahn-Hilliard
equation. It uses FFTW (version 3.0) 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.

See the README file for further details.

************************************************************/

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

#ifndef S_SPLINT_S
#include <complex.h>
#endif

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

#include "functions.h"

#define NAMELEN 50
#define LABELLEN 7

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

    int step = 0;
    int stone;
    int i1, i2;
    int half_nx, half_ny;
    double kx, ky, delta_kx, delta_ky;
    double kx2, ky2, k2, k4;
    double denom;
    double *c;
    double compos;
    double norm;

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

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

    fftw_complex *g;
    fftw_plan plan1, plan2, plan3;

    norm = 1./(n_x*n_y);

    g = fftw_malloc(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));

    /* 
     * Store initial microstructure
     */

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

    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.
     * delta_kx and delta_ky include the factor of 2PI.
     */

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

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

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

        /* calculate g */

        for (i1 = 0; i1 < n_x; ++i1) {
            for (i2 = 0; i2 < n_y; ++i2) {
                compos = creal(comp[i2 + n_y * i1]);
                g[i2 + n_y * i1] = 2.0 * A * compos * (1.0 - compos) * (1.0 - 2.0 * compos);
            }
        }

        /* Fourie transforms of g and comp */

        fftw_execute(plan1);
        fftw_execute(plan2);

        /* evolve the composition profile */

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

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

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

                denom = 1.0 + 2.0 * kappa * k4 * delta_t;
                comp[i2 + n_y * i1] = (comp[i2 + n_y * i1] - k2 * delta_t * g[i2 + n_y * i1]) / denom;

            }
        }

        /* Transform comp back to real space */

        fftw_execute(plan3);

        for (i1 = 0; i1 < n_x; ++i1) {
            for (i2 = 0; i2 < n_y; ++i2) {
                comp[i2 + n_y * i1] *= 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 (i1 = 0; i1 < n_x; ++i1) {
                for (i2 = 0; i2 < n_y; ++i2) {
                    c[i2 + n_y * i1] = creal(comp[i2 + n_y * i1]);
                }
            }

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

            generate_snapshot(name, label, n_x, n_y, c);
            printf("finished step %7d\r", step);
            fflush(stdout);
        }

    }
    printf("Done                         \n");

    free(c);

    fftw_free(g);

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