Merge branch 'feature/fieldlines' into feature/ABuffer

include/openspace/rendering/renderablefieldlines.h

@@ -52,14 +52,14 @@ private:
 	std::vector<std::string> _filenames;
 	std::vector<glm::vec3> _seedPoints;
 
-	ghoul::opengl::ProgramObject *_fieldlinesProgram, *_seedpointsProgram;
+	ghoul::opengl::ProgramObject *_fieldlinesProgram;
 	GLuint _VAO, _seedpointVAO;
 
 	ghoul::filesystem::File* _vertexSourceFile;
 	ghoul::filesystem::File* _fragmentSourceFile;
 
-	std::vector<GLint> _lineStart, _seedpointStart;
-	std::vector<GLsizei> _lineCount, _seedpointCount;
+	std::vector<GLint> _lineStart;
+	std::vector<GLsizei> _lineCount;
 
 	bool _programUpdateOnSave;
 	bool _update;

include/openspace/util/kameleonwrapper.h

@@ -25,7 +25,6 @@
 #ifndef KAMELEONWRAPPER_H_
 #define KAMELEONWRAPPER_H_
 
-//#include <glm/glm.hpp>
 #include <glm/gtx/std_based_type.hpp>
 
 namespace ccmc {
@@ -43,12 +42,12 @@ public:
 		BATSRUS		// Magnetosphere
 	};
 
