Several minor fixes.

Author: imckstewart

doc/usermanual.rst

@@ -321,6 +321,8 @@ thousands and about ten thousand.
 
 If this is left at the default value of 0, grid point sampling is performed according to the LIME<1.7 algorithm, as governed by parameter :ref:`par->sampling <par-sampling>`. If 1 is chosen, a new algorithm is employed which can quickly generate points with a distribution which accurately follows any feasible :ref:`gridDensity <grid-density>` function - including with sharp step-changes. This algorithm also incorporates a quasi-random choice of point candidates which avoids the requirement for the relatively time-consuming post-gridding smoothing phase.
 
+A user who selects par->samplingAlgorithm=1 and constructs their own :ref:`gridDensity <grid-density>` function obtains full control over the distribution of points. With this control however come some hazards. LIME still relies on 3rd-party software called qhull to triangulate the points after they are chosen, and qhull is a little flaky. It is prone to failing silently if it doesn't like the set of points one gives it. We have tried to trap these instances, to at least head off segmentation faults, but it is hard to guess all the ways in which somebody else's package may fail. If you have problems, try to smooth out any steps in your :ref:`gridDensity <grid-density>` function. If that doesn't fix things, you may have to go back to par->samplingAlgorithm=0.
+
 .. _par-sampling:
 
 .. code:: c

src/aux.c

@@ -723,6 +723,13 @@ exit(1);
 exit(1);
   }
 
+  if(!silent && !par->doMolCalcs && par->init_lte)
+    warning("Your choice of par.init_lte will have no effect.");
+  if(!silent && !par->doMolCalcs && par->lte_only)
+    warning("Your choice of par.lte_only will have no effect.");
+  if(!silent && par->nSolveIters>0 && par->lte_only)
+    warning("Requesting par.nSolveIters>0 will have no effect if LTE calculation is also requested.");
+
   /* Allocate moldata array.
   */
   if(par->nSpecies>0){

src/grid.c

@@ -419,6 +419,14 @@ Elements of structs are set as follows:
     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].numNeigh<=0){
+      if(!silent){
+        sprintf(message, "qhull failed silently, grid point %lu has 0 neighbours. Smoother gridDensity() might help.", id);
+        bail_out(message);
+      }
+exit(1);
+    }
+
     free(gp[id].neigh);
     gp[id].neigh=malloc(sizeof(struct grid *)*gp[id].numNeigh);
     for(k=0;k<gp[id].numNeigh;k++) {
@@ -498,7 +506,7 @@ Elements of structs are set as follows:
                 sprintf(message, "Something weird going on. Cannot find a cell with ID %lu", (unsigned long)(neighbor->id));
                 bail_out(message);
               }
-              exit(1);
+exit(1);
             }
           }
           i++;
