New optional parameter collPartNames added.

example/model.c

@@ -136,7 +136,7 @@ density(double x, double y, double z, double *density){
    * Define variable for radial coordinate
    */
   double r, rToUse;
-  const double rMin = 0.1*AU; /* This cutoff should be chosen smaller than par->minScale but greater than zero (to avoid a singularity at the origin). */
+  const double rMin = 0.5*AU; /* This cutoff should be chosen smaller than par->minScale but greater than zero (to avoid a singularity at the origin). */
 
   /*
    * Calculate radial distance from origin

src/aux.c

@@ -185,6 +185,8 @@ The parameters visible to the user have now been strictly confined to members of
   for(i=0;i<MAX_N_COLL_PART;i++) par->nMolWeights[i] = inpar.nMolWeights[i];
   par->dustWeights  = malloc(sizeof(double)*MAX_N_COLL_PART);
   for(i=0;i<MAX_N_COLL_PART;i++) par->dustWeights[i] = inpar.dustWeights[i];
+  par->collPartNames = malloc(sizeof(char*)*MAX_N_COLL_PART);
+  for(i=0;i<MAX_N_COLL_PART;i++) copyInparStr(inpar.collPartNames[i], &(par->collPartNames[i]));
 
   /* Copy over the grid-density-maximum data and find out how many were set:
   */
@@ -604,11 +606,33 @@ Eventually I hope readOrBuildGrid() will be unilaterally called within LIME; if
 
 void checkUserDensWeights(configInfo *par){
   /*
-This deals with three user-settable vectors: par->collPartIds, par->nMolWeights and par->dustWeights. We have to see if these (optional) parameters were set, do some basic checks on them, and make sure they have the same numbers of elements as the number of density values, which by this time should be stored in par->numDensities.
+This deals with four user-settable vectors: par->collPartIds, par->nMolWeights, par->dustWeights and par->collPartNames. We have to see if these (optional) parameters were set, do some basic checks on them, and make sure they match the number of density values, which by this time should be stored in par->numDensities.
+
+	* par->collPartIds: this list acts as a link between the N density function returns (I'm using here N as shorthand for par->numDensities) and the M collision partner ID integers found in the moldatfiles. This allows us to associate density functions with the collision partner transition rates provided in the moldatfiles.
+
+	* par->collPartNames: essentially this has only cosmetic importance since it has no effect on the functioning of LIME, only on the names of the collision partners which are printed to stdout. It's main purpose is to reassure the user who has provided transition rates for a non-LAMDA collision species in their moldatfile that they are actually getting these values and not some mysterious reversion to LAMDA.
+
+	The user can specify either, none, or both of these two parameters, with the following effects:
+
+		Ids	Names	Effect
+		----------------------
+		0	0	LAMDA coll parts are assumed and the association between the density functions and the moldatfiles is essentially guessed.
+
+		0	1	par->collPartIds is constructed to contain integers in a sequence from 1 to N. Naturally the user should write matching collision partner ID integers in their moldatfiles.
+
+		1	0	LAMDA coll parts are assumed.
+
+		1	1	User will get what they ask for. Elements of par->collPartIds are taken (after subtracting 1) to be indices to par->collPartNames; thus in this situation par->collPartNames may have a different order to the density functions, or have more values than par->numDensities. Any extra values are ignored however, and par->collPartNames must contain a value for each par->collPartIds element.
+		----------------------
+
+	* par->nMolWeights: this list gives the weights to be applied to the N density values when calculating molecular densities from abundances.
+
+	* par->dustWeights: ditto for the calculation of dust opacity.
   */
-  int i,j,numUserSetCPIds,numUserSetNMWs,numUserSetDWs;
+  int i,j,numUserSetCPIds,numUserSetNMWs,numUserSetDWs,numUserSetCPNames;
   int *uniqueCPIds=NULL;
   double sum;
+  char message[80];
 
   /* Get the number of par->collPartIds set by the user:
   */
@@ -670,7 +694,47 @@ This deals with three user-settable vectors: par->collPartIds, par->nMolWeights
     }
   }
 
