diff --git a/Makefile b/Makefile
index f6dc66d..0ac0c5d 100644
--- a/Makefile
+++ b/Makefile
@@ -50,23 +50,23 @@ endif
 
 TARGET  = lime.x 
 CC		= gcc -fopenmp
-SRCS    = src/aux.c src/messages.c src/grid.c src/LTEsolution.c   \
-		  src/main.c src/molinit.c src/photon.c src/popsin.c    \
-		  src/popsout.c src/predefgrid.c src/ratranInput.c      \
-          src/raytrace.c src/smooth.c src/sourcefunc.c          \
-		  src/stateq.c src/statistics.c src/magfieldfit.c       \
-		  src/stokesangles.c src/writefits.c src/weights.c      \
-		  src/velospline.c src/getclosest.c src/frees.c  \
-		  src/tcpsocket.c src/defaults.c src/fastexp.c
+SRCS	= src/aux.c src/messages.c src/grid.c src/LTEsolution.c	\
+	  src/main.c src/molinit.c src/photon.c src/popsin.c	\
+	  src/popsout.c src/predefgrid.c src/ratranInput.c	\
+	  src/raytrace.c src/smooth.c src/sourcefunc.c src/frees.c	\
+	  src/stateq.c src/statistics.c src/magfieldfit.c	\
+	  src/stokesangles.c src/writefits.c src/weights.c	\
+	  src/velospline.c src/getclosest.c src/raythrucells.c	\
+	  src/tcpsocket.c src/defaults.c src/fastexp.c
 MODELS  = model.c
-OBJS    = src/aux.o src/messages.o src/grid.o src/LTEsolution.o   \
-		  src/main.o src/molinit.o src/photon.o src/popsin.o    \
-		  src/popsout.o src/predefgrid.o src/raytrace.o         \
-		  src/ratranInput.o src/smooth.o src/sourcefunc.o       \
-		  src/stateq.o src/statistics.o src/magfieldfit.o       \
-		  src/stokesangles.o src/writefits.o src/weights.o      \
-		  src/velospline.o src/getclosest.o src/frees.o  \
-		  src/tcpsocket.o src/defaults.o src/fastexp.o
+OBJS    = src/aux.o src/messages.o src/grid.o src/LTEsolution.o	\
+	  src/main.o src/molinit.o src/photon.o src/popsin.o	\
+	  src/popsout.o src/predefgrid.o src/ratranInput.o	\
+	  src/raytrace.o src/smooth.o src/sourcefunc.o src/frees.o	\
+	  src/stateq.o src/statistics.o src/magfieldfit.o	\
+	  src/stokesangles.o src/writefits.o src/weights.o	\
+	  src/velospline.o src/getclosest.o src/raythrucells.o  \
+	  src/tcpsocket.o src/defaults.o src/fastexp.o
 MODELO 	= src/model.o
 
 #CCFLAGS = -O3 -falign-loops=16 -fno-strict-aliasing -DTEST
diff --git a/doc/usermanual.rst b/doc/usermanual.rst
index 3359f9a..2e2e7c1 100644
--- a/doc/usermanual.rst
+++ b/doc/usermanual.rst
@@ -435,13 +435,7 @@ file. The default is blend=0 (no line blending).
 
     (integer) par->antialias (optional)
 
-If set, LIME will anti-alias the output image. anti-alias can take the
-value of any positive integer, with the value 1 (default) being no
-anti-aliasing. Greater values correspond to stronger anti-aliasing. LIME
-uses stochastic super-sampling anti-aliasing. This is very effective in
-minimizing artifacts in the image, but it also slows down the ray-tracer
-by a factor equal to the value of antialias. This parameter is the only
-one that will not be ignored in case par->restart is set.
+This parameter is no longer used, although it is retained for the present for purposes of backward compatibility.
 
 .. code:: c
 
@@ -481,6 +475,12 @@ The final one of the density-linked parameters controls how the dust opacity is
 
 If none of the three density-linked parameters is provided, LIME will attempt to guess the information, in a manner as close as possible to the way it was done in version 1.5 and earlier. This is safe enough when a single density value is returned, and only H2 provided as collision partner in the moldata file(s), but more complicated situations can very easily result in the code guessing wrongly. For this reason we encourage users to make use of these three parameters, although in order to preserve backward compatibility with old model.c files, we have not (yet) made them mandatory.
 
+    (integer) par->traceRayAlgorithm (optional)
+
+This parameter specifies the algorithm used by LIME to solve the radiative-transfer equations during ray-tracing. The default value of zero invokes the algorithm used in LIME-1.5 and previous; a value of 1 invokes a new algorithm which is much more time-consuming but which produces much smoother images, free from step-artifacts.
+
+Note that there have been additional modifications to the raytracing algorithm which have significant effects on the output images since LIME-1.5. Image-plane interpolation is now employed in areas of the image where the grid point spacing is larger than the image pixel spacing. This leads both to a smoother image and a shorter processing time.
+
 .. _par-nthreads:
 
 .. code:: c
