Partial fix to #181 and a bug fix in photon()

example/model.c

@@ -27,6 +27,7 @@ input(inputPars *par, image *img){
   par->antialias                = 4;
   par->sampling                 = 2; // log distr. for radius, directions distr. uniformly on a sphere.
   par->nSolveIters              = 14;
+  par->resetRNG	                = 0;
 
   par->outputfile               = "populations.pop";
   par->binoutputfile            = "restart.pop";

src/aux.c

@@ -50,6 +50,7 @@ parseInput(inputPars inpar, configInfo *par, image **img, molData **m){
   par->gridInFile   = inpar.gridInFile;
   par->nSolveIters  = inpar.nSolveIters;
   par->traceRayAlgorithm = inpar.traceRayAlgorithm;
+  par->resetRNG     = inpar.resetRNG;
 
   par->gridOutFiles = malloc(sizeof(char *)*NUM_GRID_STAGES);
   for(i=0;i<NUM_GRID_STAGES;i++)
@@ -703,6 +704,7 @@ levelPops(molData *md, configInfo *par, struct grid *gp, int *popsdone, double *
   struct blendInfo blends;
   _Bool luWarningGiven=0;
   gsl_error_handler_t *defaultErrorHandler=NULL;
+  int RNG_seeds[par->nThreads];
 
   nlinetot = 0;
   for(ispec=0;ispec<par->nSpecies;ispec++)
@@ -730,7 +732,8 @@ levelPops(molData *md, configInfo *par, struct grid *gp, int *popsdone, double *
 
     for (i=0;i<par->nThreads;i++){
       threadRans[i] = gsl_rng_alloc(ranNumGenType);
-      gsl_rng_set(threadRans[i],(int)(gsl_rng_uniform(ran)*1e6));
+      if (par->resetRNG==1) RNG_seeds[i] = (int)(gsl_rng_uniform(ran)*1e6);
+      else gsl_rng_set(threadRans[i],(int)(gsl_rng_uniform(ran)*1e6));
     }
 
     calcGridCollRates(par,md,gp);
@@ -777,12 +780,14 @@ While this is off however, other gsl_* etc calls will not exit if they encounter
 
       calcGridMolSpecNumDens(par,md,gp);
 
+
       nVerticesDone=0;
       omp_set_dynamic(0);
 #pragma omp parallel private(id,ispec,threadI,nextMolWithBlend) num_threads(par->nThreads)
       {
         threadI = omp_get_thread_num();
 
+        if (par->resetRNG==1) gsl_rng_set(threadRans[threadI],RNG_seeds[threadI]);
         /* Declare and allocate thread-private variables */
         gridPointData *mp;	// Could have declared them earlier
         double *halfFirstDs;	// and included them in private() I guess.
@@ -795,7 +800,7 @@ While this is off however, other gsl_* etc calls will not exit if they encounter
         }
         halfFirstDs = malloc(sizeof(*halfFirstDs)*max_phot);
 
-#pragma omp for
+#pragma omp for schedule(dynamic)
         for(id=0;id<par->pIntensity;id++){
 #pragma omp atomic
           ++nVerticesDone;

src/inpars.h

@@ -25,6 +25,7 @@ typedef struct {
   char *gridInFile,**gridOutFiles;
   int nSolveIters;
   double (*gridDensMaxLoc)[DIM], *gridDensMaxValues;
+  _Bool resetRNG;
 } inputPars;
 
 #endif /* INPARS_H */

src/lime.h

@@ -68,8 +68,7 @@
 #define PI                      3.14159265358979323846	/* pi	*/
 #define SPI                     1.77245385091		/* sqrt(pi)	*/
 #define maxp                    0.15
-#define OtoP                    3.
-#define NITERATIONS             16
+#define NITERATIONS 		16
 #define max_phot                10000		/* don't set this value higher unless you have enough memory. */
 #define ininphot                9
 #define minpop                  1.e-6
@@ -94,9 +93,9 @@
 #define ERF_TABLE_LIMIT		6.0             /* For x>6 erf(x)-1<double precision machine epsilon, so no need to store the values for larger x. */
 #define ERF_TABLE_SIZE		6145
 #define BIN_WIDTH		(ERF_TABLE_LIMIT/(ERF_TABLE_SIZE-1.))
-#define IBIN_WIDTH 		1./BIN_WIDTH
+#define IBIN_WIDTH 		(1./BIN_WIDTH)
 #define N_VEL_SEG_PER_HALF	1
-#define NUM_VEL_COEFFS		1+2*N_VEL_SEG_PER_HALF /* This is the number of velocity samples per edge (not including the grid vertices at each end of the edge). Currently this is elsewhere hard-wired at 3, the macro just being used in the file I/O modules. Note that we want an odd number of velocity samples per edge if we want to have the ability to do 2nd-order interpolation of velocity within Delaunay tetrahedra. */
+#define NUM_VEL_COEFFS		(1+2*N_VEL_SEG_PER_HALF) /* This is the number of velocity samples per edge (not including the grid vertices at each end of the edge). Currently this is elsewhere hard-wired at 3, the macro just being used in the file I/O modules. Note that we want an odd number of velocity samples per edge if we want to have the ability to do 2nd-order interpolation of velocity within Delaunay tetrahedra. */
 
 /* Collision partner ID numbers from LAMDA */
 #define CP_H2			1
@@ -152,7 +151,7 @@ typedef struct {
   int sampling,lte_only,init_lte,antialias,polarization,nThreads,numDims;
   int nLineImages, nContImages;
   char **moldatfile;
-  _Bool writeGridAtStage[NUM_GRID_STAGES];
+  _Bool writeGridAtStage[NUM_GRID_STAGES], resetRNG;
   char *gridInFile,**gridOutFiles;
   int dataFlags,nSolveIters;
   double (*gridDensMaxLoc)[DIM], *gridDensMaxValues;

src/main.c

@@ -85,6 +85,7 @@ initParImg(inputPars *par, image **img)
   par->nThreads = NTHREADS;
   par->nSolveIters=17;
   par->traceRayAlgorithm=0;
+  par->resetRNG=0;
 
   par->gridOutFiles = malloc(sizeof(char *)*NUM_GRID_STAGES);
   for(i=0;i<NUM_GRID_STAGES;i++)

src/photon.c

@@ -335,18 +335,25 @@ photon(int id, struct grid *gp, molData *md, int iter, const gsl_rng *ran\
           /* End of line blending part */
 
 	  /* as said above, out-in split should be done also for blended lines... */
-          dtau=alpha_line_in*ds_in+alpha_line_out*ds_out+(alpha_cont+alpha_blend)*(ds_in+ds_out); 
-          if(dtau < -30) dtau = -30;
 
+	  dtau=(alpha_line_out+alpha_cont+alpha_blend)*ds_out;
+          if(dtau < -30.) dtau = -30.;
           calcSourceFn(dtau, par, &remnantSnu, &expDTau);
-          remnantSnu *= jnu_line_in*ds_in+jnu_line_out*ds_out+(jnu_cont+jnu_blend)*(ds_in+ds_out);
+          remnantSnu *= (jnu_line_out+jnu_cont+jnu_blend)*ds_out;
+          mp[molI].phot[lineI+iphot*md[molI].nline]+=expTau[iline]*remnantSnu;
+	  expTau[iline]*=expDTau;
+
+	  dtau=(alpha_line_in+alpha_cont+alpha_blend)*ds_in;
+          if(dtau < -30.) dtau = -30.;
+          calcSourceFn(dtau, par, &remnantSnu, &expDTau);
+          remnantSnu *= (jnu_line_in+jnu_cont+jnu_blend)*ds_in;
+          mp[molI].phot[lineI+iphot*md[molI].nline]+=expTau[iline]*remnantSnu;
 	  expTau[iline]*=expDTau;
 
-          if(expTau[iline] > exp(30)){
+          if(expTau[iline] > exp(30.)){
             if(!silent) warning("Maser warning: optical depth has dropped below -30");
-            expTau[iline]=exp(30);
+            expTau[iline]=exp(30.);
           }
-          mp[molI].phot[lineI+iphot*md[molI].nline]+=expTau[iline]*remnantSnu;
 
           iline++;
         } /* Next line this molecule. */