text:
#pragma once

#include "BoundingRegion.h"
#include "Ellipsoid.h"
#include "S2CellID.h"

#include <CesiumGeometry/CullingResult.h>
#include <CesiumGeometry/Plane.h>

#include <glm/vec3.hpp>
#include <gsl/span>

#include <array>
#include <string_view>

namespace CesiumGeospatial {

/**
 * A tile bounding volume specified as an S2 cell token with minimum and maximum
 * heights. The bounding volume is a k DOP. A k-DOP is the Boolean intersection
 * of extents along k directions.
 *
 * @param cellID The S2 cell ID.
 * @param minimumHeight The minimum height of the bounding volume.
 * @param maximumHeight The maximum height of the bounding volume.
 * @param ellipsoid The ellipsoid.
 */
class CESIUMGEOSPATIAL_API S2CellBoundingVolume final {
public:
  S2CellBoundingVolume(
      const S2CellID& cellID,
      double minimumHeight,
      double maximumHeight,
      const Ellipsoid& ellipsoid = Ellipsoid::WGS84);

  /**
   * @brief Gets this bounding volume's cell ID.
   */
  const S2CellID& getCellID() const { return this->_cellID; }

  /**
   * @brief Gets the minimum height of the cell.
   */
  double getMinimumHeight() const noexcept { return this->_minimumHeight; }

  /**
   * @brief Gets the maximum height of the cell.
   */
  double getMaximumHeight() const noexcept { return this->_maximumHeight; }

  /**
   * @brief Gets the center of this bounding volume in ellipsoid-fixed (ECEF)
   * coordinates.
   */
  glm::dvec3 getCenter() const noexcept;

  /**
   * @brief Gets the either corners of the bounding volume, in ellipsoid-fixed
   * (ECEF) coordinates.
   *
   * @return An array of positions with a `size()` of 8.
   */
  gsl::span<const glm::dvec3> getVertices() const noexcept;

  /**
   * @brief Determines on which side of a plane the bounding volume is located.
   *
   * @param plane The plane to test against.
   * @return The {@link CesiumGeometry::CullingResult}
   *  * `Inside` if the entire region is on the side of the plane the normal is
   * pointing.
   *  * `Outside` if the entire region is on the opposite side.
   *  * `Intersecting` if the region intersects the plane.
   */
  CesiumGeometry::CullingResult
  intersectPlane(const CesiumGeometry::Plane& plane) const noexcept;

  /**
   * @brief Computes the distance squared from a given position to the closest
   * point on this bounding volume. The position must be expressed in
   * ellipsoid-centered (ECEF) coordinates.
   *
   * @param position The position
   * @return The estimated distance squared from the bounding box to the point.
   *
   * @snippet TestOrientedBoundingBox.cpp distanceSquaredTo
   */
  double
  computeDistanceSquaredToPosition(const glm::dvec3& position) const noexcept;

  /**
   * @brief Gets the six planes that bound the volume.
   *
   * @return An array of planes with a `size()` of 6.
   */
  gsl::span<const CesiumGeometry::Plane> getBoundingPlanes() const noexcept;

  /**
   * @brief Computes the bounding begion that best fits this S2 cell volume.
   *
   * @return The bounding region.
   */
  BoundingRegion computeBoundingRegion() const noexcept;

private:
  S2CellID _cellID;
  double _minimumHeight;
  double _maximumHeight;
  glm::dvec3 _center;
  std::array<CesiumGeometry::Plane, 6> _boundingPlanes;
  std::array<glm::dvec3, 8> _vertices;
};

} // namespace CesiumGeospatial

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native


text:
#pragma once

#include <d3d11.h>
#include <gsl\gsl>
#include "Rectangle.h"

namespace Library
{
	class TextureHelper final
	{
	public:
		static Point GetTextureSize(gsl::not_null<ID3D11Texture2D*> texture);
		static Rectangle GetTextureBounds(gsl::not_null<ID3D11Texture2D*> texture);
		static std::uint32_t BitsPerPixel(const DXGI_FORMAT format);

		TextureHelper() = delete;
		TextureHelper(const TextureHelper&) = delete;
		TextureHelper& operator=(const TextureHelper&) = delete;
		TextureHelper(TextureHelper&&) = delete;
		TextureHelper& operator=(TextureHelper&&) = delete;
		~TextureHelper() = default;
	};
}

ext: h


text:
/* specfunc/test_legendre.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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 2 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, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sf.h>
#include "test_sf.h"

int test_legendre(void)
{
  gsl_sf_result r;
  double L[256];
  int s = 0;
  int sa;

  TEST_SF(s,  gsl_sf_legendre_P1_e, (-0.5, &r), -0.5, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_legendre_P1_e, ( 0.5, &r), 0.5, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_legendre_P2_e, (0.0, &r), -0.5  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_legendre_P2_e, (0.5, &r), -0.125, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_legendre_P2_e, (1.0, &r), 1.0  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_legendre_P2_e, (100.0, &r), 14999.5  , TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_legendre_P3_e, ( -0.5, &r), 0.4375, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_legendre_P3_e, (  0.5, &r), -0.4375, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_legendre_P3_e, (  1.0, &r), 1.0        , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_legendre_P3_e, (100.0, &r), 2.49985e+06, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Pl_e, (1, -0.5, &r), -0.5, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (1,  1.0e-8, &r), 1.0e-08, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (1,  0.5, &r), 0.5, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (1,  1.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Pl_e, (10, -0.5, &r), -0.1882286071777345, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (10,  1.0e-8, &r), -0.24609374999999864648, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (10,  0.5, &r), -0.18822860717773437500, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (10,  1.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Pl_e, (99, -0.5, &r), 0.08300778172138770477, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (99,  1.0e-8, &r), -7.958923738716563193e-08, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (99,  0.5, &r), -0.08300778172138770477, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (99,  0.999, &r), -0.3317727359254778874, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (99,  1.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Pl_e, (1000, -0.5, &r),   -0.019168251091650277878, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (1000,  1.0e-8, &r), 0.02522501817709828,     TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (1000,  0.5, &r),   -0.019168251091650277878, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (1000,  1.0, &r),    1.0,                     TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Pl_e, (4000, -0.5, &r), -0.009585404456573080972, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (4000,  0.5, &r), -0.009585404456573080972, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Pl_e, (4000,  1.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);

  sa = 0;
  gsl_sf_legendre_Pl_array(100, 0.5, L);
  sa += ( test_sf_frac_diff(L[0],    1.0 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(L[10],  -0.18822860717773437500 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(L[100], -0.06051802596186118687 ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_legendre_Pl_array(100)");
  s += sa;

  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 0, -0.5, &r), -0.18822860717773437500, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 0, 1.0e-08, &r), -0.24609374999999864648, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 0, 0.5, &r), -0.18822860717773437500, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 1, -0.5, &r), -2.0066877394361256516, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 1, 1.0e-08, &r), -2.7070312499999951725e-07, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 1, 0.5, &r), 2.0066877394361256516, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 5, -0.5, &r),    -30086.169706116174977,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 5, 1.0e-08, &r), -0.0025337812499999964949, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 5, 0.5, &r),      30086.169706116174977,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (10, 5, 0.999, &r),   -0.5036411489013270406,    TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Plm_e, (100, 5, -0.5, &r), -6.617107444248382171e+08, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (100, 5, 1.0e-08, &r), 817.8987598063712851, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (100, 5, 0.5, &r), 6.617107444248382171e+08, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Plm_e, (100, 5, 0.999, &r), -1.9831610803806212189e+09, TEST_TOL2, GSL_SUCCESS);

  sa = 0;
  gsl_sf_legendre_Plm_array(100, 5, 0.5, L);
  sa += ( test_sf_frac_diff(L[0],  -460.3466286991656682 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(L[10],  38852.51334152290535 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(L[95],  6.617107444248382171e+08 ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_legendre_Plm_array(100, 5, 0.5)");
  s += sa;

  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 0, -0.5, &r), -0.24332702369300133776, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 0, 0.5, &r), -0.24332702369300133776, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 0, 0.999, &r), 1.2225754122797385990, TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 5, -0.5, &r),    -0.3725739049803293972,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 5, 1.0e-08, &r), -3.1377233589376792243e-08, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 5, 0.5, &r),      0.3725739049803293972,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 5, 0.999, &r),   -6.236870674727370094e-06,  TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 10, -0.5, &r), 0.12876871185785724117, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 10, 0.5, &r), 0.12876871185785724117,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (10, 10, 0.999, &r), 1.7320802307583118647e-14, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (200, 1, -0.5, &r),   0.3302975570099492931, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (200, 1, 0.5, &r),   -0.3302975570099492931, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_sphPlm_e, (200, 1, 0.999, &r), -1.4069792055546256912, TEST_TOL2, GSL_SUCCESS);

  sa = 0;
  gsl_sf_legendre_sphPlm_array(100, 5, 0.5, L);
  sa += ( test_sf_frac_diff(L[0],   -0.22609703187800460722 ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(L[10],   0.07452710323813558940 ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(L[95],   0.25865355990880161717 ) > TEST_TOL2 );
  gsl_test(s, "  gsl_sf_legendre_sphPlm_array(100, 5, 0.5)");
  s += sa;

  TEST_SF(s, gsl_sf_conicalP_half_e, (0.0, -0.5, &r),   0.8573827581049917129, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (0.0,  0.5, &r),   0.8573827581049917129, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (0.0,  2.0, &r),   0.6062611623284649811, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (0.0,  100.0, &r), 0.07979045091636735635, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_half_e, (10.0, -0.5, &r),    5.345484922591867188e+08, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (10.0,  0.5, &r),    15137.910380385258370, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (10.0,  2.0, &r),    0.4992680691891618544, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (10.0,  100.0, &r), -0.07272008163718195685, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_half_e, (200.0, -1.0e-3, &r),  1.3347639529084185010e+136, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (200.0,  1.0e-8, &r),  1.0928098010940058507e+136, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (200.0,  0.5, &r),     3.895546021611205442e+90,   TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (200.0,  10.0, &r),   -0.04308567180833581268,     TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (200.0,  100.0, &r),  -0.04694669186576399194,     TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (200.0,  1000.0, &r),  0.023698140704121273277,    TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (200.0,  1.0e+8, &r), -0.00006790983312124277891,  TEST_TOL3, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_half_e, (1.0e+8,  1.1, &r),   1.1599311133054742944,  TEST_SQRT_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_half_e, (1.0e+8,  100.0, &r), 0.07971967557381557875, TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (0.0, -0.5, &r),  1.7956982494514644808, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (0.0,  0.5, &r),  0.8978491247257322404, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (0.0,  2.0, &r),  0.7984204253272901551, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (0.0,  100.0, &r),  0.4227531369388072584, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (10.0, -0.5, &r),  5.345484922591867181e+07, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (10.0,  0.5, &r),  1513.7910356104985334, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (10.0,  2.0, &r),  0.03439243987215615642, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (10.0,  100.0, &r),  0.003283756665952609624, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (200.0, -0.5, &r),  1.7699538115312304280e+179, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (200.0,  1.0e-8, &r),  5.464049005470029253e+133, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (200.0,  0.5, &r),  1.9477730108056027211e+88, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (200.0,  10.0, &r),  0.0012462575917716355362, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (200.0,  100.0, &r),  -0.0003225881344802625149, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (200.0,  1000.0, &r), -0.00004330652890886567623, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (200.0,  1.0e+8, &r),  2.0943091278037078483e-07, TEST_TOL3, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (1.0e+8,  1.1, &r), 2.092320445620989618e-09, 16.0*TEST_SQRT_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_mhalf_e, (1.0e+8,  100.0, &r),  -3.359967833599016923e-11, 256.0*TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_0_e, (0.0, -0.5, &r),  1.3728805006183501647, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_0_e, (0.0,  0.5, &r),  1.0731820071493643751, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_0_e, (0.0,  2.0, &r),  0.9012862993604472987, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_0_e, (0.0,  100.0, &r),  0.30091748588199264556, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_0_e, (10.0, -0.5, &r),  1.6795592815421804669e+08, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_0_e, (10.0,  0.5, &r),  4826.034132009618240,      TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_0_e, (10.0,  2.0, &r),  0.18798468917758716146,    TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_0_e, (10.0,  100.0, &r), -0.008622130749987962529, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_0_e, (200.0,  -0.5, &r), 2.502194818646823e+180, TEST_TOL4, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_0_e, (1000.0,  100.0, &r),   0.0017908817653497715844, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_0_e, (1000.0,  1000.0, &r), -0.0006566893804926284301, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_0_e, (1000.0,  1.0e+8, &r),  2.3167213561756390068e-06, TEST_TOL4, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_1_e, (0.0, -0.5, &r),    0.4939371126656998499,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (0.0,  0.5, &r),    0.14933621085538265636, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (0.0,  2.0, &r),   -0.13666874968871549533, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (0.0,  100.0, &r), -0.10544528203156629098, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_1_e, (10.0, -0.5, &r),    1.7253802958788312520e+09, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (10.0,  0.5, &r),    46781.02294059967988,      TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (10.0,  2.0, &r),    0.26613342643657444400,    TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (10.0,  100.0, &r), -0.23281959695501029796,    TEST_TOL2, GSL_SUCCESS);

  /* FIXME: Mathematica gets some brain-damaged numbers for
   * these x < 0 points. I have checked what I am doing in detail,
   * and it must be right because you can do it by summing
   * manifestly positive definite quantities.
   */
  TEST_SF(s, gsl_sf_conicalP_1_e, (200.0, -0.999, &r), 2.71635193199341135e+270, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (200.0, -0.9, &r),   4.2952493176812905e+234,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (200.0, -0.5, &r),   5.01159205956053439e+182, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (200.0,  0.999, &r), 195733.0396081538,        TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (200.0,  10.0, &r), -2.9272610662414349553,    TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_1_e, (1000.0, 100.0, &r),  -1.7783258105862399857,    TEST_TOL6, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (1000.0, 1000.0, &r),  0.4535161075156427179,    TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_1_e, (1000.0, 1.0e+8, &r),  0.0009983414549874888478, TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (2,  1.0, -0.5, &r),  1.6406279287008789526,      TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (10, 1.0, -0.5, &r),  0.000029315266725049129448, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (20, 1.0, -0.5, &r),  7.335769429462034431e-15,   TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (30, 1.0, -0.5, &r),  1.3235612394267378871e-26,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (10, 1.0, 0.5, &r),  2.7016087199857873954e-10, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (20, 1.0, 0.5, &r),  1.1782569701435933399e-24, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (30, 1.0, 0.5, &r),  3.636240588303797919e-41,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (10, 1.0, 2.0, &r),  2.4934929626284934483e-10, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (20, 1.0, 2.0, &r),  1.1284762488012616191e-24, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_sph_reg_e, (30, 100.0, 100.0, &r),  -1.6757772087159526048e-64, TEST_TOL6, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (2, 1.0, -0.5, &r),   2.2048510472375258708,       TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (10, 1.0, -0.5, &r),  0.00007335034531618655690,   TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (20, 1.0, -0.5, &r),  2.5419860619212164696e-14,   TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (30, 1.0, -0.5, &r),  5.579714972260536827e-26,    TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (10, 1.0, 0.5, &r),  1.1674078819646475282e-09,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (20, 1.0, 0.5, &r),  7.066408031229072207e-24,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (30, 1.0, 0.5, &r),  2.6541973286862588488e-40,    TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (10, 1.0, 2.0, &r),  1.0736109751890863051e-09,    TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (20, 1.0, 2.0, &r),  6.760965304863386741e-24,     TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_conicalP_cyl_reg_e, (30, 100.0, 100.0, &r), -4.268753482520651007e-63, TEST_TOL4, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (1.0e-06, 1.0e-06, &r), 0.9999999999998333333    , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (1.0, 0.0, &r), 1.0                      , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (1.0, 1.0, &r), 0.7160229153604338713    , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (1.0, 100.0, &r), -3.767437313149604566e-44 , TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (1.0, 500.0, &r), -6.665351935878582205e-218, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (100.0, 1.0, &r), -0.004308757035378200029  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (100.0, 10.0, &r), 7.508054627912986427e-07 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (1000.0, 1.0, &r), 0.0007036067909088818319 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (1.0e+08, 1.0, &r), 7.927485371429105968e-09 , TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_0_e, (1.0e+08, 100.0, &r), -3.627118904186918957e-52 , 32.0*TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1.0e-06, 1.0e-06, &r), 3.333333333334222222e-07, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1.0, 1.0e-10, &r), 4.714045207910316829e-11,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1.0, 1.0, &r), 0.3397013994799344639,           TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1.0, 100.0, &r), -7.200624449531811272e-44,     TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1.0, 500.0, &r), 4.192260336821728677e-218,     TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (100.0, 0.01, &r), 0.30117664944267412324   , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (100.0, 1.0, &r), -0.007393833425336299309  , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (100.0, 10.0, &r), -5.031062029821254982e-07 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1000.0, 0.001, &r), 0.30116875865090396421   , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1000.0, 1.0, &r), -0.0004776144516074971885 , TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1.0e+08, 1.0e-08, &r), 0.30116867893975679722   , TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1.0e+08, 1.0, &r), 3.0921097047369081582e-09, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_1_e, (1.0e+08, 100.0, &r), -6.496142701296286936e-52 , 32.0*TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0e-06, 1.0e-06, &r),  1.1544011544013627977e-32, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0, 1.0e-10, &r),      2.0224912016958766992e-52, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0, 1.0, &r),          0.011498635037491577728,   TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0, 5.0, &r),          0.0020696945662545205776,  TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0, 7.0, &r),     -0.0017555303787488993676,   TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0, 10.0, &r),     0.00008999979724504887101,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0, 100.0, &r),   -4.185397793298567945e-44,   TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0, 500.0, &r),    1.4235113901091961263e-217, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 100.0, 0.001, &r),  9.642762597222417946e-10,   TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 100.0, 0.002, &r),  3.0821201254308036109e-08,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 100.0, 0.01, &r),   0.00009281069019005840532,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 100.0, 1.0, &r),   -0.008043100696178624653,    TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 100.0, 10.0, &r),  -3.927678432813974207e-07,   TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1000.0, 0.001, &r),  0.00009256365284253254503, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1000.0, 0.01, &r),  -0.05553733815473079983, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0e+08, 1.0e-08, &r),   0.00009256115861125841299, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_H3d_e, (5, 1.0e+08, 100.0, &r),    -6.496143209092860765e-52 , 128.0*TEST_SQRT_TOL0, GSL_SUCCESS);

