From ee0f0e7e8f46bf9f25af2462c3afeac4b8b74c52 Mon Sep 17 00:00:00 2001
From: Simone Rossoni <srossoni@compile-infinix1.physnet.uni-hamburg.de>
Date: Wed, 6 Sep 2023 17:51:02 +0200
Subject: [PATCH] Modified directed emission calculation

---
 include/crpropa/Source.h | 28 ------------------
 src/Source.cpp           | 64 ++--------------------------------------
 2 files changed, 2 insertions(+), 90 deletions(-)

diff --git a/include/crpropa/Source.h b/include/crpropa/Source.h
index b1d8807ee..7cc0978b8 100644
--- a/include/crpropa/Source.h
+++ b/include/crpropa/Source.h
@@ -582,10 +582,6 @@ class SourceIsotropicEmission: public SourceFeature {
 class SourceDirectedEmission: public SourceFeature {
 	Vector3d mu; // Mean emission direction in the vMF distribution
 	double kappa; // Concentration parameter of the vMF distribution
-	double ca; // helpers for the efficient calculation of frame rotation
-	double sa;
-	double cd;
-	double sd;
 public:
 	/** Constructor
 	 @param mu	mean direction of the emission, mu should be normelized
@@ -593,30 +589,6 @@ class SourceDirectedEmission: public SourceFeature {
 	*/
 	SourceDirectedEmission(Vector3d mu, double kappa);
 	void prepareCandidate(Candidate &candidate) const;
-	/**
-	 set sampling parameter Ca
-	 @param alpha	angle between x and y component of direction. alpha = arctan(mu.y / mu.x)
-	*/
-	void setCa(double alpha);
-	/**
-	 set sampling parameter Sa
-	 @param alpha	angle between x and y component of direction. alpha = arctan(mu.y / mu.x)
-	*/
-	void setSa(double alpha);
-	/**
-	 set sampling parameter Cd
-	 @param delta	angle between mu vector and z-axis. delta = arcsin(mu.z) 
-	*/
-	void setCd(double delta);
-	/**
-	 set sampling parameter Sd
-	 @param delta	angle between mu vector and z-axis. delta = arcsin(mu.z) 
-	*/
-	void setSd(double delta);
-	double getCa() const;
-	double getSa() const;
-	double getCd() const;
-	double getSd() const;
 	void setDescription();
 };
 
diff --git a/src/Source.cpp b/src/Source.cpp
index 52cabae88..2ac28a100 100644
--- a/src/Source.cpp
+++ b/src/Source.cpp
@@ -752,38 +752,14 @@ void SourceIsotropicEmission::setDescription() {
 SourceDirectedEmission::SourceDirectedEmission(Vector3d mu, double kappa): mu(mu), kappa(kappa) {
 	if (kappa <= 0)
 		throw std::runtime_error("The concentration parameter kappa should be larger than 0.");
-	double alpha = atan2(mu.y,mu.x);
-	double delta = asin(mu.z);
-	setCa(alpha);
-	setSa(alpha);
-	setCd(delta);
-	setSd(delta);
 	setDescription();
 }
 
 void SourceDirectedEmission::prepareCandidate(Candidate &candidate) const {
 	Random &random = Random::instance();
 
-	//generate sample from von Mises Fisher distribution following
-	// http://people.csail.mit.edu/jstraub/download/straub2017vonMisesFisherInference.pdf
-	//sample normalized direction vector from unit Gaussian distribution
-	Vector3d v(random.randNorm(0., 1.),random.randNorm(0., 1.),0.);
-	v = v.getUnitVector();
-
-	//sample uniform random number
-	double xi = random.rand();
-
-	double u = 1. + 1. / kappa * log(xi + (1. - xi) * exp(-2. * kappa));
-
-	//sample vector from von-Mises distribution
-	Vector3d n = sqrt(1. - u * u) * v;
-	n.z += u;
-
-	//we are in the frame m = (0,0,1)
-	//so rotate to target frame
-	v = Vector3d(ca * sd * n.x - sa * n.y + ca * cd * n.z,
-		sa * sd * n.x + ca * n.y + sa * cd * n.z,
-		- cd * n.x + sd * n.z);
+	Vector3d muvec = mu / mu.getR();
+        Vector3d v = random.randFisherVector(muvec, kappa);
 
 	v = v.getUnitVector();
 	candidate.source.setDirection(v);
@@ -797,42 +773,6 @@ void SourceDirectedEmission::prepareCandidate(Candidate &candidate) const {
 	candidate.setWeight(weight);
 }
 
-void SourceDirectedEmission::setCa(double alpha) {
-	ca = cos(alpha);
-	return;
-}
-
-void SourceDirectedEmission::setSa(double alpha) {
-	sa = sin(alpha);
-	return;
-}
-
-void SourceDirectedEmission::setCd(double delta) {
-	cd = cos(delta);
-	return;
-}
-
-void SourceDirectedEmission::setSd(double delta) {
-	sd = sin(delta);
-	return;
-}
-
-double SourceDirectedEmission::getCa() const {
-	return ca;
-}    
-
-double SourceDirectedEmission::getSa() const {
-	return sa;
-} 
-
-double SourceDirectedEmission::getCd() const {
-	return cd;
-} 
-
-double SourceDirectedEmission::getSd() const {
-	return sd;
-} 
-
 void SourceDirectedEmission::setDescription() {
 	std::stringstream ss;
 	ss << "SourceDirectedEmission: Random directed emission following the von-Mises-Fisher distribution with mean direction ";