diff --git a/src/LTEsolution.c b/src/LTEsolution.c
index 6ea7df3..6e6c3c7 100644
--- a/src/LTEsolution.c
+++ b/src/LTEsolution.c
@@ -15,7 +15,7 @@ LTE(configInfo *par, struct grid *g, molData *m){
 
   for(id=0;id<par->pIntensity;id++){
     for(ispec=0;ispec<par->nSpecies;ispec++){
-      g[id].nmol[ispec]=g[id].abun[ispec]*g[id].dens[0];
+      g[id].mol[ispec].nmol = g[id].abun[ispec]*g[id].dens[0];
       lteOnePoint(par, m, ispec, g[id].t[0], g[id].mol[ispec].pops);
     }
   }
diff --git a/src/aux.c b/src/aux.c
index 7360302..8e244fc 100644
--- a/src/aux.c
+++ b/src/aux.c
@@ -41,6 +41,7 @@ parseInput(inputPars inpar, configInfo *par, image **img, molData **m){
   par->antialias    = inpar.antialias;
   par->polarization = inpar.polarization;
   par->nThreads     = inpar.nThreads;
+  par->traceRayAlgorithm = inpar.traceRayAlgorithm;
 
   /* Now set the additional values in par. */
   par->ncell = inpar.pIntensity + inpar.sinkPoints;
@@ -108,9 +109,9 @@ parseInput(inputPars inpar, configInfo *par, image **img, molData **m){
   while((*img)[par->nImages].filename!=NULL && par->nImages<MAX_NIMAGES)
     par->nImages++;
 
-  /* Check that the user has supplied this function (needed unless par->pregrid):
+  /* Check that the user has supplied this function (needed unless par->doPregrid):
   */
-  if(!par->doPregrid)
+  if(!par->doPregrid || par->traceRayAlgorithm==1)
     velocity(0.0,0.0,0.0, dummyVel);
 
   /*
diff --git a/src/frees.c b/src/frees.c
index 83e6ca5..84c88e0 100644
--- a/src/frees.c
+++ b/src/frees.c
@@ -10,7 +10,18 @@
 #include "lime.h"
 
 void
-freeGrid(configInfo *par, const molData *m ,struct grid *g){
+freeGAux(const unsigned long numPoints, const int numSpecies, struct gAuxType *gAux) {
+  unsigned long ppi;
+
+  if(gAux !=NULL){
+    for(ppi=0;ppi<numPoints;ppi++)
+      freePop2(numSpecies, gAux[ppi].mol);
+    free(gAux);
+  }
+}
+
+void
+freeGrid(configInfo *par, molData *m , struct grid *g){
   int i;
   if(g != NULL){
     for(i=0;i<(par->pIntensity+par->sinkPoints); i++){
@@ -23,7 +34,6 @@ freeGrid(configInfo *par, const molData *m ,struct grid *g){
       free(g[i].neigh);
       free(g[i].w);
       free(g[i].dens);
-      free(g[i].nmol);
       free(g[i].abun);
       free(g[i].ds);
       freePopulation( par, m, g[i].mol );
@@ -80,11 +90,8 @@ void freeMolsWithBlends(struct molWithBlends *mols, const int numMolsWithBlends)
   if(mols != NULL){
     for(mi=0;mi<numMolsWithBlends;mi++){
       if(mols[mi].lines != NULL){
-        for(li=0;li<mols[mi].numLinesWithBlends;li++){
-          if(mols[mi].lines[li].blends != NULL){
-            free(mols[mi].lines[li].blends);
-          }
-        }
+        for(li=0;li<mols[mi].numLinesWithBlends;li++)
+          free(mols[mi].lines[li].blends);
         free(mols[mi].lines);
       }
     }
@@ -114,7 +121,21 @@ freeParImg(const int nImages, inputPars *par, image *img){
 }
 
 void
-freePopulation(configInfo *par, const molData *m, struct populations *pop ){
+freePop2(const int numSpecies, struct pop2 *mol){
+  int i;
+
+  if(mol !=NULL){
+    for(i=0;i<numSpecies;i++){
+      free(mol[i].specNumDens);
+      free(mol[i].knu);
+      free(mol[i].dust);
+    }
+    free(mol);
+  }
+}
+
+void
+freePopulation(configInfo *par, molData *m, struct populations *pop ){
   if(pop != NULL){
     int j,k;
     for(j=0; j<par->nSpecies; j++){
diff --git a/src/grid.c b/src/grid.c
index f570f19..b5ef829 100644
--- a/src/grid.c
+++ b/src/grid.c
@@ -3,10 +3,8 @@
  *  This file is part of LIME, the versatile line modeling engine
  *
  *  Copyright (C) 2006-2014 Christian Brinch
- *  Copyright (C) 2015 The LIME development team
+ *  Copyright (C) 2016 The LIME development team
  *
-TODO:
-  - There is no need to malloc nmol if all the images are non-line.
  */
 
 #include "lime.h"
@@ -32,7 +30,6 @@ gridAlloc(configInfo *par, struct grid **g){
     (*g)[i].ds = NULL;
     (*g)[i].dens=malloc(sizeof(double)*par->numDensities);
     (*g)[i].abun=malloc(sizeof(double)*par->nSpecies);
-    (*g)[i].nmol=malloc(sizeof(double)*par->nSpecies);
     (*g)[i].t[0]=-1;
     (*g)[i].t[1]=-1;
   }
@@ -63,9 +60,9 @@ void calcGridMolDensities(configInfo *par, struct grid *g){
 
   for(id=0;id<par->ncell; id++){
     for(ispec=0;ispec<par->nSpecies;ispec++){
-      g[id].nmol[ispec] = 0.0;
+      g[id].mol[ispec].nmol = 0.0;
       for(i=0;i<par->numDensities;i++)
-        g[id].nmol[ispec] += g[id].abun[ispec]*g[id].dens[i]*par->nMolWeights[i];
+        g[id].mol[ispec].nmol += g[id].abun[ispec]*g[id].dens[i]*par->nMolWeights[i];
     }
   }
 }
@@ -127,8 +124,8 @@ void calcGridDustOpacity(configInfo *par, molData *md, struct grid *gp){
       for(di=0;di<par->numDensities;di++)
         densityForDust += gp[id].dens[di]*par->dustWeights[di];
 
+      gasIIdust(gp[id].x[0],gp[id].x[1],gp[id].x[2],&gtd);
       for(iline=0;iline<md[si].nline;iline++){
-        gasIIdust(gp[id].x[0],gp[id].x[1],gp[id].x[2],&gtd);
         gp[id].mol[si].knu[iline]=kappatab[iline]*2.4*AMU*densityForDust/gtd;
         //Check if input model supplies a dust temperature. Otherwise use the kinetic temperature
         if(gp[id].t[1]==-1) {
@@ -182,74 +179,154 @@ void calcGridCollRates(configInfo *par, molData *md, struct grid *g){
   } /* End loop over radiating molecules. */
 }
 
-void
-qhull(configInfo *par, struct grid *g){
-  int i,j,k,id;
+/*....................................................................*/
+void delaunay(const int numDims, struct grid *gp, const unsigned long numPoints\
+  , const _Bool getCells, struct cell **dc, unsigned long *numCells){
+  /*
+The principal purpose of this function is to perform a Delaunay triangulation for the set of points defined by the input argument 'g'. This is achieved via routines in the 3rd-party package qhull.
+
+A note about qhull nomenclature: a vertex is what you think it is - i.e., a point; but a facet means in this context a Delaunay triangle (in 2D) or tetrahedron (in 3D).
+
+Required elements of structs:
+	struct grid *gp:
+		.id
+		.x
+
+Elements of structs are set as follows:
+	struct grid *gp:
+		.numNeigh
+		.neigh (this is malloc'd too large and at present not realloc'd.)
+
+	cellType *dc (if getCells>0):
+		.id
+		.neigh
+		.vertx
+  */
+
+  coordT *pt_array=NULL;
+  unsigned long ppi,id,pointIdsThisFacet[numDims+1],idI,idJ,fi,ffi;
+  int i,j,k;
   char flags[255];
   boolT ismalloc = False;
-  facetT *facet;
   vertexT *vertex,**vertexp;
-  coordT *pt_array;
-  int simplex[DIM+1];
+  facetT *facet, *neighbor, **neighborp;
   int curlong, totlong;
+  _Bool neighbourNotFound;
 
-  pt_array=malloc(DIM*sizeof(coordT)*par->ncell);
-
-  for(i=0;i<par->ncell;i++) {
-    for(j=0;j<DIM;j++) {
-      pt_array[i*DIM+j]=g[i].x[j];
+  pt_array=malloc(sizeof(coordT)*numDims*numPoints);
+  for(ppi=0;ppi<numPoints;ppi++) {
+    for(j=0;j<numDims;j++) {
+      pt_array[ppi*numDims+j]=gp[ppi].x[j];
     }
   }
 
   sprintf(flags,"qhull d Qbb");
-  if (!qh_new_qhull(DIM, par->ncell, pt_array, ismalloc, flags, NULL, NULL)) {
-    /* Identify points */
-    FORALLvertices {
-      id=qh_pointid(vertex->point);
-      g[id].numNeigh=qh_setsize(vertex->neighbors);
-      if(  g[id].neigh != NULL )
-        {
-          free( g[id].neigh );
+  if (qh_new_qhull(numDims, (int)numPoints, pt_array, ismalloc, flags, NULL, NULL)) {
+    if(!silent) bail_out("Qhull failed to triangulate");
+    exit(1);
+  }
+
+  /* Identify points */
+  FORALLvertices {
+    id=(unsigned long)qh_pointid(vertex->point);
+    /* Note this is NOT the same value as vertex->id. Only the id gained via the call to qh_pointid() is the same as the index of the point in the input list. */
+
+    gp[id].numNeigh=qh_setsize(vertex->neighbors);
+    /* Note that vertex->neighbors refers to facets abutting the vertex, not other vertices. In general there seem to be more facets surrounding a point than vertices (in fact there seem to be exactly 2x as many). In any case, mallocing to N_facets gives extra room. */
+
+    if(gp[id].neigh!=NULL)
+      free( gp[id].neigh );
+    gp[id].neigh=malloc(sizeof(struct grid *)*gp[id].numNeigh);
+    for(k=0;k<gp[id].numNeigh;k++) {
+      gp[id].neigh[k]=NULL;
+    }
+  }
+
+  /* Identify the Delaunay neighbors of each point. This is a little involved, because the only direct information we have about which vertices are linked to which others is stored in qhull's facetT objects.
+  */
+  *numCells = 0;
+  FORALLfacets {
+    if (!facet->upperdelaunay) {
+      /* Store the point IDs in a list for convenience. These ID values are conveniently ordered such that qh_pointid() returns ppi for gp[ppi]. 
+      */
+      j=0;
+      FOREACHvertex_ (facet->vertices) pointIdsThisFacet[j++]=(unsigned long)qh_pointid(vertex->point);
+
+      for(i=0;i<numDims+1;i++){
+        idI = pointIdsThisFacet[i];
+        for(j=0;j<numDims+1;j++){
+          idJ = pointIdsThisFacet[j];
+          if(i!=j){
+            /* Cycle through all the non-NULL links of gp[idI], storing the link if it is new.
+            */
+            k=0;
+            while(gp[idI].neigh[k] != NULL && gp[idI].neigh[k]->id != gp[idJ].id)
+              k++;
+            gp[idI].neigh[k]=&gp[idJ];
+          }
         }
-      g[id].neigh=malloc(sizeof(struct grid *)*g[id].numNeigh);
-      for(k=0;k<g[id].numNeigh;k++) {
-        g[id].neigh[k]=NULL;
       }
+      (*numCells)++;
+    }
+  }
+
+  for(ppi=0;ppi<numPoints;ppi++){
+    j=0;
+    for(k=0;k<gp[ppi].numNeigh;k++){
+      if(gp[ppi].neigh[k] != NULL)
+        j++;
     }
-    
-    /* Identify neighbors */
+    gp[ppi].numNeigh=j;
+  }
+
+  if(getCells){
+    (*dc) = malloc(sizeof(**dc)*(*numCells));
+    fi = 0;
     FORALLfacets {
       if (!facet->upperdelaunay) {
-        j=0;
-        FOREACHvertex_ (facet->vertices) simplex[j++]=qh_pointid(vertex->point);
-        for(i=0;i<DIM+1;i++){
-          for(j=0;j<DIM+1;j++){
-            k=0;
-            if(i!=j){
-              while(g[simplex[i]].neigh[k] != NULL && g[simplex[i]].neigh[k]->id != g[simplex[j]].id) {
-                k++;
+        (*dc)[fi].id = (unsigned long)facet->id; /* Do NOT expect this to be equal to fi. */
+        fi++;
+      }
+    }
+
+    fi = 0;
+    FORALLfacets {
+      if (!facet->upperdelaunay) {
+        i = 0;
+        FOREACHneighbor_(facet) {
+          if(neighbor->upperdelaunay){
+            (*dc)[fi].neigh[i] = NULL;
+          }else{
+            /* Have to find the member of *dc with the same id as neighbour.*/
+            ffi = 0;
+            neighbourNotFound=1;
+            while(ffi<(*numCells) && neighbourNotFound){
+              if((*dc)[ffi].id==(unsigned long)neighbor->id){
+                (*dc)[fi].neigh[i] = &(*dc)[ffi];
+                neighbourNotFound = 0;
               }
-              g[simplex[i]].neigh[k]=&g[simplex[j]];
+              ffi++;
+            }
+            if(ffi>=(*numCells) && neighbourNotFound){
+              if(!silent) bail_out("Something weird going on.");
+              exit(1);
             }
           }
+          i++;
         }
-      }
-    }
-  } else {
-    if(!silent) bail_out("Qhull failed to triangulate");
-    exit(1);
-  }
 
-  for(i=0;i<par->ncell;i++){
-    j=0;
-    for(k=0;k<g[i].numNeigh;k++){
-      if(g[i].neigh[k] != NULL)
-        {
-          j++;
+        i = 0;
+        FOREACHvertex_( facet->vertices ) {
+          id = (unsigned long)qh_pointid(vertex->point);
+          (*dc)[fi].vertx[i] = &gp[id];
+          i++;
         }
+
+        fi++;
+      }
     }
-    g[i].numNeigh=j;
   }
+
   qh_freeqhull(!qh_ALL);
   qh_memfreeshort (&curlong, &totlong);
   free(pt_array);
@@ -554,6 +631,8 @@ buildGrid(configInfo *par, struct grid *g){
   double temp;
   int k=0,i,j;            /* counters									*/
   int flag;
+  struct cell *dc=NULL; /* Not used at present. */
+  unsigned long numCells;
   const int maxNumAttempts=1000;
   _Bool numRandomsThisPoint;
   int numSecondRandoms=0;
@@ -670,9 +749,9 @@ buildGrid(configInfo *par, struct grid *g){
   temperature(0.0,0.0,0.0, g[0].t);
   doppler(    0.0,0.0,0.0,&g[0].dopb);	
   abundance(  0.0,0.0,0.0, g[0].abun);
-  /* Note that velocity() is the only one of the 5 mandatory functions which is still needed (in raytrace) even if par->pregrid or par->restart. Therefore we test it already in parseInput(). */
+  /* Note that velocity() is the only one of the 5 mandatory functions which is still needed (in raytrace) unless par->doPregrid. Therefore we test it already in parseInput(). */
 
-  qhull(par, g);
+  delaunay(DIM, g, (unsigned long)par->ncell, 0, &dc, &numCells);
   distCalc(par, g);
   smooth(par,g);
 
@@ -685,10 +764,22 @@ buildGrid(configInfo *par, struct grid *g){
 
   checkGridDensities(par, g);
 
+  if(par->polarization){
+    for(i=0;i<par->pIntensity;i++)
+      magfield(g[i].x[0],g[i].x[1],g[i].x[2], g[i].B);
+  }else{
+    for(i=0;i<par->pIntensity;i++){
+      g[i].B[0]=0.0;
+      g[i].B[1]=0.0;
+      g[i].B[2]=0.0;
+    }
+  }
+
   //	getArea(par,g, ran);
   //	getMass(par,g, ran);
   getVelosplines(par,g);
   dumpGrid(par,g);
+  free(dc);
 
   gsl_rng_free(ran);
   if(!silent) done(5);
diff --git a/src/inpars.h b/src/inpars.h
index bc80ec4..dee86dc 100644
--- a/src/inpars.h
+++ b/src/inpars.h
@@ -4,7 +4,7 @@
 /* input parameters */
 typedef struct {
   double radius,minScale,tcmb,*nMolWeights,*dustWeights;
-  int sinkPoints,pIntensity,blend,*collPartIds;
+  int sinkPoints,pIntensity,blend,*collPartIds,traceRayAlgorithm;
   char *outputfile,*binoutputfile;
 //  char *inputfile; unused at present.
   char *gridfile;
diff --git a/src/lime.h b/src/lime.h
index f80f295..6053c86 100644
--- a/src/lime.h
+++ b/src/lime.h
@@ -112,7 +112,7 @@
 typedef struct {
   double radius,minScale,tcmb,*nMolWeights,*dustWeights;
   double radiusSqu,minScaleSqu,taylorCutoff;
-  int sinkPoints,pIntensity,blend,*collPartIds;
+  int sinkPoints,pIntensity,blend,*collPartIds,traceRayAlgorithm;
   int ncell,nImages,nSpecies,numDensities,doPregrid;
   char *outputfile, *binoutputfile;
   char *gridfile;
@@ -160,15 +160,14 @@ struct rates {
 
 struct populations {
   double *pops, *knu, *dust;
-  double dopb, binv;
+  double dopb, binv, nmol;
   struct rates *partner;
 };
 
 /* Grid properties */
 struct grid {
   int id;
-  double x[3];
-  double vel[3];
+  double x[DIM], vel[DIM], B[3]; /* B field only makes physical sense in 3 dimensions. */
   double *a0,*a1,*a2,*a3,*a4;
   int numNeigh;
   point *dir;
@@ -177,15 +176,16 @@ struct grid {
   int sink;
   int nphot;
   int conv;
-  double *dens,t[2],*nmol,*abun, dopb;
+  double *dens,t[2],*abun, dopb;
   double *ds;
-  struct populations* mol;
+  struct populations *mol;
 };
 
 typedef struct {
   double *intense;
   double *tau;
   double stokes[3];
+  int numRays;
 } spec;
 
 /* Image information */
@@ -205,7 +205,10 @@ typedef struct {
   double rotMat[3][3];
 } image;
 
-typedef struct {double x,y, *intensity, *tau;} rayData;
+typedef struct {
+  double x,y, *intensity, *tau;
+  unsigned int ppi;
+} rayData;
 
 struct blend{
   int molJ, lineJ;
@@ -227,6 +230,51 @@ struct blendInfo{
   struct molWithBlends *mols;
 };
 
+/* NOTE that it is assumed that vertx[i] is opposite the face that abuts with neigh[i] for all i.
+*/ 
+struct cell {
+  struct grid *vertx[DIM+1];
+  struct cell *neigh[DIM+1]; /* ==NULL flags an external face. */
+  unsigned long id;
+  double centre[DIM];
+};
+
+struct pop2 {
+  double *specNumDens, *knu, *dust;
+  double binv;
+};
+
+typedef struct{
+  double x[DIM], xCmpntRay, B[3];
+  struct pop2 *mol;
+} gridInterp;
+
+struct gAuxType{
+  struct pop2 *mol;
+};
+
+/* This struct is meant to record all relevant information about the intersection between a ray (defined by a direction unit vector 'dir' and a starting position 'r') and a face of a Delaunay cell.
+*/
+typedef struct {
+  int fi;
+  /* The index (in the range {0...DIM}) of the face (and thus of the opposite vertex, i.e. the one 'missing' from the bary[] list of this face).
+  */
+  int orientation;
+  /* >0 means the ray exits, <0 means it enters, ==0 means the face is parallel to ray.
+  */
+  double bary[DIM], dist, collPar;
+  /* 'dist' is defined via r_int = r + dist*dir. 'collPar' is a measure of how close to any edge of the face r_int lies.
+  */
+} intersectType;
+
+typedef struct {
+  double r[DIM][DIM], centre[DIM];/*, norm[3], mat[1][1], det; */
+} faceType;
+
+typedef struct {
+  double xAxis[DIM], yAxis[DIM], r[3][2];
+} triangle2D;
+
 /* Some global variables */
 int silent;
 
@@ -246,40 +294,56 @@ void	run(inputPars, image *);
 void	assignMolCollPartsToDensities(configInfo*, molData*);
 void	binpopsout(configInfo *, struct grid *, molData *);
 void	buildGrid(configInfo *, struct grid *);
+int	buildRayCellChain(double*, double*, struct grid*, struct cell*, _Bool**, unsigned long, int, int, int, const double, unsigned long**, intersectType**, int*);
 void	calcFastExpRange(const int, const int, int*, int*, int*);
 void	calcGridCollRates(configInfo*, molData*, struct grid*);
 void	calcGridDustOpacity(configInfo*, molData*, struct grid*);
 void	calcGridMolDensities(configInfo*, struct grid*);
+void	calcLineAmpInterp(const double, const double, const double, double*);
+void	calcLineAmpLinear(struct grid*, const int, const int, const double, const double, double*);
+void	calcLineAmpSample(const double x[3], const double dx[3], const double, const double, double*, const int, const double, const double, double*);
+void   	calcLineAmpSpline(struct grid*, const int, const int, const double, const double, double*);
 void	calcMolCMBs(configInfo*, molData*);
 void	calcSourceFn(double, const configInfo*, double*, double*);
 void	calcTableEntries(const int, const int);
+void	calcTriangleBaryCoords(double vertices[3][2], double, double, double barys[3]);
+triangle2D calcTriangle2D(faceType);
 void	checkGridDensities(configInfo*, struct grid*);
 void	checkUserDensWeights(configInfo*);
 void	continuumSetup(int, image*, molData*, configInfo*, struct grid*);
+void	delaunay(const int, struct grid*, const unsigned long, const _Bool, struct cell**, unsigned long*);
 void	distCalc(configInfo*, struct grid*);
+void	doBaryInterp(const intersectType, struct grid*, struct gAuxType*, double*, unsigned long*, molData*, const int, gridInterp*);
+void	doSegmentInterp(gridInterp*, const int, molData*, const int, const double, const int);
+faceType extractFace(struct grid*, struct cell*, const unsigned long, const int);
 int	factorial(const int);
 double	FastExp(const float);
 void	fit_d1fi(double, double, double*);
 void	fit_fi(double, double, double*);
 void	fit_rr(double, double, double*);
-void	freeGrid(configInfo*, const molData*, struct grid*);
+int	followRayThroughDelCells(double*, double*, struct grid*, struct cell*, const unsigned long, const double, intersectType*, unsigned long**, intersectType**, int*);
+void	freeGAux(const unsigned long, const int, struct gAuxType*);
+void	freeGrid(configInfo*, molData*, struct grid*);
 void	freeGridPointData(configInfo*, gridPointData*);
 void	freeMolData(const int, molData*);
 void	freeMolsWithBlends(struct molWithBlends*, const int);
 void	freeParImg(const int, inputPars*, image*);
-void	freePopulation(configInfo*, const molData*, struct populations*);
+void	freePopulation(configInfo*, molData*, struct populations*);
+void	freePop2(const int, struct pop2*);
 double	gaussline(double, double);
 void	getArea(configInfo *, struct grid *, const gsl_rng *);
 void	getclosest(double, double, double, long *, long *, double *, double *, double *);
 void	getjbar(int, molData*, struct grid*, const int, configInfo*, struct blendInfo, int, gridPointData*, double*);
 void	getMass(configInfo *, struct grid *, const gsl_rng *);
 void	getmatrix(int, gsl_matrix *, molData *, struct grid *, int, gridPointData *);
+int	getNewEntryFaceI(const unsigned long, const struct cell);
 int	getNextEdge(double*, int, struct grid*, const gsl_rng*);
 void	getVelosplines(configInfo *, struct grid *);
 void	getVelosplines_lin(configInfo *, struct grid *);
 void	gridAlloc(configInfo *, struct grid **);
 void	gridLineInit(configInfo*, molData*, struct grid*);
 void	input(inputPars *, image *);
+void	intersectLineTriangle(double*, double*, faceType, intersectType*);
 float	invSqrt(float);
 void	levelPops(molData *, configInfo *, struct grid *, int *);
 void	line_plane_intersect(struct grid *, double *, int , int *, double *, double *, double);
@@ -294,7 +358,6 @@ int	pointEvaluation(configInfo*, double, double, double, double);
 void	popsin(configInfo *, struct grid **, molData **, int *);
 void	popsout(configInfo *, struct grid *, molData *);
 void	predefinedGrid(configInfo *, struct grid *);
-void	qhull(configInfo *, struct grid *);
 double	ratranInput(char *, char *, double, double, double);
 void	raytrace(int, configInfo *, struct grid *, molData *, image *);
 void	readDummyCollPart(FILE*, const int);
@@ -305,18 +368,21 @@ void	setUpConfig(configInfo *, image **, molData **);
 void	setUpDensityAux(configInfo*, int*, const int);
 void	smooth(configInfo *, struct grid *);
 void	sourceFunc(double*, double*, double, molData*, double, struct grid*, int, int, int, int);
-void	sourceFunc_cont(double*, double*, struct grid*, int, int, int);
-void	sourceFunc_line(double*, double*, molData*, double, struct grid*, int, int, int);
-void	sourceFunc_pol(double*, double*, struct grid*, int, int, int, double (*rotMat)[3]);
+void    sourceFunc_line(const molData, const double, const struct populations, const int, double*, double*);
+void    sourceFunc_cont(const struct populations, const int, double*, double*);
+void    sourceFunc_line_raytrace(const molData, const double, const struct pop2, const int, double*, double*);
+void    sourceFunc_cont_raytrace(const struct pop2, const int, double*, double*);
+void	sourceFunc_pol(double*, const struct pop2, int, double (*rotMat)[3], double*, double*);
 void	stateq(int, struct grid*, molData*, const int, configInfo*, struct blendInfo, int, gridPointData*, double*, _Bool*);
 void	statistics(int, molData *, struct grid *, int *, double *, double *, int *);
-void	stokesangles(double, double, double, double (*rotMat)[3], double*);
+void	stokesangles(double*, double (*rotMat)[3], double*);
 double	taylor(const int, const float);
-void    traceray(rayData, int, int, int, configInfo*, struct grid*, molData*, image*, double);
-void	velocityspline(struct grid *, int, int, double, double, double*);
-void	velocityspline2(double *, double *, double, double, double, double*);
+void	traceray(rayData, const int, const int, const int, configInfo*, struct grid*, molData*, image*, struct gAuxType*, const int, int*, int*, const double, const int, const double);
+void	traceray_smooth(rayData, const int, const int, const int, configInfo*, struct grid*, molData*, image*, struct gAuxType*, const int, int*, int*, struct cell*, const unsigned long, const double, gridInterp*, const int, const double, const int, const double);
 double	veloproject(double *, double *);
-void	writefits(int, configInfo *, molData *, image *);
+void	write2Dfits(int, configInfo *, molData *, image *);
+void	write3Dfits(int, configInfo *, molData *, image *);
+void	writeFits(const int, configInfo*, molData*, image*);
 void	write_VTK_unstructured_Points(configInfo *, struct grid *);
 
 
diff --git a/src/main.c b/src/main.c
index d5eda72..3b4d8bf 100644
--- a/src/main.c
+++ b/src/main.c
@@ -69,6 +69,7 @@ initParImg(inputPars *par, image **img)
   par->antialias=1;
   par->polarization=0;
   par->nThreads = NTHREADS;
+  par->traceRayAlgorithm=0;
 
   /* Allocate initial space for molecular data file names */
   par->moldatfile=malloc(sizeof(char *)*MAX_NSPECIES);
@@ -190,7 +191,7 @@ run(inputPars inpars, image *img)
       else
         continuumSetup(i,img,m,&par,g);
       raytrace(i,&par,g,m,img);
-      writefits(i,&par,m,img);
+      writeFits(i,&par,m,img);
     }
   }
 
@@ -202,7 +203,7 @@ run(inputPars inpars, image *img)
   for(i=0;i<par.nImages;i++){
     if(img[i].doline){
       raytrace(i,&par,g,m,img);
-      writefits(i,&par,m,img);
+      writeFits(i,&par,m,img);
     }
   }
   
@@ -216,6 +217,16 @@ run(inputPars inpars, image *img)
   free(par.dustWeights);
 }
 
+void writeFits(const int i, configInfo *par, molData *m, image *img){
+  if(img[i].unit<5)
+    write3Dfits(i,par,m,img);
+  else if(img[i].unit==5)
+    write2Dfits(i,par,m,img);
+  else{
+    if(!silent) bail_out("Image unit number invalid");
+    exit(0);
+  }
+}
 
 int main () {
   /* Main program for stand-alone LIME */
@@ -234,3 +245,4 @@ int main () {
 
   return 0;
 }
+
diff --git a/src/photon.c b/src/photon.c
index b7494f4..07bfd60 100644
--- a/src/photon.c
+++ b/src/photon.c
@@ -3,7 +3,7 @@
  *  This file is part of LIME, the versatile line modeling engine
  *
  *  Copyright (C) 2006-2014 Christian Brinch
- *  Copyright (C) 2015 The LIME development team
+ *  Copyright (C) 2016 The LIME development team
  *
  */
 
@@ -64,8 +64,10 @@ int getNextEdge(double *inidir, int id, struct grid *g, const gsl_rng *ran){
 
 }
 
-void
-velocityspline(struct grid *g, int id, int k, double binv, double deltav, double *vfac){
+/*....................................................................*/
+void calcLineAmpSpline(struct grid *g, const int id, const int k\
+  , const double binv, const double deltav, double *vfac){
+
   int nspline,ispline,naver,iaver;
   double v1,v2,s1,s2,sd,v,vfacsub,d;
   
@@ -96,9 +98,10 @@ velocityspline(struct grid *g, int id, int k, double binv, double deltav, double
   return;
 }
 
+/*....................................................................*/
+void calcLineAmpLinear(struct grid *g, const int id, const int k\
+  , const double binv, const double deltav, double *vfac){
 
-void
-velocityspline_lin(struct grid *g, int id, int k, double binv, double deltav, double *vfac){
   int nspline,ispline,naver,iaver;
   double v1,v2,s1,s2,sd,v,vfacsub,d;
   
@@ -242,9 +245,9 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran\
 
         for(molI=0;molI<par->nSpecies;molI++){
           if(!par->doPregrid)
-            velocityspline(    g,here,neighI,g[id].mol[molI].binv,deltav,&vfac[molI]);
+            calcLineAmpSpline(g,here,neighI,g[id].mol[molI].binv,deltav,&vfac[molI]);
           else
-            velocityspline_lin(g,here,neighI,g[id].mol[molI].binv,deltav,&vfac[molI]);
+            calcLineAmpLinear(g,here,neighI,g[id].mol[molI].binv,deltav,&vfac[molI]);
           mp[molI].vfac[iphot]=vfac[molI];
         }
       } else {
@@ -252,9 +255,9 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran\
       
         for(molI=0;molI<par->nSpecies;molI++){
           if(!par->doPregrid)
-            velocityspline(    g,here,neighI,g[id].mol[molI].binv,deltav,&vfac[molI]);
+            calcLineAmpSpline(g,here,neighI,g[id].mol[molI].binv,deltav,&vfac[molI]);
           else
-            velocityspline_lin(g,here,neighI,g[id].mol[molI].binv,deltav,&vfac[molI]);
+            calcLineAmpLinear(g,here,neighI,g[id].mol[molI].binv,deltav,&vfac[molI]);
         }
       }
 
@@ -266,8 +269,8 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran\
           jnu=0.;
           alpha=0.;
 
-          sourceFunc_line(&jnu,&alpha,m,vfac[molI],g,here,molI,lineI);
-          sourceFunc_cont(&jnu,&alpha,g,here,molI,lineI);
+          sourceFunc_line(m[molI],vfac[molI],g[here].mol[molI],lineI,&jnu,&alpha);
+          sourceFunc_cont(g[here].mol[molI],lineI,&jnu,&alpha);
 
           dtau=alpha*ds;
           if(dtau < -30) dtau = -30;
@@ -295,11 +298,11 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran\
               velProj = deltav - blends.mols[nextMolWithBlend].lines[nextLineWithBlend].blends[bi].deltaV;
 
               if(!par->doPregrid)
-                velocityspline(    g,here,neighI,g[id].mol[molJ].binv,velProj,&vblend);
+                calcLineAmpSpline(g,here,neighI,g[id].mol[molJ].binv,velProj,&vblend);
               else
-                velocityspline_lin(g,here,neighI,g[id].mol[molJ].binv,velProj,&vblend);
+                calcLineAmpLinear(g,here,neighI,g[id].mol[molJ].binv,velProj,&vblend);
 
-              sourceFunc_line(&jnu,&alpha,m,vblend,g,here,molJ,lineJ);
+              sourceFunc_line(m[molJ],vblend,g[here].mol[molJ],lineJ,&jnu,&alpha);
               dtau=alpha*ds;
               if(dtau < -30) dtau = -30;
               calcSourceFn(dtau, par, &remnantSnu, &expDTau);
@@ -364,8 +367,8 @@ getjbar(int posn, molData *m, struct grid *g, const int molI\
         jnu=0.;
         alpha=0.;
 
-        sourceFunc_line(&jnu,&alpha,m,mp[molI].vfac[iphot],g,posn,molI,lineI);
-        sourceFunc_cont(&jnu,&alpha,g,posn,molI,lineI);
+        sourceFunc_line(m[molI],mp[molI].vfac[iphot],g[posn].mol[molI],lineI,&jnu,&alpha);
+        sourceFunc_cont(g[posn].mol[molI],lineI,&jnu,&alpha);
         tau=alpha*halfFirstDs[iphot];
         calcSourceFn(tau, par, &remnantSnu, &expTau);
         remnantSnu *= jnu*m[molI].norminv*halfFirstDs[iphot];
@@ -382,7 +385,7 @@ getjbar(int posn, molData *m, struct grid *g, const int molI\
             molJ  = blends.mols[nextMolWithBlend].lines[nextLineWithBlend].blends[bi].molJ;
             lineJ = blends.mols[nextMolWithBlend].lines[nextLineWithBlend].blends[bi].lineJ;
 
-            sourceFunc_line(&jnu,&alpha,m,mp[molI].vfac[iphot],g,posn,molJ,lineJ);
+            sourceFunc_line(m[molJ],mp[molI].vfac[iphot],g[posn].mol[molJ],lineJ,&jnu,&alpha);
             tau=alpha*halfFirstDs[iphot];
             calcSourceFn(tau, par, &remnantSnu, &expTau);
             remnantSnu *= jnu*m[molJ].norminv*halfFirstDs[iphot];
diff --git a/src/popsin.c b/src/popsin.c
index c5ba49a..316c769 100644
--- a/src/popsin.c
+++ b/src/popsin.c
@@ -5,6 +5,8 @@
  *  Copyright (C) 2006-2014 Christian Brinch
  *  Copyright (C) 2015 The LIME development team
  *
+***TODO:
+	- Change the definition of the file format so that nmol is now written with the other mol[] scalars.
  */
 
 #include "lime.h"
@@ -15,6 +17,8 @@ popsin(configInfo *par, struct grid **gp, molData **md, int *popsdone){
   FILE *fp;
   int i,j,k,dummyNPart,dummyNTrans;
   double dummy;
+  struct cell *dc=NULL; /* Not used at present. */
+  unsigned long numCells;
 
   if((fp=fopen(par->restart, "rb"))==NULL){
     if(!silent) bail_out("Error reading binary output populations file!");
@@ -70,7 +74,6 @@ popsin(configInfo *par, struct grid **gp, molData **md, int *popsdone){
     (*gp)[i].a4 = NULL;
     (*gp)[i].dens = NULL;
     (*gp)[i].abun = NULL;
-    (*gp)[i].nmol = NULL;
     (*gp)[i].dir = NULL;
     (*gp)[i].neigh = NULL;
     (*gp)[i].w = NULL;
@@ -79,10 +82,10 @@ popsin(configInfo *par, struct grid **gp, molData **md, int *popsdone){
     fread(&(*gp)[i].x, sizeof (*gp)[i].x, 1, fp);
     fread(&(*gp)[i].vel, sizeof (*gp)[i].vel, 1, fp);
     fread(&(*gp)[i].sink, sizeof (*gp)[i].sink, 1, fp);
-    (*gp)[i].nmol=malloc(par->nSpecies*sizeof(double));
-    for(j=0;j<par->nSpecies;j++) fread(&(*gp)[i].nmol[j], sizeof(double), 1, fp);
-    fread(&(*gp)[i].dopb, sizeof (*gp)[i].dopb, 1, fp);
     (*gp)[i].mol=malloc(par->nSpecies*sizeof(struct populations));
+    for(j=0;j<par->nSpecies;j++)
+      fread(&(*gp)[i].mol[j].nmol, sizeof(double), 1, fp);
+    fread(&(*gp)[i].dopb, sizeof (*gp)[i].dopb, 1, fp);
     for(j=0;j<par->nSpecies;j++){
       (*gp)[i].mol[j].pops=malloc(sizeof(double)*(*md)[j].nlev);
       for(k=0;k<(*md)[j].nlev;k++) fread(&(*gp)[i].mol[j].pops[k], sizeof(double), 1, fp);
@@ -100,9 +103,11 @@ popsin(configInfo *par, struct grid **gp, molData **md, int *popsdone){
   }
   fclose(fp);
 
-  qhull(par, *gp);
+  delaunay(DIM, *gp, (unsigned long)par->ncell, 0, &dc, &numCells);
   distCalc(par, *gp);
   getVelosplines(par,*gp);
   *popsdone=1;
+
+  free(dc);
 }
 
diff --git a/src/popsout.c b/src/popsout.c
index a3b42d4..f97b4dd 100644
--- a/src/popsout.c
+++ b/src/popsout.c
@@ -5,6 +5,8 @@
  *  Copyright (C) 2006-2014 Christian Brinch
  *  Copyright (C) 2015 The LIME development team
  *
+***TODO:
+	- Change the definition of the file format so that nmol is now read with the other mol[] scalars.
  */
 
 #include "lime.h"
@@ -25,7 +27,7 @@ popsout(configInfo *par, struct grid *g, molData *m){
   for(j=0;j<par->ncell-par->sinkPoints;j++){
     dens=0.;
     for(l=0;l<par->numDensities;l++) dens+=g[j].dens[l];
-    fprintf(fp,"%e %e %e %e %e %e %d ", g[j].x[0], g[j].x[1], g[j].x[2], dens, g[j].t[0], g[j].nmol[0]/dens, g[j].conv);
+    fprintf(fp,"%e %e %e %e %e %e %d ", g[j].x[0], g[j].x[1], g[j].x[2], dens, g[j].t[0], g[j].mol[0].nmol/dens, g[j].conv);
     for(k=0;k<m[0].nlev;k++) fprintf(fp,"%e ",g[j].mol[0].pops[k]);
     fprintf(fp,"\n");
     //fprintf(fp,"%i %lf %lf %lf %lf %lf %lf %lf %lf\n", g[j].id, g[j].x[0], g[j].x[1], g[j].x[2],  g[j].dens[0], g[j].t[0], g[j].vel[0], g[j].vel[1], g[j].vel[2]);
@@ -76,7 +78,8 @@ binpopsout(configInfo *par, struct grid *g, molData *m){
     fwrite(&g[i].x,  3*sizeof(double),   1, fp);
     fwrite(&g[i].vel,3*sizeof(double),   1, fp);
     fwrite(&g[i].sink, sizeof(int),      1, fp);
-    fwrite(g[i].nmol,  sizeof(double)*par->nSpecies,1, fp);
+    for(j=0;j<par->nSpecies;j++)
+      fwrite(&g[i].mol[j].nmol,  sizeof(double),           1, fp);
     fwrite(&g[i].dopb, sizeof g[i].dopb, 1, fp);
     for(j=0;j<par->nSpecies;j++){
       fwrite(g[i].mol[j].pops,  sizeof(double)*m[j].nlev, 1, fp);
@@ -96,5 +99,3 @@ binpopsout(configInfo *par, struct grid *g, molData *m){
   free(nTrans);
 }
 
-  
-  
diff --git a/src/predefgrid.c b/src/predefgrid.c
index 07bdd90..29b55b0 100644
--- a/src/predefgrid.c
+++ b/src/predefgrid.c
@@ -14,6 +14,9 @@ predefinedGrid(configInfo *par, struct grid *g){
   FILE *fp;
   int i;
   double x,y,z,scale;
+  struct cell *dc=NULL; /* Not used at present. */
+  unsigned long numCells;
+
   gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2);
 #ifdef TEST
   gsl_rng_set(ran,6611304);
@@ -39,13 +42,14 @@ predefinedGrid(configInfo *par, struct grid *g){
 
 
     g[i].sink=0;
-	g[i].t[1]=g[i].t[0];
-		
-	/* This next step needs to be done, even though it looks stupid */
-	g[i].dir=malloc(sizeof(point)*1);
-	g[i].ds =malloc(sizeof(double)*1);
-	g[i].neigh =malloc(sizeof(struct grid *)*1);
-	if(!silent) progressbar((double) i/((double)par->pIntensity-1), 4);	
+    g[i].t[1]=g[i].t[0];
+    g[i].mol[0].nmol=g[i].abun[0]*g[i].dens[0];
+
+    /* This next step needs to be done, even though it looks stupid */
+    g[i].dir=malloc(sizeof(point)*1);
+    g[i].ds =malloc(sizeof(double)*1);
+    g[i].neigh =malloc(sizeof(struct grid *)*1);
+    if(!silent) progressbar((double) i/((double)par->pIntensity-1), 4);	
   }
 
   checkGridDensities(par, g);
@@ -70,13 +74,14 @@ predefinedGrid(configInfo *par, struct grid *g){
   }
   fclose(fp);
 
-  qhull(par,g);
+  delaunay(DIM, g, (unsigned long)par->ncell, 0, &dc, &numCells);
   distCalc(par,g);
   //  getArea(par,g, ran);
   //  getMass(par,g, ran);
   getVelosplines_lin(par,g);
   if(par->gridfile) write_VTK_unstructured_Points(par, g);
   gsl_rng_free(ran);
+  free(dc);
 
   par->numDensities = 1;
 }
diff --git a/src/raythrucells.c b/src/raythrucells.c
new file mode 100644
index 0000000..4bd6724
--- /dev/null
+++ b/src/raythrucells.c
@@ -0,0 +1,526 @@
+/*
+ *  raythrucells.c
+ *  This file is part of LIME, the versatile line modeling engine
+ *
+ *  Copyright (C) 2006-2014 Christian Brinch
+ *  Copyright (C) 2015 The LIME development team
+ *
+ */
+
+#include "lime.h"
+
+/*....................................................................*/
+int followRayThroughDelCells(double x[DIM], double dir[DIM], struct grid *gp\
+  , struct cell *dc, const unsigned long numCells, const double epsilon\
+  , intersectType *entryIntcpt, unsigned long **chainOfCellIds\
+  , intersectType **cellExitIntcpts, int *lenChainPtrs){
+  /*
+The present function follows a ray through a connected, convex set of Delaunay cells and returns information about the chain of cells it passes through. If the ray is found to pass through 1 or more cells, the function returns 0, indicating success; if not, it returns a non-zero value. The chain description consists of three pieces of information: (i) intercept information for the entry face of the first cell encountered; (ii) the IDs of the cells in the chain; (iii) intercept information for the exit face of the ith cell.
+  */
+
+  const int numFaces=DIM+1, maxNumEntryFaces=100;
+  int numEntryFaces, fi, entryFis[maxNumEntryFaces], i, status;
+  faceType face;
+  unsigned long dci, entryDcis[maxNumEntryFaces];//, trialDci;
+  intersectType intcpt, entryIntcpts[maxNumEntryFaces];
+  _Bool *cellVisited=NULL;
+
+  /* Choose a set of starting faces by testing all the 'external' faces of cells which have some. */
+  numEntryFaces = 0;
+  for(dci=0;dci<numCells;dci++){
+    for(fi=0;fi<numFaces;fi++){
+      if(dc[dci].neigh[fi]==NULL){ /* means that this face lies on the outside of the model. */
+        /* Store points for this face: */
+        face = extractFace(gp, dc, dci, fi);
+
+        /* Now calculate the intercept: */
+        intersectLineTriangle(x, dir, face, &intcpt);
+        intcpt.fi = fi; /* Ultimately we need this so we can relate the bary coords for the face back to the Delaunay cell. */
+
+        if(intcpt.orientation<0){ /* it is an entry face. */
+          if(intcpt.collPar+epsilon>0.0){
+            if(numEntryFaces>maxNumEntryFaces){
+              if(!silent) bail_out("Too many entry faces.");
+              exit(1);
+            }
+
+            entryDcis[   numEntryFaces] = dci;
+            entryFis[    numEntryFaces] = fi;
+            entryIntcpts[numEntryFaces] = intcpt;
+            numEntryFaces++;
+          }
+        }
+      }
+    }
+  }
+
+  if(numEntryFaces<=0)
+    return 1;
+
+  *lenChainPtrs = 1024; /* This can be increased within followCellChain(). */
+  *chainOfCellIds  = malloc(sizeof(**chainOfCellIds) *(*lenChainPtrs));
+  *cellExitIntcpts = malloc(sizeof(**cellExitIntcpts)*(*lenChainPtrs));
+  cellVisited = malloc(sizeof(*cellVisited)*numCells);
+  for(dci=0;dci<numCells;dci++)
+    cellVisited[dci] = 0;
+
+  i = 0;
+  status = 1; /* default */
+  while(i<numEntryFaces && status>0){
+    status = buildRayCellChain(x, dir, gp, dc, &cellVisited, entryDcis[i], entryFis[i], 0, 0, epsilon, chainOfCellIds, cellExitIntcpts, lenChainPtrs);
+    i++;
+  }
+
+  if(status==0)
+    *entryIntcpt = entryIntcpts[i-1];
+  /* Note that the order of the bary coords, and the value of fi, are with reference to the vertx list of the _entered_ cell. This can't of course be any other way, because this ray enters this first cell from the exterior of the model, where there are no cells. For all the intersectType objects in the list cellExitIntcpts, the bary coords etc are with reference to the exited cell. */
+
+  free(cellVisited);
+
+  return status;
+}
+
+/*....................................................................*/
+int buildRayCellChain(double x[DIM], double dir[DIM], struct grid *gp\
+  , struct cell *dc, _Bool **cellVisited, unsigned long dci, int entryFaceI\
+  , int levelI, int nCellsInChain, const double epsilon\
+  , unsigned long **chainOfCellIds, intersectType **cellExitIntcpts, int *lenChainPtrs){
+  /*
+This function is designed to follow a ray (defined by a starting locus 'x' and a direction vector 'dir') through a convex connected set of Delaunay cells. The function returns an integer status value directly, and two lists (plus their common length) via the argument interface: chainOfCellIds and cellExitIntcpts. Taken together, these lists define a chain of Delaunay cells traversed by the ray.
+
+The task of determining which cells are traversed by the ray is simple in principle, but complications arise in computational practice due to the finite precision of floating-point calculations. Where the ray crosses a cell face near to one of its edges, numerical calculation of the 'impact parameter' may return an answer which is erroneous either way: i.e., a ray which just misses a face may be reported as hitting it, and vice versa. To deal with this, a distinction is made between impacts which are (i) 'good', that is far from any face edge; (ii) 'bad', i.e. clearly missing the face; and (iii) 'marginal', i.e. closer to the edge than some preset cutoff which is represented in the argument list by the number 'epsilon'. Note that a tally is kept of those cells which have already been tested for the current ray, and any face which abuts a neighbouring cell which has been visited already will be flagged as 'bad'.
+
+The function therefore looks at all the exit faces of the cell and sorts them into these three categories. How it proceeds then depends on the relative numbers of each, as described below.
+
+	- If there is more than 1 good exit face, an exception is generated. This is a sign that epsilon has been chosen with too small a value.
+
+	- If there is just 1 good exit face, the marginal ones are ignored. The current cell data are appended to the chain and the cell abutting the exit face becomes the new working cell.
+
+	- If there are no good exit faces, we look at the marginal ones. If there is only 1 of these, we loop as before. If there are more than 1, the function is called recursively for each cell on the far side of an exit face.
+
+Thus there are two alternate modes of operation within the function: a straightforward loop along the cells in a single chain, which will continue so long as there is only a single exit face of type either 'good' or 'marginal'; and a recursive launch of the function at a fork in the chain into each of its several possible branches.
+
+The function terminates under the following conditions:
+	- It detects that the edge of the model is reached (returns success).
+	- There are no exit faces (returns failure).
+	- There are no good exit faces and either
+		* all of the recursive calls to marginal faces have been unsuccessful (returns failure), or
+		* one of these has been successful (returns success).
+
+At a successful termination, therefore, details of all the cells to the edge of the model are correctly stored in chainOfCellIds and cellExitIntcpts, and the number of these cells is returned in lenChainPtrs.
+
+***** Note that it is assumed here that a mis-indentification of the actual cell traversed by a ray will not ultimately matter to the external calling routine. This is only reasonable if whatever function or property is being sampled by the ray does not vary in a step function across cell boundaries. *****
+  */
+
+  const int numFaces=DIM+1;
+  _Bool followingSingleChain;
+  const int bufferSize=1024;
+  int numGoodExits, numMarginalExits, fi, goodExitFis[numFaces], marginalExitFis[numFaces], exitFi, i, status, newEntryFaceI;
+  faceType face;
+  intersectType intcpt[numFaces];
+
+  followingSingleChain = 1; /* default */
+  do{ /* Follow the chain through 'good' cells, i.e. ones for which entry and exit are nicely distant from face edges. (Here we also follow marginal exits if there are no good ones, and only 1 marginal one.) */
+    (*cellVisited)[dci] = 1;
+
+    /* Store the current cell ID (we leave storing the exit face for later, when we know what it is). (If there is not enough room in chainOfCellIds and cellExitIntcpts, realloc them to new value of lenChainPtrs.) */
+    if(nCellsInChain>=(*lenChainPtrs)){
+      *lenChainPtrs += bufferSize;
+      *chainOfCellIds  = realloc(*chainOfCellIds,  sizeof(**chainOfCellIds) *(*lenChainPtrs));
+      *cellExitIntcpts = realloc(*cellExitIntcpts, sizeof(**cellExitIntcpts)*(*lenChainPtrs));
+    }
+    (*chainOfCellIds)[nCellsInChain] = dci;
+
+    /* calculate num good and bad exits */
+    numGoodExits = 0;
+    numMarginalExits = 0;
+    for(fi=0;fi<numFaces;fi++){
+      if(fi!=entryFaceI && (dc[dci].neigh[fi]==NULL || !(*cellVisited)[dc[dci].neigh[fi]->id])){
+        /* Store points for this face: */
+        face = extractFace(gp, dc, dci, fi);
+
+        /* Now calculate the intercept: */
+        intersectLineTriangle(x, dir, face, &intcpt[fi]);
+        intcpt[fi].fi = fi; /* Ultimately we need this so we can relate the bary coords for the face back to the Delaunay cell. */
+
+        if(intcpt[fi].orientation>0){ /* it is an exit face. */
+          if(intcpt[fi].collPar-epsilon>0.0){
+            goodExitFis[numGoodExits] = fi;
+            numGoodExits++;
+          }else if (intcpt[fi].collPar+epsilon>0.0){
+            marginalExitFis[numMarginalExits] = fi;
+            numMarginalExits++;
+          }
+        }
+      }
+    }
+
+    if(numGoodExits>1){
+      if(!silent) bail_out("Some sort of bug: more than 1 firm candidate found for ray exit from cell.");
+      exit(1);
+
+    }else if(numGoodExits==1 || numMarginalExits==1){
+      if(numGoodExits==1)
+        exitFi = goodExitFis[0];
+      else
+        exitFi = marginalExitFis[0];
+
+      /* Store the exit face details: */
+      (*cellExitIntcpts)[nCellsInChain] = intcpt[exitFi];
+
+      nCellsInChain++;
+
+      if(dc[dci].neigh[exitFi]==NULL){ /* Signals that we have reached the edge of the model. */
+        /* Realloc the ptrs to their final sizes: */
+        *chainOfCellIds  = realloc(*chainOfCellIds,  sizeof(**chainOfCellIds) *nCellsInChain);
+        *cellExitIntcpts = realloc(*cellExitIntcpts, sizeof(**cellExitIntcpts)*nCellsInChain);
+        *lenChainPtrs = nCellsInChain;
+
+        return 0;
+      }
+
+      entryFaceI = getNewEntryFaceI(dci, *(dc[dci].neigh[exitFi]));
+      dci = dc[dci].neigh[exitFi]->id;
+
+    } else
+      followingSingleChain = 0;
+  } while(followingSingleChain);
+
+  /* Now we have run out of good (or at least single) exit-face options, let's try the marginal ones. */
+
+  if(numMarginalExits<1)
+    return 1; /* Unsuccessful end of this chain. */
+
+  /* If we have got to this point, we must have numMarginalExits>1; thus we have a fork in the chain, and must explore each branch. We recurse here because a recursive scheme is the best way to do that. */
+
+  i = 0;
+  status = 1; /* default */
+  while(i<numMarginalExits && status>0){
+    exitFi = marginalExitFis[i];
+    (*cellExitIntcpts)[nCellsInChain] = intcpt[exitFi];
+
+    if(dc[dci].neigh[exitFi]==NULL){
+      status = 0;
+    }else{
+      newEntryFaceI = getNewEntryFaceI(dci, *(dc[dci].neigh[exitFi]));
+
+      /* Now we dive into the branch: */
+      status = buildRayCellChain(x, dir, gp, dc, cellVisited, dc[dci].neigh[exitFi]->id\
+        , newEntryFaceI, levelI+1, nCellsInChain+1, epsilon\
+        , chainOfCellIds, cellExitIntcpts, lenChainPtrs);
+    }
+    i++;
+  }
+
+  return status;
+}
+
+/*....................................................................*/
+int getNewEntryFaceI(const unsigned long dci, const struct cell newCell){
+  /* Finds the index of the old cell in the face list of the new cell. */
+
+  const int numFaces=DIM+1;
+  _Bool matchFound = 0;
+  int ffi = 0, newEntryFaceI;
+
+  while(ffi<numFaces && matchFound==0){
+    if(newCell.neigh[ffi]!=NULL && newCell.neigh[ffi]->id==dci){
+      matchFound = 1;
+      newEntryFaceI = ffi;
+    }
+    ffi++;
+  }
+
+  /* Sanity check: */
+  if(!matchFound){
+    bail_out("Cannot find old cell ID in new cell data.");
+    exit(1);
+  }
+
+  return newEntryFaceI;
+}
+
+/*....................................................................*/
+faceType extractFace(struct grid *gp, struct cell *dc, const unsigned long dci, const int fi){
+  /* Given a simplex dc[dci] and the face index (in the range {0...DIM}) fi, this returns the desired information about that face. Note that the ordering of the elements of face.r[] is the same as the ordering of the vertices of the simplex, dc[dci].vertx[]; just the vertex fi is omitted.
+
+Note that the element 'centre' of the faceType struct is mean to contain the spatial coordinates of the centre of the simplex, not of the face. This is designed to facilitate orientation of the face and thus to help determine whether rays which cross it are entering or exiting the simplex.
+ */
+
+  const int numFaces=DIM+1;
+  int vi, vvi, di;
+  struct grid point;
+  faceType face;
+
+  unsigned long gi;
+
+  vvi = 0;
+  for(vi=0;vi<numFaces;vi++){
+    if(vi!=fi){
+      gi = dc[dci].vertx[vi]->id;
+      point = gp[gi];
+      for(di=0;di<DIM;di++){
+        face.r[vvi][di] = point.x[di];
+      }
+      vvi++;
+    }
+  }
+
+  for(di=0;di<DIM;di++){
+    face.centre[di] = dc[dci].centre[di];
+  }
+
+  return face;
+}
+
+/*....................................................................*/
+void intersectLineTriangle(double x[3], double dir[3], faceType face\
+  , intersectType *intcpt){
+  /*
+This function calculates the intersection between a line in 3 space and a triangle oriented in that space. The intersection point may be expressed as
+
+	px = x + a*dir
+
+where px, x and dir are vectors and 'a' is a scalar. The routine returns the following information:
+	- The value of 'a'.
+	- The so-called barycentric coordinates (BC) of px in the triangle.
+
+This routine works best when the sides of the triangle are not too disparate in size.
+
+Notes:
+	* The line and the triangle may be parallel. In this case, the function returns 1; otherwise 0.
+
+	* There is, of course, no guarantee that the line actually intersects the triangle, even if the line and the triangular plane intersect. There are also borderline cases, i.e. where the line passes close to a vertex, in which an exact calculation would show that the intersection occurs (or doesn't occur), but the imprecise computed value claims that it doesn't (or does). Intersection may be judged via the values of the barycentric coordinates: if these all fall in the interval [0,1], the line intersects the triangle; if one of the BC is negative, it doesn't.
+
+***** NOTE ***** that the intersectType now includes a field for the index of the face. This is not known to or supplied by the present function, so it sets it to -1 just to avoid an uninitialized value.
+  */
+
+  const int numDims=3;
+  double norm[numDims], normDotDx, numerator, px2D[numDims-1], mat[numDims-1][numDims-1], vec[numDims-1], det;
+  int j, k, di;
+  double testSumForCW=0.0;
+  triangle2D face2D;
+
+  intcpt->fi = -1;
+
+  /* First calculate a normal vector to the triangle from the cross product. Since we don't know a priori whether the vertices of the face are listed CW or ACW seen from inside the cell, we will work this out by dotting the norm vector with a vector from the centre of the cell to vertex 0.
+  */
+  for(di=0;di<numDims;di++){
+    j = (di+1)%numDims;
+    k = (di+2)%numDims;
+    norm[di] = (face.r[1][j]-face.r[0][j])*(face.r[2][k]-face.r[0][k])\
+             - (face.r[1][k]-face.r[0][k])*(face.r[2][j]-face.r[0][j]);
+
+    testSumForCW += norm[di]*(face.r[0][di] - face.centre[di]);
+  }
+
+  if(testSumForCW<0.0){
+    for(di=0;di<numDims;di++)
+      norm[di] *= -1.0;
+  }
+
+  /* Calculate the scalar (or dot) product between norm and dir.
+  */
+  normDotDx = norm[0]*dir[0] + norm[1]*dir[1] + norm[2]*dir[2];
+
+  if(normDotDx>0.0){ /* it is an exit face. */
+    intcpt->orientation = 1;
+  }else if(normDotDx<0.0){ /* it is an entry face. */
+    intcpt->orientation = -1;
+  }else{ /* normDotDx==0.0, i.e. line and plane are parallel. */
+    intcpt->orientation = 0;
+    intcpt->bary[0] = 0.0;
+    intcpt->bary[1] = 0.0;
+    intcpt->bary[2] = 0.0;
+    intcpt->dist    = 0.0;
+    intcpt->collPar = 0.0;
+
+    return;
+  }
+
+  /* If we've got to here, we can be sure the line and plane are not parallel, and that we therefore expect meaningful results. */
+
+  numerator = norm[0]*(face.r[0][0] - x[0])\
+            + norm[1]*(face.r[0][1] - x[1])\
+            + norm[2]*(face.r[0][2] - x[2]);
+
+  intcpt->dist = numerator/normDotDx;
+
+  /* In order to calculate the barycentric coordinates, we need to set up a 2D coordinate basis in the plane of the triangle. I'll define the X axis as 10^, where for short I use the generic notation ji_ to mean the vector (tri[j]-tri[i]), and ji^ to mean the unit vector in the same direction. The Y axis is the unit vector in the direction (20_ - (20_.10^)*10^).
+  */
+  face2D = calcTriangle2D(face);
+
+  /* Now we want to express the intersection point in these coordinates:
+  */
+  px2D[0] = (x[0] + intcpt->dist*dir[0] - face.r[0][0])*face2D.xAxis[0]\
+          + (x[1] + intcpt->dist*dir[1] - face.r[0][1])*face2D.xAxis[1]\
+          + (x[2] + intcpt->dist*dir[2] - face.r[0][2])*face2D.xAxis[2];
+  px2D[1] = (x[0] + intcpt->dist*dir[0] - face.r[0][0])*face2D.yAxis[0]\
+          + (x[1] + intcpt->dist*dir[1] - face.r[0][1])*face2D.yAxis[1]\
+          + (x[2] + intcpt->dist*dir[2] - face.r[0][2])*face2D.yAxis[2];
+
+  /*
+The barycentric coordinates (L0,L1,L2) are given by
+
+	(tri2D[0][0]-tri2D[2][0]  tri2D[1][0]-tri2D[2][0]) (L0)   (px2D[0]-tri2D[2][0])
+	(                                                )*(  ) = (                   ),
+	(tri2D[0][1]-tri2D[2][1]  tri2D[1][1]-tri2D[2][1]) (L1)   (px2D[1]-tri2D[2][1])
+
+with L2 = 1 - L0 - L1.
+  */
+  mat[0][0] = face2D.r[0][0]-face2D.r[2][0];
+  mat[0][1] = face2D.r[1][0]-face2D.r[2][0];
+  mat[1][0] = face2D.r[0][1]-face2D.r[2][1];
+  mat[1][1] = face2D.r[1][1]-face2D.r[2][1];
+  vec[0] = px2D[0]-face2D.r[2][0];
+  vec[1] = px2D[1]-face2D.r[2][1];
+  det = mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
+  /*** We're assuming that the triangle is not pathological, i.e that det!=0. */
+  intcpt->bary[0] = ( mat[1][1]*vec[0] - mat[0][1]*vec[1])/det;
+  intcpt->bary[1] = (-mat[1][0]*vec[0] + mat[0][0]*vec[1])/det;
+  intcpt->bary[2] = 1.0 - intcpt->bary[0] - intcpt->bary[1];
+
+  di = 0;
+  if(intcpt->bary[di] < 0.5)
+    intcpt->collPar = intcpt->bary[di];
+  else
+    intcpt->collPar = 1.0 - intcpt->bary[di];
+
+  for(di=1;di<numDims;di++){
+    if(intcpt->bary[di] < 0.5){
+      if(intcpt->bary[di] < intcpt->collPar)
+        intcpt->collPar = intcpt->bary[di];
+    }else{ /* intcpt->bary[di]>=0.5 */
+      if(1.0 - intcpt->bary[di] < intcpt->collPar)
+        intcpt->collPar = 1.0 - intcpt->bary[di];
+    }
+  }
+
+  return;
+}
+
+/*....................................................................*/
+triangle2D calcTriangle2D(faceType face){
+/**** all this could be precalculated (norm and mat too). */
+
+  triangle2D face2D;
+  double lengthX, lengthY, dotResult;
+
+  face2D.xAxis[0] = face.r[1][0]-face.r[0][0];
+  face2D.xAxis[1] = face.r[1][1]-face.r[0][1];
+  face2D.xAxis[2] = face.r[1][2]-face.r[0][2];
+  lengthX = sqrt(face2D.xAxis[0]*face2D.xAxis[0]\
+               + face2D.xAxis[1]*face2D.xAxis[1]\
+               + face2D.xAxis[2]*face2D.xAxis[2]);
+  face2D.xAxis[0] /= lengthX;
+  face2D.xAxis[1] /= lengthX;
+  face2D.xAxis[2] /= lengthX;
+
+  dotResult = face2D.xAxis[0]*(face.r[2][0]-face.r[0][0])\
+            + face2D.xAxis[1]*(face.r[2][1]-face.r[0][1])\
+            + face2D.xAxis[2]*(face.r[2][2]-face.r[0][2]);
+
+  face2D.yAxis[0] = (face.r[2][0]-face.r[0][0]) - dotResult*face2D.xAxis[0];
+  face2D.yAxis[1] = (face.r[2][1]-face.r[0][1]) - dotResult*face2D.xAxis[1];
+  face2D.yAxis[2] = (face.r[2][2]-face.r[0][2]) - dotResult*face2D.xAxis[2];
+  lengthY = sqrt(face2D.yAxis[0]*face2D.yAxis[0]\
+               + face2D.yAxis[1]*face2D.yAxis[1]\
+               + face2D.yAxis[2]*face2D.yAxis[2]);
+  face2D.yAxis[0] /= lengthY;
+  face2D.yAxis[1] /= lengthY;
+  face2D.yAxis[2] /= lengthY;
+
+  /* The expressions for the triangle vertices in this basis are simple:
+  */
+  face2D.r[0][0] = 0.0;
+  face2D.r[0][1] = 0.0;
+  face2D.r[1][0] = lengthX;
+  face2D.r[1][1] = 0.0;
+  face2D.r[2][0] = dotResult;
+  face2D.r[2][1] = lengthY;
+
+  return face2D;
+}
+
+/*....................................................................*/
+void doBaryInterp(const intersectType intcpt, struct grid *gp\
+  , struct gAuxType *gAux, double xCmpntsRay[3], unsigned long gis[3]\
+  , molData *md, const int numMols, gridInterp *gip){
+
+  int di, molI, levelI, lineI;
+
+  (*gip).xCmpntRay = intcpt.bary[0]*xCmpntsRay[0]\
+                   + intcpt.bary[1]*xCmpntsRay[1]\
+                   + intcpt.bary[2]*xCmpntsRay[2];
+  for(di=0;di<DIM;di++){
+    (*gip).x[di] = intcpt.bary[0]*gp[gis[0]].x[di]\
+                 + intcpt.bary[1]*gp[gis[1]].x[di]\
+                 + intcpt.bary[2]*gp[gis[2]].x[di];
+  }
+
+  for(di=0;di<3;di++){ /* 3 not DIM, because a B field only makes sense in 3 dimensions. */
+/* ****** Maybe some test of DIM==3?? Would need to be systematic across the entire code..? */
+    (*gip).B[di] = intcpt.bary[0]*gp[gis[0]].B[di]\
+                 + intcpt.bary[1]*gp[gis[1]].B[di]\
+                 + intcpt.bary[2]*gp[gis[2]].B[di];
+  }
+
+  for(molI=0;molI<numMols;molI++){
+    (*gip).mol[molI].binv = intcpt.bary[0]*gp[gis[0]].mol[molI].binv\
+                          + intcpt.bary[1]*gp[gis[1]].mol[molI].binv\
+                          + intcpt.bary[2]*gp[gis[2]].mol[molI].binv;
+
+    for(levelI=0;levelI<md[molI].nlev;levelI++){
+      (*gip).mol[molI].specNumDens[levelI]\
+        = intcpt.bary[0]*gAux[gis[0]].mol[molI].specNumDens[levelI]\
+        + intcpt.bary[1]*gAux[gis[1]].mol[molI].specNumDens[levelI]\
+        + intcpt.bary[2]*gAux[gis[2]].mol[molI].specNumDens[levelI];
+    }
+
+    for(lineI=0;lineI<md[molI].nline;lineI++){
+      (*gip).mol[molI].dust[lineI] = intcpt.bary[0]*gp[gis[0]].mol[molI].dust[lineI]\
+                                   + intcpt.bary[1]*gp[gis[1]].mol[molI].dust[lineI]\
+                                   + intcpt.bary[2]*gp[gis[2]].mol[molI].dust[lineI];
+
+      (*gip).mol[molI].knu[lineI]  = intcpt.bary[0]*gp[gis[0]].mol[molI].knu[lineI]\
+                                   + intcpt.bary[1]*gp[gis[1]].mol[molI].knu[lineI]\
+                                   + intcpt.bary[2]*gp[gis[2]].mol[molI].knu[lineI];
+    }
+  }
+}
+
+/*....................................................................*/
+void doSegmentInterp(gridInterp gips[3], const int iA, molData *md\
+  , const int numMols, const double oneOnNumSegments, const int si){
+
+  const double fracA = (si + 0.5)*oneOnNumSegments, fracB = 1.0 - fracA;
+  const int iB = 1 - iA;
+  int di, molI, levelI, lineI;
+
+  gips[2].xCmpntRay = fracA*gips[iB].xCmpntRay + fracB*gips[iA].xCmpntRay; /* This does not seem to be used. */
+
+  for(di=0;di<DIM;di++){
+    gips[2].x[di] = fracA*gips[iB].x[di] + fracB*gips[iA].x[di];
+  }
+  for(di=0;di<3;di++){ /* 3 not DIM, because a B field only makes sense in 3 dimensions. */
+    gips[2].B[di] = fracA*gips[iB].B[di] + fracB*gips[iA].B[di];
+  }
+
+  for(molI=0;molI<numMols;molI++){
+    gips[2].mol[molI].binv = fracA*gips[iB].mol[molI].binv + fracB*gips[iA].mol[molI].binv;
+
+    for(levelI=0;levelI<md[molI].nlev;levelI++){
+      gips[2].mol[molI].specNumDens[levelI] = fracA*gips[iB].mol[molI].specNumDens[levelI]\
+                                            + fracB*gips[iA].mol[molI].specNumDens[levelI];
+    }
+
+    for(lineI=0;lineI<md[molI].nline;lineI++){
+      gips[2].mol[molI].dust[lineI] = fracA*gips[iB].mol[molI].dust[lineI]\
+                                    + fracB*gips[iA].mol[molI].dust[lineI];
+
+      gips[2].mol[molI].knu[ lineI] = fracA*gips[iB].mol[molI].knu[ lineI]\
+                                    + fracB*gips[iA].mol[molI].knu[ lineI];
+    }
+  }
+}
+
diff --git a/src/raytrace.c b/src/raytrace.c
index ae679ea..d1a37f8 100644
--- a/src/raytrace.c
+++ b/src/raytrace.c
@@ -5,24 +5,27 @@
  *  Copyright (C) 2006-2014 Christian Brinch
  *  Copyright (C) 2016 The LIME development team
  *
+TODO:
+  - In raytrace(), look at rearranging the code to do the qhull step before choosing the rays. This would allow cells with all vertices outside the image boundaries to be excluded. If the image is much smaller than the model, this could lead to significant savings in time. The only downside might be memory useage...
+  - We should not be multiplying remnantSnu by md[0].norminv, but to make it correct would be pretty fiddly.
  */
 
 #include "lime.h"
 
 
-void
-velocityspline2(double x[3], double dx[3], double ds, double binv, double deltav, double *vfac){
+/*....................................................................*/
+void calcLineAmpSample(const double x[3], const double dx[3], const double ds\
+  , const double binv, double *projVels, const int nSteps\
+  , const double oneOnNSteps, const double deltav, double *vfac){
   /*
 The bulk velocity of the model material can vary significantly with position, thus so can the value of the line-shape function at a given frequency and direction. The present function calculates 'vfac', an approximate average of the line-shape function along a path of length ds in the direction of the line of sight.
   */
-  int i,steps=10;
-  double v,d,val,vel[3];
+  int i;
+  double v,val;
 
   *vfac=0.;
-  for(i=0;i<steps;i++){
-    d=i*ds/steps;
-    velocity(x[0]+(dx[0]*d),x[1]+(dx[1]*d),x[2]+(dx[2]*d),vel);
-    v=deltav-veloproject(dx,vel); /* veloproject returns the component of the local bulk velocity in the direction of dx, whereas deltav is the recession velocity of the channel we are interested in (corrected for bulk source velocity and line displacement from the nominal frequency). Remember also that, since dx points away from the observer, positive values of the projected velocity also represent recessions. Line centre occurs when v==0, i.e. when deltav==veloproject(dx,vel). That is the reason for the subtraction here. */
+  for(i=0;i<nSteps;i++){
+    v = deltav - projVels[i]; /* projVels contains the component of the local bulk velocity in the direction of dx, whereas deltav is the recession velocity of the channel we are interested in (corrected for bulk source velocity and line displacement from the nominal frequency). Remember also that, since dx points away from the observer, positive values of the projected velocity also represent recessions. Line centre occurs when v==0, i.e. when deltav==veloproject(dx,vel). That is the reason for the subtraction here. */
     val=fabs(v)*binv;
     if(val <=  2500.){
 #ifdef FASTEXP
@@ -32,11 +35,31 @@ The bulk velocity of the model material can vary significantly with position, th
 #endif
     }
   }
-  *vfac=*vfac/steps;
+  *vfac *= oneOnNSteps;
   return;
 }
 
+/*....................................................................*/
+void calcLineAmpInterp(const double projVelRay, const double binv\
+  , const double deltav, double *vfac){
+  /*
+The bulk velocity of the model material can vary significantly with position, thus so can the value of the line-shape function at a given frequency and direction. The present function calculates 'vfac', an approximate average of the line-shape function along a path of length ds in the direction of the line of sight.
+  */
+  double v,val;
 
+  v = deltav - projVelRay; /* projVelRay is the component of the local bulk velocity in the direction of the ray, whereas deltav is the recession velocity of the channel we are interested in (corrected for bulk source velocity and line displacement from the nominal frequency). Remember also that, since the ray points away from the observer, positive values of the projected velocity also represent recessions. Line centre occurs when v==0, i.e. when deltav==projVelRay. That is the reason for the subtraction here. */
+  val = fabs(v)*binv;
+  if(val <=  2500.){
+#ifdef FASTEXP
+    *vfac = FastExp(val*val);
+#else
+    *vfac =   exp(-(val*val));
+#endif
+  }else
+    *vfac = 0.0;
+}
+
+/*....................................................................*/
 void
 line_plane_intersect(struct grid *g, double *ds, int posn, int *nposn, double *dx, double *x, double cutoff){
   /*
@@ -66,18 +89,23 @@ This function returns ds as the (always positive-valued) distance between the pr
   if(*nposn==-1) *nposn=posn;
 }
 
-
+/*....................................................................*/
 void
-traceray(rayData ray, int cmbMolI, int cmbLineI, int im, configInfo *par, struct grid *gp, molData *md, image *img, double cutoff){
+traceray(rayData ray, const int cmbMolI, const int cmbLineI, const int im\
+  , configInfo *par, struct grid *gp, molData *md, image *img, struct gAuxType *gAux\
+  , const int nlinetot, int *allLineMolIs, int *allLineLineIs, const double cutoff\
+  , const int nSteps, const double oneOnNSteps){
   /*
 For a given image pixel position, this function evaluates the intensity of the total light emitted/absorbed along that line of sight through the (possibly rotated) model. The calculation is performed for several frequencies, one per channel of the output image.
 
 Note that the algorithm employed here is similar to that employed in the function photon() which calculates the average radiant flux impinging on a grid cell: namely the notional photon is started at the side of the model near the observer and 'propagated' in the receding direction until it 'reaches' the far side. This is rather non-physical in conception but it makes the calculation easier.
   */
   const int stokesIi=0;
-  int ichan,posn,nposn,i,molI,lineI;
-  double vfac=0.,x[3],dx[3],vThisChan;
-  double deltav,ds,dist2,ndist2,xp,yp,zp,col,lineRedShift,jnu,alpha,remnantSnu,dtau,expDTau,snu_pol[3];
+  int ichan,stokesId,di,i,posn,nposn,polMolI,polLineI,iline,molI,lineI;
+  double xp,yp,zp,x[DIM],dx[DIM],dist2,ndist2,col,ds,snu_pol[3],dtau;
+  double contJnu,contAlpha,jnu,alpha,lineRedShift,vThisChan,deltav,vfac=0.;
+  double remnantSnu,expDTau,brightnessIncrement;
+  double projVels[nSteps],d,vel[DIM];
 
   for(ichan=0;ichan<img[im].nchan;ichan++){
     ray.tau[ichan]=0.0;
@@ -87,144 +115,393 @@ Note that the algorithm employed here is similar to that employed in the functio
   xp=ray.x;
   yp=ray.y;
 
-  if((xp*xp+yp*yp)/par->radiusSqu <= 1 ) {
-    zp=-sqrt(par->radiusSqu-(xp*xp+yp*yp)); /* There are two points of intersection between the line of sight and the spherical model surface; this is the Z coordinate (in the unrotated frame) of the one nearer to the observer. */
+  /* The model is circular in projection. We only follow the ray if it will intersect the model.
+  */
+  if((xp*xp+yp*yp)>par->radiusSqu)
+    return;
+
+  zp=-sqrt(par->radiusSqu-(xp*xp+yp*yp)); /* There are two points of intersection between the line of sight and the spherical model surface; this is the Z coordinate (in the unrotated frame) of the one nearer to the observer. */
 
-    /* Rotate the line of sight as desired. */
-    for(i=0;i<3;i++){
-      x[i]=xp*img[im].rotMat[i][0] + yp*img[im].rotMat[i][1] + zp*img[im].rotMat[i][2];
-      dx[i]= img[im].rotMat[i][2]; /* This points away from the observer. */
+  /* Rotate the line of sight as desired. */
+  for(di=0;di<DIM;di++){
+    x[di]=xp*img[im].rotMat[di][0] + yp*img[im].rotMat[di][1] + zp*img[im].rotMat[di][2];
+    dx[di]= img[im].rotMat[di][2]; /* This points away from the observer. */
+  }
+
+  /* Find the grid point nearest to the starting x. */
+  i=0;
+  dist2=(x[0]-gp[i].x[0])*(x[0]-gp[i].x[0]) + (x[1]-gp[i].x[1])*(x[1]-gp[i].x[1]) + (x[2]-gp[i].x[2])*(x[2]-gp[i].x[2]);
+  posn=i;
+  for(i=1;i<par->ncell;i++){
+    ndist2=(x[0]-gp[i].x[0])*(x[0]-gp[i].x[0]) + (x[1]-gp[i].x[1])*(x[1]-gp[i].x[1]) + (x[2]-gp[i].x[2])*(x[2]-gp[i].x[2]);
+    if(ndist2<dist2){
+      posn=i;
+      dist2=ndist2;
     }
+  }
+
+  col=0;
+  do{
+    ds=-2.*zp-col; /* This default value is chosen to be as large as possible given the spherical model boundary. */
+    nposn=-1;
+    line_plane_intersect(gp,&ds,posn,&nposn,dx,x,cutoff); /* Returns a new ds equal to the distance to the next Voronoi face, and nposn, the ID of the grid cell that abuts that face. */
+
+    if(par->polarization){
+      polMolI = 0; /****** Always?? */
+      polLineI = 0; /****** Always?? */
+      sourceFunc_pol(gp[posn].B, gAux[posn].mol[polMolI], polLineI, img[im].rotMat, snu_pol, &alpha);
+      dtau=alpha*ds;
+      calcSourceFn(dtau, par, &remnantSnu, &expDTau);
+      remnantSnu *= md[polMolI].norminv*ds;
+
+      for(stokesId=0;stokesId<img[im].nchan;stokesId++){ /* Loop over I, Q and U */
+#ifdef FASTEXP
+        brightnessIncrement = FastExp(ray.tau[stokesId])*remnantSnu*snu_pol[stokesId];
+#else
+        brightnessIncrement =    exp(-ray.tau[stokesId])*remnantSnu*snu_pol[stokesId];
+#endif
+        ray.intensity[stokesId] += brightnessIncrement;
+        ray.tau[stokesId]+=dtau; //**** But this will be the same for I, Q or U.
+      }
+    } else {
+      if(!par->doPregrid){
+        for(i=0;i<nSteps;i++){
+          d = i*ds*oneOnNSteps;
+          velocity(x[0]+(dx[0]*d),x[1]+(dx[1]*d),x[2]+(dx[2]*d),vel);
+          projVels[i] = veloproject(dx,vel);
+        }
+      }
+
+      /* Calculate first the continuum stuff because it is the same for all channels:
+      */
+      contJnu = 0.0;
+      contAlpha = 0.0;
+      sourceFunc_cont_raytrace(gAux[posn].mol[cmbMolI], cmbLineI, &contJnu, &contAlpha);
+
+      for(ichan=0;ichan<img[im].nchan;ichan++){
+        jnu = contJnu;
+        alpha = contAlpha;
+        vThisChan = (ichan-(img[im].nchan-1)*0.5)*img[im].velres; /* Consistent with the WCS definition in writefits(). */
+
+        if(img[im].doline){
+          for(molI=0;molI<par->nSpecies;molI++){
+            for(lineI=0;lineI<md[molI].nline;lineI++){
+              if(md[molI].freq[lineI] > img[im].freq-img[im].bandwidth*0.5\
+              && md[molI].freq[lineI] < img[im].freq+img[im].bandwidth*0.5){
+                /* Calculate the red shift of the transition wrt to the frequency specified for the image.
+                */
+                if(img[im].trans > -1){
+                  lineRedShift=(md[molI].freq[img[im].trans]-md[molI].freq[lineI])/md[molI].freq[img[im].trans]*CLIGHT;
+                } else {
+                  lineRedShift=(img[im].freq-md[molI].freq[lineI])/img[im].freq*CLIGHT;
+                }
+
+                deltav = vThisChan - img[im].source_vel - lineRedShift;
+                /* Line centre occurs when deltav = the recession velocity of the radiating material. Explanation of the signs of the 2nd and 3rd terms on the RHS: (i) A bulk source velocity (which is defined as >0 for the receding direction) should be added to the material velocity field; this is equivalent to subtracting it from deltav, as here. (ii) A positive value of lineRedShift means the line is red-shifted wrt to the frequency specified for the image. The effect is the same as if the line and image frequencies were the same, but the bulk recession velocity were higher. lineRedShift should thus be added to the recession velocity, which is equivalent to subtracting it from deltav, as here. */
+
+                /* Calculate an approximate average line-shape function at deltav within the Voronoi cell. */
+                if(!par->doPregrid)
+                  calcLineAmpSample(x,dx,ds,gp[posn].mol[molI].binv,projVels,nSteps,oneOnNSteps,deltav,&vfac);
+                else
+                  vfac = gaussline(deltav-veloproject(dx,gp[posn].vel),gp[posn].mol[molI].binv);
 
-    /* Find the grid point nearest to the starting x. */
-    i=0;
-    dist2=(x[0]-gp[i].x[0])*(x[0]-gp[i].x[0]) + (x[1]-gp[i].x[1])*(x[1]-gp[i].x[1]) + (x[2]-gp[i].x[2])*(x[2]-gp[i].x[2]);
-    posn=i;
-    for(i=1;i<par->ncell;i++){
-      ndist2=(x[0]-gp[i].x[0])*(x[0]-gp[i].x[0]) + (x[1]-gp[i].x[1])*(x[1]-gp[i].x[1]) + (x[2]-gp[i].x[2])*(x[2]-gp[i].x[2]);
-      if(ndist2<dist2){
-        posn=i;
-        dist2=ndist2;
+                /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. */
+                sourceFunc_line_raytrace(md[molI],vfac,gAux[posn].mol[molI],lineI,&jnu,&alpha);
+              }
+            }
+          }
+        }
+
+        dtau=alpha*ds;
+//???          if(dtau < -30) dtau = -30; // as in photon()?
+        calcSourceFn(dtau, par, &remnantSnu, &expDTau);
+        remnantSnu *= jnu*md[0].norminv*ds;
+#ifdef FASTEXP
+        brightnessIncrement = FastExp(ray.tau[ichan])*remnantSnu;
+#else
+        brightnessIncrement =    exp(-ray.tau[ichan])*remnantSnu;
+#endif
+        ray.intensity[ichan] += brightnessIncrement;
+        ray.tau[ichan]+=dtau;
       }
     }
 
-    col=0;
-    do{
-      ds=-2.*zp-col; /* This default value is chosen to be as large as possible given the spherical model boundary. */
-      nposn=-1;
-      line_plane_intersect(gp,&ds,posn,&nposn,dx,x,cutoff); /* Returns a new ds equal to the distance to the next Voronoi face, and nposn, the ID of the grid cell that abuts that face. */ 
+    /* Move the working point to the edge of the next Voronoi cell. */
+    for(di=0;di<DIM;di++) x[di]+=ds*dx[di];
+    col+=ds;
+    posn=nposn;
+  } while(col < 2.0*fabs(zp));
+
+  /* Add or subtract cmb. */
+  if(par->polarization){ /* just add it to Stokes I */
+#ifdef FASTEXP
+    ray.intensity[stokesIi]+=FastExp(ray.tau[stokesIi])*md[cmbMolI].local_cmb[cmbLineI];
+#else
+    ray.intensity[stokesIi]+=exp(   -ray.tau[stokesIi])*md[cmbMolI].local_cmb[cmbLineI];
+#endif
+
+  }else{
+#ifdef FASTEXP
+    for(ichan=0;ichan<img[im].nchan;ichan++){
+      ray.intensity[ichan]+=FastExp(ray.tau[ichan])*md[cmbMolI].local_cmb[cmbLineI];
+    }
+#else
+    for(ichan=0;ichan<img[im].nchan;ichan++){
+      ray.intensity[ichan]+=exp(-ray.tau[ichan])*md[cmbMolI].local_cmb[cmbLineI];
+    }
+#endif
+  }
+}
+
+/*....................................................................*/
+void traceray_smooth(rayData ray, const int cmbMolI, const int cmbLineI, const int im\
+  , configInfo *par, struct grid *gp, molData *md, image *img, struct gAuxType *gAux\
+  , const int nlinetot, int *allLineMolIs, int *allLineLineIs, struct cell *dc\
+  , const unsigned long numCells, const double epsilon, gridInterp gips[3]\
+  , const int numSegments, const double oneOnNumSegments, const int nSteps, const double oneOnNSteps){
+  /*
+For a given image pixel position, this function evaluates the intensity of the total light emitted/absorbed along that line of sight through the (possibly rotated) model. The calculation is performed for several frequencies, one per channel of the output image.
+
+Note that the algorithm employed here to solve the RTE is similar to that employed in the function photon() which calculates the average radiant flux impinging on a grid cell: namely the notional photon is started at the side of the model near the observer and 'propagated' in the receding direction until it 'reaches' the far side. This is rather non-physical in conception but it makes the calculation easier.
+
+This version of traceray implements a new algorithm in which the population values are interpolated linearly from those at the vertices of the Delaunay cell which the working point falls within.
+  */
+  const int stokesIi=0;
+  const int numFaces = DIM+1, nVertPerFace=3;
+  int ichan,stokesId,di,status,lenChainPtrs,entryI,exitI,vi,vvi,ci;
+  int si, polMolI, polLineI, iline, molI, lineI;
+  double xp,yp,zp,x[DIM],dir[DIM],projVelRay,vel[DIM];
+  double xCmpntsRay[nVertPerFace], ds, snu_pol[3], dtau, contJnu, contAlpha;
+  double jnu, alpha, lineRedShift, vThisChan, deltav, vfac, remnantSnu, expDTau;
+  double brightnessIncrement;
+  intersectType entryIntcptFirstCell, *cellExitIntcpts=NULL;
+  unsigned long *chainOfCellIds=NULL, dci;
+  unsigned long gis[2][nVertPerFace];
+
+  for(ichan=0;ichan<img[im].nchan;ichan++){
+    ray.tau[ichan]=0.0;
+    ray.intensity[ichan]=0.0;
+  }
+
+  xp=ray.x;
+  yp=ray.y;
+
+  /* The model is circular in projection. We only follow the ray if it will intersect the model.
+  */
+  if((xp*xp+yp*yp)>par->radiusSqu)
+    return;
+
+  zp=-sqrt(par->radiusSqu-(xp*xp+yp*yp)); /* There are two points of intersection between the line of sight and the spherical model surface; this is the Z coordinate (in the unrotated frame) of the one nearer to the observer. */
+
+  /* Rotate the line of sight as desired. */
+  for(di=0;di<DIM;di++){
+    x[di]=xp*img[im].rotMat[di][0] + yp*img[im].rotMat[di][1] + zp*img[im].rotMat[di][2];
+    dir[di]= img[im].rotMat[di][2]; /* This points away from the observer. */
+  }
+
+  /* Find the chain of cells the ray passes through.
+  */
+  status = followRayThroughDelCells(x, dir, gp, dc, numCells, epsilon\
+    , &entryIntcptFirstCell, &chainOfCellIds, &cellExitIntcpts, &lenChainPtrs);//, 0);
+
+  if(status!=0){
+    free(chainOfCellIds);
+    free(cellExitIntcpts);
+    return;
+  }
+
+  entryI = 0;
+  exitI  = 1;
+  dci = chainOfCellIds[0];
+
+  /* Obtain the indices of the grid points on the vertices of the entry face.
+  */
+  vvi = 0;
+  for(vi=0;vi<numFaces;vi++){
+    if(vi!=entryIntcptFirstCell.fi){
+      gis[entryI][vvi++] = dc[dci].vertx[vi]->id;
+    }
+  }
+
+  /* Calculate, for each of the 3 vertices of the entry face, the displacement components in the direction of 'dir'. *** NOTE *** that if all the rays are parallel, we could precalculate these for all the vertices.
+  */
+  for(vi=0;vi<nVertPerFace;vi++)
+    xCmpntsRay[vi] = veloproject(dir, gp[gis[entryI][vi]].x);
+
+  doBaryInterp(entryIntcptFirstCell, gp, gAux, xCmpntsRay, gis[entryI]\
+    , md, par->nSpecies, &gips[entryI]);
+
+  for(ci=0;ci<lenChainPtrs;ci++){
+    /* For each cell we have 2 data structures which give information about respectively the entry and exit points of the ray, including the barycentric coordinates of the intersection point between the ray and the appropriate face of the cell. (If we follow rays in 3D space then the cells will be tetrahedra and the faces triangles.) If we know the value of a quantity Q for each of the vertices, then the linear interpolation of the Q values for any face is (for a 3D space) bary[0]*Q[0] + bary[1]*Q[1] + bary[2]*Q[2], where the indices are taken to run over the vertices of that face. Thus we can calculate the interpolated values Q_entry and Q_exit. Further linear interpolation along the path between entry and exit is straightforward.
+    */
+
+    dci = chainOfCellIds[ci];
+
+    /* Obtain the indices of the grid points on the vertices of the exit face. */
+    vvi = 0;
+    for(vi=0;vi<numFaces;vi++){
+      if(vi!=cellExitIntcpts[ci].fi){
+        gis[exitI][vvi++] = dc[dci].vertx[vi]->id;
+      }
+    }
+
+    /* Calculate, for each of the 3 vertices of the exit face, the displacement components in the direction of 'dir'. *** NOTE *** that if all the rays are parallel, we could precalculate these for all the vertices.
+    */
+    for(vi=0;vi<nVertPerFace;vi++)
+      xCmpntsRay[vi] = veloproject(dir, gp[gis[exitI][vi]].x);
+
+    doBaryInterp(cellExitIntcpts[ci], gp, gAux, xCmpntsRay, gis[exitI]\
+      , md, par->nSpecies, &gips[exitI]);
+
+    /* At this point we have interpolated all the values of interest to both the entry and exit points of the cell. Now we break the path between entry and exit into several segments and calculate all these values at the midpoint of each segment.
+
+At the moment I will fix the number of segments, but it might possibly be faster to rather have a fixed segment length (in barycentric coordinates) and vary the number depending on how many of these lengths fit in the path between entry and exit.
+    */
+    ds = (gips[exitI].xCmpntRay - gips[entryI].xCmpntRay)*oneOnNumSegments;
+
+    for(si=0;si<numSegments;si++){
+      doSegmentInterp(gips, entryI, md, par->nSpecies, oneOnNumSegments, si);
+
       if(par->polarization){
-        sourceFunc_pol(snu_pol,&alpha,gp,posn,0,0,img[im].rotMat);
+        polMolI = 0; /****** Always?? */
+        polLineI = 0; /****** Always?? */
+        sourceFunc_pol(gips[2].B, gips[2].mol[polMolI], polLineI, img[im].rotMat, snu_pol, &alpha);
         dtau=alpha*ds;
         calcSourceFn(dtau, par, &remnantSnu, &expDTau);
-        remnantSnu *= md[0].norminv*ds;
+        remnantSnu *= md[polMolI].norminv*ds;
 
-        for(ichan=0;ichan<img[im].nchan;ichan++){ /* Loop over I, Q and U */
+        for(stokesId=0;stokesId<img[im].nchan;stokesId++){ /* Loop over I, Q and U */
 #ifdef FASTEXP
-          ray.intensity[ichan]+=FastExp(ray.tau[ichan])*remnantSnu*snu_pol[ichan];
+          brightnessIncrement = FastExp(ray.tau[stokesId])*remnantSnu*snu_pol[stokesId];
 #else
-          ray.intensity[ichan]+=   exp(-ray.tau[ichan])*remnantSnu*snu_pol[ichan];
+          brightnessIncrement =    exp(-ray.tau[stokesId])*remnantSnu*snu_pol[stokesId];
 #endif
-          ray.tau[ichan]+=dtau; //**** But this will be the same for I, Q or U.
+          ray.intensity[stokesId] += brightnessIncrement;
+          ray.tau[stokesId]+=dtau; //**** But this will be the same for I, Q or U.
         }
       } else {
+        /* It appears to be necessary to sample the velocity function in the following way rather than interpolating it from the vertices of the Delaunay cell in the same way as with all the other quantities of interest. Velocity varies too much across the cells, and in a nonlinear way, for linear interpolation to yield a totally satisfactory result.
+        */
+        velocity(gips[2].x[0], gips[2].x[1], gips[2].x[2], vel);
+        projVelRay = veloproject(dir, vel);
+
+        /* Calculate first the continuum stuff because it is the same for all channels:
+        */
+        contJnu = 0.0;
+        contAlpha = 0.0;
+        sourceFunc_cont_raytrace(gips[2].mol[cmbMolI], cmbLineI, &contJnu, &contAlpha);
+
         for(ichan=0;ichan<img[im].nchan;ichan++){
-          jnu=.0;
-          alpha=0.;
+          jnu = contJnu;
+          alpha = contAlpha;
+          vThisChan=(ichan-(img[im].nchan-1)*0.5)*img[im].velres; /* Consistent with the WCS definition in writefits(). */
 
           if(img[im].doline){
             for(molI=0;molI<par->nSpecies;molI++){
               for(lineI=0;lineI<md[molI].nline;lineI++){
-                if(md[molI].freq[lineI] > img[im].freq-img[im].bandwidth*0.5
+                if(md[molI].freq[lineI] > img[im].freq-img[im].bandwidth*0.5\
                 && md[molI].freq[lineI] < img[im].freq+img[im].bandwidth*0.5){
-                  /* Calculate the red shift of the transition wrt to the frequency specified for the image. */
+                  /* Calculate the red shift of the transition wrt to the frequency specified for the image.
+                  */
                   if(img[im].trans > -1){
                     lineRedShift=(md[molI].freq[img[im].trans]-md[molI].freq[lineI])/md[molI].freq[img[im].trans]*CLIGHT;
                   } else {
                     lineRedShift=(img[im].freq-md[molI].freq[lineI])/img[im].freq*CLIGHT;
                   }
 
-                  vThisChan=(ichan-(img[im].nchan-1)/2.)*img[im].velres; /* Consistent with the WCS definition in writefits(). */
                   deltav = vThisChan - img[im].source_vel - lineRedShift;
                   /* Line centre occurs when deltav = the recession velocity of the radiating material. Explanation of the signs of the 2nd and 3rd terms on the RHS: (i) A bulk source velocity (which is defined as >0 for the receding direction) should be added to the material velocity field; this is equivalent to subtracting it from deltav, as here. (ii) A positive value of lineRedShift means the line is red-shifted wrt to the frequency specified for the image. The effect is the same as if the line and image frequencies were the same, but the bulk recession velocity were higher. lineRedShift should thus be added to the recession velocity, which is equivalent to subtracting it from deltav, as here. */
 
-                  /* Calculate an approximate average line-shape function at deltav within the Voronoi cell. */
-                  if(!par->pregrid) velocityspline2(x,dx,ds,gp[posn].mol[molI].binv,deltav,&vfac);
-                  else vfac=gaussline(deltav-veloproject(dx,gp[posn].vel),gp[posn].mol[molI].binv);
+                  calcLineAmpInterp(projVelRay, gips[2].mol[molI].binv, deltav, &vfac);
 
-                  /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. */
-                  sourceFunc_line(&jnu,&alpha,md,vfac,gp,posn,molI,lineI);
-                }
-              }
-            }
-          }
+                  /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line.
+                  */
+                  sourceFunc_line_raytrace(md[molI], vfac, gips[2].mol[molI], lineI, &jnu, &alpha);
+                } /* end if within freq range. */
+              } /* end loop over lines this mol. */
+            } /* end loop over all mols. */
+          } /* end if doLine. */
 
-          sourceFunc_cont(&jnu,&alpha,gp,posn,cmbMolI,cmbLineI);
-          dtau=alpha*ds;
+          dtau = alpha*ds;
+//???          if(dtau < -30) dtau = -30; // as in photon()?
           calcSourceFn(dtau, par, &remnantSnu, &expDTau);
           remnantSnu *= jnu*md[0].norminv*ds;
 #ifdef FASTEXP
-          ray.intensity[ichan]+=FastExp(ray.tau[ichan])*remnantSnu;
+          brightnessIncrement = FastExp(ray.tau[ichan])*remnantSnu;
 #else
-          ray.intensity[ichan]+=   exp(-ray.tau[ichan])*remnantSnu;
+          brightnessIncrement =    exp(-ray.tau[ichan])*remnantSnu;
 #endif
-          ray.tau[ichan]+=dtau;
-        }
-      }
-
-      /* Move the working point to the edge of the next Voronoi cell. */
-      for(i=0;i<3;i++) x[i]+=ds*dx[i];
-      col+=ds;
-      posn=nposn;
-    } while(col < 2.0*fabs(zp));
-
-    /* Add or subtract cmb. */
-    if(par->polarization){ /* just add it to Stokes I */
+          ray.intensity[ichan] += brightnessIncrement;
+          ray.tau[ichan] += dtau;
+        } /* End loop over channels. */
+      } /* End if(par->polarization). */
+    } /* End loop over segments within cell. */
+
+    entryI = exitI;
+    exitI = 1 - exitI;
+  } /* End loop over cells in the chain traversed by the ray. */
+
+  /* Add or subtract cmb. */
+  if(par->polarization){ /* just add it to Stokes I */
 #ifdef FASTEXP
-      ray.intensity[stokesIi]+=FastExp(ray.tau[stokesIi])*md[cmbMolI].local_cmb[cmbLineI];
+    ray.intensity[stokesIi]+=FastExp(ray.tau[stokesIi])*md[cmbMolI].local_cmb[cmbLineI];
 #else
-      ray.intensity[stokesIi]+=exp(   -ray.tau[stokesIi])*md[cmbMolI].local_cmb[cmbLineI];
+    ray.intensity[stokesIi]+=exp(   -ray.tau[stokesIi])*md[cmbMolI].local_cmb[cmbLineI];
 #endif
 
-    }else{
+  }else{
 #ifdef FASTEXP
-      for(ichan=0;ichan<img[im].nchan;ichan++){
-        ray.intensity[ichan]+=FastExp(ray.tau[ichan])*md[cmbMolI].local_cmb[cmbLineI];
-      }
+    for(ichan=0;ichan<img[im].nchan;ichan++){
+      ray.intensity[ichan]+=FastExp(ray.tau[ichan])*md[cmbMolI].local_cmb[cmbLineI];
+    }
 #else
-      for(ichan=0;ichan<img[im].nchan;ichan++){
-        ray.intensity[ichan]+=exp(-ray.tau[ichan])*md[cmbMolI].local_cmb[cmbLineI];
-      }
-#endif
+    for(ichan=0;ichan<img[im].nchan;ichan++){
+      ray.intensity[ichan]+=exp(-ray.tau[ichan])*md[cmbMolI].local_cmb[cmbLineI];
     }
+#endif
   }
-}
 
+  free(chainOfCellIds);
+  free(cellExitIntcpts);
+}
 
+/*....................................................................*/
 void
 raytrace(int im, configInfo *par, struct grid *gp, molData *md, image *img){
-  int aa,ichan,px,iline,cmbMolI,cmbLineI,i,threadI,nRaysDone,molI,lineI;
-  double size,minfreq,absDeltaFreq,totalNumPixelsMinus1=(double)(img[im].pxls*img[im].pxls-1);
-  double cutoff;
-  const gsl_rng_type *ranNumGenType = gsl_rng_ranlxs2;
+  /*
+This function constructs an image cube by following sets of rays (at least 1 per image pixel) through the model, solving the radiative transfer equations as appropriate for each ray. The ray locations within each pixel are chosen randomly within the pixel, but the number of rays per pixel is set equal to the number of projected model grid points falling within that pixel, down to a minimum equal to par->alias.
+  */
+  const int maxNumRaysPerPixel=20; /**** Arbitrary - could make this a global, or an argument. Set it to zero to indicate there is no maximum. */
+  const double cutoff = par->minScale*1.0e-7;
+  const int numFaces=1+DIM, numInterpPoints=3, numSegments=5, minNumRaysForAverage=2;
+  const double oneOnNFaces=1.0/(double)numFaces, oneOnNumSegments = 1.0/(double)numSegments;
+  const double epsilon = 1.0e-6; // Needs thinking about. Double precision is much smaller than this.
+  const int nStepsThruCell=10;
+  const double oneOnNSteps=1.0/(double)nStepsThruCell;
+
+  double size,oneOnNumActiveRaysMinus1,imgCentreXPixels,imgCentreYPixels,minfreq,absDeltaFreq,x[2],sum,oneOnNumRays;//,oneOnTotalNumPixelsMinus1
+  unsigned int totalNumImagePixels,ppi,numPixelsForInterp;
+  int gi,molI,lineI,ei,li,nlinetot,iline,tmptrans,ichan,numActiveRays,i,di,xi,yi,ri,c,id,ids[3],vi;
+  int *allLineMolIs,*allLineLineIs,cmbMolI,cmbLineI;
+  rayData *rays;
+  struct cell *dc=NULL;
+  unsigned long numCells,dci;
+  struct gAuxType *gAux=NULL; /* This will hold some precalculated values for the grid points. */
+  coordT *pt_array, point[3];
+  char flags[255];
+  boolT ismalloc = False,isoutside;
+  realT bestdist;
+  facetT *facet;//, *neighbor, **neighborp;
+  vertexT *vertex,**vertexp;
+  int curlong, totlong;
+  double triangle[3][2],barys[3];
+  _Bool isOutsideImage;
   gsl_error_handler_t *defaultErrorHandler=NULL;
 
-  gsl_rng *ran = gsl_rng_alloc(ranNumGenType);	/* Random number generator */
-#ifdef TEST
-  gsl_rng_set(ran,178490);
-#else
-  gsl_rng_set(ran,time(0));
-#endif
-
-  gsl_rng **threadRans;
-  threadRans = malloc(sizeof(gsl_rng *)*par->nThreads);
-
-  for (i=0;i<par->nThreads;i++){
-    threadRans[i] = gsl_rng_alloc(ranNumGenType);
-    gsl_rng_set(threadRans[i],(int)(gsl_rng_uniform(ran)*1e6));
-  }
-
-  size=img[im].distance*img[im].imgres;
+  size = img[im].distance*img[im].imgres;
+  totalNumImagePixels = img[im].pxls*img[im].pxls;
+  imgCentreXPixels = img[im].pxls/2.0;
+  imgCentreYPixels = img[im].pxls/2.0;
 
   if(img[im].doline){
     /* The user may have set img.trans/img.molI but not img.freq. If so, we calculate freq now.
@@ -272,15 +549,114 @@ For a line image however, we may have a problem, because of the pair of quantiti
     cmbLineI = 0;
   }
 
-  cutoff = par->minScale*1.0e-7;
-
-  for(px=0;px<(img[im].pxls*img[im].pxls);px++){
+  for(ppi=0;ppi<totalNumImagePixels;ppi++){
     for(ichan=0;ichan<img[im].nchan;ichan++){
-      img[im].pixel[px].intense[ichan]=0.0;
-      img[im].pixel[px].tau[ichan]=0.0;
+      img[im].pixel[ppi].intense[ichan] = 0.0;
+      img[im].pixel[ppi].tau[    ichan] = 0.0;
     }
   }
 
+  for(ppi=0;ppi<totalNumImagePixels;ppi++)
+    img[im].pixel[ppi].numRays = 0;
+
+  /* The following is the first of the 3 main loops in raytrace. Here we loop over the (internal or non-sink) grid points. We're doing 2 things: loading the rotated, projected coordinates into the rays list, and counting the rays per image pixel.
+  */
+  rays = malloc(sizeof(rayData)*par->pIntensity); /* We may need to reallocate this later. */
+  numActiveRays = 0;
+  for(gi=0;gi<par->pIntensity;gi++){
+    /* Apply the inverse (i.e. transpose) rotation matrix. (We use the inverse matrix here because here we want to rotate grid coordinates to the observer frame, whereas inside traceray() we rotate observer coordinates to the grid frame.)
+    */
+    for(i=0;i<2;i++){
+      x[i]=0.0;
+      for(di=0;di<DIM;di++){
+        x[i] += gp[gi].x[di]*img[im].rotMat[di][i];
+      }
+    }
+
+    /* Calculate which pixel the projected position (x[0],x[1]) falls within.
+    */
+    xi = floor(x[0]/size + imgCentreXPixels);
+    yi = floor(x[1]/size + imgCentreYPixels);
+    if(xi<0 || xi>=img[im].pxls || yi<0 || yi>=img[im].pxls){
+      isOutsideImage = 1;
+      ppi = 0; /* Under these circumstances it ought never to be accessed, but it is not good to leave it without a value at all. */
+    }else{
+      isOutsideImage = 0;
+      ppi = (unsigned int)yi*(unsigned int)img[im].pxls + (unsigned int)xi;
+    }
+
+    /* See if we want to keep the ray. For the time being we will include those outside the image bounds, but a cleverer algorithm would exclude some of them. Note that maxNumRaysPerPixel<1 is used to flag that there is no upper limit to the number of rays per pixel.
+    */
+    if(isOutsideImage || maxNumRaysPerPixel<1 || img[im].pixel[ppi].numRays<maxNumRaysPerPixel){
+      if(!isOutsideImage)
+        img[im].pixel[ppi].numRays++;
+
+      rays[numActiveRays].ppi = ppi;
+      rays[numActiveRays].x = x[0];
+      rays[numActiveRays].y = x[1];
+      rays[numActiveRays].tau       = malloc(sizeof(double)*img[im].nchan);
+      rays[numActiveRays].intensity = malloc(sizeof(double)*img[im].nchan);
+      for(ichan=0;ichan<img[im].nchan;ichan++) {
+        rays[numActiveRays].tau[ichan] = 0.0;
+        rays[numActiveRays].intensity[ichan] = 0.0;
+      }
+
+      numActiveRays++;
+    }
+  } /* End loop 1, over grid points. */
+
+  oneOnNumActiveRaysMinus1 = 1.0/(double)(numActiveRays-1);
+
+  if(numActiveRays<par->pIntensity)
+    rays = realloc(rays, sizeof(rayData)*numActiveRays);
+
+  if(par->traceRayAlgorithm==1){
+    delaunay(DIM, gp, (unsigned long)par->ncell, 1, &dc, &numCells); /* mallocs dc if getCells==T */
+
+    /* We need to process the list of cells a bit further - calculate their centres, and reset the id values to be the same as the index of the cell in the list. (This last because we are going to construct other lists to indicate which cells have been visited etc.)
+    */
+    for(dci=0;dci<numCells;dci++){
+      for(di=0;di<DIM;di++){
+        sum = 0.0;
+        for(vi=0;vi<numFaces;vi++){
+          sum += dc[dci].vertx[vi]->x[di];
+        }
+        dc[dci].centre[di] = sum*oneOnNFaces;
+      }
+
+      dc[dci].id = dci;
+    }
+
+  }else if(par->traceRayAlgorithm!=0){
+    if(!silent) bail_out("Unrecognized value of par.traceRayAlgorithm");
+    exit(1);
+  }
+
+  /* Precalculate binv*nmol*pops for all grid points.
+  */
+  gAux = malloc(sizeof(*gAux)*par->ncell);
+  for(gi=0;gi<par->ncell;gi++){
+    gAux[gi].mol = malloc(sizeof(*(gAux[gi].mol))*par->nSpecies);
+    for(molI=0;molI<par->nSpecies;molI++){
+      gAux[gi].mol[molI].specNumDens = malloc(sizeof(*(gAux[gi].mol[molI].specNumDens))*md[molI].nlev);
+      gAux[gi].mol[molI].dust        = malloc(sizeof(*(gAux[gi].mol[molI].dust))       *md[molI].nline);
+      gAux[gi].mol[molI].knu         = malloc(sizeof(*(gAux[gi].mol[molI].knu))        *md[molI].nline);
+
+      for(ei=0;ei<md[molI].nlev;ei++)
+        gAux[gi].mol[molI].specNumDens[ei] = gp[gi].mol[molI].binv\
+          *gp[gi].mol[molI].nmol*gp[gi].mol[molI].pops[ei];
+
+      /* This next is repetition. I do it in order to be able to use the same sourcefunc.c functions for the interpolated grid values as for the 'standard' algorithm. With a sensible arrangement of memory for the grid values, this would be unnecessary.
+      */
+      for(li=0;li<md[molI].nline;li++){
+        gAux[gi].mol[molI].dust[li] = gp[gi].mol[molI].dust[li];
+        gAux[gi].mol[molI].knu[ li] = gp[gi].mol[molI].knu[ li];
+      }
+    }
+  }
+
+  /* This is the start of loop 2/3, which loops over the rays. We trace each ray, then load into the image cube those for which the number of rays per pixel exceeds a minimum. The remaining image pixels we handle via an interpolation algorithm in loop 3.
+  */
   defaultErrorHandler = gsl_set_error_handler_off();
   /*
 The GSL documentation does not recommend leaving the error handler at the default within multi-threaded code.
@@ -288,52 +664,165 @@ The GSL documentation does not recommend leaving the error handler at the defaul
 While this is off however, gsl_* calls will not exit if they encounter a problem. We may need to pay some attention to trapping their errors.
   */
 
-  nRaysDone=0;
   omp_set_dynamic(0);
-  #pragma omp parallel private(px,aa,threadI) num_threads(par->nThreads)
+  #pragma omp parallel num_threads(par->nThreads)
   {
-    threadI = omp_get_thread_num();
+    /* Declaration of thread-private pointers.
+    */
+    int threadI = omp_get_thread_num();
+    int ii, si, ri;
+    gridInterp gips[numInterpPoints];
+
+    if(par->traceRayAlgorithm==1){
+      /* Allocate memory for the interpolation points:
+      */
+      for(ii=0;ii<numInterpPoints;ii++){
+        gips[ii].mol = malloc(sizeof(*(gips[ii].mol))*par->nSpecies);
+        for(si=0;si<par->nSpecies;si++){
+          gips[ii].mol[si].specNumDens = malloc(sizeof(*(gips[ii].mol[si].specNumDens))*md[si].nlev);
+          gips[ii].mol[si].dust        = malloc(sizeof(*(gips[ii].mol[si].dust))       *md[si].nline);
+          gips[ii].mol[si].knu         = malloc(sizeof(*(gips[ii].mol[si].knu))        *md[si].nline);
+        }
+      }
+    }
 
-    /* Declaration of thread-private pointers. */
-    rayData ray;
-    ray.intensity=malloc(sizeof(double) * img[im].nchan);
-    ray.tau=malloc(sizeof(double) * img[im].nchan);
+    #pragma omp for schedule(dynamic)
+    for(ri=0;ri<numActiveRays;ri++){
+      if(par->traceRayAlgorithm==0)
+        traceray(rays[ri], cmbMolI, cmbLineI, im, par, gp, md, img, gAux\
+          , nlinetot, allLineMolIs, allLineLineIs, cutoff, nStepsThruCell, oneOnNSteps);
+      else if(par->traceRayAlgorithm==1)
+        traceray_smooth(rays[ri], cmbMolI, cmbLineI, im, par, gp, md, img, gAux\
+          , nlinetot, allLineMolIs, allLineLineIs, dc, numCells, epsilon, gips\
+          , numSegments, oneOnNumSegments, nStepsThruCell, oneOnNSteps);
 
-    #pragma omp for
-    /* Main loop through pixel grid. */
-    for(px=0;px<(img[im].pxls*img[im].pxls);px++){
-      #pragma omp atomic
-      ++nRaysDone;
+      if (threadI == 0){ /* i.e., is master thread */
+        if(!silent) progressbar((double)(ri)*oneOnNumActiveRaysMinus1, 13);
+      }
+    }
 
-      for(aa=0;aa<par->antialias;aa++){
-        ray.x = -size*(gsl_rng_uniform(threadRans[threadI]) + px%img[im].pxls - 0.5*img[im].pxls);
-        ray.y =  size*(gsl_rng_uniform(threadRans[threadI]) + px/img[im].pxls - 0.5*img[im].pxls);
+    if(par->traceRayAlgorithm==1){
+      for(ii=0;ii<numInterpPoints;ii++)
+        freePop2(par->nSpecies, gips[ii].mol);
+    }
+  } /* End of parallel block. */
 
-        traceray(ray, cmbMolI, cmbLineI, im, par, gp, md, img, cutoff);
+  gsl_set_error_handler(defaultErrorHandler);
 
-        #pragma omp critical
-        {
-          for(ichan=0;ichan<img[im].nchan;ichan++){
-            img[im].pixel[px].intense[ichan] += ray.intensity[ichan]/(double) par->antialias;
-            img[im].pixel[px].tau[ichan] += ray.tau[ichan]/(double) par->antialias;
-          }
-        }
+  /* For pixels with more than a cutoff number of rays, just average those rays into the pixel:
+  */
+  for(ri=0;ri<numActiveRays;ri++){
+    if(img[im].pixel[rays[ri].ppi].numRays >= minNumRaysForAverage){
+      for(ichan=0;ichan<img[im].nchan;ichan++){
+        img[im].pixel[rays[ri].ppi].intense[ichan] += rays[ri].intensity[ichan];
+        img[im].pixel[rays[ri].ppi].tau[    ichan] += rays[ri].tau[      ichan];
       }
-      if (threadI == 0){ /* i.e., is master thread */
-        if(!silent) progressbar((double)(nRaysDone)/totalNumPixelsMinus1, 13);
+    }
+  }
+  numPixelsForInterp = 0;
+  for(ppi=0;ppi<totalNumImagePixels;ppi++){
+    if(img[im].pixel[ppi].numRays >= minNumRaysForAverage){
+      oneOnNumRays = 1.0/(double)img[im].pixel[ppi].numRays;
+      for(ichan=0;ichan<img[im].nchan;ichan++){
+        img[im].pixel[ppi].intense[ichan] *= oneOnNumRays;
+        img[im].pixel[ppi].tau[    ichan] *= oneOnNumRays;
       }
+    }else
+      numPixelsForInterp++;
+  }
+
+  if(numPixelsForInterp>0){
+    /* Now we enter main loop 3/3, in which we loop over image pixels, and for any we need to interpolate, we do so. But first we need to invoke qhull to get a Delaunay triangulation of the projected points.
+    */
+    pt_array=malloc(sizeof(coordT)*2*numActiveRays);
+
+    for(ri=0;ri<numActiveRays;ri++) {
+      pt_array[ri*2+0] = rays[ri].x;
+      pt_array[ri*2+1] = rays[ri].y;
     }
 
-    free(ray.tau);
-    free(ray.intensity);
-  } /* End of parallel block. */
+    sprintf(flags,"qhull d Qbb");
+    if(qh_new_qhull(2, numActiveRays, pt_array, ismalloc, flags, NULL, NULL)) {
+      if(!silent) bail_out("Qhull failed to triangulate");
+      exit(1);
+    }
 
-  gsl_set_error_handler(defaultErrorHandler);
+    point[2]=0.;
+    for(ppi=0;ppi<totalNumImagePixels;ppi++){
+      if(img[im].pixel[ppi].numRays < minNumRaysForAverage){
+        xi = (int)(ppi%(unsigned int)img[im].pxls);
+        yi = floor(ppi/(double)img[im].pxls);
+        point[0] = size*(0.5 + xi - imgCentreXPixels);
+        point[1] = size*(0.5 + yi - imgCentreYPixels);
+
+        qh_setdelaunay (3, 1, point);
+        facet = qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside);
+        if(isoutside){
+          c=0;
+          FOREACHvertex_( facet->vertices ) {
+            id = qh_pointid(vertex->point);
+            triangle[c][0] = rays[id].x;
+            triangle[c][1] = rays[id].y;
+            ids[c]=id;
+            c++;
+          }
 
-  for (i=0;i<par->nThreads;i++){
-    gsl_rng_free(threadRans[i]);
+          calcTriangleBaryCoords(triangle, (double)point[0], (double)point[1], barys);
+
+          /* Interpolate: */
+          for(ichan=0;ichan<img[im].nchan;ichan++){
+            img[im].pixel[ppi].intense[ichan] += barys[0]*rays[ids[0]].intensity[ichan]\
+                                               + barys[1]*rays[ids[1]].intensity[ichan]\
+                                               + barys[2]*rays[ids[2]].intensity[ichan];
+            img[im].pixel[ppi].tau[    ichan] += barys[0]*rays[ids[0]].tau[ichan]\
+                                               + barys[1]*rays[ids[1]].tau[ichan]\
+                                               + barys[2]*rays[ids[2]].tau[ichan];
+          }
+        } /* end if !isoutside */
+      } /* end if(img[im].pixel[ppi].numRays < minNumRaysForAverage) */
+    } /* end loop over image pixels */
+
+    qh_freeqhull(!qh_ALL);
+    qh_memfreeshort (&curlong, &totlong);
+    free(pt_array);
+  } /* end if(numPixelsForInterp>0) */
+
+  img[im].trans=tmptrans;
+
+  freeGAux((unsigned long)par->ncell, par->nSpecies, gAux);
+  if(par->traceRayAlgorithm==1)
+    free(dc);
+  for(ri=0;ri<numActiveRays;ri++){
+    free(rays[ri].tau);
+    free(rays[ri].intensity);
   }
-  free(threadRans);
-  gsl_rng_free(ran);
+  free(rays);
+  free(allLineMolIs);
+  free(allLineLineIs);
+}
+
+/*....................................................................*/
+void calcTriangleBaryCoords(double vertices[3][2], double x, double y, double barys[3]){
+  double mat[2][2], vec[2], det;
+  /*
+The barycentric coordinates (L0,L1,L2) of a point r[] in a triangle v[] are given by
+
+	(v[0][0]-v[2][0]  v[1][0]-v[2][0]) (L0)   (r[0]-v[2][0])
+	(                                )*(  ) = (            ),
+	(v[0][1]-v[2][1]  v[1][1]-v[2][1]) (L1)   (r[1]-v[2][1])
+
+with L2 = 1 - L0 - L1.
+  */
+  mat[0][0] = vertices[0][0] - vertices[2][0];
+  mat[0][1] = vertices[1][0] - vertices[2][0];
+  mat[1][0] = vertices[0][1] - vertices[2][1];
+  mat[1][1] = vertices[1][1] - vertices[2][1];
+  vec[0] = x - vertices[2][0];
+  vec[1] = y - vertices[2][1];
+  det = mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0];
+  /*** We're assuming that the triangle is not pathological, i.e that det!=0. */
+  barys[0] = ( mat[1][1]*vec[0] - mat[0][1]*vec[1])/det;
+  barys[1] = (-mat[1][0]*vec[0] + mat[0][0]*vec[1])/det;
+  barys[2] = 1.0 - barys[0] - barys[1];
 }
 
diff --git a/src/smooth.c b/src/smooth.c
index 1c31d70..fa90db5 100644
--- a/src/smooth.c
+++ b/src/smooth.c
@@ -19,7 +19,9 @@ smooth(configInfo *par, struct grid *g){
   int cn;
   double move[3];	/* Auxillary array for smoothing the grid	*/
   double dist;		/* Distance to a neighbor					*/
-  	
+  struct cell *dc=NULL; /* Not used at present. */
+  unsigned long numCells;
+
   for(sg=0;sg<N_SMOOTH_ITERS;sg++){
     for(i=0;i<par->ncell && !g[i].sink;i++){
       mindist=1e30;
@@ -67,9 +69,10 @@ smooth(configInfo *par, struct grid *g){
       }	
     }
 		
-    qhull(par, g);	
+    delaunay(DIM, g, (unsigned long)par->ncell, 0, &dc, &numCells);
     distCalc(par, g);	    
     if(!silent) progressbar((double)(sg+1)/(double)N_SMOOTH_ITERS, 5);	
+    free(dc);
   }	
 }
 
diff --git a/src/sourcefunc.c b/src/sourcefunc.c
index 910e5ae..f31ae82 100644
--- a/src/sourcefunc.c
+++ b/src/sourcefunc.c
@@ -5,6 +5,8 @@
  *  Copyright (C) 2006-2014 Christian Brinch
  *  Copyright (C) 2015 The LIME development team
  *
+TODO:
+  - Merge sourceFunc_*_raytrace() and sourceFunc_*() after changes to the way grid data is stored makes this possible.
  */
 
 #include "lime.h"
@@ -18,7 +20,7 @@ sourceFunc(double *snu, double *dtau, double ds, molData *m,double vfac,struct g
   jnu	 = g[pos].mol[ispec].dust[iline]*g[pos].mol[ispec].knu[iline];
   
   /* Line part:		j_nu = v*consts*1/b*rho*n_i*A_ij */
-  if(doline) jnu	+= vfac*HPIP*g[pos].mol[ispec].binv*g[pos].nmol[ispec]*g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].aeinst[iline];
+  if(doline) jnu	+= vfac*HPIP*g[pos].mol[ispec].binv*g[pos].mol[ispec].nmol*g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].aeinst[iline];
   
   
   
@@ -28,8 +30,10 @@ sourceFunc(double *snu, double *dtau, double ds, molData *m,double vfac,struct g
   
   
   /* Line part: alpha_nu = v*const*1/b*rho*(n_j*B_ij-n_i*B_ji) */
-  if(doline) alpha += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].nmol[ispec]*(g[pos].mol[ispec].pops[m[ispec].lal[iline]]*m[ispec].beinstl[iline]
-                                                                          -g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].beinstu[iline]);
+  if(doline)\
+    alpha += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].mol[ispec].nmol\
+    *(g[pos].mol[ispec].pops[m[ispec].lal[iline]]*m[ispec].beinstl[iline]\
+    - g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].beinstu[iline]);
   
   
   /* Calculate source function and tau */
@@ -43,34 +47,72 @@ sourceFunc(double *snu, double *dtau, double ds, molData *m,double vfac,struct g
 }
 
 
+/*....................................................................*/
 void
-sourceFunc_line(double *jnu, double *alpha, molData *m,double vfac,struct grid *g,int pos,int ispec, int iline){
-  
+sourceFunc_line(const molData md, const double vfac, const struct populations gm\
+  , const int lineI, double *jnu, double *alpha){
+
   /* Line part:		j_nu = v*consts*1/b*rho*n_i*A_ij */
-  *jnu   += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].nmol[ispec]*g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].aeinst[iline];
-  
+  *jnu   += vfac*HPIP*gm.binv*gm.nmol*gm.pops[md.lau[lineI]]*md.aeinst[lineI];
+
   /* Line part: alpha_nu = v*const*1/b*rho*(n_j*B_ij-n_i*B_ji) */
-  *alpha += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].nmol[ispec]*(g[pos].mol[ispec].pops[m[ispec].lal[iline]]*m[ispec].beinstl[iline]
-                                                                -g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].beinstu[iline]);
-  
+  *alpha += vfac*HPIP*gm.binv*gm.nmol*(gm.pops[md.lal[lineI]]*md.beinstl[lineI]
+                                      -gm.pops[md.lau[lineI]]*md.beinstu[lineI]);
+
   return;
 }
 
+/*....................................................................*/
 void
-sourceFunc_cont(double *jnu, double *alpha,struct grid *g,int pos,int ispec, int iline){
-  
+sourceFunc_line_raytrace(const molData md, const double vfac\
+  , const struct pop2 gm, const int lineI, double *jnu, double *alpha){
+
+  /* Line part:		j_nu = v*consts*1/b*rho*n_i*A_ij */
+  *jnu   += vfac*HPIP*gm.specNumDens[md.lau[lineI]]*md.aeinst[lineI];
+
+  /* Line part: alpha_nu = v*const*1/b*rho*(n_j*B_ij-n_i*B_ji) */
+  *alpha += vfac*HPIP*(gm.specNumDens[md.lal[lineI]]*md.beinstl[lineI]
+                      -gm.specNumDens[md.lau[lineI]]*md.beinstu[lineI]);
+
+  return;
+}
+
+/*....................................................................*/
+void
+sourceFunc_cont(const struct populations gm, const int lineI, double *jnu\
+  , double *alpha){
+
   /* Emission */
   /* Continuum part:	j_nu = T_dust * kappa_nu */
-  *jnu   += g[pos].mol[ispec].dust[iline]*g[pos].mol[ispec].knu[iline];
-  
+  *jnu   += gm.dust[lineI]*gm.knu[lineI];
+
   /* Absorption */
   /* Continuum part: Dust opacity */
-  *alpha += g[pos].mol[ispec].knu[iline];
-  
+  *alpha += gm.knu[lineI];
+
   return;
 }
 
-void sourceFunc_pol(double *snu, double *alpha, struct grid *g, int pos, int ispec, int iline, double (*rotMat)[3]){
+/*....................................................................*/
+void
+sourceFunc_cont_raytrace(const struct pop2 gm, const int lineI, double *jnu\
+  , double *alpha){
+
+  /* Emission */
+  /* Continuum part:	j_nu = T_dust * kappa_nu */
+  *jnu   += gm.dust[lineI]*gm.knu[lineI];
+
+  /* Absorption */
+  /* Continuum part: Dust opacity */
+  *alpha += gm.knu[lineI];
+
+  return;
+}
+
+/*....................................................................*/
+void
+sourceFunc_pol(double B[3], const struct pop2 gm, int iline\
+  , double (*rotMat)[3], double *snu, double *alpha){
   /*
 The theory behind this was originally drawn from
 
@@ -78,20 +120,20 @@ The theory behind this was originally drawn from
 
 and references therein. However, as pointed out in Ade, P. A. R. et al, A&A 576, A105 (2015), Padovani's expression for sigma2 is too small by a factor of 2. This correction has been propagated here.
   */
+
   double jnu, trigFuncs[3];
-  
-  stokesangles(g[pos].x[0], g[pos].x[1], g[pos].x[2], rotMat, trigFuncs);
-  
+
+  stokesangles(B, rotMat, trigFuncs);
+
   /* Emission */
   /* Continuum part:	j_nu = rho_dust * kappa_nu */
-  jnu = g[pos].mol[ispec].dust[iline]*g[pos].mol[ispec].knu[iline];
+  jnu = gm.dust[iline]*gm.knu[iline];
   snu[0] = jnu*(1.0 - maxp*(trigFuncs[0] - 2.0/3.0));
   snu[1] = jnu*maxp*trigFuncs[1]*trigFuncs[0];
   snu[2] = jnu*maxp*trigFuncs[2]*trigFuncs[0];
   
   /* Absorption */
   /* Continuum part: Dust opacity */
-  *alpha = g[pos].mol[ispec].knu[iline];
+  *alpha = gm.knu[iline];
 }
 
-
diff --git a/src/stokesangles.c b/src/stokesangles.c
index c66f046..fe77d10 100644
--- a/src/stokesangles.c
+++ b/src/stokesangles.c
@@ -3,13 +3,13 @@
  *  This file is part of LIME, the versatile line modeling engine
  *
  *  Copyright (C) 2006-2014 Christian Brinch
- *  Copyright (C) 2015 The LIME development team
+ *  Copyright (C) 2016 The LIME development team
  *
  */
 
 #include "lime.h"
 
-void stokesangles(double x, double y, double z, double (*rotMat)[3], double *trigFuncs){
+void stokesangles(double B[3], double (*rotMat)[3], double *trigFuncs){
   /*
 This function rotates the B-field vector from the model frame to the observer frame, then calculates and returns some useful values which will in function sourceFunc_pol() make it easy to obtain the Stokes parameters of polarized submillimetre dust emission. (For an explanation of the reasons for choosing the particular quantities we do, see the comment in that function.)
 
@@ -42,12 +42,10 @@ A vector defined in the LIME model basis can be converted to the observer basis
 where ^T denotes transpose.
   */
   const int nDim=3;
-  double B[nDim],Bp[nDim];
+  double Bp[nDim];
   int i, j;
   double BxySquared, BSquared;
-
-  magfield(x,y,z,B);
-
+	
   /* Rotate B to the observer frame.
   */
   for(i=0;i<nDim;i++){
diff --git a/src/velospline.c b/src/velospline.c
index 0aea5ec..0fd9ea8 100644
--- a/src/velospline.c
+++ b/src/velospline.c
@@ -53,12 +53,14 @@ getVelosplines(configInfo *par, struct grid *g){
   }
   
   for(i=par->pIntensity;i<par->ncell;i++){
+    /* Set velocity values also for sink points (otherwise Delaunay ray-tracing has problems) */
+    velocity(g[i].x[0],g[i].x[1],g[i].x[2],g[i].vel);
+
     g[i].a0=malloc(g[i].numNeigh*sizeof(double));
     g[i].a1=malloc(g[i].numNeigh*sizeof(double));
     g[i].a2=malloc(g[i].numNeigh*sizeof(double));
     g[i].a3=malloc(g[i].numNeigh*sizeof(double));
     g[i].a4=malloc(g[i].numNeigh*sizeof(double));
-    for(j=0;j<3;j++) g[i].vel[j]=0.;
     for(j=0;j<g[i].numNeigh;j++){
       g[i].a0[j]=0.;
       g[i].a1[j]=0.;
diff --git a/src/writefits.c b/src/writefits.c
index 0c77fe8..2701d55 100644
--- a/src/writefits.c
+++ b/src/writefits.c
@@ -3,17 +3,17 @@
  *  This file is part of LIME, the versatile line modeling engine
  *
  *  Copyright (C) 2006-2014 Christian Brinch
- *  Copyright (C) 2015 The LIME development team
+ *  Copyright (C) 2016 The LIME development team
  *
  */
 
 #include "lime.h"
 
 void 
-writefits(int im, configInfo *par, molData *m, image *img){
+write3Dfits(int im, configInfo *par, molData *m, image *img){
   double bscale,bzero,epoch,lonpole,equinox,restfreq;
   double cdelt1,crpix1,crval1,cdelt2,crpix2,crval2;
-  double cdelt3,crpix3,crval3,ru3;
+  double cdelt3,crpix3,crval3,ru3,scale;
   int velref;
   float *row;
   int px,py,ichan;
@@ -23,6 +23,7 @@ writefits(int im, configInfo *par, molData *m, image *img){
   long naxes[3];
   long int fpixels[3],lpixels[3];
   char negfile[100]="! ";
+  unsigned long ppi;
 
   row = malloc(sizeof(*row)*img[im].pxls);
 
@@ -30,7 +31,7 @@ writefits(int im, configInfo *par, molData *m, image *img){
   naxes[1]=img[im].pxls;
   if(img[im].doline==1) naxes[2]=img[im].nchan;
   else if(img[im].doline==0 && par->polarization) naxes[2]=3;
-  else naxes[2]=1;
+  else naxes[2]=1;//********** should call write2Dfits for this.
 
   fits_create_file(&fptr, img[im].filename, &status);
 
@@ -88,20 +89,33 @@ writefits(int im, configInfo *par, molData *m, image *img){
   if(img[im].unit==1) fits_write_key(fptr, TSTRING, "BUNIT", &"JY/PIXEL", "", &status);
   if(img[im].unit==2) fits_write_key(fptr, TSTRING, "BUNIT", &"WM2HZSR ", "", &status);
   if(img[im].unit==3) fits_write_key(fptr, TSTRING, "BUNIT", &"Lsun/PX ", "", &status);
-  if(img[im].unit==3) fits_write_key(fptr, TSTRING, "BUNIT", &"        ", "", &status);
+  if(img[im].unit==4) fits_write_key(fptr, TSTRING, "BUNIT", &"        ", "", &status);
+
+  if(     img[im].unit==0)
+    scale=(CLIGHT/img[im].freq)*(CLIGHT/img[im].freq)/2./KBOLTZ*m[0].norm; 
+  else if(img[im].unit==1)
+    scale=1e26*img[im].imgres*img[im].imgres*m[0].norm;
+  else if(img[im].unit==2)
+    scale=m[0].norm;
+  else if(img[im].unit==3) {
+    ru3 = img[im].distance/1.975e13;
+    scale=4.*PI*ru3*ru3*img[im].freq*img[im].imgres*img[im].imgres*m[0].norm;
+  }
+  else if(img[im].unit!=4) {
+    if(!silent) bail_out("Image unit number invalid");
+    exit(0);
+  }
 
   /* Write FITS data */
-  for(py=0;py<img[im].pxls;py++){
-    for(ichan=0;ichan<img[im].nchan;ichan++){
+  for(ichan=0;ichan<img[im].nchan;ichan++){
+    for(py=0;py<img[im].pxls;py++){
       for(px=0;px<img[im].pxls;px++){
-        if(img[im].unit==0) row[px]=(float) img[im].pixel[px+py*img[im].pxls].intense[ichan]*(CLIGHT/img[im].freq)*(CLIGHT/img[im].freq)/2./KBOLTZ*m[0].norm; 
-        else if(img[im].unit==1) row[px]=(float) img[im].pixel[px+py*img[im].pxls].intense[ichan]*1e26*img[im].imgres*img[im].imgres*m[0].norm;
-        else if(img[im].unit==2) row[px]=(float) img[im].pixel[px+py*img[im].pxls].intense[ichan]*m[0].norm;
-        else if(img[im].unit==3) {
-          ru3 = img[im].distance/1.975e13;
-          row[px]=(float) img[im].pixel[px+py*img[im].pxls].intense[ichan]*4.*PI*ru3*ru3*img[im].freq*img[im].imgres*img[im].imgres*m[0].norm;
-        }
-        else if(img[im].unit==4) row[px]=(float) img[im].pixel[px+py*img[im].pxls].tau[ichan];
+        ppi = py*img[im].pxls + px;
+
+        if(img[im].unit>-1 && img[im].unit<4)
+          row[px]=(float) img[im].pixel[ppi].intense[ichan]*scale;
+        else if(img[im].unit==4)
+          row[px]=(float) img[im].pixel[ppi].tau[ichan];
         else {
           if(!silent) bail_out("Image unit number invalid");
           exit(0);
@@ -123,3 +137,131 @@ writefits(int im, configInfo *par, molData *m, image *img){
 
   if(!silent) done(13);
 }
+
+void 
+write2Dfits(int im, configInfo *par, molData *m, image *img){
+  double bscale,bzero,epoch,lonpole,equinox,restfreq;
+  double cdelt1,crpix1,crval1,cdelt2,crpix2,crval2;
+  double ru3,scale;
+  int velref;
+  float *row,minVal;
+  int px,py;
+  fitsfile *fptr;
+  int status = 0;
+  int naxis=2, bitpix=-32;
+  long naxes[2];
+  long int fpixels[2],lpixels[2];
+  char negfile[100]="! ";
+  unsigned long ppi;
+
+  row = malloc(sizeof(*row)*img[im].pxls);
+
+  naxes[0]=img[im].pxls;
+  naxes[1]=img[im].pxls;
+  if(img[im].unit!=5 && (img[im].doline==1 || (img[im].doline==0 && par->polarization))){
+    if(!silent) bail_out("You need to write a 3D FITS output in this case");
+    exit(0);
+  }
+
+  fits_create_file(&fptr, img[im].filename, &status);
+
+  if(status!=0){
+    if(!silent) warning("Overwriting existing fits file                   ");
+    status=0;
+    strcat(negfile,img[im].filename);
+    fits_create_file(&fptr, negfile, &status);
+  }
+
+  /* Write FITS header */ 
+  fits_create_img(fptr, bitpix, naxis, naxes, &status);
+  epoch   =2.0e3;
+  lonpole =1.8e2;
+  equinox =2.0e3;
+  restfreq=img[im].freq;
+  velref  =257;
+  cdelt1  =-1.8e2*img[im].imgres/PI;
+  crpix1  =(double) img[im].pxls/2+0.5;
+  crval1  =0.0e0;
+  cdelt2  =1.8e2*img[im].imgres/PI;
+  crpix2  =(double) img[im].pxls/2+0.5;
+  crval2  =0.0e0;
+  bscale  =1.0e0;
+  bzero   =0.0e0;
+
+  fits_write_key(fptr, TSTRING, "OBJECT  ", &"LIMEMDL ",    "", &status);
+  fits_write_key(fptr, TDOUBLE, "EPOCH   ", &epoch,         "", &status);
+  fits_write_key(fptr, TDOUBLE, "LONPOLE ", &lonpole,       "", &status);
+  fits_write_key(fptr, TDOUBLE, "EQUINOX ", &equinox,       "", &status);
+  fits_write_key(fptr, TSTRING, "SPECSYS ", &"LSRK    ",    "", &status);
+  fits_write_key(fptr, TDOUBLE, "RESTFREQ", &restfreq,      "", &status);
+  fits_write_key(fptr, TINT,    "VELREF  ", &velref,        "", &status);
+  fits_write_key(fptr, TSTRING, "CTYPE1  ", &"RA---SIN",    "", &status);
+  fits_write_key(fptr, TDOUBLE, "CDELT1  ", &cdelt1,        "", &status);
+  fits_write_key(fptr, TDOUBLE, "CRPIX1  ", &crpix1,        "", &status);
+  fits_write_key(fptr, TDOUBLE, "CRVAL1  ", &crval1,        "", &status);
+  fits_write_key(fptr, TSTRING, "CUNIT1  ", &"DEG     ",    "", &status);	
+  fits_write_key(fptr, TSTRING, "CTYPE2  ", &"DEC--SIN",    "", &status);
+  fits_write_key(fptr, TDOUBLE, "CDELT2  ", &cdelt2,        "", &status);
+  fits_write_key(fptr, TDOUBLE, "CRPIX2  ", &crpix2,        "", &status);  	
+  fits_write_key(fptr, TDOUBLE, "CRVAL2  ", &crval2,        "", &status);
+  fits_write_key(fptr, TSTRING, "CUNIT2  ", &"DEG     ",    "", &status);
+  fits_write_key(fptr, TDOUBLE, "BSCALE  ", &bscale,        "", &status);
+  fits_write_key(fptr, TDOUBLE, "BZERO   ", &bzero,         "", &status);
+  if(img[im].unit==0) fits_write_key(fptr, TSTRING, "BUNIT", &"K       ", "", &status);
+  if(img[im].unit==1) fits_write_key(fptr, TSTRING, "BUNIT", &"JY/PIXEL", "", &status);
+  if(img[im].unit==2) fits_write_key(fptr, TSTRING, "BUNIT", &"WM2HZSR ", "", &status);
+  if(img[im].unit==3) fits_write_key(fptr, TSTRING, "BUNIT", &"Lsun/PX ", "", &status);
+  if(img[im].unit==4) fits_write_key(fptr, TSTRING, "BUNIT", &"        ", "", &status);
+  if(img[im].unit==5) fits_write_key(fptr, TSTRING, "BUNIT", &"N_RAYS  ", "", &status);
+
+  if(     img[im].unit==0)
+    scale=(CLIGHT/img[im].freq)*(CLIGHT/img[im].freq)/2./KBOLTZ*m[0].norm; 
+  else if(img[im].unit==1)
+    scale=1e26*img[im].imgres*img[im].imgres*m[0].norm;
+  else if(img[im].unit==2)
+    scale=m[0].norm;
+  else if(img[im].unit==3) {
+    ru3 = img[im].distance/1.975e13;
+    scale=4.*PI*ru3*ru3*img[im].freq*img[im].imgres*img[im].imgres*m[0].norm;
+  }
+  else if(img[im].unit!=4 && img[im].unit!=5) {
+    if(!silent) bail_out("Image unit number invalid");
+    exit(0);
+  }
+
+  if(img[im].unit<5)
+    minVal = eps;
+  else
+    minVal = 0.0;
+
+  /* Write FITS data */
+  for(py=0;py<img[im].pxls;py++){
+    for(px=0;px<img[im].pxls;px++){
+      ppi = py*img[im].pxls + px;
+
+      if(img[im].unit>-1 && img[im].unit<4)
+        row[px]=(float) img[im].pixel[ppi].intense[0]*scale;
+      else if(img[im].unit==4)
+        row[px]=(float) img[im].pixel[ppi].tau[0];
+      else if(img[im].unit==5)
+        row[px]=(float) img[im].pixel[ppi].numRays;
+      else {
+        if(!silent) bail_out("Image unit number invalid");
+        exit(0);
+      }
+      if (fabs(row[px])<minVal) row[px]=minVal;
+
+      fpixels[0]=1;
+      fpixels[1]=py+1;
+      lpixels[0]=img[im].pxls;
+      lpixels[1]=py+1;
+      fits_write_subset(fptr, TFLOAT, fpixels, lpixels, row, &status);
+    }
+  }
+  fits_close_file(fptr, &status);
+
+  free(row);
+
+  if(!silent) done(13);
+}
+