#if 0
  sa = 0;
  gsl_sf_legendre_H3d_array(100, 1.0, 3.0, L);
  sa += ( test_sf_frac_diff(L[  0], gsl_sf_legendre_H3d(  0, 1.0, 3.0)) > 1.0e-12 );
  sa += ( test_sf_frac_diff(L[  1], gsl_sf_legendre_H3d(  1, 1.0, 3.0)) > 1.0e-12 );
  sa += ( test_sf_frac_diff(L[ 10], gsl_sf_legendre_H3d( 10, 1.0, 3.0)) > 1.0e-12 );
  sa += ( test_sf_frac_diff(L[100], gsl_sf_legendre_H3d(100, 1.0, 3.0)) > 1.0e-12 );
  gsl_test(sa, "  gsl_sf_legendre_H3d_array(100, 1.0, 3.0)");
  s += sa;
#endif

  /* x = -1 + 2^-16 */
  TEST_SF(s, gsl_sf_legendre_Q0_e, (-0.9999847412109375, &r), -5.8917472200477175158028143531855, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, (-0.5, &r), -0.5493061443340548457, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, (-1e-10, &r), -1.000000000000000000e-10, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, (0.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, (1e-10, &r), 1.000000000000000000e-10, TEST_TOL0, GSL_SUCCESS);
  /* x = 1 - 2^-16 */
  TEST_SF(s, gsl_sf_legendre_Q0_e, (0.9999847412109375, &r), 5.8917472200477175158028143531855, TEST_TOL4, GSL_SUCCESS);
  /* x = 1 + 2^-16 */
  TEST_SF(s, gsl_sf_legendre_Q0_e, ( 1.0000152587890625, &r), 5.8917548494422489138325509750429, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, ( 1.5, &r), 0.8047189562170501873, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, ( 9.99, &r), 0.1004364599660005447, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, ( 10.0, &r), 0.1003353477310755806, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, ( 10.01, &r), 0.1002344395571710243, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, ( 100, &r), 0.010000333353334762015, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q0_e, ( 1e10, &r), 1.000000000000000000e-10, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Q1_e, (-0.9999847412109375, &r), 4.8916573191196772369, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (-0.5, &r), -0.7253469278329725772, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (-0.01, &r), -0.9998999966664666524, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (-1e-10, &r), -0.999999999999999999, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (0.0, &r), -1.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (1e-10, &r), -0.999999999999999999, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (0.0001, &r), -0.9999999899999999667, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (0.01, &r), -0.9998999966664666524, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (0.5, &r), -0.7253469278329725772, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (0.9999847412109375, &r), 4.8916573191196772369, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, (1.0000152587890625, &r), 4.8918447504867045145, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, ( 1.5, &r), 0.20707843432557528095, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, ( 9.99, &r), 3.360235060345441639e-3, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, ( 10.0, &r), 3.353477310755806357e-3, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, ( 10.01, &r), 3.346739967281953346e-3, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, ( 100.0, &r), 3.333533347620158821e-5, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Q1_e, ( 1e10, &r), 3.333333333333333333e-21, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Ql_e, (10, -0.5, &r), -0.29165813966586752393,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Ql_e, (10,  0.5, &r), 0.29165813966586752393,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Ql_e, (10,  1.5, &r), 0.000014714232718207477406, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Ql_e, (100, -0.5, &r), -0.09492507395207282096,   TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Ql_e, (100,  0.5, &r), 0.09492507395207282096,    TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Ql_e, (100,  1.5, &r), 1.1628163435044121988e-43, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_legendre_Ql_e, (1000, -0.5, &r), -0.030105074974005303500, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Ql_e, (1000,  0.5, &r), 0.030105074974005303500,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_legendre_Ql_e, (1000,  1.1, &r), 1.0757258447825356443e-194, TEST_TOL3, GSL_SUCCESS);

  return s;
}

ext: c


text:
/* Copyright (c) 2016 PaddlePaddle Authors. All Rights Reserved.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License. */

#pragma once
#ifdef PADDLE_WITH_MKLML
#include "paddle/fluid/platform/dynload/mklml.h"
#endif

#ifdef PADDLE_USE_OPENBLAS
#include <cblas.h>
// remove typedef in openblas
#undef FLOAT
#undef INT
#undef SIZE
#endif

#include <cmath>
#include <vector>

#include "paddle/fluid/framework/eigen.h"
#include "paddle/fluid/framework/operator.h"
#include "paddle/fluid/framework/tensor.h"
#include "paddle/fluid/framework/tensor_util.h"
#include "paddle/fluid/platform/device_context.h"
#include "paddle/fluid/platform/enforce.h"

namespace paddle {
namespace operators {
namespace math {
template <typename DeviceContext, typename T, int Rank>
struct Transpose {
  void operator()(const DeviceContext& context, const framework::Tensor& in,
                  framework::Tensor* out, const std::vector<int>& axis);
};

template <typename DeviceContext, typename T>
struct SetConstant {
  void operator()(const DeviceContext& context, framework::Tensor* tensor,
                  T num);
};

template <typename Place>
void set_constant_with_place(const platform::DeviceContext& context,
                             framework::Tensor* tensor, float value);

void set_constant(const platform::DeviceContext& context,
                  framework::Tensor* tensor, float value);

template <typename DeviceContext, typename T>
struct RowwiseAdd {
  void operator()(const DeviceContext& context, const framework::Tensor& input,
                  const framework::Tensor& vec, framework::Tensor* output);
};

template <typename DeviceContext, typename T>
struct ColwiseSum {
  void operator()(const DeviceContext& context, const framework::Tensor& input,
                  framework::Tensor* vec);
};

template <typename DeviceContext, typename T>
struct RowwiseSum {
  void operator()(const DeviceContext& context, const framework::Tensor& input,
                  framework::Tensor* vec);
};

template <typename DeviceContext, typename T>
struct RowwiseMean {
  void operator()(const DeviceContext& context, const framework::Tensor& input,
                  framework::Tensor* vec);
};

}  // namespace math
}  // namespace operators
}  // namespace paddle

