#include "macrodef.h"
#include "share_func.h"
#include <cstddef>

/*
 * fderivs_shc — shell first derivatives converted to Cartesian via chain rule.
 *
 * Calls fderivs_sh internally, then:
 *   fx = drhodx * df/drho + dsigmadx * df/dsigma + dRdx * df/dR
 *   fy = drhody * df/drho + dsigmady * df/dsigma + dRdy * df/dR
 *   fz = drhodz * df/drho + dsigmadz * df/dsigma + dRdz * df/dR
 */

// Forward declaration (defined in fderivs_sh_c.C with extern "C" name fderivs_sh_)
extern "C" {
void fderivs_sh_(const int ex[3], const double *f,
                 double *fx, double *fy, double *fz,
                 const double *X, const double *Y, const double *Z,
                 double SYM1, double SYM2, double SYM3,
                 int Symmetry, int onoff, int sst);

void fderivs_shc_(int *ex,
                   double *f,
                   double *fx, double *fy, double *fz,
                   double *crho, double *sigma, double *R,
                   double &SYM1, double &SYM2, double &SYM3,
                   int &Symmetry, int &Lev, int &sst,
                   double *drhodx, double *drhody, double *drhodz,
                   double *dsigmadx, double *dsigmady, double *dsigmadz,
                   double *dRdx, double *dRdy, double *dRdz)
{
    const int ex3[3] = { ex[0], ex[1], ex[2] };
    const size_t n = (size_t)ex[0] * (size_t)ex[1] * (size_t)ex[2];

    // Temporary shell-coordinate derivatives
    double *gx = (double*)malloc(n * sizeof(double));
    double *gy = (double*)malloc(n * sizeof(double));
    double *gz = (double*)malloc(n * sizeof(double));
    if (!gx || !gy || !gz) { free(gx); free(gy); free(gz); return; }

    // Compute shell-coordinate derivatives
    fderivs_sh_(ex3, f, gx, gy, gz, crho, sigma, R, SYM1, SYM2, SYM3, Symmetry, Lev, sst);

    // Chain rule to Cartesian
    for (size_t i = 0; i < n; ++i) {
        fx[i] = drhodx[i]   * gx[i] + dsigmadx[i] * gy[i] + dRdx[i] * gz[i];
        fy[i] = drhody[i]   * gx[i] + dsigmady[i] * gy[i] + dRdy[i] * gz[i];
        fz[i] = drhodz[i]   * gx[i] + dsigmadz[i] * gy[i] + dRdz[i] * gz[i];
    }

    free(gx); free(gy); free(gz);
}

} // extern "C"