-  /* Check if we have either 0 par->collPartIds or the same number as the number of density values.
+  /* Get the number of par->collPartNames set by the user:
+  */
+  i = 0;
+  while(i<MAX_N_COLL_PART && par->collPartNames[i]!=NULL) i++;
+  numUserSetCPNames = i;
+
+  /*
+Re the interaction between par->collPartIds and par->collPartNames: there are three possible scenarios as follows.
+	- User wishes to use only collision partners from the set in the LAMDA database, in which case there is no need to set par->collPartNames. Note that the integer ID values for the collision partners in the standard LAMDA moldatfiles begin at 1.
+
+	- User wishes to use some non-LAMDA collision partners: it is perfectly possible to replace the values in a moldatfile for any LAMDA species with values calculated for another; LIME has no internal knowledge of whether a collision partner is LAMDA or not, it just uses the transition rates you give it. The only use it makes of collision partner names is in the the initial greeting message it prints. If the user wishes to have the correct name(s) printed, this can be transmitted to LIME via the par->collPartNames parameter.
+
+	For added convenience, it is allowed in this case that the user can supply just par->collPartNames and not par->collPartIds.
+		* If this is done, the number of elements of par->collPartNames is required to be the same as the number of density values. A list of par->collPartIds of the same length will be constructed, the elements being {1,2,3...} etc. These ID integers must match the ones supplied by the user in their moldatfiles.
+
+		* If the user supplies both par->collPartIds and par->collPartNames, the values in par->collPartIds must (after subtraction of 1) be the indices of that collision partner in par->collPartNames, as well of course as matching the value the user has supplied in the moldatfiles.
+  */
+  if(numUserSetCPNames>0){
+    if(numUserSetCPNames <= par->numDensities){
+      if(!silent) bail_out("There must be at least 1 value of par.collPartNames for each density() return.");
+      exit(1);
+    }
+
+    if(numUserSetCPIds<=0){
+      for(i=0;i<par->numDensities;i++)
+        par->collPartIds[i] = i+1;
+      numUserSetCPIds = par->numDensities;
+    }else{
+      for(i=0;i<numUserSetCPIds;i++){
+        if(par->collPartIds[i]<1 || par->collPartIds[i]>numUserSetCPNames){
+          if(!silent){
+            sprintf(message, "par.collPartIds[%d] value %d is out of range (<1 or >%d)", i, par->collPartIds[i], numUserSetCPNames);
+            bail_out(message);
+          }
+          exit(1);
+        }
+      }
+    }
+  }
+
+  /* Check that we have either 0 par->collPartIds or the same number as the number of density values. If not, the par->collPartIds values the user set will be thrown away, the pointer will be reallocated, and new values will be written to it in setUpDensityAux(), taken from the values in the moldatfiles.
   */
   if(numUserSetCPIds != par->numDensities){
     free(par->collPartIds);

src/frees.c

@@ -16,6 +16,11 @@ freeConfigInfo(configInfo par){
   free(par.nMolWeights);
   free(par.dustWeights);
   free(par.collPartIds);
+  if(par.collPartNames!= NULL){
+    for(i=0;i<MAX_N_COLL_PART;i++)
+      free(par.collPartNames[i]);
+    free(par.collPartNames);
+  }
 
   free(par.outputfile);
   free(par.binoutputfile);
@@ -34,6 +39,18 @@ freeConfigInfo(configInfo par){
   free(par.gridDensMaxLoc);
 }
 
+void
+freeInputPars(inputPars par){
+  free(par.collPartIds);
+  free(par.nMolWeights);
+  free(par.dustWeights);
+  free(par.moldatfile);
+  free(par.collPartNames);
+  free(par.gridOutFiles);
+  free(par.gridDensMaxValues);
+  free(par.gridDensMaxLoc);
+}
+
 void
 freeGrid(const unsigned int numPoints, const unsigned short numSpecies\
   , struct grid *gp){

src/grid.c

@@ -623,8 +623,14 @@ write_VTK_unstructured_Points(configInfo *par, struct grid *g){
   }
   fprintf(fp,"SCALARS Mol_density float 1\n");
   fprintf(fp,"LOOKUP_TABLE default\n");
-  for(i=0;i<par->ncell;i++){
-    fprintf(fp, "%e\n", g[i].abun[0]*g[i].dens[0]);
+  if(par->nSpecies>0){
+    for(i=0;i<par->ncell;i++){
+      fprintf(fp, "%e\n", g[i].abun[0]*g[i].dens[0]);
+    }
+  }else{
+    for(i=0;i<par->ncell;i++){
+      fprintf(fp, "%e\n", 0.0);
+    }
   }
   fprintf(fp,"SCALARS Gas_temperature float 1\n");
   fprintf(fp,"LOOKUP_TABLE default\n");
@@ -991,7 +997,7 @@ readOrBuildGrid(configInfo *par, struct grid **gp){
   int i,j,k,di,si,levelI=0,status=0,numCollPartRead;
   double theta,semiradius,z,dummyScalar;
   double *outRandDensities=NULL,*dummyPointer=NULL,x[DIM];
-  double (*outRandLocations)[DIM]=NULL,dummyTemp[2],dummyVel[DIM];
+  double (*outRandLocations)[DIM]=NULL,dummyTemp[]={-1.0,-1.0},dummyVel[DIM];
   treeRandConstType rinc;
   treeRandVarType rinv;
   struct cell *dc=NULL; /* Not used at present. */

src/inpars.h

@@ -21,7 +21,7 @@ typedef struct {
   char *restart;
   char *dust;
   int sampling,lte_only,init_lte,antialias,polarization,nThreads;
-  char **moldatfile;
+  char **moldatfile,**collPartNames;
   char *gridInFile,**gridOutFiles;
   int nSolveIters;
   double (*gridDensMaxLoc)[DIM], *gridDensMaxValues;

src/lime.h

@@ -86,7 +86,7 @@
 #define FAST_EXP_MAX_TAYLOR     3
 #define FAST_EXP_NUM_BITS       8
 #define NUM_GRID_STAGES         4
-#define MAX_N_COLL_PART         7
+#define MAX_N_COLL_PART         20
 #define N_SMOOTH_ITERS          20
 #define TYPICAL_ISM_DENS        1000.0
 #define STR_LEN_0               80
@@ -153,16 +153,17 @@ typedef struct {
   char *dust;
   int sampling,lte_only,init_lte,antialias,polarization,nThreads,numDims;
   int nLineImages,nContImages;
-  char **moldatfile;
+  char **moldatfile,**collPartNames;
   _Bool writeGridAtStage[NUM_GRID_STAGES],resetRNG,doInterpolateVels;
   char *gridInFile,**gridOutFiles;
   int dataFlags,nSolveIters;
-  double (*gridDensMaxLoc)[DIM], *gridDensMaxValues;
+  double (*gridDensMaxLoc)[DIM],*gridDensMaxValues;
 } configInfo;
 
 struct cpData {
   double *down,*temp;
   int collPartId,ntemp,ntrans,*lcl,*lcu,densityIndex;
+  char *name;
 };
 
 /* Molecular data: shared attributes */
@@ -338,6 +339,7 @@ void	freeConfigInfo(configInfo par);
 void	freeGrid(const unsigned int, const unsigned short, struct grid*);
 void	freeGridPointData(const int, gridPointData*);
 void	freeImgInfo(const int, imageInfo*);
+void	freeInputPars(inputPars par);
 void	freeMolData(const int, molData*);
 void	freeMolsWithBlends(struct molWithBlends*, const int);
 void	freePopulation(const unsigned short, struct populations*);
@@ -402,7 +404,7 @@ void	write_VTK_unstructured_Points(configInfo*, struct grid*);
 void	bail_out(char*);
 void	casaStyleProgressBar(const int, int);
 void	collpartmesg(char*, int);
-void	collpartmesg2(char*, int);
+void	collpartmesg2(char*);
 void	collpartmesg3(int, int);
 void	goodnight(int, char*);
 void	greetings(void);

src/main.c

@@ -98,6 +98,12 @@ initParImg(inputPars *par, image **img)
     par->moldatfile[id]=NULL;
   }
 
+  /* Allocate initial space for (non-LAMDA) collision partner names */
+  par->collPartNames=malloc(sizeof(char *)*MAX_N_COLL_PART);
+  for(i=0;i<MAX_N_COLL_PART;i++){
+    par->collPartNames[i]=NULL;
+  }
+
   /* Allocate initial space for output fits images */
   (*img)=malloc(sizeof(**img)*MAX_NIMAGES);
   for(i=0;i<MAX_NIMAGES;i++){
@@ -292,13 +298,7 @@ int main() {
   run(par, img, nImages);
 
   free(img);
-  free(par.collPartIds);
-  free(par.nMolWeights);
-  free(par.dustWeights);
-  free(par.moldatfile);
-  free(par.gridOutFiles);
-  free(par.gridDensMaxValues);
-  free(par.gridDensMaxLoc);
+  freeInputPars(par);
 
   return 0;
 }

src/messages.c

@@ -315,7 +315,7 @@ collpartmesg(char molecule[80], int partners){//, int specnumber){
 }
 
 void
-collpartmesg2(char name[10], int partner){
+collpartmesg2(char name[10]){
 #ifdef NO_NCURSES
   printf("      %s\n ", name);
 #else