max_stars_repo_name: SnailTowardThesun/Paddle

max_issues_repo_name: SnailTowardThesun/Paddle

max_forks_repo_name: SnailTowardThesun/Paddle


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_expint.h>
#include "ccl.h"

static double einasto_norm_integrand(double x, void *params)
{
  double alpha = *((double *)(params));
  return x*x*exp(-2*(pow(x,alpha)-1)/alpha);
}

void ccl_einasto_norm_integral(int n_m, double *r_s, double *r_delta, double *alpha,
			       double *norm_out,int *status)
{
#pragma omp parallel default(none)			\
  shared(n_m, r_s, r_delta, alpha, norm_out, status)
  {
    int ii;
    int status_this=0;
    gsl_function F;
    gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);

    if (w == NULL)
      status_this = CCL_ERROR_MEMORY;

    if(status_this == 0) {
#pragma omp for
      for(ii=0;ii<n_m;ii++) {
	int qagstatus;
	double result, eresult;
	double x_max = r_delta[ii]/r_s[ii];
	F.function = &einasto_norm_integrand;
	F.params = &(alpha[ii]);
	qagstatus = gsl_integration_qag(&F, 0, x_max, 0, 1E-4,
					1000, GSL_INTEG_GAUSS31,
					w, &result, &eresult);
	if(qagstatus != GSL_SUCCESS) {
	  ccl_raise_gsl_warning(qagstatus, "ccl_haloprofile.c: ccl_einasto_norm_integral():");
	  status_this = CCL_ERROR_INTEG;
	  result = NAN;
	}
	norm_out[ii] = 4 * M_PI * r_s[ii] * r_s[ii] * r_s[ii] * result;
      }
    } //end omp for

    gsl_integration_workspace_free(w);
    if(status_this) {
      #pragma omp atomic write
      *status = status_this;
    }
  } //end omp parallel
}

static double hernquist_norm_integrand(double x, void *params)
{
  double opx=1+x;
  return x*x/(x*opx*opx*opx);
}

void ccl_hernquist_norm_integral(int n_m, double *r_s, double *r_delta,
			       double *norm_out,int *status)
{
#pragma omp parallel default(none)		\
  shared(n_m, r_s, r_delta, norm_out, status)
  {
    int ii;
    int status_this=0;
    gsl_function F;
    gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);

    if (w == NULL)
      status_this = CCL_ERROR_MEMORY;

    if(status_this == 0) {
#pragma omp for
      for(ii=0;ii<n_m;ii++) {
	int qagstatus;
	double result, eresult;
	double x_max = r_delta[ii]/r_s[ii];
	F.function = &hernquist_norm_integrand;
	F.params = NULL;
	qagstatus = gsl_integration_qag(&F, 0, x_max, 0, 1E-4,
					1000, GSL_INTEG_GAUSS31,
					w, &result, &eresult);
	if(qagstatus != GSL_SUCCESS) {
	  ccl_raise_gsl_warning(qagstatus, "ccl_haloprofile.c: ccl_hernquist_norm_integral():");
	  status_this = CCL_ERROR_INTEG;
	  result = NAN;
	}
	norm_out[ii] = 4 * M_PI * r_s[ii] * r_s[ii] * r_s[ii] * result;
      }
    } //end omp for

    gsl_integration_workspace_free(w);
    if(status_this) {
      #pragma omp atomic write
      *status = status_this;
    }
  } //end omp parallel
}

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
#ifndef _HistWidget_H
#define _HistWidget_H

#include <QWidget>
#include <gsl/gsl_histogram.h>

#include "src_config.h"

class QwtPlot;

class _HIST_EXPORT_ HistWidget : public QWidget
{
public:
    HistWidget (double xmin, double xmax, double ymin, double ymax, QWidget * parent=0, Qt::WindowFlags flags=0);
    virtual ~HistWidget (void);

    void setHistogram (const gsl_histogram * hist);

private:
    //
    // Variables
    //
    QwtPlot *plot;
    gsl_histogram * w_hist;
private:
    Q_OBJECT

};

#endif

max_stars_repo_name: YuriyRusinov/Histogram

max_issues_repo_name: YuriyRusinov/Histogram

max_forks_repo_name: YuriyRusinov/Histogram


text:
/*
Performing approximate Bayesian computation sequential Monte Carlo (Toni et al.
2009) for linear regression.

Synthetic data is generated by `ground_truth_and_analysis.ipynb`.

This script writes particles.csv, where each row corresponds to a particle and
each column corresponds to a round of SMC.

Defining #DEBUG_MODE will silence all writing to stdout. One may then add
printf statements in the code, and perhaps write the output to file as:
`./run.sh > output.txt`

Parameter ordering convention:
0 - gradient
1 - intercept
2 - standard deviation

Author: Juvid Aryaman
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sort_double.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_fit.h>

#define N_DATA 30
#define N_PARAMETERS 3

#define N_PARTICLES 20000

#define RND gsl_rng_uniform(r)
#define SEED 1

#define X_DATA_FILENAME "x.csv"
#define Y_DATA_FILENAME "y.csv"

// Global variables
/*Define the distance threshold for every round of SMC*/
double distance_threshold_schedule[] = {7.0, 6.375, 5.75, 5.125, 4.5, 3.875,
	3.25, 2.625, 2.0};
int N_ROUNDS_SMC = (int)(sizeof(distance_threshold_schedule) / sizeof(double));

#include "smc.h"
#include "lin_reg.h"

//#define DEBUG_MODE

int main(int argc, char *argv[]) {

#ifndef DEBUG_MODE
	printf("Threshold schedule:\n");
	print_double_array(distance_threshold_schedule, N_ROUNDS_SMC);
	printf("\n");
#endif

/* set up GSL RNG */
gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
/* end of GSL setup */

gsl_rng_set(r, SEED);

/////////////////////////
/*Read data*/
/////////////////////////

FILE *data_pointer_x, *data_pointer_y;

data_pointer_x = fopen(X_DATA_FILENAME, "r");
data_pointer_y = fopen(Y_DATA_FILENAME, "r");

double data_x[N_DATA];
double data_y[N_DATA];
int i, j, read_error_status_x, read_error_status_y;
for (i=0; i < N_DATA; i++){
	read_error_status_x = fscanf(data_pointer_x, "%lf\n", &data_x[i]);
	read_error_status_y = fscanf(data_pointer_y, "%lf\n", &data_y[i]);
}
if (read_error_status_x != 1){printf("Error reading X data\n"); return 0;}
if (read_error_status_y != 1){printf("Error reading Y data\n"); return 0;}

/////////////////////////
/*Initialise variables*/
/////////////////////////

/*Fit a linear model to the data, which will be used as summary statistics of
the data*/
double gradient_fit_data, intercept_fit_data, sigma_fit_data, cov00, cov01, cov11, sumsq;
int gsl_fit_return_value;
gsl_fit_return_value = gsl_fit_linear(data_x, 1, data_y, 1, N_DATA,
																			&intercept_fit_data, &gradient_fit_data,
																			&cov00, &cov01, &cov11, &sumsq);
if (gsl_fit_return_value != 0) {printf("Fit failed.\n"); return -1;}
sigma_fit_data = sqrt(sumsq/(N_DATA-2));

#ifndef DEBUG_MODE
	printf("gradient ML = %.8f\n", gradient_fit_data);
	printf("intercept ML = %.8f\n", intercept_fit_data);
	printf("sigma ML = %.8f\n", sigma_fit_data);
#endif

/*Make a (N_PARAMETERS X N_ROUNDS_SMC X N_PARTICLES) array to store all
particles at all rounds of SMC*/
double*** theta_particle;
theta_particle = (double***) malloc(N_PARAMETERS * sizeof(double**));
for (i = 0; i < N_PARAMETERS; i++){
	theta_particle[i] = (double**) malloc(N_ROUNDS_SMC * sizeof(double*));
}
for (i = 0; i < N_PARAMETERS; i++){
	for (j = 0; j < N_ROUNDS_SMC; j++){
		theta_particle[i][j] = (double*) malloc(N_PARTICLES * sizeof(double));
	}
}

double *simulated_data = (double*) malloc(N_DATA * sizeof(double));
double *distance = malloc(N_PARTICLES * sizeof(double));
double **weight;
weight = (double**) malloc(N_ROUNDS_SMC * sizeof(double*));
for (i = 0; i < N_ROUNDS_SMC; i++) {
	weight[i] = (double*) malloc(N_PARTICLES * sizeof(double));
}

double weight_normalizer = 0.0;

int time_smc=0; // an index of each round of SMC
int param_index_chosen, prior_violated;
int particle_index;
/////////////////////////
/*Perform ABC SMC*/
/////////////////////////

/*For every round of SMC*/
for (time_smc = 0; time_smc < N_ROUNDS_SMC; time_smc++) {
	#ifndef DEBUG_MODE
		printf("Round %d of SMC\n", time_smc);
	#endif

	/*Draw or perturb a particle and compute distance*/
	for (particle_index = 0; particle_index < N_PARTICLES; particle_index++) {
		// reset distance of particle
		distance[particle_index] = distance_threshold_schedule[time_smc] + 1.0;

		while (distance[particle_index] > distance_threshold_schedule[time_smc]) {
			if (time_smc == 0) {
				// Sample from the prior
				sample_prior(r, theta_particle, particle_index);
			}
			else{
				/*Sample from the old weights*/
				param_index_chosen = weighted_choice(r, weight[time_smc-1]);
				if ((param_index_chosen < 0)||(param_index_chosen >= N_PARTICLES)) {
					printf("Error in param_index_chosen\n");
					printf("time_smc = %d\n", time_smc);
					return -1;
				}

				perturb_particle(r, theta_particle, time_smc, param_index_chosen,
												 particle_index);

				// Check if prior support is 0
				prior_violated = check_prior_violated(theta_particle, time_smc,
																							particle_index);
				if(prior_violated == 1) continue;
			}

			simulate_dataset(r, theta_particle, data_x, simulated_data, time_smc,
											 particle_index);

			distance[particle_index] = distance_metric_sum_abs_res(simulated_data, data_y);
		}
	}

	#ifndef DEBUG_MODE
		printf("Particles sampled.\n");
	#endif

	/*Compute & Normalise weights. For uniform priors, and a uniform perturbation
	kernel, all surviving particles are weighted identically as 1/N_PARTICLES.
	Whilst seemingly inefficient, I keep this code here for clarity/generality.
	The bottleneck in computation time is the while loop above.*/
	if (time_smc==0){ for (i = 0; i < N_PARTICLES; i++) weight[0][i] = 1.0;}
	else{
		weight_normalizer = 0.0;
		for (particle_index = 0; particle_index < N_PARTICLES; particle_index++) {
			// kernel_pdf is uniform, so is independent of the parameters
			weight_normalizer +=  weight[time_smc-1][particle_index]*kernel_pdf();
		}
		for (particle_index = 0; particle_index < N_PARTICLES; particle_index++) {
			weight[time_smc][particle_index] = prior_pdf(theta_particle, time_smc,
																				 particle_index)/weight_normalizer;
		}
	}
	weight_normalizer = 0.0;
	for (i = 0; i < N_PARTICLES; i++)	weight_normalizer += weight[time_smc][i];
	for (i = 0; i < N_PARTICLES; i++){
		weight[time_smc][i] = weight[time_smc][i]/weight_normalizer;
	}

}

#ifndef DEBUG_MODE
	printf("Writing particles to file\n");
#endif
	write_particles_to_csv(theta_particle);

char weight_filename[] = "weights.csv";
write_2d_double_array_to_csv(weight, N_ROUNDS_SMC, N_PARTICLES, weight_filename);

#ifndef DEBUG_MODE
	printf("Done!\n");
#endif

return 0; //return from main
} //close main


