/***********************************************************************************
                 OpenCS Project: www.daetools.com
                 Copyright (C) Dragan Nikolic
************************************************************************************
OpenCS is free software; you can redistribute it and/or modify it under the terms
of the GNU Lesser General Public License version 3 as published by the Free Software
Foundation. OpenCS is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License along with
the OpenCS software; if not, see <http://www.gnu.org/licenses/>.
***********************************************************************************/
#ifndef DIURNAL_KINETICS_KERNELS_MODEL_H
#define DIURNAL_KINETICS_KERNELS_MODEL_H

#include <string>
#include <vector>
#include <map>
#include <set>
#include <vector>
#include <iostream>
#include <math.h>
#include <OpenCS/models/cs_number.h>
#include <OpenCS/models/cs_kernel.h>
using namespace cs;

// An auxiliary class to handle the centered finite difference scheme on a 2D domain
// Homogenous Neumann BCs are assumed at all four edges: dudn|boundary = 0.0.
class DiurnalKineticsGroups_2D
{
public:
    DiurnalKineticsGroups_2D(int nx, int ny):
        Nx(nx), Ny(ny)
    {
        x0 =  0.0;
        x1 = 20.0;
        y0 = 30.0;
        y1 = 50.0;
        dx = (x1-x0) / (Nx-1);
        dy = (y1-y0) / (Ny-1);

        x_domain.resize(Nx);
        y_domain.resize(Ny);
        for(int x = 0; x < Nx; x++)
            x_domain[x] = x0 + x * dx;
        for(int y = 0; y < Ny; y++)
            y_domain[y] = y0 + y * dy;

        C1_start_index = 0*Nx*Ny;
        C2_start_index = 1*Nx*Ny;

        C1_values = NULL;
        C2_values = NULL;

        Nequations = 2*Nx*Ny;

        V   =  1.00E-03;
        Kh  =  4.00E-06;
        Kv0 =  1.00E-08;
        q1  =  1.63E-16;
        q2  =  4.66E-16;
        C3  =  3.70E+16;
        a3  = 22.62;
        a4  =  7.601;
    }

    void SetInitialConditions(std::vector<real_t>& C0)
    {
        int ix, iy;
        C0.assign(2*Nx*Ny, 0.0);

        real_t* C1_0 = &C0[C1_start_index];
        real_t* C2_0 = &C0[C2_start_index];

        for(ix = 0; ix < Nx; ix++)
        {
            for(iy = 0; iy < Ny; iy++)
            {
                int index = getIndex(ix,iy);

                real_t x = x_domain[ix];
                real_t y = y_domain[iy];

                C1_0[index] = 1E6  * alfa(x) * beta(y);
                C2_0[index] = 1E12 * alfa(x) * beta(y);
            }
        }
    }

    void GetVariableNames(std::vector<std::string>& names)
    {
        const int bsize = 32;
        char buffer[bsize];
        int index = 0;

        names.resize(Nequations);
        for(int x = 0; x < Nx; x++)
        {
            for(int y = 0; y < Ny; y++)
            {
                std::snprintf(buffer, bsize, "%s(%d,%d)", "C1", x, y);
                names[index] = buffer;
                index++;
            }
        }
        for(int x = 0; x < Nx; x++)
        {
            for(int y = 0; y < Ny; y++)
            {
                std::snprintf(buffer, bsize, "%s(%d,%d)", "C2", x, y);
                names[index] = buffer;
                index++;
            }
        }
    }

