diff --git a/doc/usermanual.rst b/doc/usermanual.rst index e6c9578..4267600 100644 --- a/doc/usermanual.rst +++ b/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 diff --git a/src/aux.c b/src/aux.c index ce7f383..aa53415 100644 --- a/src/aux.c +++ b/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); diff --git a/src/lime.h b/src/lime.h index 2a5353b..e5e933c 100644 --- a/src/lime.h +++ b/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; diff --git a/src/photon.c b/src/photon.c index 4533870..8991e40 100644 --- a/src/photon.c +++ b/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;lnSpecies;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;lnSpecies;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;jlinedoPregrid) 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; diff --git a/src/predefgrid.c b/src/predefgrid.c index 4aed807..57a95b8 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -25,6 +25,13 @@ predefinedGrid(inputPars *par, struct grid *g){ par->ncell=par->pIntensity+par->sinkPoints; for(i=0;ipIntensity;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;incell;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); } diff --git a/src/raytrace.c b/src/raytrace.c index 8c6131e..dc330a6 100644 --- a/src/raytrace.c +++ b/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]);