Combined raytrace

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

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

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);
     }
   }

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);
 
   /*

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++){