Inpars

.gitignore

@@ -5,6 +5,8 @@
 .DS_Store
 
 # Ignore example files
-example/grid.vtk
+example/*.vtk
 example/*.pop
-example/*.fits
+example/*.fits*
+example/*.dat
+

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)

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;
 

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);

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								*/

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 */