diff --git a/.gitignore b/.gitignore
index 182c93e..d929ee0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,6 +5,8 @@
 .DS_Store
 
 # Ignore example files
-example/grid.vtk
+example/*.vtk
 example/*.pop
-example/*.fits
+example/*.fits*
+example/*.dat
+
diff --git a/doc/usermanual.rst b/doc/usermanual.rst
index ada64a6..a25a18f 100644
--- a/doc/usermanual.rst
+++ b/doc/usermanual.rst
@@ -420,13 +420,6 @@ abundance, velocity etc. There is no default value.
 If set, LIME performs an LTE calculation. Useful for quick checks. The
 default lte\_only=0, i.e., full non-LTE calculation.
 
-.. code:: c
-
-    (integer) par->init_lte (optional)
-
-If set, LIME use LTE approximation as initial one for subsequent non-LTE calculations. The
-default init\_lte=0, i.e., the code will use constant value for level populations as initial solution.
-
 .. code:: c
 
     (integer) par->blend (optional)
diff --git a/src/LTEsolution.c b/src/LTEsolution.c
index e81fa13..6ea7df3 100644
--- a/src/LTEsolution.c
+++ b/src/LTEsolution.c
@@ -10,7 +10,7 @@
 #include "lime.h"
 
 void
-LTE(inputPars *par, struct grid *g, molData *m){
+LTE(configInfo *par, struct grid *g, molData *m){
   int id,ispec;
 
   for(id=0;id<par->pIntensity;id++){
@@ -22,7 +22,7 @@ LTE(inputPars *par, struct grid *g, molData *m){
   if(par->outputfile) popsout(par,g,m);
 }
 
-void lteOnePoint(inputPars *par, molData *m, const int ispec, const double temp, double *pops){
+void lteOnePoint(configInfo *par, molData *m, const int ispec, const double temp, double *pops){
   int ilev;
   double sum;
 
diff --git a/src/aux.c b/src/aux.c
index 2b98adb..4b7c070 100644
--- a/src/aux.c
+++ b/src/aux.c
@@ -14,15 +14,71 @@
 
 
 void
-parseInput(inputPars *par, image **img, molData **m){
+parseInput(inputPars inpar, configInfo *par, image **img, molData **m){
   int i,id, ispec;
   double BB[3];
   double cosPhi,sinPhi,cosTheta,sinTheta,dummyVel[DIM];
+  FILE *fp;
+
+  /* Copy over user-set parameters to the configInfo versions. (This seems like duplicated effort but it is a good principle to separate the two structs, for several reasons, as follows. (i) We will usually want more config parameters than user-settable ones. The separation leaves it clearer which things the user needs to (or can) set. (ii) The separation allows checking and screening out of impossible combinations of parameters. (iii) We can adopt new names (for clarity) for config parameters without bothering the user with a changed interface.) */
+  par->radius       = inpar.radius;
+  par->minScale     = inpar.minScale;
+  par->pIntensity   = inpar.pIntensity;
+  par->sinkPoints   = inpar.sinkPoints;
+  par->sampling     = inpar.sampling;
+  par->tcmb         = inpar.tcmb;
+  par->dust         = inpar.dust;
+  par->outputfile   = inpar.outputfile;
+  par->binoutputfile= inpar.binoutputfile;
+  par->restart      = inpar.restart;
+  par->gridfile     = inpar.gridfile;
+  par->pregrid      = inpar.pregrid;
+  par->lte_only     = inpar.lte_only;
+  par->init_lte     = inpar.init_lte;
+  par->blend        = inpar.blend;
+  par->antialias    = inpar.antialias;
+  par->polarization = inpar.polarization;
+  par->nThreads     = inpar.nThreads;
+
+  /* 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;
+  while(inpar.moldatfile[par->nSpecies]!=NULL && par->nSpecies<MAX_NSPECIES)
+    par->nSpecies++;
+
+  /* Copy over the moldatfiles.
+  */
+  if(par->nSpecies == 0){
+    par->nSpecies = 1;
+    par->moldatfile = NULL;
+
+  } else {
+    par->moldatfile=malloc(sizeof(char *)*par->nSpecies);
+    for(id=0;id<par->nSpecies;id++){
+      par->moldatfile[id] = inpar.moldatfile[id];
+    }
+
+    /* Check if files exist. */
+    for(id=0;id<par->nSpecies;id++){
+      if((fp=fopen(par->moldatfile[id], "r"))==NULL) {
+        openSocket(par->moldatfile[id]);
+      } else {
+        fclose(fp);
+      }
+    }
+  }
+
+  /* If the user has provided a list of image filenames, the corresponding elements of (*img).filename will be non-NULL. Thus we can deduce the number of images from the number of non-NULL elements.
+  */
+  par->nImages=0;
+  while((*img)[par->nImages].filename!=NULL && par->nImages<MAX_NIMAGES)
+    par->nImages++;
 
-  par->ncell=par->pIntensity+par->sinkPoints;
-  par->radiusSqu=par->radius*par->radius;
-  par->minScaleSqu=par->minScale*par->minScale;
-  if(par->pregrid!=NULL) par->doPregrid=1;
+  /* 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;
 
   /* Check that the user has supplied this function (needed unless par->pregrid):
   */
@@ -52,29 +108,29 @@ The cutoff will be the value of abs(x) for which the error in the exact expressi
 
       /* Check for polarization */
       BB[0]=0.;
-      magfield(par->minScale,par->minScale,par->minScale,BB);
+      magfield(inpar.minScale,inpar.minScale,inpar.minScale,BB);
       if(fabs(BB[0]) > 0.) par->polarization=1;
 
       if(par->polarization) (*img)[i].nchan=3;
       else (*img)[i].nchan=1;
