diff --git a/src/densities.cc b/src/densities.cc index 56f7445..f4aee69 100644 --- a/src/densities.cc +++ b/src/densities.cc @@ -17,6 +17,8 @@ #include +#include "math/special.hh" + #include "densities.hh" #include "random.hh" #include "convolution_kernel.hh" @@ -24,20 +26,6 @@ //TODO: this should be a larger number by default, just to maintain consistency with old default #define DEF_RAN_CUBE_SIZE 32 -double Meyer_scaling_function( double k, double kmax ) -{ - constexpr double twopithirds{2.0*M_PI/3.0}; - constexpr double fourpithirds{4.0*M_PI/3.0}; - auto nu = []( double x ){ return x<0.0?0.0:(x<1.0?x:1.0); }; - - k = std::abs(k)/kmax * fourpithirds; - - if( k < twopithirds ) return 1.0; - else if( k< fourpithirds ){ - return std::cos( 0.5*M_PI * nu(3*k/(2*M_PI)-1.0) ); - } - return 0.0; -} template void fft_coarsen(m1 &v, m2 &V) @@ -131,6 +119,7 @@ void fft_interpolate(m1 &V, m2 &v, int margin, bool from_basegrid = false) size_t nxf = v.size(0), nyf = v.size(1), nzf = v.size(2), nzfp = nzf + 2; size_t mxf = v.margin(0), myf = v.margin(1), mzf = v.margin(2); + // adjust offsets to respect margins, all grids have 'margins' except basegrid (which is periodic) if (!from_basegrid) { oxf += mxf/2; @@ -217,22 +206,21 @@ void fft_interpolate(m1 &V, m2 &v, int margin, bool from_basegrid = false) double kz = (k <= (int)nzc / 2) ? (double)k : (double)(k - (int)nzc); double phase = -0.5 * M_PI * (kx / nxc + ky / nyc + kz / nzc); + // double phase = -0.5 * M_PI * (kx / nxc + ky / nyc + kz / nzc); std::complex val_phas(cos(phase), sin(phase)); std::complex val(RE(ccoarse[qc]), IM(ccoarse[qc])); val *= val_phas * 8.0; - if(i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){ - double blend_coarse_x = Meyer_scaling_function(kx, nxc / 2); - double blend_coarse_y = Meyer_scaling_function(ky, nyc / 2); - double blend_coarse_z = Meyer_scaling_function(kz, nzc / 2); - double blend_coarse = blend_coarse_x*blend_coarse_y*blend_coarse_z; - double blend_fine = 1.0-blend_coarse; + double blend_coarse_x = Meyer_scaling_function(kx, nxc / 4); + double blend_coarse_y = Meyer_scaling_function(ky, nyc / 4); + double blend_coarse_z = Meyer_scaling_function(kz, nzc / 4); + double blend_coarse = blend_coarse_x*blend_coarse_y*blend_coarse_z; + double blend_fine = 1.0-blend_coarse; - RE(cfine[qf]) = blend_fine * RE(cfine[qf]) + blend_coarse * val.real(); - IM(cfine[qf]) = blend_fine * IM(cfine[qf]) + blend_coarse * val.imag(); - } + RE(cfine[qf]) = blend_fine * RE(cfine[qf]) + blend_coarse * val.real(); + IM(cfine[qf]) = blend_fine * IM(cfine[qf]) + blend_coarse * val.imag(); } delete[] rcoarse; diff --git a/src/math/special.hh b/src/math/special.hh new file mode 100644 index 0000000..59b0044 --- /dev/null +++ b/src/math/special.hh @@ -0,0 +1,42 @@ +// This file is part of monofonIC (MUSIC2) +// A software package to generate ICs for cosmological simulations +// Copyright (C) 2020 by Oliver Hahn +// +// monofonIC 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. +// +// monofonIC 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 . +#pragma once + +#include + +inline double Meyer_scaling_function( double k, double kmax ) +{ + constexpr double twopithirds{2.0*M_PI/3.0}; + constexpr double fourpithirds{4.0*M_PI/3.0}; + auto nu = []( double x ){ return x<0.0?0.0:(x<1.0?x:1.0); }; + + // k = std::abs(k)/kmax * fourpithirds; + k = std::abs(k)/kmax * 2 * M_PI; + + if( k < twopithirds ) return 1.0; + else if( k< fourpithirds ){ + return std::cos( 0.5*M_PI * nu(3*k/(2*M_PI)-1.0) ); + } + // if( k < fourpithirds ) return 1.0; + return 0.0; +} + +inline double Shannon_scaling_function( double k, double kmax ) +{ + if( std::abs(k) < kmax ) return 1.0; + return 0.0; +} \ No newline at end of file diff --git a/src/plugins/random_music_wnoise_generator.cc b/src/plugins/random_music_wnoise_generator.cc index 6d7f405..2eed5f7 100644 --- a/src/plugins/random_music_wnoise_generator.cc +++ b/src/plugins/random_music_wnoise_generator.cc @@ -4,6 +4,8 @@ #include #include +#include "math/special.hh" + #include "random.hh" #include "random_music_wnoise_generator.hh" @@ -627,18 +629,33 @@ music_wnoise_generator::music_wnoise_generator(music_wnoise_generator &rc, val *= val_phas * sqrt8; - // if (x0_ == NULL || lx_ == NULL){ - if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2) - { - RE(cfine[qf]) = val.real(); - IM(cfine[qf]) = val.imag(); - } - else - { - RE(cfine[qf]) = val.real(); - IM(cfine[qf]) = 0.0; - } - // } + // // if (x0_ == NULL || lx_ == NULL){ + // if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2) + // { + // RE(cfine[qf]) = val.real(); + // IM(cfine[qf]) = val.imag(); + // } + // else + // { + // RE(cfine[qf]) = val.real(); + // IM(cfine[qf]) = 0.0; + // } + // // } + if(i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){ + double blend_coarse_x = Meyer_scaling_function(kx, nxc / 2); + double blend_coarse_y = Meyer_scaling_function(ky, nyc / 2); + double blend_coarse_z = Meyer_scaling_function(kz, nzc / 2); + + // double blend_coarse_x = Shannon_scaling_function(kx, nxc / 2); + // double blend_coarse_y = Shannon_scaling_function(ky, nyc / 2); + // double blend_coarse_z = Shannon_scaling_function(kz, nzc / 2); + + double blend_coarse = blend_coarse_x*blend_coarse_y*blend_coarse_z; + double blend_fine = std::sqrt(1.0-blend_coarse*blend_coarse); + + RE(cfine[qf]) = blend_fine * RE(cfine[qf]) + blend_coarse * val.real(); + IM(cfine[qf]) = blend_fine * IM(cfine[qf]) + blend_coarse * val.imag(); + } } delete[] rcoarse;