Addition of facility to read and write the entire grid to FITS file.

.gitignore

@@ -8,9 +8,9 @@
 example/*.vtk
 example/*.pop
 example/*.fits*
+example/lime_*.x
 example/*.dat
 example/*_archive
-example/lime_*.x
 example/vlog.txt
 
 # Ignore html doc files

Makefile

@@ -49,26 +49,26 @@ endif
 ##
 
 TARGET  = lime.x 
-#CC	= gcc -fopenmp -g
+#CC	= gcc -fopenmp -g -Wunused -Wno-unused-result
 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/tree_random.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/tree_random.c		\
+	  src/velospline.c src/getclosest.c src/grid2fits.c		\
+	  src/tcpsocket.c src/defaults.c src/fastexp.c src/gridio.c	\
 	  src/raythrucells.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/tree_random.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/raytrace.o			\
+	  src/ratranInput.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/tree_random.o		\
+	  src/velospline.o src/getclosest.o src/grid2fits.o		\
+	  src/tcpsocket.o src/defaults.o src/fastexp.o src/gridio.o	\
 	  src/raythrucells.o
 MODELO 	= src/model.o
 

example/model.c

@@ -26,6 +26,7 @@ input(inputPars *par, image *img){
   par->moldatfile[0]            = "hco+@xpol.dat";
   par->antialias                = 4;
   par->sampling                 = 2; // log distr. for radius, directions distr. uniformly on a sphere.
+  par->nSolveIters              = 14;
 
   par->outputfile               = "populations.pop";
   par->binoutputfile            = "restart.pop";
@@ -81,6 +82,17 @@ input(inputPars *par, image *img){
   par->nMolWeights[0]           = 1.0;
   par->dustWeights[0]           = 1.0;
 
+  /* Set one or more of the following parameters for full output of the grid-specific data at any of 4 stages during the processing. (See the header of gridio.c for information about the stages.)
+  par->gridOutFiles[0] = "grid_1.ds";
+  par->gridOutFiles[1] = "grid_2.ds";
+  par->gridOutFiles[2] = "grid_3.ds";
+  par->gridOutFiles[3] = "grid_4.ds";
+  */
+
+  /* You can also optionally read in a FITS file stored via the previous parameters, or prepared externally. See the header of grid2fits.c for information about the correct file format. LIME can cope with almost any sensible subset of the recognized columns; it will use the file values if they are present, then calculate the missing ones.
+  par->gridInFile = "grid_4.ds";
+  */
+
   /*
    * Definitions for image #0. Add blocks with successive values of i for additional images.
    */

src/aux.c

