Modification of predefgrid

doc/usermanual.rst

@@ -409,6 +409,11 @@ If this option is used, LIME will not generate its own grid, but rather
 use the grid defined in this file. The file needs to contain all
 physical properties of the grid points, i.e., density, temperature,
 abundance, velocity etc. There is no default value.
+An experienced user can change what physical properties will be taken 
+from the input file. This can be done modifying the predefinedGrid 
+function in the predefgrid.c file. The physical properties that will not 
+be taken from the input file will be computed with the model 
+functions (see below).
 
 .. code:: c
 

src/aux.c

@@ -39,6 +39,7 @@ parseInput(inputPars *par, image **img, molData **m){
   par->sinkPoints=0;
   par->doPregrid=0;
   par->nThreads=0;
+  par->pregridVel=0;
 
   /* Allocate space for output fits images */
   (*img)=malloc(sizeof(image)*MAX_NSPECIES);

src/lime.h

@@ -87,7 +87,7 @@ typedef struct {
   char *pregrid;
   char *restart;
   char *dust;
-  int sampling,collPart,lte_only,antialias,polarization,doPregrid,nThreads;
+  int sampling,collPart,lte_only,antialias,polarization,doPregrid,nThreads,pregridVel;
   char **moldatfile;
 } inputPars;
 

src/photon.c

@@ -217,7 +217,7 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar
         ds=g[here].ds[dir]/2.;
         halfFirstDs[iphot]=ds;
         for(l=0;l<par->nSpecies;l++){
-          if(!par->doPregrid) velocityspline(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]);
+          if(!par->doPregrid  || !par->pregridVel) velocityspline(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]);
           else velocityspline_lin(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]);
           mp[l].vfac[iphot]=vfac[0];
         }
@@ -228,7 +228,7 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar
       }
 
       for(l=0;l<par->nSpecies;l++){
-        if(!par->doPregrid) velocityspline(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]);
+        if(!par->doPregrid || !par->pregridVel) velocityspline(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]);
         else velocityspline_lin(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]);
       }
 