    void CreateEquations(const std::vector<csNumber_t>& C_values,
                         const csNumber_t& time,
                         std::vector<csEquation_t>& equations)
    {
        if(C_values.size() != Nequations)
            throw std::runtime_error("Invalid size of data arrays");

        C1_values = &C_values[C1_start_index];
        C2_values = &C_values[C2_start_index];

        group_C1 = csGroupPtr(new csGroup_t("DiurnalKinetics_C1", 1)); // the group for C1-component (group 1)
        group_C2 = csGroupPtr(new csGroup_t("DiurnalKinetics_C2", 2)); // the group for C2-component (group 2)

        csEquation_t equation_C1(group_C1.get()); // Belongs to the group "DiurnalKinetics_C1"
        csEquation_t equation_C2(group_C2.get()); // Belongs to the group "DiurnalKinetics_C2"

        for(int x = 0; x < Nx; x++)
        {
            for(int y = 0; y < Ny; y++)
            {
                /* Component 1 */
                csNumber_t dC1_dt = V  * dC1_dx(x,y) +                              /* x-axis convection term */
                                    Kh * d2C1_dx2(x,y) +                            /* x-axis diffusion term  */
                                    Kv(y) * (0.2 * dC1_dy(x,y) + d2C1_dy2(x,y)) +   /* y-axis diffusion term  */
                                    R1(C1(x,y), C2(x,y), time);                     /* generation term        */

                equation_C1[ C1(x,y) ] = dC1_dt;
                equations.push_back(equation_C1);
            }
        }
        for(int x = 0; x < Nx; x++)
        {
            for(int y = 0; y < Ny; y++)
            {
                /* Component 2 */
                csNumber_t dC2_dt = V  * dC2_dx(x,y) +                              /* x-axis convection term */
                                    Kh * d2C2_dx2(x,y) +                            /* x-axis diffusion term  */
                                    Kv(y) * (0.2 * dC2_dy(x,y) + d2C2_dy2(x,y)) +   /* y-axis diffusion term  */
                                    R2(C1(x,y), C2(x,y), time);                     /* generation term        */

                equation_C2[ C2(x,y) ] = dC2_dt;
                equations.push_back(equation_C2);
            }
        }
    }

protected:
    int getIndex(int x, int y)
    {
        if(x < 0 || x >= Nx)
            throw std::runtime_error("Invalid x index");
        if(y < 0 || y >= Ny)
            throw std::runtime_error("Invalid y index");
        return Ny*x + y;
    }

    real_t alfa(real_t x)
    {
        real_t xmid = (x1+x0)/2;
        real_t cx = 0.1 * (x - xmid);
        real_t cx2 = cx*cx;
        return 1 - cx2 + 0.5*cx2*cx2;
    }
    real_t beta(real_t y)
    {
        real_t ymid = (y1+y0)/2;
        real_t cy = 0.1 * (y - ymid);
        real_t cy2 = cy*cy;
        return 1 - cy2 + 0.5*cy2*cy2;
    }

    csNumber_t Kv(int y)
    {
        return Kv0 * std::exp(y_domain[y] / 5.0);
    }

    csNumber_t q3(const csNumber_t& time)
    {
        real_t w = std::acos(-1.0) / 43200.0;
        csNumber_t sinwt = cs::fmax(1E-30, cs::sin(w*time));
        return cs::exp(-a3 / sinwt);
    }

    csNumber_t q4(const csNumber_t& time)
    {
        real_t w = std::acos(-1.0) / 43200.0;
        csNumber_t sinwt = cs::fmax(1E-30, cs::sin(w*time));
        return cs::exp(-a4 / sinwt);
    }

    csNumber_t R1(const csNumber_t& c1, const csNumber_t& c2, const csNumber_t& time)
    {
        return -q1*c1*C3 - q2*c1*c2 + 2*q3(time)*C3 + q4(time)*c2;
    }

    csNumber_t R2(const csNumber_t& c1, const csNumber_t& c2, const csNumber_t& time)
    {
        return q1*c1*C3 - q2*c1*c2 - q4(time)*c2;
    }

    const csNumber_t& C1(int x, int y)
    {
        int index = getIndex(x, y);
        return C1_values[index];
    }
    const csNumber_t& C2(int x, int y)
    {
        int index = getIndex(x, y);
        return C2_values[index];
    }