@@ -47,14 +47,24 @@ parseInput(inputPars inpar, configInfo *par, image **img, molData **m){
   par->antialias    = inpar.antialias;
   par->polarization = inpar.polarization;
   par->nThreads     = inpar.nThreads;
+  par->gridInFile   = inpar.gridInFile;
+  par->nSolveIters  = inpar.nSolveIters;
   par->traceRayAlgorithm = inpar.traceRayAlgorithm;
 
+  par->gridOutFiles = malloc(sizeof(char *)*NUM_GRID_STAGES);
+  for(i=0;i<NUM_GRID_STAGES;i++)
+    par->gridOutFiles[i] = inpar.gridOutFiles[i];
+
   /* Now set the additional values in par. */
   par->ncell = inpar.pIntensity + inpar.sinkPoints;
   par->radiusSqu = inpar.radius*inpar.radius;
   par->minScaleSqu=inpar.minScale*inpar.minScale;
   par->doPregrid = (inpar.pregrid==NULL)?0:1;
 
+  for(i=0;i<NUM_GRID_STAGES;i++)
+    par->writeGridAtStage[i] = 0;
+  par->dataFlags = 0;
+
   /* If the user has provided a list of moldatfile names, the corresponding elements of par->moldatfile will be non-NULL. Thus we can deduce the number of files (species) from the number of non-NULL elements.
   */
   par->nSpecies=0;
@@ -104,6 +114,11 @@ parseInput(inputPars inpar, configInfo *par, image **img, molData **m){
   while(i<MAX_N_HIGH && par->gridDensMaxValues[i]>=0) i++;
   par->numGridDensMaxima = i;
 
+  /* Check that the user has supplied the velocity function (needed in raytracing unless par->doPregrid). Note that the other previously mandatory functions (density, abundance, doppler and temperature) may not be necessary if the user reads in the appropriate values from a file. This is tested at the appropriate place in readOrBuildGrid().
+  */
+  if(!par->doPregrid || par->traceRayAlgorithm==1)
+    velocity(0.0,0.0,0.0, dummyVel);
+
   /* Calculate par->numDensities.
   */
   if(!(par->doPregrid || par->restart)){ /* These switches cause par->numDensities to be set in routines they call. */
@@ -143,10 +158,10 @@ 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->doPregrid):
-  */
-  if(!par->doPregrid || par->traceRayAlgorithm==1)
-    velocity(0.0,0.0,0.0, dummyVel);
+  for(i=0;i<NUM_GRID_STAGES;i++){
+    if(par->gridOutFiles[i] != NULL)
+      par->writeGridAtStage[i] = 1;
+  };
 
   /*
 Now we need to calculate the cutoff value used in calcSourceFn(). The issue is to decide between
@@ -680,7 +695,7 @@ Pointers are indicated by a * before the attribute name and an arrow to the memo
 
 void
 levelPops(molData *md, configInfo *par, struct grid *gp, int *popsdone, double *lamtab, double *kaptab, const int nEntries){
-  int id,conv=0,iter,ilev,prog=0,ispec,c=0,n,i,threadI,nVerticesDone,nlinetot;
+  int id,iter,ilev,ispec,c=0,n,i,threadI,nVerticesDone,nItersDone,nlinetot;//,conv=0
   double percent=0.,*median,result1=0,result2=0,snr,delta_pop;
   int nextMolWithBlend;
   struct statistics { double *pop, *ave, *sigma; } *stat;
@@ -749,8 +764,9 @@ This is done to allow proper handling of errors which may arise in the LU solver
 While this is off however, other gsl_* etc calls will not exit if they encounter a problem. We may need to pay some attention to trapping their errors.
     */
 
-    do{
-      if(!silent) progressbar2(0, prog++, 0, result1, result2);
+    nItersDone=0;
+    while(nItersDone < par->nSolveIters){ /* Not a 'for' loop because we will probably later want to add a convergence criterion. */
+      if(!silent) progressbar2(par, 0, nItersDone, 0, result1, result2);
 
       for(id=0;id<par->pIntensity;id++){
         for(ilev=0;ilev<md[0].nlev;ilev++) {
@@ -840,15 +856,16 @@ While this is off however, other gsl_* etc calls will not exit if they encounter
       }
 
       gsl_sort(median, 1, c);
-      if(conv>1){
+      if(nItersDone>1){
         result1=median[0];
         result2 =gsl_stats_median_from_sorted_data(median, 1, c);
       }
       free(median);
 
-      if(!silent) progressbar2(1, prog, percent, result1, result2);
-      if(par->outputfile) popsout(par,gp,md);
-    } while(conv++<NITERATIONS);
+      if(!silent) progressbar2(par, 1, nItersDone, percent, result1, result2);
+      if(par->outputfile != NULL) popsout(par,gp,md);
+      nItersDone++;
+    }
     gsl_set_error_handler(defaultErrorHandler);
 
     freeMolsWithBlends(blends.mols, blends.numMolsWithBlends);
@@ -867,8 +884,46 @@ While this is off however, other gsl_* etc calls will not exit if they encounter
     free(stat);
   }
 
-  if(par->binoutputfile) binpopsout(par,gp,md);
+  par->dataFlags |= (1 << DS_bit_populations);
+
+  if(par->binoutputfile != NULL) binpopsout(par,gp,md);
 
   *popsdone=1;
 }
 
+_Bool allBitsSet(const int flags, const int mask){
+  /* Returns true only if all the masked bits of flags are set. */
+
+  if(~flags & mask)
+    return 0;
+  else
+    return 1;
+}
+
+_Bool anyBitSet(const int flags, const int mask){
+  /* Returns true if any of the masked bits of flags are set. */
+
+  if(flags & mask)
+    return 1;
+  else
+    return 0;
+}
+
+_Bool bitIsSet(const int flags, const int bitI){
+  /* Returns true if the designated bit of flags is set. */
+
+  if(flags & (1 << bitI))
+    return 1;
+  else
+    return 0;
+}
+
+_Bool onlyBitsSet(const int flags, const int mask){
+  /* Returns true if flags has no bits set apart from those which are true in mask. */
+
+  if(flags & ~mask)
+    return 0;
+  else
+    return 1;
+}
+

src/defaults.c

@@ -72,7 +72,6 @@ The grid points within the model are chosen randomly via the rejection method wi
 
       density(r[0],r[1],r[2],val);
       for (i=0;i<par->numDensities;i++) totalDensity += val[i];
-//      fracDensity = pow(totalDensity/par->gridDensGlobalMax,DENSITY_POWER);
       fracDensity = pow(totalDensity,DENSITY_POWER)/par->gridDensGlobalMax;
 
       return fracDensity;

src/frees.c

@@ -25,6 +25,9 @@ void freeGrid(const unsigned int numPoints, const unsigned short numSpecies\
 
   if(gp != NULL){
     for(i_u=0;i_u<numPoints;i_u++){
+      free(gp[i_u].v1);
+      free(gp[i_u].v2);
+      free(gp[i_u].v3);
       free(gp[i_u].dir);
       free(gp[i_u].neigh);
       free(gp[i_u].w);
@@ -115,6 +118,7 @@ freeParImg(const int nImages, inputPars *par, image *img){
   free(par->collPartIds);
   free(par->nMolWeights);
   free(par->dustWeights);
+  free(par->gridOutFiles);
 }
 
 void freePopulation(const unsigned short numSpecies, struct populations *pop){