Fixes #212

src/raythrucells.c

@@ -16,6 +16,21 @@ error(int errCode, char *message){
   exit(1);
 }
 
+/*....................................................................*/
+intersectType initializeIntersect(const int numDims){
+  int di;
+  intersectType intcpt;
+
+  intcpt.fi = -1;
+  intcpt.orientation = 0;
+  for(di=0;di<numDims;di++)
+    intcpt.bary[di] = 0.0;
+  intcpt.dist = 0.0;
+  intcpt.collPar = 0.0;
+
+  return intcpt;
+}
+
 /*....................................................................*/
 faceType
 extractFace(const int numDims, double *vertexCoords, struct simplex *dc\
@@ -569,11 +584,11 @@ followRayThroughCells(const int numDims, double *x, double *dir, double *vertexC
   /*
 The present function follows a ray through a connected, convex set of cells (assumed to be simplices) 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.
 
-The calling routine should free chainOfCellIds, cellExitIntcpts & cellVisited after use.
+The pointers *chainOfCellIds and *cellExitIntcpts should be freed after the function is called.
 
 The argument facePtrs may be set to NULL, in which case the function will construct each face from the list of cells etc as it needs it. This saves on memory but takes more time. If the calling routine supplies these values it needs to do something like as follows:
 
-	faceType face,*facePtrs[N_DIMS+1]=malloc(sizeof(*(*facePtrs[N_DIMS+1]))*numFaces); // numFaces must of course be calculated beforehand.
+	faceType face,*facePtrs[N_DIMS+1]=malloc(sizeof(*(*facePtrs[N_DIMS+1]))*numFaces);
 	for(i=0;i<numFaces;i++){
 	  for(j=0;j<numDims+1;j++){
 	    face = extractFace(numDims, vertexCoords, dc, i, j);
@@ -634,8 +649,14 @@ and filled as
     }
   }
 
-  if(numEntryFaces<=0)
+  if(numEntryFaces<=0){
+    *entryIntcpt = initializeIntersect(numDims);
+    *lenChainPtrs=0;
+    *chainOfCellIds=NULL;
+    *cellExitIntcpts=NULL;
+
     return 2;
+  }
 
   *lenChainPtrs = 1024; /* This can be increased within followCellChain(). */
   *chainOfCellIds  = malloc(sizeof(**chainOfCellIds) *(*lenChainPtrs));
@@ -654,9 +675,16 @@ and filled as
 
   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. */
-
-//*** this is not too good because *entryIntcpt, *chainOfCellIds, *cellExitIntcpts and lenChainPtrs are left at unsuitable values if status!=0.
+    /* 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. */
+
+  else{
+    *entryIntcpt = initializeIntersect(numDims);
+    free(*chainOfCellIds);
+    *chainOfCellIds=NULL;
+    free(*cellExitIntcpts);
+    *cellExitIntcpts=NULL;
+    *lenChainPtrs=0;
+  }
 
   free(cellVisited);
 

src/raytrace.c

@@ -313,7 +313,7 @@ traceray_smooth(rayData ray, const int im\
   , configInfo *par, struct grid *gp, double *vertexCoords, molData *md\
   , imageInfo *img, struct simplex *dc, const unsigned long numCells\
   , const double epsilon, gridInterp gips[3], const int numSegments\
-  , const double oneOnNumSegments, const int nSteps, const double oneOnNSteps){
+  , const double oneOnNumSegments){
   /*
 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.
 
@@ -322,7 +322,7 @@ Note that the algorithm employed here to solve the RTE is similar to that employ
 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 numFaces = DIM+1, nVertPerFace=3;
-  int ichan,stokesId,di,status,lenChainPtrs,entryI,exitI,vi,vvi,ci;
+  int ichan,stokesId,di,status,lenChainPtrs=0,entryI,exitI,vi,vvi,ci;
   int si,molI,lineI;
   double xp,yp,zp,x[DIM],dir[DIM],projVelRay,vel[DIM];
   double xCmpntsRay[nVertPerFace], ds, snu_pol[3], dtau, contJnu, contAlpha;
@@ -891,7 +891,7 @@ How to calculate this distance? Well if we have N points randomly but evenly dis
       locateRayOnImage(xs, pixelSize, imgCentreXPixels, imgCentreYPixels, img, im, maxNumRaysPerPixel, rays, &numActiveRays);
     }
   }
-  oneOnNumActiveRaysMinus1 = 1.0/(double)(numActiveRays-1);
+  oneOnNumActiveRaysMinus1 = 1.0/(double)(numActiveRaysInternal-1);
 
   if(numActiveRays<par->pIntensity+numCircleRays)
     rays = realloc(rays, sizeof(rayData)*numActiveRays);
@@ -970,7 +970,7 @@ While this is off however, gsl_* calls will not exit if they encounter a problem
       else if(par->traceRayAlgorithm==1)
         traceray_smooth(rays[ri], im, par, gp, vertexCoords, md, img\
           , cells, numCells, epsilon, gips\
-          , numSegments, oneOnNumSegments, nStepsThruCell, oneOnNSteps);
+          , numSegments, oneOnNumSegments);
 
       if (threadI == 0){ /* i.e., is master thread */
         if(!silent) progressbar((double)(ri)*oneOnNumActiveRaysMinus1, 13);
@@ -1029,7 +1029,7 @@ While this is off however, gsl_* calls will not exit if they encounter a problem
     unsigned long num2DCells,gis[3];
     intersectType entryIntcptFirstCell,*cellExitIntcpts=NULL;
     unsigned long *chainOfCellIds=NULL,*rasterCellIDs=NULL;
-    int lenChainPtrs,status=0,startYi,si;
+    int lenChainPtrs=0,status=0,startYi,si;
     double triangle[3][2],barys[3],y,deltaY;
     _Bool *rasterPixelIsInCells=NULL;