text:
#pragma once

#ifndef __linux__
#error We only support Linux
#endif

#ifndef _GNU_SOURCE
#define _GNU_SOURCE
#endif

#ifdef __FAAS_CPP_WORKER
#if !defined(__FAAS_USED_IN_BINDING) && !defined(__FAAS_CPP_WORKER_SRC)
#error Need the source file to have __FAAS_USED_IN_BINDING defined
#endif
#endif

#ifdef __FAAS_NODE_ADDON
#if !defined(__FAAS_USED_IN_BINDING) && !defined(__FAAS_NODE_ADDON_SRC)
#error Need the source file to have __FAAS_USED_IN_BINDING defined
#endif
#define __FAAS_CXX_NO_EXCEPTIONS
#endif

#ifdef __FAAS_PYTHON_BINDING
#if !defined(__FAAS_USED_IN_BINDING) && !defined(__FAAS_PYTHON_BINDING_SRC)
#error Need the source file to have __FAAS_USED_IN_BINDING defined
#endif
#include <Python.h>
#if PY_VERSION_HEX < 0x03070000
#error FaaS Python binding requires Python 3.7+
#endif
#endif

// C includes
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stdint.h>
#include <time.h>
#include <unistd.h>

// C++ includes
#include <limits>
#include <atomic>
#include <memory>
#include <utility>
#include <optional>
#include <sstream>
#include <string>
#include <string_view>
#include <functional>
#include <algorithm>

// STL containers
#include <vector>
#include <queue>
#include <map>
#include <set>
#include <unordered_map>
#include <unordered_set>
#ifdef __FAAS_CPP_WORKER
#include <mutex>
#endif

// fmtlib
#define FMT_HEADER_ONLY
#include <fmt/core.h>
#include <fmt/format.h>

// Guidelines Support Library (GSL)
#include <gsl/gsl>

#include "base/diagnostic.h"

#ifdef __FAAS_SRC
#define __FAAS_HAVE_ABSL
#endif

#if defined(__FAAS_HAVE_ABSL) && !defined(__FAAS_USED_IN_BINDING)

// Will not include common absl headers in source files
// with __FAAS_USED_IN_BINDING defined

__BEGIN_THIRD_PARTY_HEADERS

#include <absl/base/call_once.h>
#include <absl/flags/flag.h>
#include <absl/time/time.h>
#include <absl/time/clock.h>
#include <absl/strings/str_cat.h>
#include <absl/strings/match.h>
#include <absl/strings/strip.h>
#include <absl/strings/str_split.h>
#include <absl/strings/numbers.h>
#include <absl/strings/ascii.h>
#include <absl/random/random.h>
#include <absl/random/distributions.h>
#include <absl/container/flat_hash_map.h>
#include <absl/container/flat_hash_set.h>
#include <absl/container/fixed_array.h>
#include <absl/container/inlined_vector.h>
#include <absl/synchronization/mutex.h>
#include <absl/synchronization/notification.h>
#include <absl/functional/bind_front.h>
#include <absl/algorithm/container.h>

__END_THIRD_PARTY_HEADERS

#endif  // defined(__FAAS_HAVE_ABSL) && !defined(__FAAS_USED_IN_BINDING)

#include "base/macro.h"
#include "base/logging.h"
#include "base/std_span.h"

max_stars_repo_name: jhweintraub/boki

max_issues_repo_name: jhweintraub/boki

max_forks_repo_name: jhweintraub/boki


text:
/**
 * @file spectralestimator.h
 * @brief Speech recognition front end.
 * @author John McDonough, Matthias Woelfel, Kenichi Kumatani
 */

#ifndef SPECTRALESTIMATOR_H
#define SPECTRALESTIMATOR_H

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <math.h>
#include <algorithm>

#include "matrix/gslmatrix.h"
#include "stream/stream.h"
#include "common/mlist.h"

#include <stdio.h>
#include <iostream>

/**
   @brief calculate the auto-correlation vector, r, with the time-discrete signal, x[n].
          Each component, r[k], represents the auto-correlation function of the input sample for lag, k.
          Accordingly, r[k] can be expressed as \cE{ x[n] x[n-k] }.
   @param const gsl_vector *samples[in] time-discrete signal which corresponds to x[n].
   @param gsl_vector *autoCorrelationVector[out] auto-correlation vector which corresponds to r.
   @note  the size of "autoCorrelationVector" should be "order+1" in order to calculate the LP model of the order, "order".
 */
void calcAutoCorrelationVector( const gsl_vector_float *samples, gsl_vector *autoCorrelationVector );

/**
   @brief calculate the linear prediction coefficients (LPC) with the Levinson-Durbin Algorithm.
   @param gsl_vector_float *autoCorrelationVector[in]
   @param gsl_vector_float *lpcs[out] lpcs[0,...,p-1] where p is the order of the model.
   @return prediction error
   @notice The estimated signal is defined as
           \hat{x}[n] = sum_{k=1}^p lpcs[k-1] * x[n-k].
           That is different from the definition in MATLAB lpc()
	   where it is \hat{x}(n) = 1.0 - sum_{k=1}^p lpcs(k+1) * x(n-k).
 */
float calcLPCswithLevinsonDurbin( gsl_vector_float *autoCorrelationVector, gsl_vector_float *lpcs );

// ----- definition for class `LPCSpectrumEstimator' -----
//

/**
   @class calculate LPCs with the Levinson-Durbin algorithm, and estimate the LPC envelope.
   @note  You can feed samples to this object with calcLPCs( const gsl_vector_float* samples ),
	  and the LPCs are then computed.
 */
class LPCSpectrumEstimator : public VectorFloatFeatureStream {
public:
  LPCSpectrumEstimator( const VectorFloatFeatureStreamPtr& samp, unsigned order, unsigned fftLen, const String& nm = "LPCSpectrumEstimator");
  ~LPCSpectrumEstimator();

  const gsl_vector_float* next(int frame_no = -5);
  void reset();

  /**
     @brief get the auto-correlation vector, r[0,...,order-1],
            where r[k] represents the auto-correlation function of the input sample for lag, k+1.
     @note you have to call calcLPCs() before you call this method.
   */
  const gsl_vector_float* getAutoCorrelationVector(){
    return (const gsl_vector_float*)_r;
  }
  /**
     @brief get the prediction error of the LP model.
     @note you have to call calcLPCs() before you call this method.
   */
  float getPredictionError(){
    return _e;
  }
  /**
     @brief get the LPCs.
     @note you have to call calcLPCs() before you call this method.
   */
  const gsl_vector_float* getLPCs(){
    return _lpcs; //_vector;
  }

  gsl_vector_complex* getFTLPCs(){
    return _ftlpcs;
  }

  const gsl_vector_float* calcLPCs( const gsl_vector_float* samples );
  void calcLPEnvelope();

private:
  void complexPackedArrayUnpack(const double* halfcomplex_coefficient, gsl_vector_complex* complex_coefficient) const;
private:
  static const double             _epsilon;
  VectorFloatFeatureStreamPtr        _samp;
  unsigned                          _order; /* the LP model of the order */
  unsigned                         _fftLen;
  unsigned                        _fftLen2;
  float                                 _e; /* prediction error */
  gsl_vector_float*                     _r; /* _order x 1 auto-correlation vector _r[1,...,p]; See [1] */
  gsl_vector*             _autoCorrelation; /* (_order+1) x 1 auto-correlation vector, where the first element corresponds to the power */
  gsl_vector_float*                  _lpcs; /* LPCs */
  double*                          _wslpcs; /* workspace for the LPCs */
  gsl_vector_complex*              _ftlpcs;
};

typedef Inherit<LPCSpectrumEstimator, VectorFloatFeatureStreamPtr> LPCSpectrumEstimatorPtr;

// ----- definition for class `CepstralSpectrumEstimator' -----
//
class CepstralSpectrumEstimator : public VectorFloatFeatureStream {
public:
  CepstralSpectrumEstimator(const VectorComplexFeatureStreamPtr& source, unsigned order = 14, unsigned fftLen = 256,
			    double logPadding = 1.0, const String& nm = "CepstralSpectrumEstimator");
  ~CepstralSpectrumEstimator();

  const gsl_vector_float* next(int frame_no = -5);
  void reset() { _source->reset();  VectorFloatFeatureStream::reset(); }

private:
  VectorComplexFeatureStreamPtr			_source;
  const unsigned				_order;
  const unsigned				_fftLen;
  const double					_logPadding;

  double*					_samples;
  gsl_vector_complex*				_spectrum;
  gsl_vector_complex*				_cepstrum;
};

typedef Inherit<CepstralSpectrumEstimator, VectorFloatFeatureStreamPtr> CepstralSpectrumEstimatorPtr;

// ----- definition for class `SEMNB' -----
//

/**
   @class
 */

class SEMNB : public Countable {
public:
  SEMNB( unsigned order, unsigned fftLen, const String& nm = "SEMNB");
  ~SEMNB();

  const gsl_vector*	calcDerivativeOfDeviation( LPCSpectrumEstimatorPtr &lpcSEPtr );
  virtual void		reset();
  const gsl_vector*	getLPEnvelope() const { return _lpEnvelope; }

private:
  void			_calcEigen();
  void			_calcLambda_bar();

  double		_dLambdak_dPm(unsigned k, unsigned m);
  void			_dLambdabar_dPm(unsigned m);
  void			_da_dPm(unsigned m);
  float			_dSigma_dP(const gsl_vector_float* lpcs, double epsilonp, unsigned m);

  double		_dSp_dan(const gsl_vector_float* lpcs, unsigned m, unsigned n);
  double		_dSp_dPm(const gsl_vector_float* lpcs, unsigned m);
  double		_depsilonp_dPm(const gsl_vector_float* lpcs, unsigned m);
  void			_calc_dr_dPm(unsigned m);
  double		_dPhat_dP(const gsl_vector_float* lpcs, double epsilonp, unsigned m);

  const unsigned				_order;
  const unsigned				_fftLen;
  const unsigned				_fftLen2;

  gsl_vector*					_derivative;
  gsl_matrix*					_R;
  gsl_matrix*					_Rcopy;
  gsl_eigen_symmv_workspace*			_workSpace;
  gsl_vector*					_evalues;
  gsl_matrix*					_evectors;
  gsl_vector*	                                _lpEnvelope;
  gsl_vector_complex*                           _ftlpcs;

  gsl_vector*					_r;
  gsl_vector*					_dr_dPm;
  gsl_vector*					_dak_dPm;
  gsl_matrix*					_Lambda_bar;
  gsl_matrix*					_dLambdamatrix;

  gsl_matrix*					_C;
  gsl_matrix*					_D;
};

typedef refcountable_ptr<SEMNB> SEMNBPtr;

#endif


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

#include "allvars.h"
#include "proto.h"
#ifdef COSMIC_RAYS
#include "cosmic_rays.h"
#endif
#ifdef CS_MODEL
#include "cs_metals.h"
#endif

/*! Structure for communication during the density computation. Holds data that is sent to other processors.
 */
static struct densdata_in
{
  MyDouble Pos[3];
  MyFloat Vel[3];
  MyFloat Hsml;
#ifdef VOLUME_CORRECTION
  MyFloat DensityOld;
#endif
#ifdef WINDS
  MyFloat DelayTime;
#endif
  int NodeList[NODELISTLENGTH];
#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS) || defined(TRACEDIVB)
  MyFloat BPred[3];
#endif
#ifdef EULERPOTENTIALS
  MyFloat EulerA, EulerB;
#endif
#ifdef CS_MODEL
  MyFloat DensityOld, Entropy;
#endif
}
 *DensDataIn, *DensDataGet;