-      if((*img)[i].trans>-1 || (*img)[i].bandwidth>-1. || (*img)[i].freq==0 || par->dust==NULL){
-        if(!silent) bail_out("Error: Image keywords are ambiguous");
+      if((*img)[i].trans>-1 || (*img)[i].bandwidth>-1. || (*img)[i].freq==0 || inpar.dust==NULL){
+        if(!silent) bail_out("Image keywords are ambiguous");
         exit(1);
       }
       (*img)[i].doline=0;
     } else if (((*img)[i].nchan>0 || (*img)[i].velres > 0)){
       /* Assume line image */
       par->polarization=0;
-      if(par->moldatfile==NULL){
-        if(!silent) bail_out("Error: No data file is specified for line image.");
+      if(inpar.moldatfile==NULL){
+        if(!silent) bail_out("No data file is specified for line image.");
         exit(1);
       }
       if(((*img)[i].trans>-1 && (*img)[i].freq>-1) || ((*img)[i].trans<0 && (*img)[i].freq<0)){
-        if(!silent) bail_out("Error: Specify either frequency or transition ");
+        if(!silent) bail_out("Specify either frequency or transition ");
         exit(1);
       }
       if(((*img)[i].nchan==0 && (*img)[i].bandwidth<0) || ((*img)[i].bandwidth<0 && (*img)[i].velres<0)){
-        if(!silent) bail_out("Error: Image keywords are not set properly");
+        if(!silent) bail_out("Image keywords are not set properly");
         exit(1);
       }
       (*img)[i].doline=1;
@@ -119,33 +175,32 @@ The cutoff will be the value of abs(x) for which the error in the exact expressi
 
   /* Allocate moldata array */
   (*m)=malloc(sizeof(molData)*par->nSpecies);
-  for( i=0; i<par->nSpecies; i++ )
-    {
-      (*m)[i].ntrans = NULL;
-      (*m)[i].lal = NULL;
-      (*m)[i].lau = NULL;
-      (*m)[i].lcl = NULL;
-      (*m)[i].lcu = NULL;
-      (*m)[i].aeinst = NULL;
-      (*m)[i].freq = NULL;
-      (*m)[i].beinstu = NULL;
-      (*m)[i].beinstl = NULL;
-      (*m)[i].down = NULL;
-      (*m)[i].ntemp = NULL;
-      (*m)[i].eterm = NULL;
-      (*m)[i].gstat = NULL;
-      (*m)[i].cmb = NULL;
-      (*m)[i].local_cmb = NULL;
-    }
+  for( i=0; i<par->nSpecies; i++ ){
+    (*m)[i].ntrans = NULL;
+    (*m)[i].lal = NULL;
+    (*m)[i].lau = NULL;
+    (*m)[i].lcl = NULL;
+    (*m)[i].lcu = NULL;
+    (*m)[i].aeinst = NULL;
+    (*m)[i].freq = NULL;
+    (*m)[i].beinstu = NULL;
+    (*m)[i].beinstl = NULL;
+    (*m)[i].down = NULL;
+    (*m)[i].eterm = NULL;
+    (*m)[i].gstat = NULL;
+    (*m)[i].cmb = NULL;
+    (*m)[i].local_cmb = NULL;
+  }
 }
 
+
 void
-freeMoldata( inputPars *par, molData* mol )
+freeMoldata(const int nSpecies, molData *mol)
 {
   int i;
   if( mol!= 0 )
     {
-      for( i=0; i<par->nSpecies; i++ )
+      for( i=0; i<nSpecies; i++ )
         {
           if( mol[i].ntrans != NULL )
             {
@@ -214,7 +269,7 @@ freeMoldata( inputPars *par, molData* mol )
 }
 
 void
-freeGridPointData(inputPars *par, gridPointData *mol){
+freeGridPointData(configInfo *par, gridPointData *mol){
   int i;
   if (mol!= 0){
     for (i=0;i<par->nSpecies;i++){
@@ -245,7 +300,7 @@ invSqrt(float x){
   return x;
 }
 
-void checkGridDensities(inputPars *par, struct grid *g){
+void checkGridDensities(configInfo *par, struct grid *g){
   int i;
   static _Bool warningAlreadyIssued=0;
   char errStr[80];
@@ -265,7 +320,7 @@ void checkGridDensities(inputPars *par, struct grid *g){
 }
 
 void
-continuumSetup(int im, image *img, molData *m, inputPars *par, struct grid *g){
+continuumSetup(int im, image *img, molData *m, configInfo *par, struct grid *g){
   int id;
   img[im].trans=0;
   m[0].nline=1;
@@ -306,7 +361,7 @@ lineCount(int n,molData *m,int **counta,int **countb,int *nlinetot){
 }
 
 void
-lineBlend(molData *m, inputPars *par, blend **matrix){
+lineBlend(molData *m, configInfo *par, blend **matrix){
   int iline, jline, nlinetot=0,c;
   int *counta,*countb;
 
@@ -346,7 +401,7 @@ lineBlend(molData *m, inputPars *par, blend **matrix){
 }
 
 void
-levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){
+levelPops(molData *m, configInfo *par, struct grid *g, int *popsdone){
   int id,conv=0,iter,ilev,prog=0,ispec,c=0,n,i,threadI,nVerticesDone;
   double percent=0.,*median,result1=0,result2=0,snr,delta_pop;
   blend *matrix;
@@ -392,7 +447,7 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){
   /* Check for blended lines */
   lineBlend(m,par,&matrix);
 
-  if(par->lte_only || par->init_lte) LTE(par,g,m);
+  if(par->lte_only!=0) LTE(par,g,m);
 
   for(id=0;id<par->pIntensity;id++){
     stat[id].pop=malloc(sizeof(double)*m[0].nlev*5);
diff --git a/src/grid.c b/src/grid.c
index 3eb9fd4..61b916d 100644
--- a/src/grid.c
+++ b/src/grid.c
@@ -11,7 +11,7 @@
 
 
 void
-gridAlloc(inputPars *par, struct grid **g){
+gridAlloc(configInfo *par, struct grid **g){
   int i;
   double temp[99];
 
@@ -47,7 +47,7 @@ gridAlloc(inputPars *par, struct grid **g){
 }
 
 void
-freePopulation(const inputPars *par, const molData* m, struct populations* pop ) {
+freePopulation(const configInfo *par, const molData* m, struct populations* pop ) {
   if( pop !=NULL )
     {
       int j,k;
@@ -74,7 +74,7 @@ freePopulation(const inputPars *par, const molData* m, struct populations* pop )
     }
 }
 void
-freeGrid(const inputPars *par, const molData* m ,struct grid* g){
+freeGrid(const configInfo *par, const molData* m ,struct grid* g){
   int i;
   if( g != NULL )
     {
@@ -137,7 +137,7 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){
 }
 
 void
-qhull(inputPars *par, struct grid *g){
+qhull(configInfo *par, struct grid *g){
   int i,j,k,id;
   char flags[255];
   boolT ismalloc = False;
@@ -210,7 +210,7 @@ qhull(inputPars *par, struct grid *g){
 }
 
 void
-distCalc(inputPars *par, struct grid *g){
+distCalc(configInfo *par, struct grid *g){
   int i,k,l;
 
   for(i=0;i<par->ncell;i++){
@@ -237,7 +237,7 @@ distCalc(inputPars *par, struct grid *g){
 
 
 void
-write_VTK_unstructured_Points(inputPars *par, struct grid *g){
+write_VTK_unstructured_Points(configInfo *par, struct grid *g){
   FILE *fp;
   double length;
   int i,j,l=0;
@@ -324,12 +324,12 @@ write_VTK_unstructured_Points(inputPars *par, struct grid *g){
 }
 
 void
-dumpGrid(inputPars *par, struct grid *g){
+dumpGrid(configInfo *par, struct grid *g){
   if(par->gridfile) write_VTK_unstructured_Points(par, g);
 }
 
 void
-getArea(inputPars *par, struct grid *g, const gsl_rng *ran){
+getArea(configInfo *par, struct grid *g, const gsl_rng *ran){
   int i,j,k,b;//=-1;
   double *angle,best;
   /*	double wsum; */
@@ -392,7 +392,7 @@ getArea(inputPars *par, struct grid *g, const gsl_rng *ran){
 
 
 void
-getMass(inputPars *par, struct grid *g, const gsl_rng *ran){
+getMass(configInfo *par, struct grid *g, const gsl_rng *ran){
   double mass=0.,dist;
   double vol=0.,dp,dpbest,*farea,suma;
   int i,k,j,best=-1;
@@ -501,7 +501,7 @@ getMass(inputPars *par, struct grid *g, const gsl_rng *ran){
 
 
 void
-buildGrid(inputPars *par, struct grid *g){
+buildGrid(configInfo *par, struct grid *g){
   double lograd;		/* The logarithm of the model radius		*/
   double logmin;	    /* Logarithm of par->minScale				*/
   double r,theta,phi,sinPhi,x,y,z,semiradius;	/* Coordinates								*/
diff --git a/src/inpars.h b/src/inpars.h
new file mode 100644
index 0000000..1aec912
--- /dev/null
+++ b/src/inpars.h
@@ -0,0 +1,18 @@
+#ifndef INPARS_H
+#define INPARS_H
+
+/* input parameters */
+typedef struct {
+  double radius,minScale,tcmb;
+  int sinkPoints,pIntensity,blend;
+  char *outputfile, *binoutputfile;
+//  char *inputfile; unused at present.
+  char *gridfile;
+  char *pregrid;
+  char *restart;
+  char *dust;
+  int sampling,lte_only,init_lte,antialias,polarization,nThreads;
+  char **moldatfile;
+} inputPars;
+
+#endif /* INPARS_H */
diff --git a/src/lime.h b/src/lime.h
index d92e40f..08726e2 100644
--- a/src/lime.h
+++ b/src/lime.h
@@ -10,6 +10,8 @@
 #ifndef LIME_H
 #define LIME_H
 
+#include "inpars.h"
+
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
@@ -80,6 +82,7 @@
 #define fixset                  1e-6
 #define blendmask               1.e4
 #define MAX_NSPECIES            100
+#define MAX_NIMAGES             100
 #define N_RAN_PER_SEGMENT       3
 #define FAST_EXP_MAX_TAYLOR     3
 #define FAST_EXP_NUM_BITS       8
@@ -96,18 +99,19 @@
 #define CP_Hplus		7
 
 
-/* input parameters */
 typedef struct {
-  double radius,radiusSqu,minScale,minScaleSqu,tcmb,taylorCutoff;
-  int ncell,sinkPoints,pIntensity,nImages,nSpecies,blend;
-  char *outputfile, *binoutputfile, *inputfile;
+  double radius,minScale,tcmb;
+  double radiusSqu,minScaleSqu,taylorCutoff;
+  int sinkPoints,pIntensity,blend;
+  int ncell,nImages,nSpecies,collPart,doPregrid;
+  char *outputfile, *binoutputfile;
   char *gridfile;
   char *pregrid;
   char *restart;
   char *dust;
-  int sampling,collPart,lte_only,init_lte,antialias,polarization,doPregrid,nThreads;
+  int sampling,lte_only,init_lte,antialias,polarization,nThreads;
   char **moldatfile;
-} inputPars;
+} configInfo;
 
 /* Molecular data: shared attributes */
 typedef struct {
@@ -196,7 +200,7 @@ typedef struct {double x,y, *intensity, *tau;} rayData;
 /* Some global variables */
 int silent;
 
-/* Some functions */
+/* User-specifiable functions */
 void density(double,double,double,double *);
 void temperature(double,double,double,double *);
 void abundance(double,double,double,double *);
@@ -206,73 +210,73 @@ void magfield(double,double,double,double *);
 void gasIIdust(double,double,double,double *);
 
 /* More functions */
-void    run(inputPars *, image *);
+void    run(inputPars, image *);
 
-void   	binpopsout(inputPars *, struct grid *, molData *);
-void   	buildGrid(inputPars *, struct grid *);
+void   	binpopsout(configInfo *, struct grid *, molData *);
+void   	buildGrid(configInfo *, struct grid *);
 void	calcFastExpRange(const int, const int, int*, int*, int*);
-void    calcSourceFn(double, const inputPars*, double*, double*);
+void    calcSourceFn(double, const configInfo*, double*, double*);
 void	calcTableEntries(const int, const int);
-void	checkGridDensities(inputPars*, struct grid*);
-void	continuumSetup(int, image*, molData*, inputPars*, struct grid*);
-void	distCalc(inputPars*, struct grid*);
+void	checkGridDensities(configInfo*, struct grid*);
+void	continuumSetup(int, image*, molData*, configInfo*, struct grid*);
+void	distCalc(configInfo*, struct grid*);
 int	factorial(const int);
 double	FastExp(const float);
 void	fit_d1fi(double, double, double*);
 void    fit_fi(double, double, double*);
 void    fit_rr(double, double, double*);
-void   	input(inputPars *, image *);
-float  	invSqrt(float);
-void    freeMoldata(inputPars *, molData* m );
-void   	freeGrid(const inputPars * par, const molData* m, struct grid * g);
-void   	freePopulation(const inputPars * par, const molData* m, struct populations * pop);
+void    freeMoldata(const int, molData*);
+void   	freeGrid(const configInfo * par, const molData* m, struct grid * g);
+void   	freePopulation(const configInfo * par, const molData* m, struct populations * pop);
 double 	gaussline(double, double);
-void    getArea(inputPars *, struct grid *, const gsl_rng *);
+void    getArea(configInfo *, struct grid *, const gsl_rng *);
 void	getclosest(double, double, double, long *, long *, double *, double *, double *);
-void    getjbar(int, molData*, struct grid*, inputPars*, gridPointData*, double*);
-void    getMass(inputPars *, struct grid *, const gsl_rng *);
+void    getjbar(int, molData*, struct grid*, configInfo*, gridPointData*, double*);
+void    getMass(configInfo *, struct grid *, const gsl_rng *);
 void   	getmatrix(int, gsl_matrix *, molData *, struct grid *, int, gridPointData *);
-void	getVelosplines(inputPars *, struct grid *);
-void	getVelosplines_lin(inputPars *, struct grid *);
-void	gridAlloc(inputPars *, struct grid **);
+void	getVelosplines(configInfo *, struct grid *);
+void	getVelosplines_lin(configInfo *, struct grid *);
+void	gridAlloc(configInfo *, struct grid **);
 void   	input(inputPars *, image *);
 float  	invSqrt(float);
-void   	kappa(molData *, struct grid *, inputPars *,int);
-void	levelPops(molData *, inputPars *, struct grid *, int *);
+void   	kappa(molData *, struct grid *, configInfo *,int);
+void	levelPops(molData *, configInfo *, struct grid *, int *);
 void	line_plane_intersect(struct grid *, double *, int , int *, double *, double *, double);
-void	lineBlend(molData *, inputPars *, blend **);
+void	lineBlend(molData *, configInfo *, blend **);
 void    lineCount(int,molData *,int **, int **, int *);
-void	LTE(inputPars *, struct grid *, molData *);
-void	lteOnePoint(inputPars*, molData*, const int, const double, double*);
-void   	molinit(molData *, inputPars *, struct grid *,int);
-void    openSocket(inputPars *par, int);
-void	parseInput(inputPars *, image **, molData **);
-void  	photon(int, struct grid*, molData*, int, const gsl_rng*, inputPars*, blend*, gridPointData*, double*);
+void	LTE(configInfo *, struct grid *, molData *);
+void	lteOnePoint(configInfo*, molData*, const int, const double, double*);
+void   	molinit(molData *, configInfo *, struct grid *,int);
+void	openSocket(char *);
+void	parseInput(inputPars, configInfo*, image**, molData**);
+void  	photon(int, struct grid*, molData*, int, const gsl_rng*, configInfo*, blend*, gridPointData*, double*);
 double 	planckfunc(int, double, molData *, int);
-int     pointEvaluation(inputPars*, double, double, double, double);
-void   	popsin(inputPars *, struct grid **, molData **, int *);
-void   	popsout(inputPars *, struct grid *, molData *);
-void	predefinedGrid(inputPars *, struct grid *);
-void	qhull(inputPars *, struct grid *);
+int     pointEvaluation(configInfo*, double, double, double, double);
+void   	popsin(configInfo *, struct grid **, molData **, int *);
+void   	popsout(configInfo *, struct grid *, molData *);
+void	predefinedGrid(configInfo *, struct grid *);
+void	qhull(configInfo *, struct grid *);
 double 	ratranInput(char *, char *, double, double, double);
-void   	raytrace(int, inputPars *, struct grid *, molData *, image *);
-void	report(int, inputPars *, struct grid *);
-void	smooth(inputPars *, struct grid *);
+void   	raytrace(int, configInfo *, struct grid *, molData *, image *);
+void	readUserInput(inputPars *, image **, int *, int *);
+void	report(int, configInfo *, struct grid *);
+void	setUpConfig(configInfo *, image **, molData **);
+void	smooth(configInfo *, struct grid *);
 int     sortangles(double *, int, struct grid *, const gsl_rng *);
 void	sourceFunc(double*, double*, double, molData*, double, struct grid*, int, int, int, int);
 void    sourceFunc_cont(double*, double*, struct grid*, int, int, int);
 void    sourceFunc_line(double*, double*, molData*, double, struct grid*, int, int, int);
 void    sourceFunc_pol(double*, double*, double, molData*, double, struct grid*, int, int, int, double);
-void   	stateq(int, struct grid*, molData*, int, inputPars*, gridPointData*, double*, _Bool*);
+void   	stateq(int, struct grid*, molData*, int, configInfo*, gridPointData*, double*, _Bool*);
 void	statistics(int, molData *, struct grid *, int *, double *, double *, int *);
 void    stokesangles(double, double, double, double, double *);
 double	taylor(const int, const float);
-void    traceray(rayData, int, int, inputPars *, struct grid *, molData *, image *, int, int *, int *, double);
+void    traceray(rayData, int, int, configInfo *, struct grid *, molData *, image *, int, int *, int *, double);
 void   	velocityspline(struct grid *, int, int, double, double, double*);
 void   	velocityspline2(double *, double *, double, double, double, double*);
 double 	veloproject(double *, double *);
-void	writefits(int, inputPars *, molData *, image *);
-void    write_VTK_unstructured_Points(inputPars *, struct grid *);
+void	writefits(int, configInfo *, molData *, image *);
+void    write_VTK_unstructured_Points(configInfo *, struct grid *);
 
 
 /* Curses functions */
@@ -292,4 +296,5 @@ void    collpartmesg(char *, int);
 void    collpartmesg2(char *, int);
 void    collpartmesg3(int, int);
 
-#endif  /* LIME_H */
+#endif /* LIME_H */
+
diff --git a/src/main.c b/src/main.c
index 971f143..f88d84a 100644
--- a/src/main.c
+++ b/src/main.c
@@ -12,8 +12,8 @@
 
 
 /* Forward declaration of functions only used in this file */
-void initParImg(inputPars *par, image **img);
-void freeParImg(inputPars *par, image *img);
+int initParImg(inputPars *par, image **img);
+void freeParImg(const int nImages, inputPars *par, image *img);
 int main ();
 
 
@@ -29,21 +29,26 @@ double EXP_TABLE_2D[1][1]; // nominal definitions so the fastexp.c module will c
 double EXP_TABLE_3D[1][1][1];
 #endif
 
-void
+int
 initParImg(inputPars *par, image **img)
 {
   /* Initialize par with default values, allocate space for the
      output fits images, initialize the images with default values,
-     and finally call the input routine from model.c to set both the
+     and finally call the input() routine from model.c to set both the
      par and image values.
   */
 
-  int i, id;
+  int i,id,nImages;
   FILE *fp;
 
-  /* Set default values */
+  /* Set 'impossible' default values for mandatory parameters */
+  par->radius    =-1;
+  par->minScale  =-1;
+  par->pIntensity=-1;
+  par->sinkPoints=-1;
+
+  /* Set default values for optional parameters */
   par->dust  	    = NULL;
-  par->inputfile    = NULL;
   par->outputfile   = NULL;
   par->binoutputfile= NULL;
   par->gridfile     = NULL;
@@ -57,65 +62,53 @@ initParImg(inputPars *par, image **img)
   par->blend=0;
   par->antialias=1;
   par->polarization=0;
-  par->pIntensity=0;
-  par->sinkPoints=0;
-  par->doPregrid=0;
-
   par->nThreads = NTHREADS;
 
-  /* Allocate space for output fits images */
-  (*img)=malloc(sizeof(image)*MAX_NSPECIES);
-  par->moldatfile=malloc(sizeof(char *) * MAX_NSPECIES);
+  /* Allocate initial space for molecular data file names */
+  par->moldatfile=malloc(sizeof(char *)*MAX_NSPECIES);
   for(id=0;id<MAX_NSPECIES;id++){
-    (*img)[id].filename=NULL;
     par->moldatfile[id]=NULL;
   }
+
+  /* Allocate initial space for output fits images */
+  (*img)=malloc(sizeof(image)*MAX_NIMAGES);
+  for(i=0;i<MAX_NIMAGES;i++){
+    (*img)[i].filename=NULL;
+  }
+
+  /* First-pass reading of the user-set parameters */
   input(par, *img);
-  id=-1;
-  while((*img)[++id].filename!=NULL);
-  par->nImages=id;
-  if(par->nImages==0) {
-    if(!silent) bail_out("Error: no images defined");
+
+  /* Check that the mandatory parameters now have 'sensible' settings (i.e., that they have been set at all). Raise an exception if not. */
+  if (par->radius<=0){
+    if(!silent) bail_out("You must define the radius parameter.");
+    exit(1);
+  }
+  if (par->minScale<=0){
+    if(!silent) bail_out("You must define the minScale parameter.");
+    exit(1);
+  }
+  if (par->pIntensity<=0){
+    if(!silent) bail_out("You must define the pIntensity parameter.");
+    exit(1);
+  }
+  if (par->sinkPoints<=0){
+    if(!silent) bail_out("You must define the sinkPoints parameter.");
     exit(1);
   }
 
-  *img=realloc(*img, sizeof(image)*par->nImages);
-
-  id=-1;
-  while(par->moldatfile[++id]!=NULL);
-  par->nSpecies=id;
-  if( par->nSpecies <= 0 )
-    {
-      par->nSpecies = 1;
-      free(par->moldatfile);
-      par->moldatfile = NULL;
-    }
-  else
-    {
-      par->moldatfile=realloc(par->moldatfile, sizeof(char *)*par->nSpecies);
-      /* Check if files exists */
-      for(id=0;id<par->nSpecies;id++){
-        if((fp=fopen(par->moldatfile[id], "r"))==NULL) {
-          openSocket(par, id);
-        }
-        else {
-          fclose(fp);
-        }
-      }
-    }
+  /* If the user has provided a list of image filenames, the corresponding elements of (*img).filename will be non-NULL. Thus we can deduce the number of images from the number of non-NULL elements. */
+  nImages=0;
+  while((*img)[nImages].filename!=NULL && nImages<MAX_NIMAGES)
+    nImages++;
 
-  if(par->dust != NULL){
-    if((fp=fopen(par->dust, "r"))==NULL){
-      if(!silent) bail_out("Error opening dust opacity data file!");
-      exit(1);
-    }
-    else  {
-      fclose(fp);
-    }
+  if(nImages==0) {
+    if(!silent) bail_out("No images defined (or you haven't set the 1st filename).");
+    exit(1);
   }
 
-  /* Set defaults and read inputPars and img[] */
-  for(i=0;i<par->nImages;i++) {
+  /* Set img defaults. */
+  for(i=0;i<nImages;i++) {
     (*img)[i].source_vel=0.0;
     (*img)[i].phi=0.0;
     (*img)[i].nchan=0;
@@ -124,18 +117,22 @@ initParImg(inputPars *par, image **img)
     (*img)[i].freq=-1.;
     (*img)[i].bandwidth=-1.;
   }
+
+  /* Second-pass reading of the user-set parameters (this time just to read the par->moldatfile and img stuff). */
   input(par,*img);
+
+  return nImages;
 }
 
 
 void
-freeParImg(inputPars *par, image *img)
+freeParImg(const int nImages, inputPars *par, image *img)
 {
   /* Release memory allocated for the output fits images
      and for par->moldatfile
   */
   int i, id;
-  for(i=0;i<par->nImages;i++){
+  for(i=0;i<nImages;i++){
     for(id=0;id<(img[i].pxls*img[i].pxls);id++){
       free( img[i].pixel[id].intense );
       free( img[i].pixel[id].tau );
@@ -154,18 +151,19 @@ freeParImg(inputPars *par, image *img)
 
 
 void
-run(inputPars *par, image *img)
+run(inputPars inpars, image *img)
 {
-  /* Run LIME with par and the output fits files specified.
+  /* Run LIME with inpars and the output fits files specified.
 
      This routine may be used as an interface to LIME from external
-     programs. In this case, par and img must be specified by the
+     programs. In this case, inpars and img must be specified by the
      external program.
   */
   int i,nLineImages;
   int initime=time(0);
   int popsdone=0;
   molData*     m = NULL;
+  configInfo  par;
   struct grid* g = NULL;
 
   /*Set locale to avoid trouble when reading files*/
@@ -178,52 +176,53 @@ run(inputPars *par, image *img)
   calcTableEntries(FAST_EXP_MAX_TAYLOR, FAST_EXP_NUM_BITS);
 #endif
 
-  parseInput(par, &img, &m);
+  parseInput(inpars, &par, &img, &m);
 
-  if(par->doPregrid)
+  if(par.doPregrid)
     {
-      gridAlloc(par,&g);
-      predefinedGrid(par,g);
+      gridAlloc(&par,&g);
+      predefinedGrid(&par,g);
     }
-  else if(par->restart)
+  else if(par.restart)
     {
-      popsin(par,&g,&m,&popsdone);
+      popsin(&par,&g,&m,&popsdone);
     }
   else
     {
-      gridAlloc(par,&g);
-      buildGrid(par,g);
+      gridAlloc(&par,&g);
+      buildGrid(&par,g);
     }
 
   /* Make all the continuum images, and count the non-continuum images at the same time:
   */
   nLineImages = 0;
-  for(i=0;i<par->nImages;i++){
+  for(i=0;i<par.nImages;i++){
     if(img[i].doline)
       nLineImages++;
     else{
-      continuumSetup(i,img,m,par,g);
-      raytrace(i,par,g,m,img);
-      writefits(i,par,m,img);
+      continuumSetup(i,img,m,&par,g);
+      raytrace(i,&par,g,m,img);
+      writefits(i,&par,m,img);
     }
   }
 
   if(nLineImages>0 && !popsdone)
-    levelPops(m,par,g,&popsdone);
+    levelPops(m,&par,g,&popsdone);
 
   /* Now make the line images.
   */
-  for(i=0;i<par->nImages;i++){
+  for(i=0;i<par.nImages;i++){
     if(img[i].doline){
-      raytrace(i,par,g,m,img);
-      writefits(i,par,m,img);
+      raytrace(i,&par,g,m,img);
+      writefits(i,&par,m,img);
     }
   }
   
   if(!silent) goodnight(initime,img[0].filename);
 
-  freeGrid( par, m, g);
-  freeMoldata(par, m);
+  freeGrid( &par, m, g);
+  freeMoldata(par.nSpecies, m);
+  free(par.moldatfile);
 }
 
 
@@ -232,14 +231,15 @@ int main () {
 
   inputPars par;
   image	*img = NULL;
+  int nImages;
 
   silent = 0;
 
-  initParImg(&par, &img);
+  nImages = initParImg(&par, &img);
 
-  run(&par, img);
+  run(par, img);
 
-  freeParImg(&par, img);
+  freeParImg(nImages, &par, img);
 
   return 0;
 }
diff --git a/src/molinit.c b/src/molinit.c
index a9044c9..ec98fb2 100644
--- a/src/molinit.c
+++ b/src/molinit.c
@@ -10,7 +10,7 @@
 #include "lime.h"
 
 void
-kappa(molData *m, struct grid *g, inputPars *par, int s){
+kappa(molData *m, struct grid *g, configInfo *par, int s){
   FILE *fp;
   char string[80];
   int i=0,k,j,iline,id;
@@ -102,7 +102,7 @@ planckfunc(int iline, double temp, molData *m,int s){
 }
 
 void
-molinit(molData *m, inputPars *par, struct grid *g,int i){
+molinit(molData *m, configInfo *par, struct grid *g,int i){
   int id, ilev, iline, itrans, ispec, itemp, *ntemp, tnint=-1, idummy, ipart, *collPartIDs,flag=0;
   char *collpartnames[] = {"H2","p-H2","o-H2","electrons","H","He","H+"}; /* definition from LAMDA */
   double fac, uprate, downrate=0, dummy, amass;
diff --git a/src/photon.c b/src/photon.c
index b24783b..a112ce3 100644
--- a/src/photon.c
+++ b/src/photon.c
@@ -131,14 +131,14 @@ double gaussline(double v, double oneOnSigma){
 }
 
 
-void calcSourceFn(double dTau, const inputPars *par, double *remnantSnu, double *expDTau){
+void calcSourceFn(double dTau, const configInfo *par, double *remnantSnu, double *expDTau){
   /*
   The source function S is defined as j_nu/alpha, which is clearly not
   defined for alpha==0. However S is used in the algorithm only in the
   term (1-exp[-alpha*ds])*S, which is defined for all values of alpha.
   The present function calculates this term and returns it in the
   argument remnantSnu. For values of abs(alpha*ds) less than a pre-
-  calculated cutoff supplied in inputPars, a Taylor approximation is
+  calculated cutoff supplied in configInfo, a Taylor approximation is
   used.
 
   Note that the same cutoff condition holds for replacement of
@@ -165,7 +165,9 @@ void calcSourceFn(double dTau, const inputPars *par, double *remnantSnu, double
 
 
 void
-photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPars *par,blend *matrix, gridPointData *mp, double *halfFirstDs){
+photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran\
+  , configInfo *par, blend *matrix, gridPointData *mp, double *halfFirstDs){
+
   int iphot,iline,jline,here,there,firststep,neighI,np_per_line,ip_at_line,l;
   int *counta, *countb,nlinetot;
   double deltav,segment,vblend,dtau,expDTau,jnu,alpha,ds,vfac[par->nSpecies],pt_theta,pt_z,semiradius;
@@ -298,7 +300,7 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar
 }
 
 void
-getjbar(int posn, molData *m, struct grid *g, inputPars *par, gridPointData *mp, double *halfFirstDs){
+getjbar(int posn, molData *m, struct grid *g, configInfo *par, gridPointData *mp, double *halfFirstDs){
   int iline,iphot;
   double tau, expTau, remnantSnu, vsum=0., jnu, alpha;
   int *counta, *countb,nlinetot;
diff --git a/src/popsin.c b/src/popsin.c
index 9c83047..4cb9b3a 100644
--- a/src/popsin.c
+++ b/src/popsin.c
@@ -11,7 +11,7 @@
 
 
 void
-popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){
+popsin(configInfo *par, struct grid **g, molData **m, int *popsdone){
   FILE *fp;
   int i,j,k;
   double dummy;
diff --git a/src/popsout.c b/src/popsout.c
index 93ff3e3..7152564 100644
--- a/src/popsout.c
+++ b/src/popsout.c
@@ -10,7 +10,7 @@
 #include "lime.h"
 
 void
-popsout(inputPars *par, struct grid *g, molData *m){
+popsout(configInfo *par, struct grid *g, molData *m){
   FILE *fp;
   int j,k,l;
   double dens;
@@ -35,7 +35,7 @@ popsout(inputPars *par, struct grid *g, molData *m){
 
 
 void
-binpopsout(inputPars *par, struct grid *g, molData *m){
+binpopsout(configInfo *par, struct grid *g, molData *m){
   FILE *fp;
   int i,j;
   
diff --git a/src/predefgrid.c b/src/predefgrid.c
index 1d9543f..82f6452 100644
--- a/src/predefgrid.c
+++ b/src/predefgrid.c
@@ -10,7 +10,7 @@
 #include "lime.h"
 
 void
-predefinedGrid(inputPars *par, struct grid *g){
+predefinedGrid(configInfo *par, struct grid *g){
   FILE *fp;
   int i;
   double x,y,z,scale;
diff --git a/src/raytrace.c b/src/raytrace.c
index f367670..a6f0a52 100644
--- a/src/raytrace.c
+++ b/src/raytrace.c
@@ -68,7 +68,7 @@ This function returns ds as the (always positive-valued) distance between the pr
 
 
 void
-traceray(rayData ray, int tmptrans, int im, inputPars *par, struct grid *g, molData *m, image *img, int nlinetot, int *counta, int *countb, double cutoff){
+traceray(rayData ray, int tmptrans, int im, configInfo *par, struct grid *g, molData *m, image *img, int nlinetot, int *counta, int *countb, double cutoff){
   /*
 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.
 
@@ -188,7 +188,7 @@ Note that the algorithm employed here is similar to that employed in the functio
 
 
 void
-raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){
+raytrace(int im, configInfo *par, struct grid *g, molData *m, image *img){
   int *counta, *countb,nlinetot,aa;
   int ichan,px,iline,tmptrans,i,threadI,nRaysDone;
   double size,minfreq,absDeltaFreq,totalNumPixelsMinus1=(double)(img[im].pxls*img[im].pxls-1);
@@ -308,7 +308,7 @@ While this is off however, gsl_* calls will not exit if they encounter a problem
 
 
 void
-raytrace_1_4(int im, inputPars *par, struct grid *g, molData *m, image *img){
+raytrace_1_4(int im, configInfo *par, struct grid *g, molData *m, image *img){
   /*
 This is an alternative raytracing algorithm which was implemented by
 C Brinch in version 1.4 (the original parallelized version) of LIME.
diff --git a/src/report.c b/src/report.c
index bd8c6e0..46211ac 100644
--- a/src/report.c
+++ b/src/report.c
@@ -13,7 +13,7 @@
 
 
 void
-report(int i, inputPars *par, struct grid *g){
+report(int i, configInfo *par, struct grid *g){
 	FILE *fp;
 	int j,k,p,q;
 	const int bins = 100;
diff --git a/src/smooth.c b/src/smooth.c
index ebbf533..1c31d70 100644
--- a/src/smooth.c
+++ b/src/smooth.c
@@ -12,7 +12,7 @@
 
 /* Based on Lloyds Algorithm (Lloyd, S. IEEE, 1982) */	
 void
-smooth(inputPars *par, struct grid *g){
+smooth(configInfo *par, struct grid *g){
   double mindist;	/* Distance to closest neighbor				*/
   int k=0,j,i;		/* counters									*/
   int sg;		/* counter for smoothing the grid			*/
diff --git a/src/stateq.c b/src/stateq.c
index c00ec37..1d19f71 100644
--- a/src/stateq.c
+++ b/src/stateq.c
@@ -14,7 +14,7 @@
 
 
 void
-stateq(int id, struct grid *g, molData *m, int ispec, inputPars *par\
+stateq(int id, struct grid *g, molData *m, int ispec, configInfo *par\
   , gridPointData *mp, double *halfFirstDs, _Bool *luWarningGiven){
 
   int t,s,iter,status;
diff --git a/src/tcpsocket.c b/src/tcpsocket.c
index b1a151e..87fa138 100644
--- a/src/tcpsocket.c
+++ b/src/tcpsocket.c
@@ -16,7 +16,7 @@
 int close(int);
 
 void
-openSocket(inputPars *par, int id){
+openSocket(char *moldatfile){
   struct sockaddr_in *remote;
   int tmpres;
   int sock;
@@ -30,20 +30,20 @@ openSocket(inputPars *par, int id){
   FILE *fp;
 
   // Check if moldatfile contains .dat
-  if(strstr(par->moldatfile[id], ".dat") == NULL){
-    size_t length = strlen(par->moldatfile[id]);
+  if(strstr(moldatfile, ".dat") == NULL){
+    size_t length = strlen(moldatfile);
     s = (char*)malloc(sizeof(char) * (length + 5));
-    strcpy(s,par->moldatfile[id]);
+    strcpy(s,moldatfile);
     strcat(s, ".dat");
     s[length+5]='\0';
-    par->moldatfile[id]=s;
+    moldatfile=s;
   }
 
   size_t length1 = strlen(page);
-  size_t length2 = strlen(par->moldatfile[id]);
+  size_t length2 = strlen(moldatfile);
   t = (char*)malloc(sizeof(char) * (length1 + length2 + 1));
   strcpy(t,page);
-  strcat(t, par->moldatfile[id]);
+  strcat(t, moldatfile);
   page=t;
 
 
@@ -88,7 +88,7 @@ openSocket(inputPars *par, int id){
   }
 
   memset(buf, 0, sizeof(buf));
-  if((fp=fopen(par->moldatfile[id], "w"))==NULL) {
+  if((fp=fopen(moldatfile, "w"))==NULL) {
     if(!silent) bail_out("Failed to write moldata!");
     exit(1);
   }
diff --git a/src/velospline.c b/src/velospline.c
index 3e063bd..0aea5ec 100644
--- a/src/velospline.c
+++ b/src/velospline.c
@@ -10,7 +10,7 @@
 #include "lime.h"
 
 void
-getVelosplines(inputPars *par, struct grid *g){
+getVelosplines(configInfo *par, struct grid *g){
   int i,k,j,l,s;
   double v[5], vel[3], x[3], d;
   gsl_matrix *matrix = gsl_matrix_alloc(5,5);
@@ -76,7 +76,7 @@ getVelosplines(inputPars *par, struct grid *g){
 
 
 void
-getVelosplines_lin(inputPars *par, struct grid *g){
+getVelosplines_lin(configInfo *par, struct grid *g){
   int i,k,j;
   double v[2];
   
diff --git a/src/weights.c b/src/weights.c
index f4de61a..c0be334 100644
--- a/src/weights.c
+++ b/src/weights.c
@@ -10,7 +10,7 @@
 #include "lime.h"
 
 int
-pointEvaluation(inputPars *par,double ran, double x, double y, double z){
+pointEvaluation(configInfo *par,double ran, double x, double y, double z){
   double weight1, weight2, val[99],normalizer=0.0,totalDensity=0.0;
   int i;
 
diff --git a/src/writefits.c b/src/writefits.c
index 2ba9971..0c77fe8 100644
--- a/src/writefits.c
+++ b/src/writefits.c
@@ -10,7 +10,7 @@
 #include "lime.h"
 
 void 
-writefits(int im, inputPars *par, molData *m, image *img){
+writefits(int im, configInfo *par, molData *m, image *img){
   double bscale,bzero,epoch,lonpole,equinox,restfreq;
   double cdelt1,crpix1,crval1,cdelt2,crpix2,crval2;
   double cdelt3,crpix3,crval3,ru3;