/* * Copyright 2017 Anton Khirnov * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ #include #include "common.h" #include "pssolve.h" #include "td_constraints.h" #define I_C sqrt(225.0 / (8.0 * M_PI)) static double I_m2(double x, double a) { return a * I_C * exp(-SQR(x)); } static double I_m1(double x, double a) { return a * I_C * (-2.0 * x * exp(-SQR(x))); } static double I_0(double x, double a) { return a * I_C * (-2.0 + 4.0 * SQR(x)) * exp(-SQR(x)); } static double I_p1(double x, double a) { return a * I_C * (12.0 * x - 8.0 * x * SQR(x)) * exp(-SQR(x)); } static double I_p2(double x, double a) { return a * I_C * (12.0 - 48.0 * SQR(x) + 16.0 * SQR(SQR(x))) * exp(-SQR(x)); } static double I_p3(double x, double a) { return a * I_C * (-120.0 * x + 144.0 * x * SQR(x) - 32.0 * x * SQR(SQR(x))) * exp(-SQR(x)); } static double K_rtheta(double a, double r, double theta) { double r2 = r * r; double r3 = r2 * r; // equal to the below return a * I_C * exp(-r2) * sin(2 * theta) * (32.0 * r3 - 48.0 * r); //return ((I_p2(-r, a) + I_p2(r, a)) * r3 + 3.0 * (I_p1(-r, a) - I_p1(r, a)) * r2 + // 6.0 * (I_0(-r, a) + I_0(r, a)) * r + 6.0 * (I_m1(-r, a) - I_m1(r, a))) * sin(2.0 * theta) / r4; } static double K_rtheta_dtheta(double a, double r, double theta) { return 2.0 * cos(2.0 * theta) * K_rtheta(a, r, theta) / sin(2.0 * theta); } static double K_rtheta_dr(double a, double r, double theta) { double r2 = r * r; return -2 * r * K_rtheta(a, r, theta) + a * I_C * exp(-r2) * sin(2 * theta) * (3.0 * 32.0 * r2 - 48.0); //return ((-I_p3(-r, a) + I_p3(r, a)) / r + 3.0 * (-I_p2(-r, a) - I_p2(r, a)) / (SQR(r)) + // 6.0 * (-I_p1(-r, a) + I_p1(r, a)) / (r * SQR(r)) + 6.0 * (-I_0(-r, a) - I_0(r, a)) / (SQR(SQR(r))) + // (I_p2(-r, a) + I_p2(r, a)) / (-SQR(r)) + 3.0 * (I_p1(-r, a) - I_p1(r, a)) * (-2.0 / (r * SQR(r))) + // 6.0 * (I_0(-r, a) + I_0(r, a)) * (-3.0 / (SQR(SQR(r)))) + 6.0 * (I_m1(-r, a) - I_m1(r, a)) * (-4.0 / (r * SQR(SQR(r))))) * sin(2.0 * theta); } #include "eval_constraints.c" #include "eval_k_rtheta.c" //typedef void (*eval_metric_func)(double a, double L, double r_0, double x_0, double x_1, double psi_val_00, double psi_val_01, double psi_val_02, double psi_val_10, double psi_val_11, double psi_val_20, double *res); // //static const eval_metric_func eval_dmetric[3] = { // eval_dmetric_0, // eval_dmetric_1, // eval_dmetric_2, //}; // //static const eval_metric_func eval_d2metric[3][3] = { // { // eval_d2metric_00, // eval_d2metric_01, // eval_d2metric_02, // }, // { // eval_d2metric_10, // eval_d2metric_11, // eval_d2metric_12, // }, // { // eval_d2metric_20, // eval_d2metric_21, // eval_d2metric_22, // }, //}; #define EQNAME constraint_eq_ham #define EQUATION TD_CONSTRAINT_EQ_HAM #include "td_constraints_template.c" #undef EQNAME #undef EQUATION #define EQNAME constraint_eq_mom_0 #define EQUATION TD_CONSTRAINT_EQ_MOM_0 #include "td_constraints_template.c" #undef EQNAME #undef EQUATION #define EQNAME constraint_eq_mom_1 #define EQUATION TD_CONSTRAINT_EQ_MOM_1 #include "td_constraints_template.c" #undef EQNAME #undef EQUATION void (*constraint_funcs[TD_CONSTRAINT_EQ_NB])(const double a, unsigned int nb_coords[2], double *coords[2], double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], double *out) = { [TD_CONSTRAINT_EQ_HAM] = constraint_eq_ham, [TD_CONSTRAINT_EQ_MOM_0] = constraint_eq_mom_0, [TD_CONSTRAINT_EQ_MOM_1] = constraint_eq_mom_1, }; static double q_gundlach(const double a, double rho, double z) { return a * SQR(rho) * exp(-(SQR(rho) + SQR(z))); } static double dq_rho_gundlach(const double a, double rho, double z) { return a * 2 * rho * exp(-SQR(rho) - SQR(z)) * (1.0 - rho * (rho)); } static double d2q_rho_gundlach(const double a, double rho, double z) { double rho2 = SQR(rho); return a * 2 * exp(-SQR(rho) - SQR(z)) * (1.0 - 4.0 * rho * (rho) - rho2 + 2.0 * rho2 * SQR(rho)); } static double d2q_rho_z_gundlach(const double a, double rho, double z) { return a * 4 * z * rho * exp(-SQR(rho) - SQR(z)) * (rho * (rho) - 1.0); } static double dq_z_gundlach(const double a, double rho, double z) { return -a * 2 * z * SQR(rho) * exp(-SQR(rho) - SQR(z)); } static double d2q_z_gundlach(const double a, double rho, double z) { return a * 2 * SQR(rho) * exp(-SQR(rho) - SQR(z)) * (2 * SQR(z) - 1.0); } int tdi_constraints_eval(enum TDConstraintEq eq, const double a, unsigned int nb_coords[2], double *coords[2], double *vars[TD_CONSTRAINT_VAR_NB][PSSOLVE_DIFF_ORDER_NB], double *out) { //for (int idx1 = 0; idx1 < nb_coords[1]; idx1++) { // const double theta = coords[1][idx1]; // for (int idx0 = 0; idx0 < nb_coords[0]; idx0++) { // const double r = coords[0][idx0]; // const double rho = r * sin(theta); // const double z = r * cos(theta); // const double ct = cos(theta); // const double st = sin(theta); // const int idx = idx1 * nb_coords[0] + idx0; // const double laplace = vars[eq][PSSOLVE_DIFF_ORDER_20][idx] + // vars[eq][PSSOLVE_DIFF_ORDER_10][idx] * (2.0 / r) + // vars[eq][PSSOLVE_DIFF_ORDER_02][idx] * (1.0 / SQR(r)) + // vars[eq][PSSOLVE_DIFF_ORDER_01][idx] * (ct / (SQR(r) * st)); // out[idx] = laplace + 0.25 * (d2q_z_gundlach(a, rho, z) + d2q_rho_gundlach(a, rho, z)) * (1.0 + vars[eq][PSSOLVE_DIFF_ORDER_00][idx]); // } //} constraint_funcs[eq](a, nb_coords, coords, vars, out); return 0; }