@@ -259,7 +259,7 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar
           alpha=0.;
           for(jline=0;jline<sizeof(matrix)/sizeof(blend);jline++){
             if(matrix[jline].line1 == jline || matrix[jline].line2 == jline){	
-              if(!par->doPregrid) velocityspline(g,here,dir,g[id].mol[counta[jline]].binv,deltav-matrix[jline].deltav,&vblend);
+              if(!par->doPregrid || !par->pregridVel) velocityspline(g,here,dir,g[id].mol[counta[jline]].binv,deltav-matrix[jline].deltav,&vblend);
               else velocityspline_lin(g,here,dir,g[id].mol[counta[jline]].binv,deltav-matrix[jline].deltav,&vblend);	
               sourceFunc_line(&jnu,&alpha,m,vblend,g,here,counta[jline],countb[jline]);
               dtau=alpha*ds;

src/predefgrid.c

@@ -25,6 +25,13 @@ predefinedGrid(inputPars *par, struct grid *g){
   par->ncell=par->pIntensity+par->sinkPoints;
 
   for(i=0;i<par->pIntensity;i++){
+    g[i].id = -1;
+    g[i].dens[0] = -1.0;
+    g[i].t[0] = -1.0;
+    g[i].dopb = -1.0;
+    g[i].abun[0] = -1.0;
+    g[i].vel[0] = 1.e40; g[i].vel[1] = 1.e40; g[i].vel[2] = 1.e40;
+    //    fscanf(fp,"%lf %lf %lf\n", &g[i].x[0], &g[i].x[1], &g[i].x[2]);
     //    fscanf(fp,"%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", &g[i].id, &g[i].x[0], &g[i].x[1], &g[i].x[2],  &g[i].dens[0], &g[i].t[0], &abun, &g[i].dopb, &g[i].vel[0], &g[i].vel[1], &g[i].vel[2]);
     //    fscanf(fp,"%d %lf %lf %lf %lf %lf %lf %lf\n", &g[i].id, &g[i].x[0], &g[i].x[1], &g[i].x[2],  &g[i].dens[0], &g[i].t[0], &abun, &g[i].dopb);
     int nRead = fscanf(fp,"%d %lf %lf %lf %lf %lf %lf %lf %lf\n", &g[i].id, &g[i].x[0], &g[i].x[1], &g[i].x[2],  &g[i].dens[0], &g[i].t[0], &g[i].vel[0], &g[i].vel[1], &g[i].vel[2]);
@@ -34,19 +41,22 @@ predefinedGrid(inputPars *par, struct grid *g){
         exit(0);
       }
 
-    g[i].dopb=200;
-    g[i].abun[0]=1e-9;
-
+    if(g[i].id==-1) g[i].id = i;
+    if(g[i].dens[0]==-1.0) density(g[i].x[0],g[i].x[1],g[i].x[2],g[i].dens);
+    if(g[i].t[0]==-1.0) temperature(g[i].x[0],g[i].x[1],g[i].x[2],g[i].t);
+    if(g[i].dopb==-1.0) doppler(g[i].x[0],g[i].x[1],g[i].x[2],&g[i].dopb);
+    if(g[i].abun[0]==-1.0) abundance(g[i].x[0],g[i].x[1],g[i].x[2],g[i].abun);
+    if(g[i].vel[0]!=1.e40 && g[i].vel[1]!=1.e40 && g[i].vel[2]!=1.e40) par->pregridVel = 1;
 
     g[i].sink=0;
-	g[i].t[1]=g[i].t[0];
-	g[i].nmol[0]=g[i].abun[0]*g[i].dens[0];
+    g[i].t[1]=g[i].t[0];
+    g[i].nmol[0]=g[i].abun[0]*g[i].dens[0];
 
-	/* This next step needs to be done, even though it looks stupid */
-	g[i].dir=malloc(sizeof(point)*1);
-	g[i].ds =malloc(sizeof(double)*1);
-	g[i].neigh =malloc(sizeof(struct grid *)*1);
-	if(!silent) progressbar((double) i/((double)par->pIntensity-1), 4);	
+    /* This next step needs to be done, even though it looks stupid */
+    g[i].dir=malloc(sizeof(point)*1);
+    g[i].ds =malloc(sizeof(double)*1);
+    g[i].neigh =malloc(sizeof(struct grid *)*1);
+    if(!silent) progressbar((double) i/((double)par->pIntensity-1), 4);	
   }
 
   for(i=par->pIntensity;i<par->ncell;i++){
@@ -73,7 +83,8 @@ predefinedGrid(inputPars *par, struct grid *g){
   distCalc(par,g);
   //  getArea(par,g, ran);
   //  getMass(par,g, ran);
-  getVelosplines_lin(par,g);
+  if(par->pregridVel) getVelosplines_lin(par,g);
+  else getVelosplines(par,g);
   if(par->gridfile) write_VTK_unstructured_Points(par, g);
   gsl_rng_free(ran);
 }

src/raytrace.c

@@ -125,7 +125,7 @@ traceray(rayData ray, int tmptrans, int im, inputPars *par, struct grid *g, molD
               vThisChan=(ichan-(img[im].nchan-1)/2.)*img[im].velres; // Consistent with the WCS definition in writefits().
               deltav = vThisChan - img[im].source_vel + shift;
 
-              if(!par->pregrid) velocityspline2(x,dx,ds,g[posn].mol[counta[iline]].binv,deltav,&vfac);
+              if(!par->pregrid || !par->pregridVel) velocityspline2(x,dx,ds,g[posn].mol[counta[iline]].binv,deltav,&vfac);
               else vfac=gaussline(deltav+veloproject(dx,g[posn].vel),g[posn].mol[counta[iline]].binv);
 
               sourceFunc_line(&jnu,&alpha,m,vfac,g,posn,counta[iline],countb[iline]);