/*
! Copyright:
!    This software was developed within the context of
!    the EUMETSAT Satellite Application Facility on
!    Numerical Weather Prediction (NWP SAF), under the
!    Cooperation Agreement dated 7 September 2021, between
!    EUMETSAT and the Met Office, UK, by one or more partners
!    within the NWP SAF. The partners in the NWP SAF are
!    the Met Office, ECMWF, DWD and MeteoFrance.
!
!    Copyright 2024, EUMETSAT, All Rights Reserved.
*/

///
/// @file RttovSafe.cpp
///
/// @brief Wrapper class for RTTOV
///

#include "RttovSafe.h"

namespace rttov {

///@brief RttovSafe class constructor method
RttovSafe::RttovSafe() : Rttov::Rttov(), profilesP(NULL) {

}

/// @brief Deallocate internal arrays and drop the coefficients
RttovSafe::~RttovSafe() {
    delete profiles;
    delete [] p_half;
    delete [] p;
    delete [] t;
    delete [] datetimes;
    delete [] angles;
    delete [] surface_fraction;
    delete [] surftype;
    delete [] skin;
    delete [] near_surface;
    delete [] surfgeom;
    delete [] simplecloud;
    delete [] clwde_param;
    delete [] icede_param;
    delete [] hydro_frac_eff;
    delete [] zeeman;
    delete [] gases;
    delete [] gas_id;
}


///@brief Return the associated profiles
std::vector<Profile>* RttovSafe::getTheProfiles() const {
    if (isProfileSet()) {
        return profilesP;
    } else {
        return NULL;
    }
}

/// @brief Associate a vector of Profile objects with this RttovSafe object; 
///        carries out checks on profiles before calling RTTOV to help prevent errors:
///        all profiles must be have the same number of levels with the same content (gases, hydrometeors, aerosols)
///        and have the same gas_units
void RttovSafe::setTheProfiles(std::vector<Profile>& profiles) {

    profileSet = false;

    profilesP = &profiles;

    nprofiles = profiles.size();
    if (nprofiles == 0) throw std::logic_error("Error: no Profile instances in the vector of profiles");

    bool isOK = profiles[0].check();
    if (! isOK ) throw std::logic_error("Error: profile 0 is not OK");

    nlevels = profiles[0].getNlevels();
    nlayers = nlevels - 1;
    nsurfaces = profiles[0].getNsurfaces();

    std::vector <itemIdType> Profilecontents = profiles[0].getItemContents();
    ngases = profiles[0].getItemContents().size();
    if (debug) std::cerr << "ngases = " << ngases << std::endl;

    gas_units = profiles[0].getGasUnits();
    mmr_hydro = profiles[0].isMmrHydro();
    mmr_aer = profiles[0].isMmrAer();

    for (int i = 1; i < nprofiles; i++) {
        if (profiles[i].getNlevels() != nlevels)
            throw std::logic_error("Error: Profiles must all have the same number of levels");
        if (profiles[i].getNsurfaces() != nsurfaces)
            throw std::logic_error("Error: Profiles must all have the same number of surfaces");
        if (profiles[i].getItemContents() != Profilecontents)
            throw std::logic_error("Error: Profiles must all have the same gas/hydrometeor/aerosol data defined");
        if (! profiles[i].check() )
            throw std::logic_error("Error: Profile check failed");
        if (profiles[i].getGasUnits() != gas_units)
            throw std::logic_error("Error: Profiles must all have the same gas_units");
        if (profiles[i].isMmrHydro() != mmr_hydro)
            throw std::logic_error("Error: Profiles must all have the same mmr_hydro");
        if (profiles[i].isMmrAer() != mmr_aer)
            throw std::logic_error("Error: Profiles must all have the same mmr_aer");
    }
    if ( gas_units == unknown ) {
        std::cerr << "Profile warning: gas_units not specified, assuming kg/kg wet" << std::endl;
        gas_units = kg_per_kg;
    }

    // convert to pointers to doubles for the wrapper
    std::vector <double> vOfd;
    std::vector <int>  vOfi;
    delete [] p_half;
    p_half = new double[nprofiles*nlevels]();
    delete [] p;
    p = new double[nprofiles*nlayers]();
    delete [] t;
    t = new double[nprofiles*nlayers]();
    delete [] t;
    t = new double[nprofiles*nlayers]();
    delete [] datetimes;
    datetimes = new int[nprofiles*6]();
    delete [] angles;
    angles = new double[nprofiles*4]();
    delete [] surfgeom;
    surfgeom = new double[nprofiles*3]();
    delete [] surface_fraction;
    surface_fraction = new double[nprofiles*(nsurfaces-1)]();
    delete [] surftype;
    surftype = new int[nprofiles*nsurfaces*2]();
    delete [] skin;
    skin = new double[nprofiles*nsurfaces*9]();
    delete [] near_surface;
    near_surface = new double[nprofiles*nsurfaces*5]();
    delete [] simplecloud;
    simplecloud = new double[nprofiles*2]();
    delete [] clwde_param;
    clwde_param = new int[nprofiles]();
    delete [] icede_param;
    icede_param = new int[nprofiles]();
    delete [] hydro_frac_eff;
    hydro_frac_eff = new double[nprofiles]();
    delete [] zeeman;
    zeeman = new double[nprofiles*2]();
    delete [] gases;
    gases = new double[ngases*nprofiles*nlayers]();
    delete [] gas_id;
    gas_id = new int[ngases]();

    for (int i = 0; i < nprofiles; i++) {

        vOfd = profiles[i].getPHalf();
        for (int j = 0; j < vOfd.size(); j++) {
            p_half[nlevels*i+j] = vOfd[j];
        }

        vOfd = profiles[i].getP();
        for (int j = 0; j < vOfd.size(); j++) {
            p[nlayers*i+j] = vOfd[j];
        }

        vOfd = profiles[i].getT();
        for (int j = 0; j < vOfd.size(); j++) {
            t[nlayers*i+j] = vOfd[j];
        }

        vOfi=(profiles[i].getDateTimes());
        for (int j = 0; j < vOfi.size(); j++) {
            datetimes[6*i+j] = vOfi[j];
        }

        vOfd = profiles[i].getAngles();
        for (int j = 0; j < vOfd.size(); j++) {
            angles[4*i+j] = vOfd[j];
        }

        vOfd = profiles[i].getSurfGeom();
        for (int j = 0; j < vOfd.size(); j++) {
            surfgeom[3*i+j] = vOfd[j];
        }

        vOfd = profiles[i].getSurfaceFraction();
        for (int j = 0; j < vOfd.size(); j++) {
            surface_fraction[(nsurfaces-1)*i+j] = vOfd[j];
        }

        for (int isurf = 0; isurf < nsurfaces; isurf++) {
            vOfi = profiles[i].getSurfType(isurf);
            for (int j = 0; j < vOfi.size(); j++) {
                surftype[nsurfaces*2*i+2*isurf+j] = vOfi[j];
            }

            vOfd = profiles[i].getSkin(isurf);
            for (int j = 0; j < vOfd.size(); j++) {
                skin[nsurfaces*9*i+9*isurf+j] = vOfd[j];
            }

            vOfd = profiles[i].getNearSurface(isurf);
            for (int j = 0; j < vOfd.size(); j++) {
                near_surface[nsurfaces*5*i+5*isurf+j] = vOfd[j];
            }
        }

        vOfd = profiles[i].getSimpleCloud();
        for (int j = 0; j < vOfd.size(); j++) {
            simplecloud[2*i+j] = vOfd[j];
        }

        clwde_param[i] = profiles[i].getClwdeParam();
        icede_param[i] = profiles[i].getIcedeParam();
        hydro_frac_eff[i] = profiles[i].getHydroFracEff();

        vOfd = profiles[i].getZeeman();
        for (int j = 0; j < vOfd.size(); j++) {
            zeeman[2*i+j] = vOfd[j];
        }
    }

    delete this->profiles;
    this->profiles = new Profiles(nprofiles, nlevels, 1);

    ItemIdPointerMap myItems;
    itemIdType key;
    int i = 0;
    for (int prof = 0; prof < nprofiles; prof++) {
        i = 0;
        if (debug) std::cerr << "Profile " << prof << std::endl;
        myItems = profiles[prof].getItems();
        for (ItemIdPointerMap::iterator iter = myItems.begin(); iter != myItems.end(); ++iter) {
            key = iter->first;
            if (myItems.count(key) > 0) {
                if (prof == 0) gas_id[i] = key;
                for (int lay = 0; lay < nlayers; lay++) {
                    gases[i*nprofiles*nlayers + nlayers*prof + lay] = myItems[key][lay];
                    // Use the Rttov->profiles object gas_index
                    this->profiles->gas_index[key] = i;
                }
            }
            i++;
        }
    }

    profileSet = true;

    if (debug) {
        for (int i = 0; i < ngases; i++) std::cerr << "gas_id : " << gas_id[i] << std::endl;
    }
}

} /* namespace rttov */