static struct densdata_out
{
  MyLongDouble Rho;
  MyLongDouble DhsmlDensity;
  MyLongDouble Ngb;
#ifndef NAVIERSTOKES
  MyLongDouble Div, Rot[3];
#else
  MyFloat DV[3][3];
#endif
#ifdef MAGNETIC
#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS)
  MyFloat RotB[3];
#endif
#ifdef TRACEDIVB
  MyFloat divB;
#endif
#endif

#ifdef RADTRANSFER_FLUXLIMITER
  MyFloat Grad_ngamma[3];
#endif

#if defined(BLACK_HOLES)
  MyLongDouble SmoothedEntr;
#endif
#ifdef CONDUCTION_SATURATION
  MyFloat GradEntr[3];
#endif

#ifdef BLACK_HOLES
  MyLongDouble GasVel[3];
#endif

#ifdef HYDRO_COST_FACTOR
  int Ninteractions;
#endif

#ifdef VOLUME_CORRECTION
  MyFloat DensityStd;
#endif

#ifdef EULERPOTENTIALS
  MyFloat dEulerA[3], dEulerB[3];
#endif

}
 *DensDataResult, *DensDataOut;

/*! \file density.c
 *  \brief SPH density computation and smoothing length determination
 *
 *  This file contains the "first SPH loop", where the SPH densities and some
 *  auxiliary quantities are computed.  There is also functionality that
 *  corrects the smoothing length if needed.
 */

/*! This function computes the local density for each active SPH particle, the
 * number of neighbours in the current smoothing radius, and the divergence
 * and rotation of the velocity field.  The pressure is updated as well.  If a
 * particle with its smoothing region is fully inside the local domain, it is
 * not exported to the other processors. The function also detects particles
 * that have a number of neighbours outside the allowed tolerance range. For
 * these particles, the smoothing length is adjusted accordingly, and the
 * density() computation is called again.  Note that the smoothing length is
 * not allowed to fall below the lower bound set by MinGasHsml (this may mean
 * that one has to deal with substantially more than normal number of
 * neighbours.)
 */
