#include "cs_defs.h"

/*----------------------------------------------------------------------------
* Standard C library headers
*----------------------------------------------------------------------------*/

#include <assert.h>
#include <math.h>
#include <stdio.h>  // Added for debugging

/*----------------------------------------------------------------------------
* PLE library headers
*----------------------------------------------------------------------------*/

#include <ple_coupling.h>

/*----------------------------------------------------------------------------
* Local headers
*----------------------------------------------------------------------------*/

#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)
{
   /* Field structure */
   const cs_field_t  *f = cs_field_by_id(f_id);
   cs_real_3_t *_st_exp = (cs_real_3_t *)st_exp;

   /* Check if the field is temperature */
   if (f != CS_F_(t)) {
       return; /* Apply only for the temperature field */
   }

   /* Mesh and material properties */
   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 cs_real_t *cvar_temp = CS_F_(t)->val;
   const cs_real_3_t *cvar_vel = (const cs_real_3_t *)CS_F_(vel)->val;

   const cs_real_t cp = 1005.0;             /* Specific heat capacity [J/(kg.K)] */
   const cs_real_t L = 1.6;                 /* Characteristic length [m] (example) */
   const cs_real_t rho = 1.0;               /* Constant density [kg/m^3] (example) */

   /* Computation of mass flow rate and surface */
   const int kimasf = cs_field_key_id("inner_mass_flux_id");
   const cs_real_t *i_massflux = cs_field_by_id(cs_field_get_key_int(f, kimasf))->val;

   cs_real_t mflow = 0;
   cs_real_t tot_sur = 0;

   for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
       const cs_real_t x_abs = i_fac_cog[face_id][0];
       if (x_abs < 0.0001) {
           const cs_real_t x_s = i_face_surf[face_id];
           const cs_real_t x_n = i_face_normal[face_id][0];
           const cs_real_t x_mf = i_massflux[face_id];

           tot_sur += x_s;
           mflow   += x_mf * x_n / x_s;
       }
   }

   cs_real_t bbuf[2] = {tot_sur, mflow};
   cs_parall_sum(2, CS_REAL_TYPE, bbuf);
   tot_sur = bbuf[0];
   mflow = bbuf[1];

   /* Open log file */
   FILE *log_file = fopen("source_term_log.txt", "a");
   if (log_file == NULL) {
       printf("Error: Could not open log file for writing.\n");
       return;
   }

   fprintf(log_file, "Computed mflow: %e\n", mflow);

   if (mflow == 0.0) {
       printf("Warning: Mass flow rate is zero, setting a_T to zero.\n");
   }

   /* Compute the thermal fluxes at selected boundary faces */
   const int location_id = CS_MESH_LOCATION_BOUNDARY_FACES;
   const cs_lnum_t n_elts = cs_mesh_location_get_n_elts(location_id)[0];
   cs_real_t *boundary_flux = NULL;
   BFT_MALLOC(boundary_flux, n_elts, cs_real_t);

/*   const cs_field_t *fl = CS_F_(t); // Use temperature field directly*/
   const cs_field_t *fl = cs_thermal_model_field();
   cs_post_boundary_flux(fl->name, n_elts, NULL, boundary_flux);

   const cs_real_t *b_face_surf = cs_glob_mesh_quantities->b_face_surf;
   const cs_lnum_t (*b_face_cog)[3] = cs_glob_mesh_quantities->b_face_cog; // Get boundary zone

   /* Identify the bottom boundary zone ID */
   /*const int bottom_zone_id = cs_boundary_zone_by_name("GROUND");  // Use known name if available
   if (bottom_zone_id == -1) {
       fprintf(stderr, "Error: Bottom boundary not found.\n");
       return;
   }

   /* Computation of phi * S */
   cs_real_t integral_flux = 0.0;
   cs_real_t total_surface_area = 0.0;
   for (cs_lnum_t j = 0; j < n_elts; j++) {
       if (b_face_cog[j][2] < 0.00001) {
           cs_real_t surface_area = b_face_surf[j];
           cs_real_t boundary_flux_cell = boundary_flux[j];
           integral_flux += boundary_flux_cell * surface_area;
           total_surface_area += surface_area;
       }
   }

   cs_real_t bbuf2[2] = {integral_flux,total_surface_area};
   cs_parall_sum(2, CS_REAL_TYPE, bbuf2);
   integral_flux = bbuf2[0];
   total_surface_area = bbuf2[1];

   BFT_FREE(boundary_flux);

   /* Compute a_T = φS / (mflow C_p L) */
   cs_real_t a_T = 0.0;
   if (mflow != 0.0) {
       a_T = integral_flux / (mflow * cp * L);
   }
   fprintf(log_file, "Computed a_T: %e\n", a_T);
   fprintf(log_file, "Computed integral_flux: %e\n", integral_flux);
   fprintf(log_file, "Computed total_surface_area: %e\n", total_surface_area);


   /* Loop through all cells to apply the source term */
   for (cs_lnum_t i = 0; i < n_cells; i++) {
       double u_x = cvar_vel[i][0];  // Corrected indexing

   /* Print values for debugging */
       if (isnan(a_T) || isnan(rho) || isnan(u_x) || isnan(cp)) {
           printf("Error: NaN detected in source term at cell %ld\n", i);
           continue;
       }

       /* Compute the volumetric source term: a_T * ρ * u_x * c_p */
       cs_real_t source_term = a_T * rho * u_x * cp;

       /* Explicit source term */
       st_exp[i] += source_term;

       /* Implicit source term is zero for this example */
       st_imp[i] = 0.0;
   }
   /* Close log file */
   fclose(log_file);
}

END_C_DECLS