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

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.

See the README file for further details.

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

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

/* 
 * The following ifndef is for the sake of splint 
 */

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

/* 
 * The Gnu Scientific Library (GSL) header file (for generating random
 * numbers) and FFTW header file (for doing the Fourier transforms)
 */

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

/*
 * local function declarations 
 */

#include "functions.h"

int main(void)
{

    FILE *fpr, *fpw;            /* File pointers to read and write data */
    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 */
    double kappa;               /* The gradient energy coefficient */
    double A;                   /* The chemical free energy constant */
    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 ave_comp;            /* The actual alloy composition  */
    double noise_str;           /* The strength of initial noise */
    gsl_rng *ran_num;           /* The random number */
    int i1, i2;                 /* Indices */

    /*
     * The input values used for any particular run is written on to the file
     * "README" in the directory named "output/".  I open it here.
     */

    if ((fpw = fopen("output/README", "w")) == NULL) {
        printf("Unable to open output/README. Exiting");
        exit(0);
    } 

    /*
     * Let us read the system size and grid size
     */

    if ((fpr = fopen("input/system_size", "r")) == NULL) {
        printf("Unable to open input/system_size. Exiting");
        exit(0);
    } 

    (void) fscanf(fpr, "%d%d%le%le", &n_x, &n_y, &delta_x, &delta_y);
    (void) fclose(fpr);

    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);

    /*
     * Let us read the chemical free energy related constants
     */

    if ((fpr = fopen("input/constants", "r")) == NULL) {
        printf("Unable to oepn input/constants. Exiting");
        exit(0);
    } 

    (void) fscanf(fpr, "%le%le", &kappa, &A);
    (void) fclose(fpr);

    fprintf(fpw, "kappa = %le\nA = %le\n", kappa, A);

    /*
     * Let us read the time step and number of time steps
     */

    if ((fpr = fopen("input/time_information", "r")) == NULL) {
        printf("Unable to open input/time_information. Exiting");
        exit(0);
    }

    (void) fscanf(fpr, "%le%d", &delta_t, &time_steps);
    (void) fclose(fpr);

    fprintf(fpw, "delta_t = %le\n", delta_t);
    fprintf(fpw, "time_steps = %d\n", time_steps);

    /*
     * Let us read the average composition and the strength of initial noise
     */

    if ((fpr = fopen("input/composition_information", "r")) == NULL) {
        printf("Unable to open input/composition_information. Exiting");
        exit(0);
    }

    (void) fscanf(fpr, "%le%le", &c_zero, &noise_str);
    (void) fclose(fpr);

    fprintf(fpw, "c_zero = %le\n", c_zero);
    fprintf(fpw, "noise_str = %le\n", noise_str);

    (void) fclose(fpw);

    /*
     * Let us check the input data
     */

    test_input_data(n_x, n_y, delta_x, delta_y, kappa, A, delta_t,
                    time_steps, c_zero, noise_str);

    /*
     * Memory allocation to the variable comp. See the FFTW manual (version
     * 3.0) for further details
     */

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

    /*
     * Setting up the GSL random number generator. See the GSL manual for
     * further details
     */

    unsigned long int s = 2UL;
    ran_num = gsl_rng_alloc (gsl_rng_taus);
    gsl_rng_set (ran_num, s);
    /*(void) gsl_rng_env_setup();
    ran_num = gsl_rng_alloc(gsl_rng_taus);*/

    /*
     * Generating the initial microstructure
     */

    ave_comp = 0.0;
    for (i1 = 0; i1 < n_x; i1++) {
        for (i2 = 0; i2 < n_y; i2++) {
            /*
             * creal(comp[i2+n_y*i1]) = c_zero + noise_str*(0.5-gsl_rng_uniform(ran_num));
             * cimag(comp[i2+n_y*i1]) = 0.0;
             */
            comp[i2 + n_y * i1] = (fftw_complex) (c_zero + noise_str * (0.5 - gsl_rng_uniform(ran_num)));
            ave_comp = ave_comp + creal(comp[i2 + n_y * i1]);
        }
    }
    ave_comp = ave_comp / (n_x * n_y);

    /*
     * The user inputs the nominal alloy composition. However, the actual alloy 
     * composition used in the calculation will be different. To be precise,
     * if "c_zero" is the nominal alloy composition, and if "delta_c" is the
     * strength of the initial noise, the actual alloy composition would be 
     * "c_zero plus_or_minus delta_c". Here I write the actual alloy
     * composition used in the calculation to the "output/README" file.
     */

    if ((fpw = fopen("output/README", "a")) == NULL) {
        printf("Unable to open output/README. Exiting");
        exit(0);
    }

    fprintf(fpw, "ave_comp = %le\n", ave_comp);
    (void) fclose(fpw);

    /*
     * Let us evolve the microstructure
     */

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


    /*
     * Release the allocated memory
     */

    fftw_free(comp);
    gsl_rng_free(ran_num);

    return(0);
}