    // First order partial derivative per x.
    csNumber_t dC1_dx(int x, int y)
    {
        // If called for x == 0 or x == Nx-1 return 0.0 (zero flux through boundaries).
        if(x == 0 || x == Nx-1)
            return csNumber_t(0.0);

        const csNumber_t& ci1 = C1(x+1, y);
        const csNumber_t& ci2 = C1(x-1, y);
        return (ci1 - ci2) / (2*dx);
    }
    csNumber_t dC2_dx(int x, int y)
    {
        // If called for x == 0 or x == Nx-1 return 0.0 (zero flux through boundaries).
        if(x == 0 || x == Nx-1)
            return csNumber_t(0.0);

        const csNumber_t& ci1 = C2(x+1, y);
        const csNumber_t& ci2 = C2(x-1, y);
        return (ci1 - ci2) / (2*dx);
    }

    // First order partial derivative per y.
    csNumber_t dC1_dy(int x, int y)
    {
        // If called for y == 0 or y == Ny-1 return 0.0 (zero flux through boundaries).
        if(y == 0 || y == Ny-1)
            return csNumber_t(0.0);

        const csNumber_t& ci1 = C1(x, y+1);
        const csNumber_t& ci2 = C1(x, y-1);
        return (ci1 - ci2) / (2*dy);
    }
    csNumber_t dC2_dy(int x, int y)
    {
        // If called for y == 0 or y == Ny-1  return 0.0 (zero flux through boundaries).
        if(y == 0 || y == Ny-1)
            return csNumber_t(0.0);

        const csNumber_t& ci1 = C2(x, y+1);
        const csNumber_t& ci2 = C2(x, y-1);
        return (ci1 - ci2) / (2*dy);
    }

    // Second order partial derivative per x.
    csNumber_t d2C1_dx2(int x, int y)
    {
        // If called for x == 0 or x == Nx-1 return 0.0 (no diffusion through boundaries).
        if(x == 0 || x == Nx-1)
            return csNumber_t(0.0);

        const csNumber_t& ci1 = C1(x+1, y);
        const csNumber_t& ci  = C1(x,   y);
        const csNumber_t& ci2 = C1(x-1, y);
        return (ci1 - 2*ci + ci2) / (dx*dx);
    }
    csNumber_t d2C2_dx2(int x, int y)
    {
        // If called for x == 0 or x == Nx-1 return 0.0 (no diffusion through boundaries).
        if(x == 0 || x == Nx-1)
            return csNumber_t(0.0);

        const csNumber_t& ci1 = C2(x+1, y);
        const csNumber_t& ci  = C2(x,   y);
        const csNumber_t& ci2 = C2(x-1, y);
        return (ci1 - 2*ci + ci2) / (dx*dx);
    }

    // Second order partial derivative per y.
    csNumber_t d2C1_dy2(int x, int y)
    {
        // If called for y == 0 or y == Ny-1 return 0.0 (no diffusion through boundaries).
        if(y == 0 || y == Ny-1)
            return csNumber_t(0.0);

        const csNumber_t& ci1 = C1(x, y+1);
        const csNumber_t& ci  = C1(x,   y);
        const csNumber_t& ci2 = C1(x, y-1);
        return (ci1 - 2*ci + ci2) / (dy*dy);
    }
    csNumber_t d2C2_dy2(int x, int y)
    {
        // If called for y == 0 or y == Ny-1 return 0.0 (no diffusion through boundaries).
        if(y == 0 || y == Ny-1)
            return csNumber_t(0.0);

        const csNumber_t& ci1 = C2(x, y+1);
        const csNumber_t& ci  = C2(x,   y);
        const csNumber_t& ci2 = C2(x, y-1);
        return (ci1 - 2*ci + ci2) / (dy*dy);
    }

public:
    int    Nequations;
    int    Nx;
    int    Ny;
    real_t V;
    real_t Kh;
    real_t Kv0;
    real_t q1;
    real_t q2;
    real_t C3;
    real_t a3;
    real_t a4;

protected:
    int    C1_start_index;
    int    C2_start_index;
    real_t x0, x1;
    real_t y0, y1;
    real_t dx;
    real_t dy;
    std::vector<real_t> x_domain;
    std::vector<real_t> y_domain;
    const csNumber_t* C1_values;
    const csNumber_t* C2_values;

    csGroupPtr group_C1;
    csGroupPtr group_C2;

};

#endif