-	enum TraceDirection {
+	enum class TraceDirection {
 		FORWARD = 1,
 		BACK 	= -1
 	};
 
-	enum FieldlineEnd {
+	enum class FieldlineEnd {
 		NORTH,
 		SOUTH,
 		OUT
@@ -60,15 +59,25 @@ public:
 	float* getUniformSampledVectorValues(const std::string& xVar, const std::string& yVar,
 			const std::string& zVar, glm::size3_t outDimensions);
 
-	std::vector<std::vector<glm::vec3> > getFieldLines(const std::string& xVar,
+	std::vector<std::vector<glm::vec3> > getClassifiedFieldLines(const std::string& xVar,
 			const std::string& yVar, const std::string& zVar,
 			std::vector<glm::vec3> seedPoints, float stepSize);
 
+	std::vector<std::vector<glm::vec3> > getFieldLines(const std::string& xVar,
+				const std::string& yVar, const std::string& zVar,
+				std::vector<glm::vec3> seedPoints, float stepSize, glm::vec3 color);
+
+	std::vector<std::vector<glm::vec3> > getLorentzTrajectories(std::vector<glm::vec3> seedPoints,
+			glm::vec3 color, float stepsize);
+
 private:
 	std::vector<glm::vec3> traceCartesianFieldline(const std::string& xVar,
 			const std::string& yVar, const std::string& zVar, glm::vec3 seedPoint,
 			float stepSize, TraceDirection direction, FieldlineEnd& end);
 
+	std::vector<glm::vec3> traceLorentzTrajectory(glm::vec3 seedPoint,
+			float stepsize, float eCharge);
+
 	void getGridVariables(std::string& x, std::string& y, std::string& z);
 	void progressBar(int current, int end);
 	glm::vec3 classifyFieldline(FieldlineEnd fEnd, FieldlineEnd bEnd);

src/rendering/renderablefieldlines.cpp

@@ -92,15 +92,7 @@ RenderableFieldlines::RenderableFieldlines(const ghoul::Dictionary& dictionary)
 			_fieldlinesProgram->attachObject(fragmentShader);
 		}
 	}
-/*
-	_seedpointsProgram = new ghoul::opengl::ProgramObject("SeedpointsProgram");
-	ghoul::opengl::ShaderObject* seedpointVertexShader = new ghoul::opengl::ShaderObject(ghoul::opengl::ShaderObject::ShaderTypeVertex,
-			"/home/hhellteg/openspace/openspace-data/scene/fieldlines/seedPoints.vert");
-	ghoul::opengl::ShaderObject* seedpointFragmentShader = new ghoul::opengl::ShaderObject(ghoul::opengl::ShaderObject::ShaderTypeFragment,
-			"/home/hhellteg/openspace/openspace-data/scene/fieldlines/seedPoints.frag");
-	_seedpointsProgram->attachObject(seedpointVertexShader);
-	_seedpointsProgram->attachObject(seedpointFragmentShader);
-*/
+
 	if(dictionary.hasKey("UpdateOnSave")) {
 		dictionary.getValue("UpdateOnSave", _programUpdateOnSave);
 	}
@@ -128,14 +120,17 @@ bool RenderableFieldlines::initialize() {
 			_lineCount.push_back(fieldlinesData[j].size()/2.0);
 			prevEnd = prevEnd + fieldlinesData[j].size()/2.0;
 
-			_seedpointStart.push_back(j);
-			_seedpointCount.push_back(1);
-
 			vertexData.insert( vertexData.end(), fieldlinesData[j].begin(), fieldlinesData[j].end());
 		}
 	}
+	LDEBUG("Number of vertices : " << vertexData.size()/2.0);
 
-	LDEBUG("Vertex orginal : " << vertexData.size()/2.0);
+	// Give seedpoints a color for visualizing as GL_POINTS
+	std::vector<glm::vec3> seedPointsData;
+	for (int i = 0; i < _seedPoints.size(); ++i) {
+		seedPointsData.push_back(_seedPoints[i]);
+		seedPointsData.push_back(glm::vec3(1.0, 0.5, 0.0));
+	}
 
 	//	------ FIELDLINES -----------------
 	GLuint vertexPositionBuffer;
@@ -158,27 +153,31 @@ bool RenderableFieldlines::initialize() {
 	glBindBuffer(GL_ARRAY_BUFFER, 0); //unbind buffer
 	glBindVertexArray(0); //unbind array
 
-//	//	------ SEEDPOINTS -----------------
+	//	------ SEEDPOINTS -----------------
 	GLuint seedpointPositionBuffer;
 	glGenVertexArrays(1, &_seedpointVAO); // generate array
 	glBindVertexArray(_seedpointVAO); // bind array
 	glGenBuffers(1, &seedpointPositionBuffer); // generate buffer
 	glBindBuffer(GL_ARRAY_BUFFER, seedpointPositionBuffer); // bind buffer
-	glBufferData(GL_ARRAY_BUFFER, _seedPoints.size()*sizeof(glm::vec3), &_seedPoints.front(), GL_STATIC_DRAW);
+	glBufferData(GL_ARRAY_BUFFER, seedPointsData.size()*sizeof(glm::vec3), &seedPointsData.front(), GL_STATIC_DRAW);
 
 	// Vertex positions
-	GLuint seedpointLocation = 0;
-	glEnableVertexAttribArray(seedpointLocation);
-	glVertexAttribPointer(seedpointLocation, 3, GL_FLOAT, GL_FALSE, sizeof(glm::vec3), reinterpret_cast<void*>(0));
+	glEnableVertexAttribArray(vertexLocation);
+	glVertexAttribPointer(vertexLocation, 3, GL_FLOAT, GL_FALSE, 2*sizeof(glm::vec3), reinterpret_cast<void*>(0));
+
+	// Texture coordinates
+	glEnableVertexAttribArray(texcoordLocation);
+	glVertexAttribPointer(texcoordLocation, 3, GL_FLOAT, GL_FALSE, 2*sizeof(glm::vec3), (void*)(sizeof(glm::vec3)));
 
 	glBindBuffer(GL_ARRAY_BUFFER, 0); //unbind buffer
 	glBindVertexArray(0); //unbind array
 
+	glPointSize(5); // size of seedpoints
+
 	//	------ SETUP SHADERS -----------------
 	auto privateCallback = [this](const ghoul::filesystem::File& file) {
 		_update = true;
 	};
-
 	if(_programUpdateOnSave) {
 		_vertexSourceFile->setCallback(privateCallback);
 		_fragmentSourceFile->setCallback(privateCallback);
@@ -208,8 +207,8 @@ void RenderableFieldlines::render(const Camera* camera,	const psc& thisPosition)
 	transform = transform*camTransform;
 	//transform = glm::translate(transform, relative.vec3());
 	transform = glm::mat4(1.0);
-	transform = glm::rotate(transform, -90.0f, glm::vec3(1.0, 0.0, 0.0)); // Model has positive Z as up
-	transform = glm::scale(transform, glm::vec3(0.036*0.5*0.5));
+//	transform = glm::scale(transform, glm::vec3(0.036*0.5*0.5));
+	transform = glm::scale(transform, glm::vec3(0.01));
 	//transform = glm::scale(transform, glm::vec3(0.1)); // Scale to avoid depth buffer problems
 
 	psc currentPosition = thisPosition;
@@ -231,14 +230,11 @@ void RenderableFieldlines::render(const Camera* camera,	const psc& thisPosition)
 	//	------ FIELDLINES -----------------
 	glBindVertexArray(_VAO);
 	glMultiDrawArrays(GL_LINE_STRIP, &_lineStart[0], &_lineCount[0], _lineStart.size());
-	glBindVertexArray(0);
 
 	//	------ SEEDPOINTS -----------------
-	// glBindVertexArray(_seedpointVAO);
-	// glPointSize(5);
-	// glMultiDrawArrays(GL_POINTS, &_lineStart[0], &_lineCount[0], _seedPoints.size());
-	// glBindVertexArray(0);
-
+	glBindVertexArray(_seedpointVAO);
+	glMultiDrawArrays(GL_POINTS, &_lineStart[0], &_lineCount[0], _seedPoints.size());
+	glBindVertexArray(0);
 
 	// Deactivate shader
 	_fieldlinesProgram->deactivate();
@@ -260,7 +256,11 @@ std::vector<std::vector<glm::vec3> > RenderableFieldlines::getFieldlinesData(std
 	KameleonWrapper::Model model;
 	std::vector<std::vector<glm::vec3> > fieldlinesData;
 
+	bool classification = false, lorentz = false;
+	glm::vec3 fieldlineColor = glm::vec3(1.0, 1.0, 1.0); // default color if no color or classification is specified
+
 	if (hintsDictionary.hasKey("Model") && hintsDictionary.getValue("Model", modelString)) {
+		//	------ MODEL -----------------
 		if (modelString == "BATSRUS") {
 			model = KameleonWrapper::Model::BATSRUS;
 		} else if (modelString == "ENLIL") {
@@ -271,25 +271,35 @@ std::vector<std::vector<glm::vec3> > RenderableFieldlines::getFieldlinesData(std
 			return fieldlinesData;
 		}
 
+		//	------ VARIBLES / LORENTZ -----------------
 		if (hintsDictionary.hasKey("Variables")) {
 			bool xVar, yVar, zVar;
 			xVar = hintsDictionary.getValue("Variables.1", xVariable);
-			yVar = hintsDictionary.getValue("Variables.2", yVariable);
-			zVar = hintsDictionary.getValue("Variables.3", zVariable);
+			if (xVar && xVariable == "Lorentz") {
+				lorentz = true;
+			} else {
 
-			if (!xVar || !yVar || !zVar) {
-				LWARNING("Error reading variables! Must be 3 and must exist in CDF data");
-				return fieldlinesData;
+				yVar = hintsDictionary.getValue("Variables.2", yVariable);
+				zVar = hintsDictionary.getValue("Variables.3", zVariable);
+
+				if (!xVar || !yVar || !zVar) {
+					LWARNING("Error reading variables! Must be 3 and must exist in CDF data");
+					return fieldlinesData;
+				}
 			}
 		} else {
 			LWARNING("Hints does not specify  valid 'Variables'");
 			return fieldlinesData;
 		}
 
+		//	------ STEPSIZE -----------------
 		if (!hintsDictionary.hasKey("Stepsize") || !hintsDictionary.getValue("Stepsize", stepSize)) {
 			LDEBUG("No stepsize set for fieldlines. Setting to default value (" << stepSize << ")");
 		}
+
+		//	------ SEEDPOINTS ---------------
 		ghoul::Dictionary seedpointsDictionary;
+		_seedPoints.clear();
 		if (hintsDictionary.hasKey("Seedpoints") && hintsDictionary.getValue("Seedpoints", seedpointsDictionary)) {
 			glm::vec3 seedPos;
 			for (auto index : seedpointsDictionary.keys()) {
@@ -297,10 +307,22 @@ std::vector<std::vector<glm::vec3> > RenderableFieldlines::getFieldlinesData(std
 				_seedPoints.push_back(seedPos);
 			}
 		}
+
+		//	------ CLASSIFICATION & COLOR -----------
+		hintsDictionary.getValue("Color", fieldlineColor);
+		hintsDictionary.getValue("Classification", classification);
 	}
 
 	KameleonWrapper kw(filename, model);
-	fieldlinesData = kw.getFieldLines(xVariable, yVariable, zVariable, _seedPoints, stepSize);
+	if (lorentz) {
+		fieldlinesData = kw.getLorentzTrajectories(_seedPoints, fieldlineColor, stepSize);
+	} else {
+		if (classification)
+			fieldlinesData = kw.getClassifiedFieldLines(xVariable, yVariable, zVariable, _seedPoints, stepSize);
+		else
+			fieldlinesData = kw.getFieldLines(xVariable, yVariable, zVariable, _seedPoints, stepSize, fieldlineColor);
+	}
+
 	return fieldlinesData;
 }
 

src/util/kameleonwrapper.cpp

@@ -36,6 +36,8 @@
 #include <iomanip>
 #include <stdlib.h>
 
+#include <glm/gtx/rotate_vector.hpp>
+
 namespace openspace {
 
 std::string _loggerCat = "KameleonWrapper";
@@ -122,7 +124,8 @@ float* KameleonWrapper::getUniformSampledValues(const std::string& var, glm::siz
                     double zPos = _zMin + stepZ*z;
                     
                     // get interpolated data value for (xPos, yPos, zPos)
-                    double value = _interpolator->interpolate(var, xPos, yPos, zPos);
+                    // swap yPos and zPos because model has Z as up
+                    double value = _interpolator->interpolate(var, xPos, zPos, yPos);
                     
                     // scale to [0,1]
                     //doubleData[index] = (value-varMin)/(varMax-varMin);
@@ -266,7 +269,7 @@ float* KameleonWrapper::getUniformSampledVectorValues(const std::string& xVar, c
 	return data;
 }
 
-std::vector<std::vector<glm::vec3> > KameleonWrapper::getFieldLines(
+std::vector<std::vector<glm::vec3> > KameleonWrapper::getClassifiedFieldLines(
 		const std::string& xVar, const std::string& yVar,
 		const std::string& zVar, std::vector<glm::vec3> seedPoints,
 		float stepSize ) {
@@ -305,6 +308,66 @@ std::vector<std::vector<glm::vec3> > KameleonWrapper::getFieldLines(
 	return fieldLines;
 }
 
+std::vector<std::vector<glm::vec3> > KameleonWrapper::getFieldLines(
+		const std::string& xVar, const std::string& yVar,
+		const std::string& zVar, std::vector<glm::vec3> seedPoints,
+		float stepSize, glm::vec3 color ) {
+	assert(_model && _interpolator);
+	assert(_type == Model::ENLIL || _type == Model::BATSRUS);
+	LINFO("Creating " << seedPoints.size() << " fieldlines from variables " << xVar << " " << yVar << " " << zVar);
+
+	std::vector<glm::vec3> fLine, bLine;
+	std::vector<std::vector<glm::vec3> > fieldLines;
+	FieldlineEnd forwardEnd, backEnd;
+
+	if (_type == Model::BATSRUS) {
+		for (glm::vec3 seedPoint : seedPoints) {
+			fLine = traceCartesianFieldline(xVar, yVar, zVar, seedPoint, stepSize, TraceDirection::FORWARD, forwardEnd);
+			bLine = traceCartesianFieldline(xVar, yVar, zVar, seedPoint, stepSize, TraceDirection::BACK, backEnd);
+
+			bLine.insert(bLine.begin(), fLine.rbegin(), fLine.rend());
+
+			// write colors
+			std::vector<glm::vec3> line;
+			for (glm::vec3 position : bLine) {
+				line.push_back(position);
+				line.push_back(color);
+			}
+
+			fieldLines.push_back(line);
+		}
+	} else {
+		LERROR("Fieldlines are only supported for BATSRUS model");
+	}
+
+	return fieldLines;
+}
+
+std::vector<std::vector<glm::vec3> > KameleonWrapper::getLorentzTrajectories(
+		std::vector<glm::vec3> seedPoints, glm::vec3 color, float stepsize) {
+	LINFO("Creating " << seedPoints.size() << " Lorentz force trajectories");
+
+	std::vector<std::vector<glm::vec3> > trajectories;
+	std::vector<glm::vec3> plusTraj, minusTraj;
+
+	for (auto seedPoint : seedPoints) {
+		plusTraj = traceLorentzTrajectory(seedPoint, stepsize, 1.0);
+		minusTraj = traceLorentzTrajectory(seedPoint, stepsize, -1.0);
+
+		minusTraj.insert(minusTraj.begin(), plusTraj.rbegin(), plusTraj.rend());
+
+		// write colors
+		std::vector<glm::vec3> trajectory;
+		for (glm::vec3 position : minusTraj) {
+			trajectory.push_back(position);
+			trajectory.push_back(color);
+		}
+		trajectories.push_back(trajectory);
+	}
+
+	return trajectories;
+}
+
 std::vector<glm::vec3> KameleonWrapper::traceCartesianFieldline(
 		const std::string& xVar, const std::string& yVar,
 		const std::string& zVar, glm::vec3 seedPoint,
@@ -313,8 +376,7 @@ std::vector<glm::vec3> KameleonWrapper::traceCartesianFieldline(
 	glm::vec3 color, pos, k1, k2, k3, k4;
 	std::vector<glm::vec3> line;
 	float stepX, stepY, stepZ;
-	int numSteps = 0;
-	int maxSteps = 5000;
+	int numSteps = 0, maxSteps = 5000;
 	pos = seedPoint;
 
 	_model->loadVariable(xVar);
@@ -329,15 +391,15 @@ std::vector<glm::vec3> KameleonWrapper::traceCartesianFieldline(
 	while ((pos.x < _xMax && pos.x > _xMin && pos.y < _yMax && pos.y > _yMin &&
 			pos.z < _zMax && pos.z > _zMin) && !(pos.x*pos.x + pos.y*pos.y + pos.z*pos.z < 1.0)) {
 
-		// Save position
-		line.push_back(pos);
+		// Save position. Model has +Z as up
+		line.push_back(glm::vec3(pos.x, pos.z, pos.y));
 
 		// Calculate new position with Runge-Kutta 4th order
 		k1.x = _interpolator->interpolate(xID, pos.x, pos.y, pos.z, stepX, stepY, stepZ);
 		k1.y = _interpolator->interpolate(yID, pos.x, pos.y, pos.z);
 		k1.z = _interpolator->interpolate(zID, pos.x, pos.y, pos.z);
 		k1 = (float)direction*glm::normalize(k1);
-//		stepX*=stepSize, stepY*=stepSize, stepZ*=stepSize;
+		stepX=stepX*stepSize, stepY=stepY*stepSize, stepZ=stepZ*stepSize;
 		k2.x = _interpolator->interpolate(xID, pos.x+(stepX/2.0)*k1.x, pos.y+(stepY/2.0)*k1.y, pos.z+(stepZ/2.0)*k1.z);
 		k2.y = _interpolator->interpolate(yID, pos.x+(stepX/2.0)*k1.x, pos.y+(stepY/2.0)*k1.y, pos.z+(stepZ/2.0)*k1.z);
 		k2.z = _interpolator->interpolate(zID, pos.x+(stepX/2.0)*k1.x, pos.y+(stepY/2.0)*k1.y, pos.z+(stepZ/2.0)*k1.z);
@@ -360,7 +422,8 @@ std::vector<glm::vec3> KameleonWrapper::traceCartesianFieldline(
 			break;
 		}
 	}
-	line.push_back(pos);
+	// Save last position. Model has +Z as up
+	line.push_back(glm::vec3(pos.x, pos.z, pos.y));
 
 	if (pos.z > 0.0 && (pos.x*pos.x + pos.y*pos.y + pos.z*pos.z < 1.0))
 		end = FieldlineEnd::NORTH;
@@ -372,6 +435,98 @@ std::vector<glm::vec3> KameleonWrapper::traceCartesianFieldline(
 	return line;
 }
 
+std::vector<glm::vec3> KameleonWrapper::traceLorentzTrajectory(glm::vec3 seedPoint,
+		float stepsize, float eCharge) {
+	glm::vec3 B, E, v0, k1, k2, k3, k4, sPos, tmpV;
+	float stepX = stepsize, stepY = stepsize, stepZ = stepsize;
+
+	long int bxID = _model->getVariableID("bx");
+	long int byID = _model->getVariableID("by");
+	long int bzID = _model->getVariableID("bz");
+	long int jxID = _model->getVariableID("jx");
+	long int jyID = _model->getVariableID("jy");
+	long int jzID = _model->getVariableID("jz");
+
+	std::vector<glm::vec3> trajectory;
+	glm::vec3 pos = seedPoint;
+	int numSteps = 0, maxSteps = 5000;
+	v0.x = _interpolator->interpolate("ux", pos.x, pos.y, pos.z);
+	v0.y = _interpolator->interpolate("uy", pos.x, pos.y, pos.z);
+	v0.z = _interpolator->interpolate("uz", pos.x, pos.y, pos.z);
+	v0 = glm::normalize(v0);
+
+	// While we are inside the models boundries and not inside earth
+	while ((pos.x < _xMax && pos.x > _xMin && pos.y < _yMax && pos.y > _yMin &&
+			pos.z < _zMax && pos.z > _zMin) && !(pos.x*pos.x + pos.y*pos.y + pos.z*pos.z < 1.0)) {
+
+		// Save position. Model has +Z as up
+		trajectory.push_back(glm::vec3(pos.x, pos.z, pos.y));
+
+		// Calculate new position with Lorentz force quation and Runge-Kutta 4th order
+		B.x = _interpolator->interpolate(bxID, pos.x, pos.y, pos.z);
+		B.y = _interpolator->interpolate(byID, pos.x, pos.y, pos.z);
+		B.z = _interpolator->interpolate(bzID, pos.x, pos.y, pos.z);
+		E.x = _interpolator->interpolate(jxID, pos.x, pos.y, pos.z);
+		E.y = _interpolator->interpolate(jyID, pos.x, pos.y, pos.z);
+		E.z = _interpolator->interpolate(jzID, pos.x, pos.y, pos.z);
+		k1 = eCharge*(E + glm::cross(v0, B));
+		k1 = glm::normalize(k1);
+
+		sPos = glm::vec3(	pos.x+(stepX/2.0)*v0.x+(stepX*stepX/8.0)*k1.x,
+							pos.y+(stepY/2.0)*v0.y+(stepY*stepY/8.0)*k1.y,
+							pos.z+(stepZ/2.0)*v0.z+(stepZ*stepZ/8.0)*k1.z);
+		B.x = _interpolator->interpolate(bxID, sPos.x, sPos.y, sPos.z);
+		B.y = _interpolator->interpolate(byID, sPos.x, sPos.y, sPos.z);
+		B.z = _interpolator->interpolate(bzID, sPos.x, sPos.y, sPos.z);
+		E.x = _interpolator->interpolate(jxID, sPos.x, sPos.y, sPos.z);
+		E.y = _interpolator->interpolate(jyID, sPos.x, sPos.y, sPos.z);
+		E.z = _interpolator->interpolate(jzID, sPos.x, sPos.y, sPos.z);
+		tmpV = v0+(stepX/2.0f)*k1;
+		k2 = eCharge*(E + glm::cross(tmpV, B));
+		k2 = glm::normalize(k2);
+
+		B.x = _interpolator->interpolate(bxID, sPos.x, sPos.y, sPos.z);
+		B.y = _interpolator->interpolate(byID, sPos.x, sPos.y, sPos.z);
+		B.z = _interpolator->interpolate(bzID, sPos.x, sPos.y, sPos.z);
+		E.x = _interpolator->interpolate(jxID, sPos.x, sPos.y, sPos.z);
+		E.y = _interpolator->interpolate(jyID, sPos.x, sPos.y, sPos.z);
+		E.z = _interpolator->interpolate(jzID, sPos.x, sPos.y, sPos.z);
+		tmpV = v0+(stepX/2.0f)*k2;
+		k3 = eCharge*(E + glm::cross(tmpV, B));
+		k3 = glm::normalize(k3);
+
+		sPos = glm::vec3(	pos.x+stepX*v0.x+(stepX*stepX/2.0)*k1.x,
+							pos.y+stepY*v0.y+(stepY*stepY/2.0)*k1.y,
+							pos.z+stepZ*v0.z+(stepZ*stepZ/2.0)*k1.z);
+		B.x = _interpolator->interpolate(bxID, sPos.x, sPos.y, sPos.z);
+		B.y = _interpolator->interpolate(byID, sPos.x, sPos.y, sPos.z);
+		B.z = _interpolator->interpolate(bzID, sPos.x, sPos.y, sPos.z);
+		E.x = _interpolator->interpolate(jxID, sPos.x, sPos.y, sPos.z);
+		E.y = _interpolator->interpolate(jyID, sPos.x, sPos.y, sPos.z);
+		E.z = _interpolator->interpolate(jzID, sPos.x, sPos.y, sPos.z);
+		tmpV = v0+stepX*k3;
+		k4 = eCharge*(E + glm::cross(tmpV, B));
+		k4 = glm::normalize(k4);
+
+		pos.x = pos.x + stepX*v0.x + (stepX*stepX/6.0)*(k1.x + k2.x + k3.x);
+		pos.y = pos.y + stepY*v0.y + (stepY*stepY/6.0)*(k1.y + k2.y + k3.y);
+		pos.z = pos.z + stepZ*v0.z + (stepZ*stepZ/6.0)*(k1.z + k2.z + k3.z);
+
+		v0.x = v0.x + (stepX/6.0)*(k1.x + 2.0*k2.x + 2.0*k3.x + k4.z);
+		v0.y = v0.y + (stepY/6.0)*(k1.y + 2.0*k2.y + 2.0*k3.y + k4.y);
+		v0.z = v0.z + (stepZ/6.0)*(k1.z + 2.0*k2.z + 2.0*k3.z + k4.z);
+
+		++numSteps;
+		if (numSteps > maxSteps) {
+			LDEBUG("Max number of steps taken (" << maxSteps <<")");
+			break;
+		}
+	}
+	// Save last position. Model has +Z as up
+	trajectory.push_back(glm::vec3(pos.x, pos.z, pos.y));
+	return trajectory;
+}
+
 void KameleonWrapper::getGridVariables(std::string& x, std::string& y, std::string& z) {
 	// get the grid system string
 	std::string gridSystem = _model->getGlobalAttribute("grid_system_1").getAttributeString();