diff options
Diffstat (limited to 'src/patch/patch.cc')
| -rw-r--r-- | src/patch/patch.cc | 57 |
1 files changed, 54 insertions, 3 deletions
diff --git a/src/patch/patch.cc b/src/patch/patch.cc index f4305f6..6b355fd 100644 --- a/src/patch/patch.cc +++ b/src/patch/patch.cc @@ -43,6 +43,7 @@ #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; #include "coords.hh" #include "grid.hh" @@ -55,7 +56,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** //****************************************************************************** @@ -300,7 +300,7 @@ else error_exit(PANIC_EXIT, " which is a multiple of 90 degrees!\n" " %s patch: [min,max]_dsigma()=[%g,%g]\n" , - name(), min_dsigma(), max_dsigma()); + name(), double(min_dsigma()), double(max_dsigma())); const fp sigma = sigma_of_dsigma(dsigma); const int isigma = isigma_of_sigma(sigma); @@ -371,7 +371,7 @@ else error_exit(PANIC_EXIT, " which is a multiple of 90 degrees!\n" " %s patch: [min,max]_drho()=[%g,%g]\n" , - name(), min_drho(), max_drho()); + name(), double(min_drho()), double(max_drho())); const fp rho = rho_of_drho(drho); const int irho = irho_of_rho(rho); @@ -611,6 +611,57 @@ fp sum = 0.0; return delta_rho() * delta_sigma() * sum; } +fp patch::integrate_gridpoint(int unknown_src_gfn, + int ghosted_radius_gfn, + int g_xx_gfn, int g_xy_gfn, int g_xz_gfn, + int g_yy_gfn, int g_yz_gfn, + int g_zz_gfn, + enum integration_method method, + int irho, int isigma) + const +{ +const bool src_is_ghosted = is_valid_ghosted_gfn(unknown_src_gfn); + +const fp fn = unknown_gridfn(src_is_ghosted, + unknown_src_gfn, irho,isigma); + +const fp rho = rho_of_irho (irho); +const fp sigma = sigma_of_isigma(isigma); +const fp r = ghosted_gridfn(ghosted_radius_gfn, irho,isigma); +const fp partial_surface_r_wrt_rho + = partial_rho (ghosted_radius_gfn, irho,isigma); +const fp partial_surface_r_wrt_sigma + = partial_sigma(ghosted_radius_gfn, irho,isigma); + +const fp g_xx = gridfn(g_xx_gfn, irho,isigma); +const fp g_xy = gridfn(g_xy_gfn, irho,isigma); +const fp g_xz = gridfn(g_xz_gfn, irho,isigma); +const fp g_yy = gridfn(g_yy_gfn, irho,isigma); +const fp g_yz = gridfn(g_yz_gfn, irho,isigma); +const fp g_zz = gridfn(g_zz_gfn, irho,isigma); + +fp g_rho_rho, g_rho_sigma, g_sigma_sigma; +const fp Jac = rho_sigma_metric(r, rho, sigma, + partial_surface_r_wrt_rho, + partial_surface_r_wrt_sigma, + g_xx, g_xy, g_xz, + g_yy, g_yz, + g_zz, + g_rho_rho, g_rho_sigma, + g_sigma_sigma); + +const fp coeff_rho = integration_coeff(method, + max_irho()-min_irho(), + irho -min_irho()); +const fp coeff_sigma = integration_coeff(method, + max_isigma()-min_isigma(), + isigma -min_isigma()); + +const fp val = coeff_rho*coeff_sigma * fn * sqrt(jtutil::abs(Jac)); + +return delta_rho() * delta_sigma() * val; +} + //****************************************************************************** // |