@@ -1080,7 +1088,7 @@ exit(1);
 exit(1);
   }
   if(gridInfoRead.nSpecies > 0){
-    if((int)gridInfoRead.nSpecies!=par->nSpecies){
+    if((int)gridInfoRead.nSpecies!=par->nSpecies && par->doMolCalcs){
       if(!silent){
         sprintf(message, "Grid file had %d species but you have provided moldata files for %d."\
           , (int)gridInfoRead.nSpecies, par->nSpecies);
@@ -1098,6 +1106,9 @@ exit(1);
     }
 exit(1);
   }
+
+  if(!silent && par->nSolveItersDone>0 && (par->init_lte || par->lte_only))
+    warning("Your choice of LTE calculation will erase the RTE solution you read from file.");
 }
 
 /*....................................................................*/
@@ -1155,10 +1166,10 @@ readOrBuildGrid(configInfo *par, struct grid **gp){
     status = readGrid(par->gridInFile, &gridInfoRead, desiredKwds\
       , numDesiredKwds, gp, &collPartNames, &numCollPartRead, &(par->dataFlags));
 
-    sanityCheckOfRead(status, par, gridInfoRead);
-
     par->nSolveItersDone = desiredKwds[1].intValue;
 
+    sanityCheckOfRead(status, par, gridInfoRead);
+
     freeKeywords(desiredKwds, numDesiredKwds);
     freeGridInfo(&gridInfoRead);
 
@@ -1373,7 +1384,8 @@ Generate the remaining values if needed.
     par->dataFlags |= DS_mask_density;
   }
 
-  checkGridDensities(par, *gp); /* Check that none of the density samples is too small. */
+  if(par->doMolCalcs)
+    checkGridDensities(par, *gp); /* Check that none of the density samples is too small. */
 
   if(!allBitsSet(par->dataFlags, DS_mask_temperatures)){
     for(i=0;i<par->pIntensity;i++)

src/lime.h

@@ -258,11 +258,6 @@ typedef struct {
   _Bool doInterpolateVels;
 } imageInfo;
 
-typedef struct {
-  double x,y, *intensity, *tau;
-  unsigned int ppi;
-} rayData;
-
 struct blend{
   int molJ, lineJ;
   double deltaV;
@@ -364,7 +359,7 @@ void	lineBlend(molData*, configInfo*, struct blendInfo*);
 void	LTE(configInfo*, struct grid*, molData*);
 void	lteOnePoint(molData*, const int, const double, double*);
 void	mallocAndSetDefaultGrid(struct grid**, const size_t, const size_t);
-void	molInit(configInfo*, molData*, int*, const int);
+void	molInit(configInfo*, molData*);
 void	openSocket(char*);
 void	parseInput(inputPars, image*, const int, configInfo*, imageInfo**, molData**);
 void	photon(int, struct grid*, molData*, const gsl_rng*, configInfo*, const int, struct blendInfo, gridPointData*, double*);
@@ -377,7 +372,6 @@ void	processFitsError(int);
 double	ratranInput(char*, char*, double, double, double);
 void	raytrace(int, configInfo*, struct grid*, molData*, imageInfo*, double*, double*, const int);
 void	readDustFile(char*, double**, double**, int*);
-void	readMolData(configInfo*, molData*, int**, int*);
 void	readOrBuildGrid(configInfo*, struct grid**);
 void	readUserInput(inputPars*, imageInfo**, int*, int*);
 unsigned long reorderGrid(const unsigned long, struct grid*);

src/main.c

@@ -169,8 +169,6 @@ run(inputPars inpars, image *inimg, const int nImages){
   char message[80];
   int nEntries=0;
   double *lamtab=NULL,*kaptab=NULL;
-  int *allUniqueCollPartIds=NULL;
-  int numUniqueCollPartsFound;
 
   /*Set locale to avoid trouble when reading files*/
   setlocale(LC_ALL, "C");
@@ -215,16 +213,13 @@ run(inputPars inpars, image *inimg, const int nImages){
     }
   }
 
-  if(!popsdone)
-    readMolData(&par, md, &allUniqueCollPartIds, &numUniqueCollPartsFound);
-
   if(par.doMolCalcs){
     if(!popsdone){
-      molInit(&par, md, allUniqueCollPartIds, numUniqueCollPartsFound);
+      molInit(&par, md);
       calcGridMolDoppler(&par, md, gp);
     }
     if(par.useAbun)
-      calcGridMolDensities(&par,&gp);
+      calcGridMolDensities(&par, &gp);
 
     for(gi=0;gi<par.ncell;gi++){
       for(si=0;si<par.nSpecies;si++)
@@ -240,8 +235,6 @@ run(inputPars inpars, image *inimg, const int nImages){
     par.nSolveItersDone = par.nSolveIters;
   }
 
-  free(allUniqueCollPartIds);
-
   if(onlyBitsSet(par.dataFlags & DS_mask_all_but_mag, DS_mask_3))
     writeGridIfRequired(&par, gp, NULL, 3);
   else if(onlyBitsSet(par.dataFlags & DS_mask_all_but_mag, DS_mask_5)){

src/molinit.c

@@ -309,11 +309,15 @@ void setUpGir(configInfo *par, molData *md){
 }
 
 /*....................................................................*/
-void molInit(configInfo *par, molData *md, int *allUniqueCollPartIds, const int numUniqueCollPartsFound){
+void molInit(configInfo *par, molData *md){
   int i,j,jStart,numActiveCollParts;
   char partstr[90];
+  int *allUniqueCollPartIds=NULL;
+  int numUniqueCollPartsFound;
 
+  readMolData(par, md, &allUniqueCollPartIds, &numUniqueCollPartsFound);
   setUpDensityAux(par, allUniqueCollPartIds, numUniqueCollPartsFound);
+
   if(par->girdatfile!=NULL){
     setUpGir(par, md);
   }
@@ -358,5 +362,6 @@ void molInit(configInfo *par, molData *md, int *allUniqueCollPartIds, const int
   }
 
   calcMolCMBs(par, md);
+  free(allUniqueCollPartIds);
 }
 

src/predefgrid.c

@@ -57,7 +57,11 @@ predefinedGrid(configInfo *par, struct grid *gp){
     if(!silent) progressbar((double) i/((double)par->pIntensity-1), 4);	
   }
 
-  checkGridDensities(par, gp);
+  fclose(fp);
+  if(!silent) printDone(4);
+
+  if(par->doMolCalcs)
+    checkGridDensities(par, gp);
 
   for(i=par->pIntensity;i<par->ncell;i++){
     x=2*gsl_rng_uniform(ran)-1.;
@@ -82,7 +86,6 @@ predefinedGrid(configInfo *par, struct grid *gp){
       for(j=0;j<DIM;j++) gp[i].vel[j]=0.;
     } else i--;
   }
-  fclose(fp);
 
   delaunay(DIM, gp, (unsigned long)par->ncell, 0, 1, &dc, &numCells);
 

src/raytrace.c

@@ -13,6 +13,12 @@
 #include "lime.h"
 #include "raythrucells.h"
 
+typedef struct {
+  double x,y, *intensity, *tau;
+  unsigned int ppi;
+  _Bool isInsideImage;
+} rayData;
+
 struct baryVelBuffType {
   int numVertices,numEdges,(*edgeVertexIndices)[2];
   double **vertexVels,**edgeVels,*entryCellBary,*midCellBary,*exitCellBary,*shapeFns;
@@ -844,8 +850,31 @@ At the moment I will fix the number of segments, but it might possibly be faster
 }
 
 /*....................................................................*/
-void
+_Bool
 locateRayOnImage(double x[2], const double size, const double imgCentreXPixels\
+  , const double imgCentreYPixels, imageInfo *img, const int im\
+  , int *xi, int *yi, unsigned int *ppi){
+
+  _Bool isInsideImage;
+
+  /* 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){
+    isInsideImage = 0;
+    *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{
+    isInsideImage = 1;
+    *ppi = (unsigned int)(*yi)*(unsigned int)img[im].pxls + (unsigned int)(*xi);
+  }
+
+  return isInsideImage;
+}
+
+/*....................................................................*/
+void
+assignRayOnImage(double x[2], const double size, const double imgCentreXPixels\
   , const double imgCentreYPixels, imageInfo *img, const int im\
   , const int maxNumRaysPerPixel, rayData *rays, int *numActiveRays){
   /*
@@ -863,27 +892,19 @@ Returned information is thus:
   */
 
   int xi,yi,ichan;
-  _Bool isOutsideImage;
+  _Bool isInsideImage;
   unsigned int ppi;
 
-  /* 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;
-  }
+  isInsideImage = locateRayOnImage(x, size, imgCentreXPixels\
+    , imgCentreYPixels, img, im, &xi, &yi, &ppi);
 
-  /* 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.
+  /* See if we want to keep the ray. For the time being we will include those outside the image bounds, because we need to interpolate between outside-image rays and inside-image ones, but a cleverer algorithm would exclude all but those which are affected by this (i.e. that immediately abut the image bounds). 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)
+  if(!isInsideImage || maxNumRaysPerPixel<1 || img[im].pixel[ppi].numRays<maxNumRaysPerPixel){
+    if(isInsideImage)
       img[im].pixel[ppi].numRays++;
 
+    rays[*numActiveRays].isInsideImage = isInsideImage;
     rays[*numActiveRays].ppi = ppi;
     rays[*numActiveRays].x = x[0];
     rays[*numActiveRays].y = x[1];
@@ -1221,7 +1242,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, &numActiveRaysInternal);
+    assignRayOnImage(xs, pixelSize, imgCentreXPixels, imgCentreYPixels, img, im, maxNumRaysPerPixel, rays, &numActiveRaysInternal);
   } /* End loop 1, over grid points. */
 
   /* Add the circle rays:
@@ -1233,7 +1254,7 @@ How to calculate this distance? Well if we have N points randomly but evenly dis
       angle = i*scale;
       xs[0] = par->radius*cos(angle);
       xs[1] = par->radius*sin(angle);
-      locateRayOnImage(xs, pixelSize, imgCentreXPixels, imgCentreYPixels, img, im, maxNumRaysPerPixel, rays, &numActiveRays);
+      assignRayOnImage(xs, pixelSize, imgCentreXPixels, imgCentreYPixels, img, im, maxNumRaysPerPixel, rays, &numActiveRays);
     }
   }
   oneOnNumActiveRaysMinus1 = 1.0/(double)(numActiveRaysInternal-1);
@@ -1362,6 +1383,7 @@ While this is off however, gsl_* calls will not exit if they encounter a problem
   } /* End of parallel block. */
 
   gsl_set_error_handler(defaultErrorHandler);
+  if(!silent) printDone(13);
 
   if(par->traceRayAlgorithm==1){
     free(cells);
@@ -1393,7 +1415,7 @@ While this is off however, gsl_* calls will not exit if they encounter a problem
   /* 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){
+    if(rays[ri].isInsideImage && 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];