Fixed various bugs. Updated README and docs.

README.md

@@ -53,7 +53,7 @@ Documentation of SOFT is available in two forms:
 SOFT is developed by the [Plasma Theory group](http://ft.nephy.chalmers.se/) at the Department of Physics, Chalmers University of Technology.
 
 ## How to cite
-When refering to SOFT, please cite the [SOFT paper](https://arxiv.org/abs/1709.00674):
+When refering to SOFT, please cite the [SOFT paper](https://doi.org/10.1088/1741-4326/aa9abb):
 ```
-[1] M. Hoppe, O. Embréus, R. A. Tinguely, R. S. Granetz, A. Stahl and T. Fülöp, "SOFT: a synthetic synchrotron diagnostic for runaway electrons", Accepted for publication in Nucl. Fusion.
+[1] M. Hoppe, O. Embréus, R. A. Tinguely, R. S. Granetz, A. Stahl and T. Fülöp, "SOFT: a synthetic synchrotron diagnostic for runaway electrons", Nucl. Fusion **58** 026032 [https://doi.org/10.1088/1741-4326/aa9abb].
 ```

docs/sphinx/conf.py

@@ -47,7 +47,7 @@
 
 # General information about the project.
 project = 'SOFT'
-copyright = '2017, Mathias Hoppe'
+copyright = '2017-2018, Mathias Hoppe'
 author = 'Mathias Hoppe'
 
 # The version info for the project you're documenting, acts as replacement for

docs/sphinx/index.rst

@@ -23,6 +23,7 @@ mathematical details of how SOFT is implemented, you should look into the
   distributions
   geomkern
   polarization
+  space3d
   paramref
   troubleshooting
 

docs/sphinx/paramref.rst

@@ -504,7 +504,7 @@ The ``image`` sycout generates a camera image.
 space3d
 ^^^^^^^
 The ``space3d`` can be used to store 3D data about the points of space that contribute
-to an image. *Will be described in more detail in a separate section*
+to an image. A description about how to use it can be found in `space3d <space3d.html>`__.
 
 .. option:: output
 

docs/sphinx/space3d.rst

@@ -0,0 +1,62 @@
+.. _space3d:
+
+3D emission maps
+================
+Due to the highly anisotropic nature of bremsstrahlung and synchrotron
+radiation combined with the fact that radiation is only detected if it's
+emitted directly at the detector, a given detector can only measure radiation
+from particles in a certain regions of space. It can be shown that these
+regions of space all satisfy (approximately) the condition
+
+.. math::
+   \hat{\boldsymbol{b}}(\boldsymbol{x})\cdot\frac{\boldsymbol{x}-\boldsymbol{X}}{\left|\boldsymbol{x}-\boldsymbol{X}\right|} = \cos\theta_{\mathrm{p}},
+   :label: bdotr
+
+where :math:`\hat{\boldsymbol{b}}` is the magnetic field unit vector,
+:math:`\boldsymbol{x}` is the particle's position, :math:`\boldsymbol{X}` is
+the detector's position and :math:`\theta_{\mathrm{p}}` denotes the particle's
+pitch angle (note that the pitch angle also varies as the particle moves in
+the inhomogeneous magnetic field, and therefore picks up a dependence on
+:math:`\boldsymbol{x}`). The solution to this equation, i.e. the points
+:math:`\boldsymbol{x}` satisfying it, trace out a surface in real space which
+we refer to as the `surface-of-visibility`. When the detector is located in
+the midplane, this surface typically takes the shape of a twisted cylinder.
+
+Solving for surface-of-visibility
+---------------------------------
+Using SOFT it is possible to solve :eq:`bdotr`, accounting for the finite
+detector size. This is done by adding the sycout ``space3d`` (see
+`space3d <paramref.html#space3d>`__ for a parameter reference) to your SOFT
+runscript. One example definition of the sycout is as follows ::
+
+  sycout space3d {
+      output=outfile.mat
+      type=pixels
+      pixels=200
+      point0=-0.5,-0.25,-0.25
+      point1=0.5,0.25,0.25
+  }
+
+The ``output`` parameter specifies the name of the output file, and the
+``type`` parameter specifies the algorithm to use for storing 3D information.
+Setting ``type=pixels`` means SOFT will divide the space into :math:`N^3`
+cells, where :math:`N` is the value assigned to the ``pixels`` parameter,
+between the two edge points ``point0`` and ``point1`` (see the figure below;
+the red dots indicate the locations of the edge points). During the SOFT run,
+each cell records the radiation being emitted from the box and accumulates it. 
+
+.. image:: box.png
+   :width: 500 px
+   :align: center
+
+The other value available for ``type`` is ``real``, which stores the exact
+coordinates of each point that contributes to the final image. This means that
+the output will be more detailed, but it will also grow with each particle.
+
+Visualizing
+-----------
+Visualization of space3d files is complicatd by the fact that each point
+represents emitted light, which adds together along lines of sights. A simple
+C program has been written by Mathias for generating sequences of PNG images
+from S3D output files. The program is available on GitHub:
+`s3dvid <https://github.com/hoppe93/s3dvid>`__.

src/magnetic/magnetic_circ.c

@@ -64,8 +64,8 @@ void magnetic_circ_init(struct general_settings *set) {
 	double *r = malloc(sizeof(double)*n),
 		   *z = malloc(sizeof(double)*n);
 	for (i = 0; i < n; i++) {
-		r[i] = mcpar->r * cos(2.0*PI * i / (double)n);
-		z[i] = mcpar->r * sin(2.0*PI * i / (double)n);
+		r[i] = mcpar->Rm + mcpar->r * cos(2.0*PI * i / (double)n);
+		z[i] = mcpar->Rm + mcpar->r * sin(2.0*PI * i / (double)n);
 	}
 
 	domain_set(r, z, n);

src/main.c

@@ -7,6 +7,7 @@
 #include <sys/time.h>
 
 #include "config.h"
+#include "constants.h"
 #include "ctsv.h"
 //#include "diag.h"
 #include "domain.h"
@@ -285,7 +286,7 @@ double main_solve(particle *p, equation *eq, magnetic_handler *mh, tool *usetool
 		r = R[1] = hypot(x, y); Z[1] = z;
 		if (domain_check(R, Z) == DOMAIN_OUTSIDE) {
 			double v0 = hypot(p->v0[0], hypot(p->v0[1], p->v0[2])),
-				   gm = 1 / sqrt(1 - v0*v0),
+				   gm = 1 / sqrt(1 - v0*v0/(LIGHTSPEED*LIGHTSPEED)),
 				   p0 = gm * p->mass * v0,
 				   xi0= p->vpar / v0;
 			printf("Particle collided with device wall!  r0 = %e, p0 = %e, xi0 = %e\n", r0, p0, xi0);

src/particles.c

@@ -486,9 +486,9 @@ particle *particles_generate_distributed(void) {
 	}
 
 	/* Last particle for this thread? */
-	if (particles_r == particles_rend &&
-		particles_param1 == particles_param1end &&
-		particles_param2 == particles_param2end)
+	if (particles_r == particles_rend-1 &&
+		particles_param1 == particles_param1end-1 &&
+		particles_param2 == particles_param2end-1)
 		particles_done = 1;
 
 	/* Move on to next particle */
@@ -555,9 +555,9 @@ particle *particles_generate_queue(void) {
 		}
 
 		/* Last particle for this thread? */
-		if (particles_g_r == particles_rn &&
-			particles_g_param1 == particles_param1n &&
-			particles_g_param2 == particles_param2n)
+		if (particles_g_r == particles_rn-1 &&
+			particles_g_param1 == particles_param1n-1 &&
+			particles_g_param2 == particles_param2n-1)
 			particles_g_done = 1;
 
 		/* Move on to next particle */
@@ -779,7 +779,7 @@ double particles_compute_X_pol(double ppar, double pperp, double r) {
 	return hypot(Xdotr, Xdotz) / Beffpar;
 }
 double particles_find_axis_r(double ppar, double pperp) {
-	double a = particles_rinner, b = particles_router, c, d, tol = 1e-4,
+	double a = particles_rinner, b = particles_router, c, d, tol = 1e-7,
 		   Xpol_c=0, Xpol_d=0,
 		   phi = 2 / (1+sqrt(5));
 
@@ -808,6 +808,14 @@ double particles_find_axis_r(double ppar, double pperp) {
 	}
 
     particles_effective_magnetic_axis = (a+b)*0.5;
+    /*
+    fprintf(
+        stderr, "%.16e,%.16e,%.16e\n",
+        ppar/(ELECTRONMASS*LIGHTSPEED), pperp/(ELECTRONMASS*LIGHTSPEED),
+        particles_effective_magnetic_axis
+    );
+    */
+
 	return particles_effective_magnetic_axis;
 }
 /**