Added Lopez-Corredoira et al. (2005) bulge model

CMakeLists.txt

@@ -110,6 +110,7 @@ cuda_add_executable(galfast
 
   src/skygen/model_expDisk.cpp
   src/skygen/model_powerLawEllipsoid.cpp
+  src/skygen/model_LCBulge.cpp
   src/skygen/os_skygen.cpp
   src/skygen/model_brokenPowerLaw.cpp
   src/skygen/model_densityCube.cpp

src/pipeline.cpp

@@ -568,6 +568,7 @@ bool os_countsMap::construct(const Config &cfg, otable &t, opipeline &pipe)
 	cfg.get(z0, "z0", z0);
 	cfg.get(z1, "z1", z1);
 	cfg.get(dz, "dz", dz);
+
 	tmp = (z1 - z0) / dz;
 	n_z = (int)(fabs(tmp - round(tmp)) < 1e-4 ? round(tmp) : ceil(tmp)) + 1;
 

src/skygen/model_LCBulge.cpp

@@ -0,0 +1,83 @@
+/***************************************************************************
+ *   Copyright (C) 2004 by Mario Juric                                     *
+ *   [email protected]                                            *
+ *                                                                         *
+ *   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.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#include "galfast_config.h"
+
+#include "model_LCBulge.h"
+#include "skyconfig_impl.h"
+#include "model_lib.h"
+#include "transform.h"
+
+#include <astro/math.h>
+
+#include <astro/useall.h>
+#include <fstream>
+
+void LCBulge::load(host_state_t &hstate, const peyton::system::Config &cfg)
+{
+	// component ID
+	int userComp = cfg.get("comp");
+	comp = componentMap.seqIdx(userComp);
+
+	// ellipsoid orientation
+	std::string ctr, orient;
+	cfg.get(ctr,    "center",      "galplane");
+	cfg.get(orient, "orientation", "galplane");
+	load_transform(&T.x, rot, ctr, orient, cfg);
+
+	// exponential model parameters
+	n = cfg.get("n");	// power law index of the exponent
+	l = cfg.get("l");	// exponential scale length
+	ba = cfg.get("b_a");	// b / a axis ratio
+	ca = cfg.get("c_a");	// c / a axis ratio
+
+	// Load and renormalize the luminosity function
+	float f;
+	cfg.get(f, "f", 1.f);
+	std::string lffile;
+	hstate.lf = load_lf_raw(*this, cfg, lffile);
+	FOREACH(hstate.lf) { *i *= f; }
+
+	//
+	if(lffile.empty()) { lffile = "-"; }
+	MLOG(verb1) << "Component " << userComp << " : " << " LCbulge {" << l << ", " << n << ", " << ca << ", " << ba << ", " << ", " << f << ", " << lffile << "}";
+}
+
+void LCBulge::prerun(host_state_t &hstate, bool draw)
+{
+	// bind the luminosity function texture to texture reference
+	LCBulgeLF.bind(hstate.lf);
+}
+
+void LCBulge::postrun(host_state_t &hstate, bool draw)
+{
+	// unbind LF texture reference
+	LCBulgeLF.unbind();
+}
+
+bool LCBulge::hint_absmag(host_state_t &hstate, float &M0, float &M1) const
+{
+	return absmag_adjust_from_lf(hstate.lf, M0, M1);
+}
+
+extern "C" skygenInterface *create_model_lcbulge()
+{
+	return new skygenHost<LCBulge>();
+}

src/skygen/model_LCBulge.h

@@ -0,0 +1,98 @@
+/***************************************************************************
+ *   Copyright (C) 2004 by Mario Juric                                     *
+ *   [email protected]                                            *
+ *                                                                         *
+ *   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.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#ifndef LCBulge_h__
+#define LCBulge_h__
+
+#include "skygen.h"
+
+#include "model_powerLawEllipsoid.h"
+
+// luminosity function texture reference
+DEFINE_TEXTURE(LCBulgeLF, float, 1, cudaReadModeElementType, false, cudaFilterModeLinear, cudaAddressModeClamp);
+
+// exponential disk model
+struct ALIGN(16) LCBulge : public modelConcept
+{
+public:
+	struct ALIGN(16) host_state_t
+	{
+		cuxTexture<float> lf;
+	};
+
+	struct state
+	{
+		float rho;
+	};
+
+public:
+	// Exponential ellipsoid parameters for
+	//
+	//	rho = LF(M) * f * exp(-t/l)
+	//	l = ( x^n + (y/ba)^n + (y/ca)^n )^(1/n)
+	//
+	// model. Note that LF(M) * f is precomputed on load.
+
+	float ba, ca;		// ellipsoid axes ratios
+	float n;		// power law index of the exponent
+	float l;		// exponential scale
+
+	// Positioning
+	float3 T;		// Translation parameter for transform() (negative distance to ellipsoid center, in Galactic coordinates)
+	float rot[3][3];	// Rotation matrix of the ellipsoid (transforms points from Galactic to ellipsoid coordinate system)
+
+public:
+	void load(host_state_t &hstate, const peyton::system::Config &cfg);
+	void prerun(host_state_t &hstate, bool draw);
+	void postrun(host_state_t &hstate, bool draw);
+	bool hint_absmag(host_state_t &hstate, float &M0, float &M1) const;
+
+public:
+	__host__ __device__ float rho(float x, float y, float z) const
+	{
+		using namespace cudacc;
+
+		// translate and rotate into ellipsoid-centric frame
+		float3 v = {x, y, z};
+		v = transform(v, T, rot);
+
+		// compute the exponent and density
+		float t = powf(powf(v.x, n) + powf(v.y / ba, n) + powf(v.z / ca, n), 1./n);
+		float rho = expf(-t / l);
+
+		return rho;
+	}
+
+public:
+	__device__ void setpos(state &s, float x, float y, float z) const
+	{
+		s.rho = rho(x, y, z);
+	}
+
+	__device__ float rho(state &s, float M) const
+	{
+		float phi = TEX1D(LCBulgeLF, M);
+		return phi * s.rho;
+	}
+};
+
+MODEL_IMPLEMENTATION(LCBulge);
+
+#endif // #ifndef LCBulge_h__

src/skygen/model_lib.h

@@ -27,20 +27,15 @@
 #include <astro/math.h>
 
 //
-// Loads the luminosity function and computes and sets the normalization based
-// on the coordinates where the LF was sampled.
+// Loads the luminosity function as-is, with no renormalization
 //
-// LFSeparableModel concept requirements:
+// SimpleLFSeparableModel concept requirements:
 //	*) \rho(x,y,z,M) density must be separable and computed as
 //	   \rho(x,y,z) * LF(M)
 //      *) \rho(x,y,z) must exist
-//      *) A non-const member LFSeparableModel::f must exist and the
-//         output of \rho(x,y,z) must scale linearly with f. This variable
-//         will be rescaled so that \rho(x,y,z) = f_user (where f_user is
-//         the "f" value set by the user in the config file).
 //
-template<typename LFSeparableModel>
-cuxTexture<float> load_lf(LFSeparableModel &m, const peyton::system::Config &cfg, std::string &lffile)
+template<typename SimpleLFSeparableModel>
+cuxTexture<float> load_lf_raw(SimpleLFSeparableModel &m, const peyton::system::Config &cfg, std::string &lffile)
 {
 	cuxTexture<float> lf;
 
@@ -56,7 +51,27 @@ cuxTexture<float> load_lf(LFSeparableModel &m, const peyton::system::Config &cfg
 		lf = load_constant_texture_1D(1.f, -100., 100.);
 	}
 
-	// location where the given luminosity function was sampled --
+	return lf;
+}
+
+//
+// Loads the luminosity function and computes and sets the normalization based
+// on the coordinates where the LF was sampled.
+//
+// LFSeparableModel concept requirements:
+//	*) Inherits SimpleLFSeparableModel's requirements
+//      *) A non-const member LFSeparableModel::f must exist and the
+//         output of \rho(x,y,z) must scale linearly with f. This variable
+//         will be rescaled so that \rho(x,y,z) = f_user (where f_user is
+//         the "f" value set by the user in the config file).
+//
+template<typename LFSeparableModel>
+cuxTexture<float> load_lf(LFSeparableModel &m, const peyton::system::Config &cfg, std::string &lffile)
+{
+	// load the raw LF
+	cuxTexture<float> lf = load_lf_raw(m, cfg, lffile);
+
+	// location where the luminosity function was sampled --
 	// default is the Solar neighborhood
 	float3 lfc;
 	cfg.get(lfc.x, "lf_l", 0.f);	// Galactic longitude

src/skygen/skygen.cu.h

@@ -662,3 +662,4 @@ void skygenHost<T>::compute(bool draw)
 #include "model_densityCube.h"
 #include "model_powerLawEllipsoid.h"
 #include "model_brokenPowerLaw.h"
+#include "model_LCBulge.h"