void density(void)
{
  MyFloat *Left, *Right;
  int i, j, ndone, ndone_flag, npleft, dt_step, dummy, iter = 0;
  int ngrp, sendTask, recvTask, place, nexport, nimport;
  long long ntot;
  double dmax1, dmax2, fac;
  double timeall = 0, timecomp1 = 0, timecomp2 = 0, timecommsumm1 = 0, timecommsumm2 = 0, timewait1 =
    0, timewait2 = 0;
  double timecomp, timecomm, timewait;
  double dt_entr, tstart, tend, t0, t1;
  double desnumngb;

#ifdef COSMIC_RAYS
  int CRpop;
#endif

#if defined(NAVIERSTOKES)
  int k;
#endif
#ifdef NAVIERSTOKES
  double dvel[3][3];
  double rotx, roty, rotz;
#endif

#if defined(SOFTEREQS)
  double a3inv;

  if(All.ComovingIntegrationOn)
    a3inv = 1 / (All.Time * All.Time * All.Time);
  else
    a3inv = 1;
#endif

#ifdef EULERPOTENTIALS
  double efak;

  if(All.ComovingIntegrationOn)
    efak = 1. / All.Time / All.HubbleParam;
  else
    efak = 1;
#endif

  CPU_Step[CPU_DENSMISC] += measure_time();

  Left = (MyFloat *) mymalloc(NumPart * sizeof(MyFloat));
  Right = (MyFloat *) mymalloc(NumPart * sizeof(MyFloat));

  for(i = FirstActiveParticle; i >= 0; i = NextActiveParticle[i])
    {
      if(density_isactive(i))
	{
	  Left[i] = Right[i] = 0;

#ifdef BLACK_HOLES
	  P[i].SwallowID = 0;
#endif
#if defined(BLACK_HOLES) && defined(FLTROUNDOFFREDUCTION)
	  if(P[i].Type == 0)
	    SphP[i].i.dInjected_BH_Energy = SphP[i].i.Injected_BH_Energy;
#endif
	}
    }

  /* allocate buffers to arrange communication */

  Ngblist = (int *) mymalloc(NumPart * sizeof(int));

  All.BunchSize =
    (int) ((All.BufferSize * 1024 * 1024) / (sizeof(struct data_index) + sizeof(struct data_nodelist) +
					     sizeof(struct densdata_in) + sizeof(struct densdata_out) +
					     sizemax(sizeof(struct densdata_in),
						     sizeof(struct densdata_out))));
  DataIndexTable = (struct data_index *) mymalloc(All.BunchSize * sizeof(struct data_index));
  DataNodeList = (struct data_nodelist *) mymalloc(All.BunchSize * sizeof(struct data_nodelist));

  t0 = second();

  desnumngb = All.DesNumNgb;

  /* we will repeat the whole thing for those particles where we didn't find enough neighbours */
  do
    {
      i = FirstActiveParticle;	/* begin with this index */

      do
	{
	  for(j = 0; j < NTask; j++)
	    {
	      Send_count[j] = 0;
	      Exportflag[j] = -1;
	    }

	  /* do local particles and prepare export list */
	  tstart = second();
	  for(nexport = 0; i >= 0; i = NextActiveParticle[i])
	    {
	      if(density_isactive(i))
		{
		  if(density_evaluate(i, 0, &nexport, Send_count) < 0)
		    break;
		}
	    }
	  tend = second();
	  timecomp1 += timediff(tstart, tend);

#ifdef MYSORT
	  mysort_dataindex(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);
#else
	  qsort(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);
#endif

	  tstart = second();

	  MPI_Allgather(Send_count, NTask, MPI_INT, Sendcount_matrix, NTask, MPI_INT, MPI_COMM_WORLD);

	  tend = second();
	  timewait1 += timediff(tstart, tend);

	  for(j = 0, nimport = 0, Recv_offset[0] = 0, Send_offset[0] = 0; j < NTask; j++)
	    {
	      Recv_count[j] = Sendcount_matrix[j * NTask + ThisTask];
	      nimport += Recv_count[j];

	      if(j > 0)
		{
		  Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
		  Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
		}
	    }

	  DensDataGet = (struct densdata_in *) mymalloc(nimport * sizeof(struct densdata_in));
	  DensDataIn = (struct densdata_in *) mymalloc(nexport * sizeof(struct densdata_in));

	  /* prepare particle data for export */
	  for(j = 0; j < nexport; j++)
	    {
	      place = DataIndexTable[j].Index;

	      DensDataIn[j].Pos[0] = P[place].Pos[0];
	      DensDataIn[j].Pos[1] = P[place].Pos[1];
	      DensDataIn[j].Pos[2] = P[place].Pos[2];
	      DensDataIn[j].Hsml = PPP[place].Hsml;

#ifdef CS_MODEL
	      DensDataIn[j].DensityOld = SphP[place].DensityOld;
	      DensDataIn[j].Entropy = SphP[place].Entropy;
#endif
	      memcpy(DensDataIn[j].NodeList,
		     DataNodeList[DataIndexTable[j].IndexGet].NodeList, NODELISTLENGTH * sizeof(int));

#if defined(BLACK_HOLES)
	      if(P[place].Type != 0)
		{
		  DensDataIn[j].Vel[0] = 0;
		  DensDataIn[j].Vel[1] = 0;
		  DensDataIn[j].Vel[2] = 0;
		}
	      else
#endif
		{
		  DensDataIn[j].Vel[0] = SphP[place].VelPred[0];
		  DensDataIn[j].Vel[1] = SphP[place].VelPred[1];
		  DensDataIn[j].Vel[2] = SphP[place].VelPred[2];
		}
#ifdef VOLUME_CORRECTION
	      DensDataIn[j].DensityOld = SphP[place].DensityOld;
#endif

#ifdef EULERPOTENTIALS
	      DensDataIn[j].EulerA = SphP[place].EulerA;
	      DensDataIn[j].EulerB = SphP[place].EulerB;
#endif

#ifdef WINDS
	      DensDataIn[j].DelayTime = SphP[place].DelayTime;
#endif

#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS) || defined(TRACEDIVB)
	      DensDataIn[j].BPred[0] = SphP[place].BPred[0];
	      DensDataIn[j].BPred[1] = SphP[place].BPred[1];
	      DensDataIn[j].BPred[2] = SphP[place].BPred[2];
#endif
	    }
	  /* exchange particle data */
	  tstart = second();
	  for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	    {
	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;

	      if(recvTask < NTask)
		{
		  if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		    {
		      /* get the particles */
		      MPI_Sendrecv(&DensDataIn[Send_offset[recvTask]],
				   Send_count[recvTask] * sizeof(struct densdata_in), MPI_BYTE,
				   recvTask, TAG_DENS_A,
				   &DensDataGet[Recv_offset[recvTask]],
				   Recv_count[recvTask] * sizeof(struct densdata_in), MPI_BYTE,
				   recvTask, TAG_DENS_A, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		    }
		}
	    }
	  tend = second();
	  timecommsumm1 += timediff(tstart, tend);

	  myfree(DensDataIn);
	  DensDataResult = (struct densdata_out *) mymalloc(nimport * sizeof(struct densdata_out));
	  DensDataOut = (struct densdata_out *) mymalloc(nexport * sizeof(struct densdata_out));

	  /* now do the particles that were sent to us */

	  tstart = second();
	  for(j = 0; j < nimport; j++)
	    density_evaluate(j, 1, &dummy, &dummy);
	  tend = second();
	  timecomp2 += timediff(tstart, tend);

	  if(i < 0)
	    ndone_flag = 1;
	  else
	    ndone_flag = 0;

	  tstart = second();
	  MPI_Allreduce(&ndone_flag, &ndone, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
	  tend = second();
	  timewait2 += timediff(tstart, tend);

	  /* get the result */
	  tstart = second();
	  for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	    {
	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;
	      if(recvTask < NTask)
		{
		  if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		    {
		      /* send the results */
		      MPI_Sendrecv(&DensDataResult[Recv_offset[recvTask]],
				   Recv_count[recvTask] * sizeof(struct densdata_out),
				   MPI_BYTE, recvTask, TAG_DENS_B,
				   &DensDataOut[Send_offset[recvTask]],
				   Send_count[recvTask] * sizeof(struct densdata_out),
				   MPI_BYTE, recvTask, TAG_DENS_B, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		    }
		}

	    }
	  tend = second();
	  timecommsumm2 += timediff(tstart, tend);

	  /* add the result to the local particles */
	  tstart = second();
	  for(j = 0; j < nexport; j++)
	    {
	      place = DataIndexTable[j].Index;

	      PPP[place].n.dNumNgb += DensDataOut[j].Ngb;
#ifdef HYDRO_COST_FACTOR
	      P[place].GravCost += HYDRO_COST_FACTOR * DensDataOut[j].Ninteractions;
#endif

	      if(P[place].Type == 0)
		{
		  SphP[place].d.dDensity += DensDataOut[j].Rho;
		  SphP[place].h.dDhsmlDensityFactor += DensDataOut[j].DhsmlDensity;

#ifndef NAVIERSTOKES
		  SphP[place].v.dDivVel += DensDataOut[j].Div;
		  SphP[place].r.dRot[0] += DensDataOut[j].Rot[0];
		  SphP[place].r.dRot[1] += DensDataOut[j].Rot[1];
		  SphP[place].r.dRot[2] += DensDataOut[j].Rot[2];
#else
		  for(k = 0; k < 3; k++)
		    {
		      SphP[place].u.DV[k][0] += DensDataOut[j].DV[k][0];
		      SphP[place].u.DV[k][1] += DensDataOut[j].DV[k][1];
		      SphP[place].u.DV[k][2] += DensDataOut[j].DV[k][2];
		    }
#endif

#ifdef VOLUME_CORRECTION
		  SphP[place].DensityStd += DensDataOut[j].DensityStd;
#endif

#ifdef CONDUCTION_SATURATION
		  SphP[place].GradEntr[0] += DensDataOut[j].GradEntr[0];
		  SphP[place].GradEntr[1] += DensDataOut[j].GradEntr[1];
		  SphP[place].GradEntr[2] += DensDataOut[j].GradEntr[2];
#endif

#ifdef RADTRANSFER_FLUXLIMITER
		  SphP[place].Grad_ngamma[0] += DensDataOut[j].Grad_ngamma[0];
		  SphP[place].Grad_ngamma[1] += DensDataOut[j].Grad_ngamma[1];
		  SphP[place].Grad_ngamma[2] += DensDataOut[j].Grad_ngamma[2];
#endif

#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS)
		  SphP[place].RotB[0] += DensDataOut[j].RotB[0];
		  SphP[place].RotB[1] += DensDataOut[j].RotB[1];
		  SphP[place].RotB[2] += DensDataOut[j].RotB[2];
#endif

#ifdef TRACEDIVB
		  SphP[place].divB += DensDataOut[j].divB;
#endif

#ifdef EULERPOTENTIALS
		  SphP[place].dEulerA[0] += DensDataOut[j].dEulerA[0];
		  SphP[place].dEulerA[1] += DensDataOut[j].dEulerA[1];
		  SphP[place].dEulerA[2] += DensDataOut[j].dEulerA[2];
		  SphP[place].dEulerB[0] += DensDataOut[j].dEulerB[0];
		  SphP[place].dEulerB[1] += DensDataOut[j].dEulerB[1];
		  SphP[place].dEulerB[2] += DensDataOut[j].dEulerB[2];
#endif
		}

#if defined(RADTRANSFER) || defined(SNIA_HEATING)
	      if(P[place].Type == 4)
		P[place].DensAroundStar += DensDataOut[j].Rho;
#endif

#ifdef BLACK_HOLES
	      if(P[place].Type == 5)
		{
		  P[place].b1.dBH_Density += DensDataOut[j].Rho;
		  P[place].b2.dBH_Entropy += DensDataOut[j].SmoothedEntr;
		  P[place].b3.dBH_SurroundingGasVel[0] += DensDataOut[j].GasVel[0];
		  P[place].b3.dBH_SurroundingGasVel[1] += DensDataOut[j].GasVel[1];
		  P[place].b3.dBH_SurroundingGasVel[2] += DensDataOut[j].GasVel[2];
		}
#endif
	    }
	  tend = second();
	  timecomp1 += timediff(tstart, tend);

	  myfree(DensDataOut);
	  myfree(DensDataResult);
	  myfree(DensDataGet);
	}
      while(ndone < NTask);

#ifdef FLTROUNDOFFREDUCTION
      for(i = FirstActiveParticle; i >= 0; i = NextActiveParticle[i])
	if(density_isactive(i))
	  {
	    PPP[i].n.NumNgb = FLT(PPP[i].n.dNumNgb);

	    if(P[i].Type == 0)
	      {
		SphP[i].d.Density = FLT(SphP[i].d.dDensity);
		SphP[i].h.DhsmlDensityFactor = FLT(SphP[i].h.dDhsmlDensityFactor);
		SphP[i].v.DivVel = FLT(SphP[i].v.dDivVel);
		for(j = 0; j < 3; j++)
		  SphP[i].r.Rot[j] = FLT(SphP[i].r.dRot[j]);
	      }

#ifdef BLACK_HOLES
	    if(P[i].Type == 5)
	      {
		P[i].b1.BH_Density = FLT(P[i].b1.dBH_Density);
		P[i].b2.BH_Entropy = FLT(P[i].b2.dBH_Entropy);
		for(j = 0; j < 3; j++)
		  P[i].b3.BH_SurroundingGasVel[j] = FLT(P[i].b3.dBH_SurroundingGasVel[j]);
	      }
#endif
	  }
#endif

      /* do final operations on results */
      tstart = second();
      for(i = FirstActiveParticle, npleft = 0; i >= 0; i = NextActiveParticle[i])
	{
	  if(density_isactive(i))
	    {
	      if(P[i].Type == 0)
		{
		  if(SphP[i].d.Density > 0)
		    {
#ifdef VOLUME_CORRECTION
		      SphP[i].DensityOld = SphP[i].DensityStd;
#endif
		      SphP[i].h.DhsmlDensityFactor *= PPP[i].Hsml / (NUMDIMS * SphP[i].d.Density);
		      if(SphP[i].h.DhsmlDensityFactor > -0.9)	/* note: this would be -1 if only a single particle at zero lag is found */
			SphP[i].h.DhsmlDensityFactor = 1 / (1 + SphP[i].h.DhsmlDensityFactor);
		      else
			SphP[i].h.DhsmlDensityFactor = 1;
#ifndef NAVIERSTOKES
		      SphP[i].r.CurlVel = sqrt(SphP[i].r.Rot[0] * SphP[i].r.Rot[0] +
					       SphP[i].r.Rot[1] * SphP[i].r.Rot[1] +
					       SphP[i].r.Rot[2] * SphP[i].r.Rot[2]) / SphP[i].d.Density;

		      SphP[i].v.DivVel /= SphP[i].d.Density;
#else
		      for(k = 0; k < 3; k++)
			{
			  dvel[k][0] = SphP[i].u.DV[k][0] / SphP[i].d.Density;
			  dvel[k][1] = SphP[i].u.DV[k][1] / SphP[i].d.Density;
			  dvel[k][2] = SphP[i].u.DV[k][2] / SphP[i].d.Density;
			}
		      SphP[i].u.s.DivVel = dvel[0][0] + dvel[1][1] + dvel[2][2];

		      SphP[i].u.s.StressDiag[0] = 2 * dvel[0][0] - 2.0 / 3 * SphP[i].u.s.DivVel;
		      SphP[i].u.s.StressDiag[1] = 2 * dvel[1][1] - 2.0 / 3 * SphP[i].u.s.DivVel;
		      SphP[i].u.s.StressDiag[2] = 2 * dvel[2][2] - 2.0 / 3 * SphP[i].u.s.DivVel;

		      SphP[i].u.s.StressOffDiag[0] = dvel[0][1] + dvel[1][0];	/* xy */
		      SphP[i].u.s.StressOffDiag[1] = dvel[0][2] + dvel[2][0];	/* xz */
		      SphP[i].u.s.StressOffDiag[2] = dvel[1][2] + dvel[2][1];	/* yz */

#ifdef NAVIERSTOKES_BULK
		      SphP[i].u.s.StressBulk = All.NavierStokes_BulkViscosity * SphP[i].u.s.DivVel;
#endif
		      rotx = dvel[1][2] - dvel[2][1];
		      roty = dvel[2][0] - dvel[0][2];
		      rotz = dvel[0][1] - dvel[1][0];
		      SphP[i].u.s.CurlVel = sqrt(rotx * rotx + roty * roty + rotz * rotz);
#endif

#ifdef CONDUCTION_SATURATION
		      SphP[i].GradEntr[0] /= SphP[i].d.Density;
		      SphP[i].GradEntr[1] /= SphP[i].d.Density;
		      SphP[i].GradEntr[2] /= SphP[i].d.Density;
#endif

#ifdef RADTRANSFER_FLUXLIMITER
		      SphP[i].Grad_ngamma[0] /= SphP[i].d.Density;
		      SphP[i].Grad_ngamma[1] /= SphP[i].d.Density;
		      SphP[i].Grad_ngamma[2] /= SphP[i].d.Density;
#endif

#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS)
		      SphP[i].RotB[0] /= SphP[i].d.Density;
		      SphP[i].RotB[1] /= SphP[i].d.Density;
		      SphP[i].RotB[2] /= SphP[i].d.Density;
#endif

#ifdef TRACEDIVB
		      SphP[i].divB /= SphP[i].d.Density;
#endif

#ifdef EULERPOTENTIALS
		      SphP[i].dEulerA[0] *= efak / SphP[i].d.Density;
		      SphP[i].dEulerA[1] *= efak / SphP[i].d.Density;
		      SphP[i].dEulerA[2] *= efak / SphP[i].d.Density;
		      SphP[i].dEulerB[0] *= efak / SphP[i].d.Density;
		      SphP[i].dEulerB[1] *= efak / SphP[i].d.Density;
		      SphP[i].dEulerB[2] *= efak / SphP[i].d.Density;
#ifdef	PERIODIC
		      SphP[i].dEulerB[2] = 1.0;
#endif
		      SphP[i].BPred[0] =
			SphP[i].dEulerA[1] * SphP[i].dEulerB[2] - SphP[i].dEulerA[2] * SphP[i].dEulerB[1];
		      SphP[i].BPred[1] =
			SphP[i].dEulerA[2] * SphP[i].dEulerB[0] - SphP[i].dEulerA[0] * SphP[i].dEulerB[2];
		      SphP[i].BPred[2] =
			SphP[i].dEulerA[0] * SphP[i].dEulerB[1] - SphP[i].dEulerA[1] * SphP[i].dEulerB[0];
#endif
		    }

#ifndef WAKEUP
		  dt_step = (P[i].TimeBin ? (1 << P[i].TimeBin) : 0);
#else
		  dt_step = P[i].dt_step;
#endif
		  dt_entr = (All.Ti_Current - (P[i].Ti_begstep + dt_step / 2)) * All.Timebase_interval;

#ifndef EOS_DEGENERATE
#ifndef MHM
#ifndef SOFTEREQS
		  SphP[i].Pressure =
		    (SphP[i].Entropy + SphP[i].e.DtEntropy * dt_entr) * pow(SphP[i].d.Density, GAMMA);
#else
		  /* use an intermediate EQS, between isothermal and the full multiphase model */
		  if(SphP[i].d.Density * a3inv >= All.PhysDensThresh)
		    SphP[i].Pressure = All.FactorForSofterEQS *
		      (SphP[i].Entropy + SphP[i].e.DtEntropy * dt_entr) * pow(SphP[i].d.Density, GAMMA) +
		      (1 - All.FactorForSofterEQS) * GAMMA_MINUS1 * SphP[i].d.Density * All.InitGasU;
		  else
		    SphP[i].Pressure =
		      (SphP[i].Entropy + SphP[i].e.DtEntropy * dt_entr) * pow(SphP[i].d.Density, GAMMA);
#endif
#else
		  /* Here we use an isothermal equation of state */
		  SphP[i].Pressure = GAMMA_MINUS1 * SphP[i].d.Density * All.InitGasU;
		  SphP[i].Entropy = SphP[i].Pressure / pow(SphP[i].d.Density, GAMMA);
#endif
#else
		  /* call tabulated eos with physical units */
#ifdef EOS_ENERGY
		  eos_calc_egiven2(SphP[i].d.Density * All.UnitDensity_in_cgs, SphP[i].xnuc,
				   SphP[i].Entropy + SphP[i].e.DtEntropy * dt_entr * All.UnitTime_in_s,
				   &SphP[i].temp, &SphP[i].Pressure);
		  SphP[i].Pressure /= All.UnitPressure_in_cgs;
#else
		  eos_calc_sgiven(SphP[i].d.Density * All.UnitDensity_in_cgs, SphP[i].xnuc,
				  SphP[i].Entropy + SphP[i].e.DtEntropy * dt_entr * All.UnitTime_in_s,
				  &SphP[i].temp, &SphP[i].Pressure, &SphP[i].u);
		  SphP[i].Pressure /= All.UnitPressure_in_cgs;
#endif
#endif

#ifdef COSMIC_RAYS
		  for(CRpop = 0; CRpop < NUMCRPOP; CRpop++)
		    {
		      CR_Particle_Update(SphP + i, CRpop);
#ifndef CR_NOPRESSURE
		      SphP[i].Pressure += CR_Comoving_Pressure(SphP + i, CRpop);
#endif
		    }
#endif
		}

#ifdef BLACK_HOLES
	      if(P[i].Type == 5)
		{
		  if(P[i].b1.BH_Density > 0)
		    {
		      P[i].b2.BH_Entropy /= P[i].b1.BH_Density;
		      P[i].b3.BH_SurroundingGasVel[0] /= P[i].b1.BH_Density;
		      P[i].b3.BH_SurroundingGasVel[1] /= P[i].b1.BH_Density;
		      P[i].b3.BH_SurroundingGasVel[2] /= P[i].b1.BH_Density;
		    }
		}
#endif

	      /* now check whether we had enough neighbours */

	      desnumngb = All.DesNumNgb;

#ifdef BLACK_HOLES
	      if(P[i].Type == 5)
		desnumngb = All.DesNumNgb * All.BlackHoleNgbFactor;
#endif

#ifdef RADTRANSFER
	      if(P[i].Type == 4)
		desnumngb = 64;//NORM_COEFF * KERNEL_COEFF_1;	/* will assign the stellar luminosity to very few (one actually) gas particles */
#endif

	      if(PPP[i].n.NumNgb < (desnumngb - All.MaxNumNgbDeviation) ||
		 (PPP[i].n.NumNgb > (desnumngb + All.MaxNumNgbDeviation)
		  && PPP[i].Hsml > (1.01 * All.MinGasHsml)))
		{
		  /* need to redo this particle */
		  npleft++;

		  if(Left[i] > 0 && Right[i] > 0)
		    if((Right[i] - Left[i]) < 1.0e-3 * Left[i])
		      {
			/* this one should be ok */
			npleft--;
			P[i].TimeBin = -P[i].TimeBin - 1;	/* Mark as inactive */
			continue;
		      }

		  if(PPP[i].n.NumNgb < (desnumngb - All.MaxNumNgbDeviation))
		    Left[i] = DMAX(PPP[i].Hsml, Left[i]);
		  else
		    {
		      if(Right[i] != 0)
			{
			  if(PPP[i].Hsml < Right[i])
			    Right[i] = PPP[i].Hsml;
			}
		      else
			Right[i] = PPP[i].Hsml;
		    }

		  if(iter >= MAXITER - 10)
		    {
		      printf
			("i=%d task=%d ID=%d Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g\n   pos=(%g|%g|%g)\n",
			 i, ThisTask, (int) P[i].ID, PPP[i].Hsml, Left[i], Right[i],
			 (float) PPP[i].n.NumNgb, Right[i] - Left[i], P[i].Pos[0], P[i].Pos[1], P[i].Pos[2]);
		      fflush(stdout);
		    }

		  if(Right[i] > 0 && Left[i] > 0)
		    PPP[i].Hsml = pow(0.5 * (pow(Left[i], 3) + pow(Right[i], 3)), 1.0 / 3);
		  else
		    {
		      if(Right[i] == 0 && Left[i] == 0)
			endrun(8188);	/* can't occur */

		      if(Right[i] == 0 && Left[i] > 0)
			{
			  if(P[i].Type == 0 && fabs(PPP[i].n.NumNgb - desnumngb) < 0.5 * desnumngb)
			    {
			      fac = 1 - (PPP[i].n.NumNgb -
					 desnumngb) / (NUMDIMS * PPP[i].n.NumNgb) *
				SphP[i].h.DhsmlDensityFactor;

			      if(fac < 1.26)
				PPP[i].Hsml *= fac;
			      else
				PPP[i].Hsml *= 1.26;
			    }
			  else
			    PPP[i].Hsml *= 1.26;
			}

		      if(Right[i] > 0 && Left[i] == 0)
			{
			  if(P[i].Type == 0 && fabs(PPP[i].n.NumNgb - desnumngb) < 0.5 * desnumngb)
			    {
			      fac = 1 - (PPP[i].n.NumNgb -
					 desnumngb) / (NUMDIMS * PPP[i].n.NumNgb) *
				SphP[i].h.DhsmlDensityFactor;

			      if(fac > 1 / 1.26)
				PPP[i].Hsml *= fac;
			      else
				PPP[i].Hsml /= 1.26;
			    }
			  else
			    PPP[i].Hsml /= 1.26;
			}
		    }

		  if(PPP[i].Hsml < All.MinGasHsml)
		    PPP[i].Hsml = All.MinGasHsml;

#ifdef BLACK_HOLES
		  if(P[i].Type == 5)
		    if(Left[i] > All.BlackHoleMaxAccretionRadius)
		      {
			/* this will stop the search for a new BH smoothing length in the next iteration */
			PPP[i].Hsml = Left[i] = Right[i] = All.BlackHoleMaxAccretionRadius;
		      }
#endif
		}
	      else
		P[i].TimeBin = -P[i].TimeBin - 1;	/* Mark as inactive */
	    }
	}
      tend = second();
      timecomp1 += timediff(tstart, tend);

      sumup_large_ints(1, &npleft, &ntot);

      if(ntot > 0)
	{
	  iter++;

	  if(iter > 0 && ThisTask == 0)
	    {
	      printf("ngb iteration %d: need to repeat for %d%09d particles.\n", iter,
		     (int) (ntot / 1000000000), (int) (ntot % 1000000000));
	      fflush(stdout);
	    }

	  if(iter > MAXITER)
	    {
	      printf("failed to converge in neighbour iteration in density()\n");
	      fflush(stdout);
	      endrun(1155);
	    }
	}
    }
  while(ntot > 0);

  myfree(DataNodeList);
  myfree(DataIndexTable);
  myfree(Ngblist);
  myfree(Right);
  myfree(Left);

#if defined(CS_MODEL) && defined(CS_FEEDBACK)
  double xhyd, yhel, ne, mu, energy, temp, a3inv;

  for(i = FirstActiveParticle; i >= 0; i = NextActiveParticle[i])
    {
      if(P[i].Type == 0 && P[i].EnergySN < 0)	/* only gas particles should enter here */
	{

	  printf("prom 2i=%d type=%d energySN%g\n", i, P[i].Type, P[i].EnergySN);
	  fflush(0);

	  P[i].EnergySN = 0;

	}

      if(All.ComovingIntegrationOn)
	a3inv = 1 / (All.Time * All.Time * All.Time);
      else
	a3inv = 1;

      if(P[i].Type == 0 && (SphP[i].TempPromotion > 0 || SphP[i].DensPromotion > 0))
	{
	  xhyd = P[i].Zm[6] / P[i].Mass;
	  yhel = (1 - xhyd) / (4. * xhyd);
	  ne = SphP[i].Ne;
	  mu = (1 + 4 * yhel) / (1 + yhel + ne);
	  energy = SphP[i].Entropy * P[i].Mass / GAMMA_MINUS1 * pow(SphP[i].d.Density * a3inv, GAMMA_MINUS1);	/* Total Energys */
	  temp = GAMMA_MINUS1 / BOLTZMANN * energy / P[i].Mass * PROTONMASS * mu;
	  temp *= All.UnitEnergy_in_cgs / All.UnitMass_in_g;	/* Temperature in Kelvin */

	  fprintf(FdPromotion, "%g %d %g %g %g %g %g %g\n", All.Time, P[i].ID, SphP[i].DensPromotion,
		  SphP[i].TempPromotion, SphP[i].d.Density, temp, SphP[i].da.DensityAvg,
		  SphP[i].ea.EntropyAvg);
	  fflush(FdPromotion);
	}
    }
#endif

  /* mark as active again */
  for(i = FirstActiveParticle; i >= 0; i = NextActiveParticle[i])
    if(P[i].TimeBin < 0)
      P[i].TimeBin = -P[i].TimeBin - 1;

  /* collect some timing information */

  t1 = WallclockTime = second();
  timeall += timediff(t0, t1);

  timecomp = timecomp1 + timecomp2;
  timewait = timewait1 + timewait2;
  timecomm = timecommsumm1 + timecommsumm2;

  CPU_Step[CPU_DENSCOMPUTE] += timecomp;
  CPU_Step[CPU_DENSWAIT] += timewait;
  CPU_Step[CPU_DENSCOMM] += timecomm;
  CPU_Step[CPU_DENSMISC] += timeall - (timecomp + timewait + timecomm);
}

/*! This function represents the core of the SPH density computation. The
 *  target particle may either be local, or reside in the communication
 *  buffer.
 */
int density_evaluate(int target, int mode, int *nexport, int *nsend_local)
{
  int j, n;
  int startnode, numngb, numngb_inbox, listindex = 0;
  double h, h2, fac, hinv, hinv3, hinv4;
  MyLongDouble rho;
  double wk, dwk;
  double dx, dy, dz, r, r2, u, mass_j;
  double dvx, dvy, dvz;
  MyLongDouble weighted_numngb;
  MyLongDouble dhsmlrho;

#ifdef HYDRO_COST_FACTOR
  int ninteractions = 0;
#endif

#ifdef CS_MODEL
  double densityold, entropy = 0;
#endif
#ifdef BLACK_HOLES
  MyLongDouble gasvel[3];
#endif
#ifndef NAVIERSTOKES
  MyLongDouble divv, rotv[3];
#else
  int k;
  double dvel[3][3];
#endif

  MyDouble *pos;
  MyFloat *vel;
  static MyFloat veldummy[3] = { 0, 0, 0 };

#ifdef EULERPOTENTIALS
  MyDouble deulera[3], deulerb[3], eulera, eulerb, dea, deb;

  deulera[0] = deulera[1] = deulera[2] = deulerb[0] = deulerb[1] = deulerb[2] = 0;
#endif

#if defined(BLACK_HOLES)
  MyLongDouble smoothentr;

  smoothentr = 0;
#endif

#ifdef WINDS
  double delaytime;
#endif

#ifdef VOLUME_CORRECTION
  double densityold, densitystd = 0;
#endif

#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS) || defined(TRACEDIVB)
  MyFloat *bflt;
  MyDouble dbx, dby, dbz;
#endif

#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS)
  MyDouble rotb[3];
#endif

#ifdef TRACEDIVB
  double divB = 0;
#endif

#ifdef CONDUCTION_SATURATION
  double gradentr[3];

  gradentr[0] = gradentr[1] = gradentr[2] = 0;
#endif

#ifdef RADTRANSFER_FLUXLIMITER
  double grad_ngamma[3];

  grad_ngamma[0] = grad_ngamma[1] = grad_ngamma[2] = 0;
#endif

#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS)
  rotb[0] = rotb[1] = rotb[2] = 0;
#endif

#ifndef NAVIERSTOKES
  divv = rotv[0] = rotv[1] = rotv[2] = 0;
#else
  for(k = 0; k < 3; k++)
    dvel[k][0] = dvel[k][1] = dvel[k][2] = 0;
#endif
#ifdef BLACK_HOLES
  gasvel[0] = gasvel[1] = gasvel[2] = 0;
#endif
  rho = weighted_numngb = dhsmlrho = 0;

  if(mode == 0)
    {
      pos = P[target].Pos;
      h = PPP[target].Hsml;
#ifdef VOLUME_CORRECTION
      densityold = SphP[target].DensityOld;
#endif
      if(P[target].Type == 0)
	{
	  vel = SphP[target].VelPred;
#ifdef WINDS
	  delaytime = SphP[target].DelayTime;
#endif
#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS) || defined(TRACEDIVB)
	  bflt = SphP[target].BPred;
#endif
#ifdef EULERPOTENTIALS
	  eulera = SphP[target].EulerA;
	  eulerb = SphP[target].EulerB;
#endif
	}
      else
	vel = veldummy;
#ifdef CS_MODEL
      if(P[target].Type == 0)
	{
	  densityold = SphP[target].DensityOld;
	  entropy = SphP[target].Entropy;
	}
      else
	densityold = 0;
#endif
    }
  else
    {
      pos = DensDataGet[target].Pos;
      vel = DensDataGet[target].Vel;
      h = DensDataGet[target].Hsml;
#ifdef WINDS
      delaytime = DensDataGet[target].DelayTime;
#endif
#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS) || defined(TRACEDIVB)
      bflt = DensDataGet[target].BPred;
#endif
#ifdef VOLUME_CORRECTION
      densityold = DensDataGet[target].DensityOld;
#endif
#ifdef EULERPOTENTIALS
      eulera = DensDataGet[target].EulerA;
      eulerb = DensDataGet[target].EulerB;
#endif
#ifdef CS_MODEL
      densityold = DensDataGet[target].DensityOld;
      entropy = DensDataGet[target].Entropy;
#endif
    }

  h2 = h * h;
  hinv = 1.0 / h;
#ifndef  TWODIMS
  hinv3 = hinv * hinv * hinv;
#else
  hinv3 = hinv * hinv / boxSize_Z;
#endif
  hinv4 = hinv3 * hinv;

  if(mode == 0)
    {
      startnode = All.MaxPart;	/* root node */
    }
  else
    {
      startnode = DensDataGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	/* open it */
    }

  numngb = 0;

  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
#ifdef CS_MODEL
	  numngb_inbox = cs_ngb_treefind_variable_decoupling(&pos[0], h, target, &startnode,
							     densityold, entropy, &vel[0], mode, nexport,
							     nsend_local);
#else
	  numngb_inbox = ngb_treefind_variable(pos, h, target, &startnode, mode, nexport, nsend_local);
#endif

	  if(numngb_inbox < 0)
	    return -1;

	  for(n = 0; n < numngb_inbox; n++)
	    {
#ifdef HYDRO_COST_FACTOR
	      ninteractions++;
#endif
	      j = Ngblist[n];
#ifdef WINDS
	      if(SphP[j].DelayTime > 0)	/* partner is a wind particle */
		if(!(delaytime > 0))	/* if I'm not wind, then ignore the wind particle */
		  continue;
#endif
#ifdef BLACK_HOLES
	      if(P[j].Mass == 0)
		continue;
#endif
	      dx = pos[0] - P[j].Pos[0];
	      dy = pos[1] - P[j].Pos[1];
	      dz = pos[2] - P[j].Pos[2];

#ifdef PERIODIC			/*  now find the closest image in the given box size  */
	      if(dx > boxHalf_X)
		dx -= boxSize_X;
	      if(dx < -boxHalf_X)
		dx += boxSize_X;
	      if(dy > boxHalf_Y)
		dy -= boxSize_Y;
	      if(dy < -boxHalf_Y)
		dy += boxSize_Y;
	      if(dz > boxHalf_Z)
		dz -= boxSize_Z;
	      if(dz < -boxHalf_Z)
		dz += boxSize_Z;
#endif
	      r2 = dx * dx + dy * dy + dz * dz;

	      if(r2 < h2)
		{
		  numngb++;

		  r = sqrt(r2);

		  u = r * hinv;

		  if(u < 0.5)
		    {
		      wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		      dwk = hinv4 * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
		    }
		  else
		    {
		      wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
		      dwk = hinv4 * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
		    }

		  mass_j = P[j].Mass;

#ifdef VOLUME_CORRECTION
		  rho += FLT(mass_j * wk * pow(densityold / SphP[j].DensityOld, VOLUME_CORRECTION));
		  densitystd += FLT(mass_j * wk);
#else
		  rho += FLT(mass_j * wk);
#endif

		  weighted_numngb += FLT(NORM_COEFF * wk / hinv3);	/* 4.0/3 * PI = 4.188790204786 */

		  dhsmlrho += FLT(-mass_j * (NUMDIMS * hinv * wk + u * dwk));

#ifdef BLACK_HOLES
		  smoothentr += FLT(mass_j * wk * SphP[j].Entropy);
		  gasvel[0] += FLT(mass_j * wk * SphP[j].VelPred[0]);
		  gasvel[1] += FLT(mass_j * wk * SphP[j].VelPred[1]);
		  gasvel[2] += FLT(mass_j * wk * SphP[j].VelPred[2]);
#endif

#ifdef CONDUCTION_SATURATION
		  if(r > 0)
		    {
		      gradentr[0] += mass_j * dwk * dx / r * SphP[j].Entropy;
		      gradentr[1] += mass_j * dwk * dy / r * SphP[j].Entropy;
		      gradentr[2] += mass_j * dwk * dz / r * SphP[j].Entropy;
		    }
#endif

#ifdef RADTRANSFER_FLUXLIMITER
		  if(r > 0)
		    {
		      grad_ngamma[0] += mass_j * dwk * dx / r * SphP[j].n_gamma;
		      grad_ngamma[1] += mass_j * dwk * dy / r * SphP[j].n_gamma;
		      grad_ngamma[2] += mass_j * dwk * dz / r * SphP[j].n_gamma;
		    }
#endif

		  if(r > 0)
		    {
		      fac = mass_j * dwk / r;

		      dvx = vel[0] - SphP[j].VelPred[0];
		      dvy = vel[1] - SphP[j].VelPred[1];
		      dvz = vel[2] - SphP[j].VelPred[2];

#ifndef NAVIERSTOKES
		      divv += FLT(-fac * (dx * dvx + dy * dvy + dz * dvz));

		      rotv[0] += FLT(fac * (dz * dvy - dy * dvz));
		      rotv[1] += FLT(fac * (dx * dvz - dz * dvx));
		      rotv[2] += FLT(fac * (dy * dvx - dx * dvy));
#else
		      dvel[0][0] -= fac * dx * dvx;
		      dvel[0][1] -= fac * dx * dvy;
		      dvel[0][2] -= fac * dx * dvz;
		      dvel[1][0] -= fac * dy * dvx;
		      dvel[1][1] -= fac * dy * dvy;
		      dvel[1][2] -= fac * dy * dvz;
		      dvel[2][0] -= fac * dz * dvx;
		      dvel[2][1] -= fac * dz * dvy;
		      dvel[2][2] -= fac * dz * dvz;
#endif
#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS) || defined(TRACEDIVB)
		      dbx = bflt[0] - SphP[j].BPred[0];
		      dby = bflt[1] - SphP[j].BPred[1];
		      dbz = bflt[2] - SphP[j].BPred[2];
#endif
#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS)
		      rotb[0] += FLT(fac * (dz * dby - dy * dbz));
		      rotb[1] += FLT(fac * (dx * dbz - dz * dbx));
		      rotb[2] += FLT(fac * (dy * dbx - dx * dby));
#endif
#ifdef TRACEDIVB
		      divB += fac * (dbx * dx + dby * dy + dbz * dz);
#endif
#ifdef EULERPOTENTIALS
		      dea = eulera - SphP[j].EulerA;
		      deb = eulerb - SphP[j].EulerB;
		      deulera[0] -= fac * dx * dea;
		      deulera[1] -= fac * dy * dea;
		      deulera[2] -= fac * dz * dea;
		      deulerb[0] -= fac * dx * deb;
		      deulerb[1] -= fac * dy * deb;
		      deulerb[2] -= fac * dz * deb;
#endif
		    }
		}
	    }
	}

      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = DensDataGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	/* open it */
	    }
	}
    }

  if(mode == 0)
    {
      PPP[target].n.dNumNgb = weighted_numngb;
#ifdef HYDRO_COST_FACTOR
      P[target].GravCost += HYDRO_COST_FACTOR * ninteractions;
#endif
      if(P[target].Type == 0)
	{
	  SphP[target].d.dDensity = rho;
#ifdef VOLUME_CORRECTION
	  SphP[target].DensityStd = densitystd;
#endif
	  SphP[target].h.dDhsmlDensityFactor = dhsmlrho;
#ifndef NAVIERSTOKES
	  SphP[target].v.dDivVel = divv;
	  SphP[target].r.dRot[0] = rotv[0];
	  SphP[target].r.dRot[1] = rotv[1];
	  SphP[target].r.dRot[2] = rotv[2];
#else
	  for(k = 0; k < 3; k++)
	    {
	      SphP[target].u.DV[k][0] = dvel[k][0];
	      SphP[target].u.DV[k][1] = dvel[k][1];
	      SphP[target].u.DV[k][2] = dvel[k][2];
	    }
#endif

#ifdef CONDUCTION_SATURATION
	  SphP[target].GradEntr[0] = gradentr[0];
	  SphP[target].GradEntr[1] = gradentr[1];
	  SphP[target].GradEntr[2] = gradentr[2];
#endif

#ifdef RADTRANSFER_FLUXLIMITER
	  SphP[target].Grad_ngamma[0] = grad_ngamma[0];
	  SphP[target].Grad_ngamma[1] = grad_ngamma[1];
	  SphP[target].Grad_ngamma[2] = grad_ngamma[2];
#endif

#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS)
	  SphP[target].RotB[0] = rotb[0];
	  SphP[target].RotB[1] = rotb[1];
	  SphP[target].RotB[2] = rotb[2];
#endif
#ifdef TRACEDIVB
	  SphP[target].divB = divB;
#endif
#ifdef EULERPOTENTIALS
	  SphP[target].dEulerA[0] = deulera[0];
	  SphP[target].dEulerA[1] = deulera[1];
	  SphP[target].dEulerA[2] = deulera[2];
	  SphP[target].dEulerB[0] = deulerb[0];
	  SphP[target].dEulerB[1] = deulerb[1];
	  SphP[target].dEulerB[2] = deulerb[2];
#endif
	}
#ifdef BLACK_HOLES
      P[target].b1.dBH_Density = rho;
      P[target].b2.dBH_Entropy = smoothentr;
      P[target].b3.dBH_SurroundingGasVel[0] = gasvel[0];
      P[target].b3.dBH_SurroundingGasVel[1] = gasvel[1];
      P[target].b3.dBH_SurroundingGasVel[2] = gasvel[2];
#endif
#if defined(RADTRANSFER) || defined(SNIA_HEATING)
      if(P[target].Type == 4)
	P[target].DensAroundStar = rho;
#endif
    }
  else
    {
#ifdef HYDRO_COST_FACTOR
      DensDataResult[target].Ninteractions = ninteractions;
#endif
      DensDataResult[target].Rho = rho;
#ifdef VOLUME_CORRECTION
      DensDataResult[target].DensityStd = densitystd;
#endif
      DensDataResult[target].Ngb = weighted_numngb;
      DensDataResult[target].DhsmlDensity = dhsmlrho;
#ifndef NAVIERSTOKES
      DensDataResult[target].Div = divv;
      DensDataResult[target].Rot[0] = rotv[0];
      DensDataResult[target].Rot[1] = rotv[1];
      DensDataResult[target].Rot[2] = rotv[2];
#else
      for(k = 0; k < 3; k++)
	{
	  DensDataResult[target].DV[k][0] = dvel[k][0];
	  DensDataResult[target].DV[k][1] = dvel[k][1];
	  DensDataResult[target].DV[k][2] = dvel[k][2];
	}
#endif

#if defined(BLACK_HOLES)
      DensDataResult[target].SmoothedEntr = smoothentr;
#endif
#ifdef CONDUCTION_SATURATION
      DensDataResult[target].GradEntr[0] = gradentr[0];
      DensDataResult[target].GradEntr[1] = gradentr[1];
      DensDataResult[target].GradEntr[2] = gradentr[2];
#endif

#ifdef RADTRANSFER_FLUXLIMITER
      DensDataResult[target].Grad_ngamma[0] = grad_ngamma[0];
      DensDataResult[target].Grad_ngamma[1] = grad_ngamma[1];
      DensDataResult[target].Grad_ngamma[2] = grad_ngamma[2];
#endif

#if defined(MAGNETIC_DIFFUSION) || defined(ROT_IN_MAG_DIS)
      DensDataResult[target].RotB[0] = rotb[0];
      DensDataResult[target].RotB[1] = rotb[1];
      DensDataResult[target].RotB[2] = rotb[2];
#endif
#ifdef TRACEDIVB
      DensDataResult[target].divB = divB;
#endif
#ifdef BLACK_HOLES
      DensDataResult[target].GasVel[0] = gasvel[0];
      DensDataResult[target].GasVel[1] = gasvel[1];
      DensDataResult[target].GasVel[2] = gasvel[2];
#endif

#ifdef EULERPOTENTIALS
      DensDataResult[target].dEulerA[0] = deulera[0];
      DensDataResult[target].dEulerA[1] = deulera[1];
      DensDataResult[target].dEulerA[2] = deulera[2];
      DensDataResult[target].dEulerB[0] = deulerb[0];
      DensDataResult[target].dEulerB[1] = deulerb[1];
      DensDataResult[target].dEulerB[2] = deulerb[2];
#endif
    }

  return 0;
}

int density_isactive(int n)
{
  if(P[n].TimeBin < 0)
    return 0;

#if defined(RADTRANSFER) || defined(SNIA_HEATING)
  if(P[n].Type == 4)
    return 1;
#endif

#ifdef BLACK_HOLES
  if(P[n].Type == 5)
    return 1;
#endif

  if(P[n].Type == 0)
    return 1;

  return 0;
}

#ifdef NAVIERSTOKES
double get_shear_viscosity(int i)
{
  return All.NavierStokes_ShearViscosity;
}
#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp