/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM 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 General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "N64.H"

#include "error.H"
#include "scalarMatrices.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

const Foam::scalar Foam::N64::a11[(M+1)*(M+2)/2] = {
        1.21327441805, -0.81508875, -0.06484549, 0.31926761, 0.35164219, -0.01484638,
        -0.02385399, -0.14337633, -0.18493828, 0.05466979, -0.0010004, 0.01936938,
        -0.0130214, 0.12179535, -0.05237791, -0.00172218, -0.00447904, 0.02211965,
        0.00964791, -0.07566612, 0.02527412
    };


const Foam::scalar Foam::N64::a22[(M+1)*(M+2)/2] = {
        1.3678302569, -8.59231154e-01, -7.02815662e-02, 2.79760234e-01,
        3.59026813e-01, -3.25082644e-01, -8.27621218e-03, -1.12825149e-01,
        -2.07675933e-01, 1.02433682e+00, -4.59772292e-04, 1.03264160e-02,
        -1.35152522e-02, 1.69188798e-01, -1.21839345e+00, -1.88957374e-03,
        -4.45369981e-03, 2.07928092e-02, 7.40762555e-03, -9.71847966e-02,
        5.07723424e-01
    };

//Foam::scalar Foam::N64::vec[(M+1)*(M+2)/2] = { 0.0 };
Foam::scalarField Foam::N64::vec((M+1), 0.0);

const Foam::scalar Foam::N64::digamma3 = 0.922784335098467139393487909917597568957840664060;

const Foam::scalar Foam::N64::trigamma3 = sqr(Ion::pi) / 6. - 5./4.;

/*
const Foam::scalar Foam::N64::digamma3 = 0.922784335098;

const Foam::scalar Foam::N64::trigamma3 = 0.394934066848;
*/

// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //


// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::N64::N64(const Ion& a, const Neutral& b)
:
    a_(a),
    b_(b),
    ksiij_(a_.alpha()/(sqr(a_.z())*sqrt(b_.alpha())*(1.0+pow(2.0*a_.alpha()/b_.alpha(),2./3.)))),
    sigmaij_(K1 * (pow(a_.alpha(), 1./3.) + pow(b_.alpha(), 1./3.)) / pow(a_.alpha()*b_.alpha()*(1. + 1./ksiij_), kappa)),
    piSigmaSqr_(Foam::Particle::pi*sqr(sigmaij_ * Foam::Particle::Angstrom)),
    eij_(5.2 * b_.alpha() * sqr(a_.z()) * (1.0 + ksiij_) / pow4(sigmaij_) / (a_.k/a_.e)),
    C6ij_(0.0),
    gammaij_(0.0),
    Wij_(1.0/(1.0/a_.W() + 1.0/b_.W())),
    b11_((M+1), 0.0),
    b22_((M+1), 0.0)
{
    scalar aij = a_.alpha() / b_.alpha();
    scalar aji = 1.0 / aij;
    C6ij_ = 2.0 * a_.C6() * b_.C6() / (aji * a_.C6() + aij * b_.C6());
    gammaij_ = ((2.0/sqr(a_.z())) * C6ij_ + b_.alphaQ()) / (b_.alpha() * sqr(sigmaij_));

    Info << a_.name() + " = " + b_.name() + " N64, gamma = " << gammaij_ << endl;

    scalarField temp((M+1)*(M+2)/2, 0.0); // temporary feature vector with fixed gammaij value

    transform(1.0, gammaij_, temp.data());

    Info << "fixed gammaij feature" << endl << temp << endl;

    scalarField t11 (temp * UList<scalar>(const_cast<scalar*>(a11), temp.size()));

    scalarField t22 (temp * UList<scalar>(const_cast<scalar*>(a22), temp.size()));

    scalarSquareMatrix mat11(M+1, 0.0);

    scalarSquareMatrix mat22(M+1, 0.0);

    label index = 0;

    for (label s = 0; s < M+1; s++)
    {
        for (label i = 0; i <= s; i++)
        {
           mat11(i, s-i) = t11[index];
           mat22(i, s-i) = t22[index];
           index++;
        }
    }


    for (label s = 0; s < M+1; s++)
    {
        b11_ += UList<scalar>(mat11[s], M+1);
        b22_ += UList<scalar>(mat22[s], M+1);
    }

    if (false) //  (a.hasRct(b.name()))
    {
        Info << "RCT" << endl;

        scalar A = a.Arct(b.name());

        scalar B = a.Brct(b.name());

        scalar C = A/B - digamma3 - loge(eij_);

        scalar D = sqr(B) / (a_.pi * sqr(sigmaij_)); // Formula in [A^2] (Angstrom)

        scalar a0 = D * (trigamma3 + sqr(C));

        scalar a1 = D * (-2.0 * C);

        scalar a2 = D;

        b11_[0] += a0;

        b11_[1] += a1;

        b11_[2] += a2;
    }

    /*

    for (label s = 0; s < M+1; s++)
    {
        Info << UList<scalar>(mat11[s], M+1) << endl;
    }
    Info << b11_ << endl;
    for (label s = 0; s < M+1; s++)
    {
        Info << UList<scalar>(mat22[s], M+1) << endl;
    }
    Info << b22_ << endl;

    */

}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::N64::~N64()
{}


// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //



// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //

void Foam::N64::operator=(const N64& rhs)
{
    // Check for assignment to self
    if (this == &rhs)
    {
        FatalErrorInFunction
            << "Attempted assignment to self"
            << abort(FatalError);
    }
}

// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //


// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //


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