diff --git a/src/equation_of_state/planetary/sesame.h b/src/equation_of_state/planetary/sesame.h index 11c16964602b28c0d1a080b6c262ff20c1f5b9cb..037de315aa7e88303f1700746c83daddb16e2949 100644 --- a/src/equation_of_state/planetary/sesame.h +++ b/src/equation_of_state/planetary/sesame.h @@ -338,15 +338,35 @@ INLINE static float SESAME_pressure_from_internal_energy( idx_u_2 = mat->num_T - 2; } - intp_rho = (log_rho - mat->table_log_rho[idx_rho]) / - (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]); - intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) / - (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] - - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]); - intp_u_2 = - (log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) / - (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] - - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]); + // Check for duplicates in SESAME tables before interpolation + if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) { + intp_rho = (log_rho - mat->table_log_rho[idx_rho]) / + (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]); + } else { + intp_rho = 1.; + } + + // Check for duplicates in SESAME tables before interpolation + if (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] != + mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) { + intp_u_1 = + (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) / + (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] - + mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]); + } else { + intp_u_1 = 1.; + } + + // Check for duplicates in SESAME tables before interpolation + if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] != + mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) { + intp_u_2 = + (log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) / + (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] - + mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]); + } else { + intp_u_2 = 1.; + } // Table values P_1 = mat->table_P_rho_T[idx_rho * mat->num_T + idx_u_1]; @@ -442,15 +462,35 @@ INLINE static float SESAME_soundspeed_from_internal_energy( idx_u_2 = mat->num_T - 2; } - intp_rho = (log_rho - mat->table_log_rho[idx_rho]) / - (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]); - intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) / - (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] - - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]); - intp_u_2 = - (log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) / - (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] - - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]); + // Check for duplicates in SESAME tables before interpolation + if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) { + intp_rho = (log_rho - mat->table_log_rho[idx_rho]) / + (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]); + } else { + intp_rho = 1.; + } + + // Check for duplicates in SESAME tables before interpolation + if (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] != + mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) { + intp_u_1 = + (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) / + (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] - + mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]); + } else { + intp_u_1 = 1.; + } + + // Check for duplicates in SESAME tables before interpolation + if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] != + mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) { + intp_u_2 = + (log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) / + (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] - + mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]); + } else { + intp_u_2 = 1.; + } // Table values c_1 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1];