#include "cs_defs.h"
#include <math.h>
#include <ple_coupling.h>
#include "cs_headers.h"

BEGIN_C_DECLS

void
cs_user_source_terms(cs_domain_t  *domain,
                     int           f_id,
                     cs_real_t    *st_exp,
                     cs_real_t    *st_imp)
{
    // Define parameters and variables
    const cs_field_t *fl = cs_field_by_id(f_id);
    cs_real_t *_st_exp = (cs_real_t *)st_exp;
    const cs_lnum_t n_cells = domain->mesh->n_cells;
    const cs_lnum_t n_i_faces = domain->mesh->n_i_faces;
    const cs_real_t *cell_vol = domain->mesh_quantities->cell_vol;
    const cs_real_t *restrict i_face_surf = domain->mesh_quantities->i_face_surf;
    const cs_real_3_t *i_fac_cog
      = (const cs_real_3_t *) domain->mesh_quantities->i_face_cog;
    const cs_real_3_t *i_face_normal
      = (const cs_real_3_t *)domain->mesh_quantities->i_face_normal;
    
    const int nt_cur = cs_glob_time_step->nt_cur;
    const int nt_max = cs_glob_time_step->nt_max;
    const int nt_prev = cs_glob_time_step->nt_prev;

    const cs_real_t *cvar_temp = CS_F_(t)->val;
    const cs_real_t *cvara_temp = CS_F_(t)->val_pre;

    cs_real_t total_enthalpy_prev = 0; // Previous total enthalpy
    cs_real_t total_enthalpy = 0; // Current total enthalpy
    cs_real_t total_volume = 0; // Total volume of selected cells

    // Define selection criteria (for example, selecting cells based on their coordinates)
    const double x_min = 0.0; // Minimum x-coordinate
    const double x_max = 1.0; // Maximum x-coordinate
    const double y_min = 0.0; // Minimum y-coordinate
    const double y_max = 0.75; // Maximum y-coordinate
    const double z_min = 0.3; // Minimum z-coordinate
    const double z_max = 1.0; // Maximum z-coordinate


    // Calculate total enthalpy and total volume in the selected part of the domain
    for (cs_lnum_t i = 0; i < n_cells; i++) {
        // Get the coordinates of the cell's centroid
        const cs_real_t temp_cell = cvar_temp[i];
        const cs_real_t temp_pre_cell = cvara_temp[i];
        const cs_real_t *centroid = &(domain->mesh_quantities->cell_cen[i*3]);
        const double x = centroid[0];
        const double y = centroid[1];
        const double z = centroid[2];
        // Check if the cell's centroid is within the selected region
        if (x >= x_min && x <= x_max && y >= y_min && y <= y_max && z >= z_min && z <= z_max) {
            total_enthalpy += temp_cell * cell_vol[i]; // Assuming temperature is stored in CS_F_(temperature)
            total_enthalpy_prev += temp_pre_cell * cell_vol[i];
            total_volume += cell_vol[i];
        }
    }
    

    if (cs_glob_rank_id >= 0) {
    	cs_parall_sum(1, CS_REAL_TYPE, &total_volume);
    	cs_parall_sum(1, CS_REAL_TYPE, &total_enthalpy);
        cs_parall_sum(1, CS_REAL_TYPE, &total_enthalpy_prev);
    }

  
    // Compute feedback term based on the change in total enthalpy per unit volume
    cs_real_t feedback_term = 0;
    if (nt_cur > 1 && total_enthalpy_prev > 0) {
        cs_real_t delta_enthalpy = total_enthalpy - total_enthalpy_prev;
        feedback_term = delta_enthalpy / total_volume;
    }
    /*total_enthalpy_prev = total_enthalpy;*/


    // Apply feedback term to the temperature field in the selected part of the domain
    if (f_id == CS_F_(t)->id) {
        for (cs_lnum_t i = 0; i < n_cells; i++) {
            const cs_real_t *centroid = &(domain->mesh_quantities->cell_cen[i*3]);
            const double x = centroid[0];
            const double y = centroid[1];
            const double z = centroid[2];
            if (x >= x_min && x <= x_max && y >= y_min && y <= y_max && z >= z_min && z <= z_max) {
                _st_exp[i] = -feedback_term; // Source or sink term applied to temperature
            }
        }
    }
    bft_printf("feedback_term: %e total_enthalpy: %e total_enthalpy_prev: %e", feedback_term, total_enthalpy, total_enthalpy_prev);

}

END_C_DECLS