/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | faSavageHutterFOAM
    \\  /    A nd           | Copyright (C) 2017 Matthias Rauter
     \\/     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/>.

Author
    Matthias Rauter matthias.rauter@uibk.ac.at

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

#include "radar.H"
#include "volFields.H"
#include "dictionary.H"
#include "foamTime.H"
#include "areaMesh.H"
#include "uniformDimensionedFields.H"
#include "edgeFields.H"
#include "cfdfunction.H"
#include "volSurfaceMapping.H"
#include "faCFD.H"

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

namespace Foam
{
    defineTypeNameAndDebug(radar, 0);
}

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



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

Foam::radar::radar
(
    const word& name,
    const objectRegistry& obr,
    const dictionary& dict,
    const bool loadFromFiles
):
    obr_(obr),
    radarFilePtr_(NULL)
{

    read(dict);


    Info << endl << "Calculating radar:" << endl
         << "   log is " << log_ << endl
         << endl << endl;


}


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

Foam::radar::~radar()
{}


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

void Foam::radar::read(const dictionary& dict)
{
    log_ = Switch(dict.lookup("log"));
    pos_ = vector(dict.lookup("position"));
    dx_ = scalar(dict.lookupOrDefault<scalar>("deltaX", 1.0));
    xmin_ = scalar(dict.lookupOrDefault<scalar>("xmin", 0));
    xmax_ = scalar(dict.lookupOrDefault<scalar>("xmax", 30));
    hmin_ = scalar(dict.lookupOrDefault<scalar>("hmin", 1e-2));
    word intensity_desc = word(dict.lookupOrDefault<word>("intensity", "sum_velocity"));
    n_gates_ = (xmax_ - xmin_)/dx_;

    if (intensity_desc == "sum_velocity")
        intensity_ = sum_velocity;
    else if (intensity_desc == "sum_volume")
        intensity_ = sum_volume;
    else if (intensity_desc == "sum_momentum")
        intensity_ = sum_momentum;
    else if (intensity_desc == "max_height")
        intensity_ = max_height;
    else if (intensity_desc == "max_velocity")
        intensity_ = max_velocity;
    else if (intensity_desc == "div_momentum")
        intensity_ = div_momentum;
    else
    {
        Warning << "unknown intensity desceription \"" << intensity_desc << "\"" << endl
                << "    taking sum_momentum." << endl;
        intensity_desc = "sum_momentum";
        intensity_ = sum_momentum;
    }

    Info << "running radar with " << endl
         << "  log = " << log_ << endl
         << "  int = " << intensity_desc << endl
         << "  pos = " << pos_ << endl
         << "  dx = "  << dx_ << endl
         << "  xmax = "  << xmin_ << endl
         << "  xmax = "  << xmax_ << endl
         << "  hmin = " << hmin_ << endl
         << "  n_gates = "  << n_gates_ << endl;

    if(Pstream::master())
    {
        radarFilePtr_.reset(new OFstream(word(dict.lookup("tableOutPut"))));
        radarFilePtr_() << "time"        << ";";

        int i;
        for(i=0; i<n_gates_-1; i++)
        {
            radarFilePtr_() << xmin_+dx_*i << ";";
        }

        radarFilePtr_() << xmin_+dx_*i++ << endl;
    }

}

void Foam::radar::execute()
{
    const areaScalarField &h = obr_.lookupObject<areaScalarField>("h");
    const areaVectorField &Us = obr_.lookupObject<areaVectorField>("Us");

    const faMesh &aMesh = h.mesh();
    Field<scalar> V;// = aMesh.S()*mag(Us.internalField());
    if (intensity_ == sum_velocity)
        V = aMesh.S()*mag(Us.internalField());
    else if (intensity_ == sum_volume)
        V = aMesh.S()*h.internalField();
    else if (intensity_ == sum_momentum)
        V = aMesh.S()*mag(Us.internalField())*h.internalField();
    else if (intensity_ == max_height)
        V = h.internalField();
    else if (intensity_ == max_velocity)
        V = mag(Us.internalField());
    else if (intensity_ == div_momentum)
        V = mag(fac::div(Us*h));


    const areaVectorField &C = aMesh.areaCentres();

    List<scalar> intensity(n_gates_);
    int i;
    for(i=0; i<n_gates_; i++)
    {
        intensity[i] = 0;
    }

    forAll(C, i)
    {
        if(h[i] > hmin_)
        {            
            scalar d = mag(C[i]-pos_);
            label gate = (d-xmin_)/dx_;
            if (gate > -1 && gate < n_gates_)
            {
                if (intensity_ == max_height || intensity_ == max_velocity)
                    intensity[gate] = max(intensity[gate], V[i]);
                else
                    intensity[gate] += V[i];

            }
        }
    }

    reduce( intensity, sumOp<List<scalar> >() );

    if(Pstream::master())
    {
        radarFilePtr_() << obr_.time().value() << ";";
        for(i=0; i<n_gates_-1; i++)
        {
            radarFilePtr_() << intensity[i] << ";";
        }

        radarFilePtr_() << intensity[i++] << endl;
    }

    if (log_)
    {

    }
}


void Foam::radar::end()
{

}


void Foam::radar::write()
{

}

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