/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "mkernel_block.h"

#ifdef GCTL_POTENTIAL_AUTODIFF

autodiff::var magobser_block_thickness_za_sig(const gctl::array<gctl::mag_block> &ele, const gctl::array<double> &sus, 
    const gctl::point3dc &obsp, const autodiff::ArrayXvar& btm);

double gctl::magobser_wrt_thickness(array<double> &btm_derivatives, loss_func &lf, 
    const array<point3dc> &obsp, const array<double> &sus, const array<mag_block> &ele, 
    magnetic_field_type_e comp_id, verbose_type_e verbose)
{
    if (ele[0].att == nullptr) throw gctl::runtime_error("[GCTL::Magkernel_Block] Magnetization parameter not set.");

	int i, j;
	int o_size = obsp.size();
	int e_size = ele.size();

    btm_derivatives.resize(e_size, 0.0);

    double tmp_val, tmp_grad;
    autodiff::var tmp_obs;
    autodiff::ArrayXvar btm_depth;
    Eigen::VectorXd btm_ders;

	gctl::progress_bar bar(o_size, "Za_wrt_thickness (openmp)");
    if (verbose == gctl::FullMsg) bar.progressed(0);
	else if (verbose == gctl::ShortMsg) bar.progressed_simple(0);

#pragma omp parallel for private (i, j, tmp_val, tmp_grad, tmp_obs, btm_depth, btm_ders) shared(e_size) schedule(guided)
    for (i = 0; i < o_size; i++)
    {
        btm_depth.resize(e_size);
        for (j = 0; j < e_size; j++)
        {
            btm_depth[j] = ele[j].dl->z;
        }

        btm_ders.setZero();
        tmp_obs = magobser_block_thickness_za_sig(ele, sus, obsp[i], btm_depth);
        btm_ders = autodiff::gradient(tmp_obs, btm_depth);
        tmp_val = autodiff::val(tmp_obs);
        tmp_grad = lf.gradient(tmp_val, i);

#pragma omp critical
{
        lf.evaluate(tmp_val, i);
        for (j = 0; j < e_size; j++)
        {
            btm_derivatives[j] += tmp_grad*btm_ders[j];
        }
}
        //if (verbose == gctl::FullMsg) bar.progressed(i);
	    //else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
    }

    if (verbose == gctl::FullMsg) bar.progressed(o_size - 1);
	else if (verbose == gctl::ShortMsg) bar.progressed_simple(o_size - 1);

    return lf.get_loss();
}

autodiff::var magobser_block_thickness_za_sig(const gctl::array<gctl::mag_block> &ele, const gctl::array<double> &sus, 
    const gctl::point3dc &obsp, const autodiff::ArrayXvar& btm)
{
    double I,A;
	double Alpha,Beta,Gamma;
    double x, y, z, su;
	double x1,x2,y1,y2,z2;
	double w111_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1;
	double w222_2,w112_2,w122_2,w212_2;
	double R222,R122,R212,R112;
	double G222,G122,G212,G112;

    autodiff::var obs = 0.0;
    autodiff::var w211_2,w121_2,w221_2,w111_2;
    autodiff::var R221,R121,R211,R111;
    autodiff::var G221,G121,G211,G111;

	gctl::magblock_para* mp;
    for (int i = 0; i < ele.size(); i++)
    {
        mp = ele[i].att;
        I = mp->inclina_deg*M_PI/180.0;
        A = (90 - mp->declina_deg)*M_PI/180.0;

        Alpha=cos(I)*cos(A);
        Beta=cos(I)*sin(A);
        Gamma=sin(I);

        x = obsp.x;
        y = obsp.y;
        z = obsp.z;
        su = sus[i];

        x1 = ele[i].dl->x; x2 = ele[i].ur->x;
        y1 = ele[i].dl->y; y2 = ele[i].ur->y;
                           z2 = ele[i].ur->z;

        R222=sqrt((x2-x)*(x2-x)+(y2-y)*(y2-y)+(z2-z)*(z2-z));
        R122=sqrt((x1-x)*(x1-x)+(y2-y)*(y2-y)+(z2-z)*(z2-z));
        R212=sqrt((x2-x)*(x2-x)+(y1-y)*(y1-y)+(z2-z)*(z2-z));
        R112=sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y)+(z2-z)*(z2-z));
        R221=sqrt((x2-x)*(x2-x)+(y2-y)*(y2-y)+(btm[i] - z)*(btm[i] - z));
        R121=sqrt((x1-x)*(x1-x)+(y2-y)*(y2-y)+(btm[i] - z)*(btm[i] - z));
        R211=sqrt((x2-x)*(x2-x)+(y1-y)*(y1-y)+(btm[i] - z)*(btm[i] - z));
        R111=sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y)+(btm[i] - z)*(btm[i] - z));

        w222_1=(x2-x)*(y2-y);
        w122_1=(x1-x)*(y2-y);
        w212_1=(x2-x)*(y1-y);
        w112_1=(x1-x)*(y1-y);
        w221_1=(x2-x)*(y2-y);
        w121_1=(x1-x)*(y2-y);
        w211_1=(x2-x)*(y1-y);
        w111_1=(x1-x)*(y1-y);

        w222_2=R222*(z2-z);
        w122_2=R122*(z2-z);
        w212_2=R212*(z2-z);
        w112_2=R112*(z2-z);
        w221_2=R221*(btm[i] - z);
        w121_2=R121*(btm[i] - z);
        w211_2=R211*(btm[i] - z);
        w111_2=R111*(btm[i] - z);
        
        G222=-Gamma*atan2(w222_1,w222_2)+Alpha*log(R222+y2-y)+Beta*log(R222+x2-x);
        G122=-Gamma*atan2(w122_1,w122_2)+Alpha*log(R122+y2-y)+Beta*log(R122+x1-x);
        G212=-Gamma*atan2(w212_1,w212_2)+Alpha*log(R212+y1-y)+Beta*log(R212+x2-x);
        G112=-Gamma*atan2(w112_1,w112_2)+Alpha*log(R112+y1-y)+Beta*log(R112+x1-x);
        G221=-Gamma*atan2(w221_1,w221_2)+Alpha*log(R221+y2-y)+Beta*log(R221+x2-x);
        G121=-Gamma*atan2(w121_1,w121_2)+Alpha*log(R121+y2-y)+Beta*log(R121+x1-x);
        G211=-Gamma*atan2(w211_1,w211_2)+Alpha*log(R211+y1-y)+Beta*log(R211+x2-x);
        G111=-Gamma*atan2(w111_1,w111_2)+Alpha*log(R111+y1-y)+Beta*log(R111+x1-x);

        obs += GCTL_T0*su/(4*M_PI)*(G222-G122-G212+G112-G221+G121+G211-G111);
    }

    return obs;
}

#endif // GCTL_POTENTIAL_AUTODIFF