diff options
Diffstat (limited to 'td_constraints_template.c')
| -rw-r--r-- | td_constraints_template.c | 43 |
1 files changed, 32 insertions, 11 deletions
diff --git a/td_constraints_template.c b/td_constraints_template.c index df09290..5770617 100644 --- a/td_constraints_template.c +++ b/td_constraints_template.c @@ -16,8 +16,8 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. */ -#define L 1.0 -#define R_0 0.0 +#define L 1 +#define R_0 5.0 static void EQNAME(double a, unsigned int nb_coords[2], double *coords[2], @@ -53,9 +53,14 @@ static void EQNAME(double a, //const double K_rtheta_val = K_rtheta(a, r, theta); //const double dK_rtheta_r_val = K_rtheta_dr(a, r, theta); //const double dK_rtheta_theta_val = K_rtheta_dtheta(a, r, theta); - const double K_rtheta_val = eval_K_rtheta(a, L, R_0, coord0, coord1); - const double dK_rtheta_r_val = eval_dK_rtheta_r(a, L, R_0, coord0, coord1); - const double dK_rtheta_theta_val = eval_dK_rtheta_theta(a, L, R_0, coord0, coord1); + double K_rtheta_val = eval_K_rtheta(a, L, R_0, coord0, coord1); + double dK_rtheta_r_val = eval_dK_rtheta_r(a, L, R_0, coord0, coord1); + double dK_rtheta_theta_val = eval_dK_rtheta_theta(a, L, R_0, coord0, coord1); + if (fabs(r - R_0) > L) { + K_rtheta_val = 0.0; + dK_rtheta_r_val = 0.0; + dK_rtheta_theta_val = 0.0; + } const double Km[3][3] = {{ K_rr, K_rtheta_val, 0 }, { K_rtheta_val / SQR(r), -(K_rr + K_phiphi), 0 }, @@ -157,8 +162,14 @@ static void EQNAME(double a, #endif if (EQUATION == TD_CONSTRAINT_EQ_HAM) { - double Rscal = eval_R_scalar(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, - K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10); + double Rscal; + if (fabs(r - R_0) < L) { + Rscal = eval_R_scalar(a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, + K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10); + } else { + Rscal = eval_R_scalar1(a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, + K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10); + } double K2 = 0.0; for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) @@ -166,11 +177,21 @@ static void EQNAME(double a, out[idx] = Rscal - K2; } else if (EQUATION == TD_CONSTRAINT_EQ_MOM_0) { - out[idx] = eval_Dcurv_m_r(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, - K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10); + if (fabs(r - R_0) < L) { + out[idx] = eval_Dcurv_m_r(a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, + K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10); + } else { + out[idx] = eval_Dcurv_m_r1(a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, + K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10); + } } else { - out[idx] = eval_Dcurv_m_t(a, L, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, - K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10); + if (fabs(r - R_0) < L) { + out[idx] = eval_Dcurv_m_t(a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, + K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10); + } else { + out[idx] = eval_Dcurv_m_t1(a, L, R_0, coord0, coord1, psi, dpsi_01, dpsi_02, dpsi_10, dpsi_11, dpsi_20, + K_rr, dK_rr_01, dK_rr_10, K_phiphi, dK_phiphi_01, dK_phiphi_10); + } } #if 0 } else { |