From e9b5c24d8bc08cf066e286af3b5438d2daf23617 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 11 Sep 2015 13:32:53 +0200 Subject: [PATCH 01/43] New member struct spec The new member 'numRays' will allow diagnostic output from raytracing. --- example/model.c | 2 +- src/lime.h | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/example/model.c b/example/model.c index 75d1e31..d3a716d 100644 --- a/example/model.c +++ b/example/model.c @@ -40,7 +40,7 @@ input(inputPars *par, image *img){ img[0].theta = 0.0; // 0: face-on, pi/2: edge-on img[0].distance = 140*PC; // source distance in m img[0].source_vel = 0; // source velocity in m/s - img[0].unit = 0; // 0:Kelvin 1:Jansky/pixel 2:SI 3:Lsun/pixel 4:tau + img[0].unit = 0; // 0:Kelvin 1:Jansky/pixel 2:SI 3:Lsun/pixel 4:tau 5:numRays img[0].filename = "image0.fits"; // Output filename } diff --git a/src/lime.h b/src/lime.h index f1d602c..45e77f0 100644 --- a/src/lime.h +++ b/src/lime.h @@ -147,6 +147,7 @@ typedef struct { double *intense; double *tau; double stokes[3]; + int numRays; } spec; /* Image information */ From 05591d037ea61d3803f6bfa1a317ae85ed0358f0 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 11 Sep 2015 13:43:59 +0200 Subject: [PATCH 02/43] Rearranged raytrace code This has been done to prepare for adding a new raytrace algorithm. --- src/raytrace.c | 68 +++++++++++++++++++++++++++++++------------------- 1 file changed, 43 insertions(+), 25 deletions(-) diff --git a/src/raytrace.c b/src/raytrace.c index a0ca190..90b1651 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -170,17 +170,20 @@ traceray(rayData ray, int tmptrans, int im, inputPars *par, struct grid *g, molD void raytrace(int im, inputPars *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); + int *counta,*countb,nlinetot,aa,numRaysThisPixel; + int ichan,gi,iline,tmptrans,i,threadI,nRaysDone; + double size,minfreq,absDeltaFreq; double cutoff; + unsigned int totalNumImagePixels,ppi; + double imgXiCentre,imgYiCentre,x[3]; + int xi,yi,ri,j; const gsl_rng_type *ranNumGenType = gsl_rng_ranlxs2; - gsl_rng *ran = gsl_rng_alloc(ranNumGenType); /* Random number generator */ + gsl_rng *randGen = gsl_rng_alloc(ranNumGenType); /* Random number generator */ #ifdef TEST - gsl_rng_set(ran,178490); + gsl_rng_set(randGen,178490); #else - gsl_rng_set(ran,time(0)); + gsl_rng_set(randGen,time(0)); #endif gsl_rng **threadRans; @@ -188,10 +191,13 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ for (i=0;inThreads;i++){ threadRans[i] = gsl_rng_alloc(ranNumGenType); - gsl_rng_set(threadRans[i],(int)gsl_rng_uniform(ran)*1e6); + gsl_rng_set(threadRans[i],(int)gsl_rng_uniform(randGen)*1e6); } size=img[im].distance*img[im].imgres; + totalNumImagePixels = img[im].pxls*img[im].pxls; + imgXiCentre = img[im].pxls/2.0; + imgYiCentre = img[im].pxls/2.0; /* Determine whether there are blended lines or not */ lineCount(par->nSpecies, m, &counta, &countb, &nlinetot); @@ -218,48 +224,61 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ } } else tmptrans=img[im].trans; - cutoff = par->minScale*1.0e-7; - - for(px=0;px<(img[im].pxls*img[im].pxls);px++){ + for(ppi=0;ppiantialias; + } + + cutoff = par->minScale*1.0e-7; + nRaysDone=0; omp_set_dynamic(0); - #pragma omp parallel private(px,aa,threadI) num_threads(par->nThreads) + #pragma omp parallel private(ppi,aa,threadI) num_threads(par->nThreads) { threadI = omp_get_thread_num(); /* Declaration of thread-private pointers */ rayData ray; - ray.intensity=malloc(sizeof(double) * img[im].nchan); - ray.tau=malloc(sizeof(double) * img[im].nchan); + ray.intensity= malloc(sizeof(double)*img[im].nchan); + ray.tau = malloc(sizeof(double)*img[im].nchan); #pragma omp for - /* Main loop through pixel grid */ - for(px=0;px<(img[im].pxls*img[im].pxls);px++){ + /* Main loop through pixel grid. ***NOTE*** though that the division of labour between threads can be uneven if par->raytraceAlgorithm==1. Probably should rearrange the code, maybe so there is only 1 'for' loop over all the rays? + */ + for(ppi=0;ppipar->antialias) + numRaysThisPixel = img[im].pixel[ppi].numRays; + else + numRaysThisPixel = par->antialias; + #pragma omp atomic ++nRaysDone; - for(aa=0;aaantialias;aa++){ - ray.x = size*(gsl_rng_uniform(threadRans[threadI])+px%img[im].pxls)-size*img[im].pxls/2.; - ray.y = size*(gsl_rng_uniform(threadRans[threadI])+px/img[im].pxls)-size*img[im].pxls/2.; + for(aa=0;aaantialias; - img[im].pixel[px].tau[ichan] += ray.tau[ichan]/(double) par->antialias; + img[im].pixel[ppi].intense[ichan] += ray.intensity[ichan]/(double)numRaysThisPixel; + img[im].pixel[ppi].tau[ ichan] += ray.tau[ ichan]/(double)numRaysThisPixel; } } } if (threadI == 0){ // i.e., is master thread - if(!silent) progressbar((double)(nRaysDone)/totalNumPixelsMinus1, 13); + if(!silent) progressbar((double)(nRaysDone)/(double)(totalNumImagePixels-1), 13); } } @@ -275,10 +294,9 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ gsl_rng_free(threadRans[i]); } free(threadRans); - gsl_rng_free(ran); + gsl_rng_free(randGen); } - void raytrace_1_4(int im, inputPars *par, struct grid *g, molData *m, image *img){ /* From 1ca8978e3f02cafd58c66323efc99c8c281d9076 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 11 Sep 2015 13:47:38 +0200 Subject: [PATCH 03/43] Added new raytrace algorithm. This algorithm addresses a shortcoming of the existing algorithm in cases where there is a large number of grid points included in a single image pixel. The previous algorithm sampled such regions poorly. The effect of the new algorithm can be seen using the standard model file: the pixels towards the centre exhibit less variation from pixel to pixel. Selection of the new algorithm is mediated by a new (optional) parameter raytraceAlgorithm. --- src/aux.c | 1 + src/lime.h | 2 +- src/raytrace.c | 37 +++++++++++++++++++++++++++++++++++-- 3 files changed, 37 insertions(+), 3 deletions(-) diff --git a/src/aux.c b/src/aux.c index f3d9d39..22d906f 100644 --- a/src/aux.c +++ b/src/aux.c @@ -40,6 +40,7 @@ parseInput(inputPars *par, image **img, molData **m){ par->sinkPoints=0; par->doPregrid=0; par->nThreads=0; + par->raytraceAlgorithm=0; /* Allocate space for output fits images */ (*img)=malloc(sizeof(image)*MAX_NSPECIES); diff --git a/src/lime.h b/src/lime.h index 45e77f0..29e7f71 100644 --- a/src/lime.h +++ b/src/lime.h @@ -81,7 +81,7 @@ /* input parameters */ typedef struct { double radius,radiusSqu,minScale,minScaleSqu,tcmb,taylorCutoff; - int ncell,sinkPoints,pIntensity,nImages,nSpecies,blend; + int ncell,sinkPoints,pIntensity,nImages,nSpecies,blend,raytraceAlgorithm; char *outputfile, *binoutputfile, *inputfile; char *gridfile; char *pregrid; diff --git a/src/raytrace.c b/src/raytrace.c index 90b1651..1da4bd0 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -231,8 +231,41 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ } } - for(ppi=0;ppiantialias; + if(par->raytraceAlgorithm==0){ + for(ppi=0;ppiantialias; + } + + } else if(par->raytraceAlgorithm==1){ + /* +In this algorithm we attempt to deal better with image pixels which cover areas of the model where grid points are packed much more densely than the pixel spacing. In algorithm 0, such regions are sampled by only at par->antialias rays per pixel, which risks missing extreme values. The approach here is to generate approximately the same number of rays per pixel as the number of projected grid points in that pixel (with par->antialias providing a minimum floor). The set of values per pixel are then averaged to obtain the image value for that pixel. + */ + + for(ppi=0;ppipIntensity;gi++){ + /* Apply the inverse (i.e. transpose) rotation matrix. (We use the inverse matrix here because here we want to rotate grid coordinates to the observer frame, whereas inside traceray() we rotate observer coordinates to the grid frame.) + */ + for(i=0;i<2;i++){ + x[i]=0.0; + for(j=0;j=0 && ppiminScale*1.0e-7; From 607709c4f11922cd7acecd6d6fdb379840132ccb Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 11 Sep 2015 13:51:49 +0200 Subject: [PATCH 04/43] Removed function raytrace_1_4() Its main idea has now been implemented in the new raytrace algorithm. --- src/raytrace.c | 241 ------------------------------------------------- 1 file changed, 241 deletions(-) diff --git a/src/raytrace.c b/src/raytrace.c index 1da4bd0..0db79e6 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -330,244 +330,3 @@ In this algorithm we attempt to deal better with image pixels which cover areas gsl_rng_free(randGen); } -void -raytrace_1_4(int im, inputPars *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. -I've adapted it slightly so it makes use of the function traceray(), -which was modified from the function tracerays() in v1.4. This algorithm -is not currently used, but may be useful as an option; that's why I have -kept it. - */ - - int *counta, *countb,nlinetot; - int ichan,i,px,iline,tmptrans,count; - double size,xp,yp,minfreq,absDeltaFreq; - double cutoff; - - gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2); /* Random number generator */ -#ifdef TEST - gsl_rng_set(ran,178490); -#else - gsl_rng_set(ran,time(0)); -#endif - rayData *rays; - - int sg,n; - double cx,cy; - - double x1,x2,x3,y1,y2,y3,z1,z2,z3,xt[3],yt[3],di,p,d1,d2,d3,temp1; - int zt[3]; - int c; - - char flags[255]; - boolT ismalloc = False; - facetT *facet, *neighbor, **neighborp;; - vertexT *vertex,**vertexp; - coordT *pt_array; - - int id; - coordT point[3]; - boolT isoutside; - realT bestdist; - - size=img[im].distance*img[im].imgres; - - /* Determine whether there are blended lines or not */ - lineCount(par->nSpecies, m, &counta, &countb, &nlinetot); - if(img[im].doline==0) nlinetot=1; - - /* Fix the image parameters */ - if(img[im].freq < 0) img[im].freq=m[0].freq[img[im].trans]; - if(img[im].nchan == 0 && img[im].bandwidth>0){ - img[im].nchan=(int) (img[im].bandwidth/(img[im].velres/CLIGHT*img[im].freq)); - } else if (img[im].velres<0 && img[im].bandwidth>0){ - img[im].velres = img[im].bandwidth*CLIGHT/img[im].freq/img[im].nchan; - } else img[im].bandwidth = img[im].nchan*img[im].velres/CLIGHT * img[im].freq; - - if(img[im].trans<0){ - iline=0; - minfreq=fabs(img[im].freq-m[0].freq[iline]); - tmptrans=iline; - for(iline=1;ilinepIntensity)); - - for(i=0;ipIntensity;i++){ - rays[i].x=g[i].x[0]; - rays[i].y=g[i].x[1]; - rays[i].tau=malloc(sizeof(double) * img[im].nchan); - rays[i].intensity=malloc(sizeof(double) * img[im].nchan); - for(ichan=0;ichanpIntensity); - - for(i=0;ipIntensity;i++) { - pt_array[i*2+0]=rays[i].x; - pt_array[i*2+1]=rays[i].y; - } - - sprintf(flags,"qhull v s Qbb T0"); - if (!qh_new_qhull(2, par->pIntensity, pt_array, ismalloc, flags, NULL, NULL)) { - - qh_setvoronoi_all(); - - FORALLvertices { - i=qh_pointid(vertex->point); - - cx=0.; - cy=0.; - n=0; - FOREACHneighbor_(vertex) { - if (!neighbor->upperdelaunay) n++; - } - if(n>0){ - - - } else { - if(!silent) bail_out("Qhull error"); - exit(0); - } - - FOREACHneighbor_(vertex) { - if (!neighbor->upperdelaunay) { - cx+=neighbor->center[0]; - cy+=neighbor->center[1]; - } - } - - rays[i].x = rays[i].x - (rays[i].x-cx/ (double) n)*0.1; - rays[i].y = rays[i].y - (rays[i].y-cy/ (double) n)*0.1; - } - } else { - printf("qhull error\n"); - } - - qh_freeqhull(!qh_ALL); - free(pt_array); - } - - cutoff = par->minScale*1.0e-7; - - /* Main loop through rays */ - count=0; - for(px=0;pxpIntensity;px++){ - traceray(rays[px], tmptrans, im, par, g, m, img, nlinetot, counta, countb, cutoff); - ++count; - if(!silent) progressbar((double)(count)/(double)(par->pIntensity-1), 13); - } - - /* Remap rays onto pixel grid */ - pt_array=malloc(2*sizeof(coordT)*par->pIntensity); - - for(i=0;ipIntensity;i++) { - pt_array[i*2+0]=rays[i].x; - pt_array[i*2+1]=rays[i].y; - } - -/* This allocation belongs to "Shepard's method" below - d=malloc(sizeof(double)*par->pIntensity); -*/ - size=img[im].distance*img[im].imgres; - - sprintf(flags,"qhull d Qbb"); - if (!qh_new_qhull(2, par->pIntensity, pt_array, ismalloc, flags, NULL, NULL)) { - for(px=0;pxpIntensity;i++){ - // d[i]=1./pow(sqrt(pow(xp-rays[i].x,2)+ pow(yp-rays[i].y,2)),8.); - temp1 = (xp-rays[i].x)*(xp-rays[i].x)+ (yp-rays[i].y)*(yp-rays[i].y) - d[i]=1./(temp1*temp1*temp1*temp1); - img[im].pixel[px].intense[ichan] += rays[i].intensity[ichan]*d[i]; -**** how to handle img[im].pixel[px].tau[ichan]? - di+=d[i]; - } - img[im].pixel[px].intense[ichan] /= di; - } -*/ - - - point[0]=xp; - point[1]=yp; - point[2]=0.; - - qh_setdelaunay (3, 1, point); - facet= qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside); - - c=0; - FOREACHvertex_( facet->vertices ) { - id=qh_pointid(vertex->point); - xt[c]=rays[id].x; yt[c]=rays[id].y; zt[c]=id; - c++; - } - - - x1=xt[0];x2=xt[1];x3=xt[2]; - y1=yt[0];y2=yt[1];y3=yt[2]; - - for(ichan=0;ichanpIntensity;i++){ - free(rays[i].tau); - free(rays[i].intensity); - } - free(rays); - free(counta); - free(countb); -} - - From 86be9af21e74df36cd6c67f6db76839ae6e774d7 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 11 Sep 2015 13:57:01 +0200 Subject: [PATCH 05/43] Rearranged and renamed writefits() This function was made a little more readable and perhaps efficient, and renamed to write3Dfits() in preparation for adding a 2D version. --- src/lime.h | 2 +- src/main.c | 2 +- src/raytrace.c | 1 + src/writefits.c | 43 +++++++++++++++++++++++++++++-------------- 4 files changed, 32 insertions(+), 16 deletions(-) diff --git a/src/lime.h b/src/lime.h index 29e7f71..3c58d77 100644 --- a/src/lime.h +++ b/src/lime.h @@ -241,7 +241,7 @@ void traceray(rayData, int, int, inputPars *, struct grid *, molData *, image 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 write3Dfits(int, inputPars *, molData *, image *); void write_VTK_unstructured_Points(inputPars *, struct grid *); int factorial(const int n); double taylor(const int maxOrder, const float x); diff --git a/src/main.c b/src/main.c index aa92ea9..1835720 100644 --- a/src/main.c +++ b/src/main.c @@ -62,7 +62,7 @@ int main () { } raytrace(i,&par,g,m,img); - writefits(i,&par,m,img); + write3Dfits(i,&par,m,img); } if(!silent) goodnight(initime,img[0].filename); diff --git a/src/raytrace.c b/src/raytrace.c index 0db79e6..9a72788 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -330,3 +330,4 @@ In this algorithm we attempt to deal better with image pixels which cover areas gsl_rng_free(randGen); } + diff --git a/src/writefits.c b/src/writefits.c index 2ba9971..f899842 100644 --- a/src/writefits.c +++ b/src/writefits.c @@ -10,10 +10,10 @@ #include "lime.h" void -writefits(int im, inputPars *par, molData *m, image *img){ +write3Dfits(int im, inputPars *par, molData *m, image *img){ double bscale,bzero,epoch,lonpole,equinox,restfreq; double cdelt1,crpix1,crval1,cdelt2,crpix2,crval2; - double cdelt3,crpix3,crval3,ru3; + double cdelt3,crpix3,crval3,ru3,scale; int velref; float *row; int px,py,ichan; @@ -23,6 +23,7 @@ writefits(int im, inputPars *par, molData *m, image *img){ long naxes[3]; long int fpixels[3],lpixels[3]; char negfile[100]="! "; + unsigned long ppi; row = malloc(sizeof(*row)*img[im].pxls); @@ -30,7 +31,7 @@ writefits(int im, inputPars *par, molData *m, image *img){ naxes[1]=img[im].pxls; if(img[im].doline==1) naxes[2]=img[im].nchan; else if(img[im].doline==0 && par->polarization) naxes[2]=3; - else naxes[2]=1; + else naxes[2]=1;//********** should call write2Dfits for this. fits_create_file(&fptr, img[im].filename, &status); @@ -88,20 +89,33 @@ writefits(int im, inputPars *par, molData *m, image *img){ if(img[im].unit==1) fits_write_key(fptr, TSTRING, "BUNIT", &"JY/PIXEL", "", &status); if(img[im].unit==2) fits_write_key(fptr, TSTRING, "BUNIT", &"WM2HZSR ", "", &status); if(img[im].unit==3) fits_write_key(fptr, TSTRING, "BUNIT", &"Lsun/PX ", "", &status); - if(img[im].unit==3) fits_write_key(fptr, TSTRING, "BUNIT", &" ", "", &status); + if(img[im].unit==4) fits_write_key(fptr, TSTRING, "BUNIT", &" ", "", &status); + + if( img[im].unit==0) + scale=(CLIGHT/img[im].freq)*(CLIGHT/img[im].freq)/2./KBOLTZ*m[0].norm; + else if(img[im].unit==1) + scale=1e26*img[im].imgres*img[im].imgres*m[0].norm; + else if(img[im].unit==2) + scale=m[0].norm; + else if(img[im].unit==3) { + ru3 = img[im].distance/1.975e13; + scale=4.*PI*ru3*ru3*img[im].freq*img[im].imgres*img[im].imgres*m[0].norm; + } + else if(img[im].unit!=4) { + if(!silent) bail_out("Image unit number invalid"); + exit(0); + } /* Write FITS data */ - for(py=0;py-1 && img[im].unit<4) + row[px]=(float) img[im].pixel[ppi].intense[ichan]*scale; + else if(img[im].unit==4) + row[px]=(float) img[im].pixel[ppi].tau[ichan]; else { if(!silent) bail_out("Image unit number invalid"); exit(0); @@ -123,3 +137,4 @@ writefits(int im, inputPars *par, molData *m, image *img){ if(!silent) done(13); } + From b7ba6f3e629083431d0e569fdd89bf50c77e6b66 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 11 Sep 2015 14:02:36 +0200 Subject: [PATCH 06/43] New function to write 2D FITS At present this is just called to write diagnostic images showing the number of rays per pixel, invoked by setting the image unit to the new value of 5. It could also be called to wrote continuum images, instead of writing a cube with just 1 layer as at present. --- src/lime.h | 1 + src/main.c | 10 +++- src/writefits.c | 127 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 137 insertions(+), 1 deletion(-) diff --git a/src/lime.h b/src/lime.h index 3c58d77..39a3000 100644 --- a/src/lime.h +++ b/src/lime.h @@ -241,6 +241,7 @@ void traceray(rayData, int, int, inputPars *, struct grid *, molData *, image void velocityspline(struct grid *, int, int, double, double, double*); void velocityspline2(double *, double *, double, double, double, double*); double veloproject(double *, double *); +void write2Dfits(int, inputPars *, molData *, image *); void write3Dfits(int, inputPars *, molData *, image *); void write_VTK_unstructured_Points(inputPars *, struct grid *); int factorial(const int n); diff --git a/src/main.c b/src/main.c index 1835720..3a4ac31 100644 --- a/src/main.c +++ b/src/main.c @@ -62,7 +62,15 @@ int main () { } raytrace(i,&par,g,m,img); - write3Dfits(i,&par,m,img); + + if(img[i].unit<5) + write3Dfits(i,&par,m,img); + else if(img[i].unit==5) + write2Dfits(i,&par,m,img); + else{ + if(!silent) bail_out("Image unit number invalid"); + exit(0); + } } if(!silent) goodnight(initime,img[0].filename); diff --git a/src/writefits.c b/src/writefits.c index f899842..4ac331e 100644 --- a/src/writefits.c +++ b/src/writefits.c @@ -138,3 +138,130 @@ write3Dfits(int im, inputPars *par, molData *m, image *img){ if(!silent) done(13); } +void +write2Dfits(int im, inputPars *par, molData *m, image *img){ + double bscale,bzero,epoch,lonpole,equinox,restfreq; + double cdelt1,crpix1,crval1,cdelt2,crpix2,crval2; + double cdelt3,crpix3,crval3,ru3,scale; + int velref; + float *row,minVal; + int px,py; + fitsfile *fptr; + int status = 0; + int naxis=2, bitpix=-32; + long naxes[2]; + long int fpixels[2],lpixels[2]; + char negfile[100]="! "; + unsigned long ppi; + + row = malloc(sizeof(*row)*img[im].pxls); + + naxes[0]=img[im].pxls; + naxes[1]=img[im].pxls; + if(img[im].unit!=5 && (img[im].doline==1 || (img[im].doline==0 && par->polarization))){ + if(!silent) bail_out("You need to write a 3D FITS output in this case"); + exit(0); + } + + fits_create_file(&fptr, img[im].filename, &status); + + if(status!=0){ + if(!silent) warning("Overwriting existing fits file "); + status=0; + strcat(negfile,img[im].filename); + fits_create_file(&fptr, negfile, &status); + } + + /* Write FITS header */ + fits_create_img(fptr, bitpix, naxis, naxes, &status); + epoch =2.0e3; + lonpole =1.8e2; + equinox =2.0e3; + restfreq=img[im].freq; + velref =257; + cdelt1 =-1.8e2*img[im].imgres/PI; + crpix1 =(double) img[im].pxls/2+0.5; + crval1 =0.0e0; + cdelt2 =1.8e2*img[im].imgres/PI; + crpix2 =(double) img[im].pxls/2+0.5; + crval2 =0.0e0; + bscale =1.0e0; + bzero =0.0e0; + + fits_write_key(fptr, TSTRING, "OBJECT ", &"LIMEMDL ", "", &status); + fits_write_key(fptr, TDOUBLE, "EPOCH ", &epoch, "", &status); + fits_write_key(fptr, TDOUBLE, "LONPOLE ", &lonpole, "", &status); + fits_write_key(fptr, TDOUBLE, "EQUINOX ", &equinox, "", &status); + fits_write_key(fptr, TSTRING, "SPECSYS ", &"LSRK ", "", &status); + fits_write_key(fptr, TDOUBLE, "RESTFREQ", &restfreq, "", &status); + fits_write_key(fptr, TINT, "VELREF ", &velref, "", &status); + fits_write_key(fptr, TSTRING, "CTYPE1 ", &"RA---SIN", "", &status); + fits_write_key(fptr, TDOUBLE, "CDELT1 ", &cdelt1, "", &status); + fits_write_key(fptr, TDOUBLE, "CRPIX1 ", &crpix1, "", &status); + fits_write_key(fptr, TDOUBLE, "CRVAL1 ", &crval1, "", &status); + fits_write_key(fptr, TSTRING, "CUNIT1 ", &"DEG ", "", &status); + fits_write_key(fptr, TSTRING, "CTYPE2 ", &"DEC--SIN", "", &status); + fits_write_key(fptr, TDOUBLE, "CDELT2 ", &cdelt2, "", &status); + fits_write_key(fptr, TDOUBLE, "CRPIX2 ", &crpix2, "", &status); + fits_write_key(fptr, TDOUBLE, "CRVAL2 ", &crval2, "", &status); + fits_write_key(fptr, TSTRING, "CUNIT2 ", &"DEG ", "", &status); + fits_write_key(fptr, TDOUBLE, "BSCALE ", &bscale, "", &status); + fits_write_key(fptr, TDOUBLE, "BZERO ", &bzero, "", &status); + if(img[im].unit==0) fits_write_key(fptr, TSTRING, "BUNIT", &"K ", "", &status); + if(img[im].unit==1) fits_write_key(fptr, TSTRING, "BUNIT", &"JY/PIXEL", "", &status); + if(img[im].unit==2) fits_write_key(fptr, TSTRING, "BUNIT", &"WM2HZSR ", "", &status); + if(img[im].unit==3) fits_write_key(fptr, TSTRING, "BUNIT", &"Lsun/PX ", "", &status); + if(img[im].unit==4) fits_write_key(fptr, TSTRING, "BUNIT", &" ", "", &status); + if(img[im].unit==5) fits_write_key(fptr, TSTRING, "BUNIT", &"N_RAYS ", "", &status); + + if( img[im].unit==0) + scale=(CLIGHT/img[im].freq)*(CLIGHT/img[im].freq)/2./KBOLTZ*m[0].norm; + else if(img[im].unit==1) + scale=1e26*img[im].imgres*img[im].imgres*m[0].norm; + else if(img[im].unit==2) + scale=m[0].norm; + else if(img[im].unit==3) { + ru3 = img[im].distance/1.975e13; + scale=4.*PI*ru3*ru3*img[im].freq*img[im].imgres*img[im].imgres*m[0].norm; + } + else if(img[im].unit!=4 && img[im].unit!=5) { + if(!silent) bail_out("Image unit number invalid"); + exit(0); + } + + if(img[im].unit<5) + minVal = eps; + else + minVal = 0.0; + + /* Write FITS data */ + for(py=0;py-1 && img[im].unit<4) + row[px]=(float) img[im].pixel[ppi].intense[0]*scale; + else if(img[im].unit==4) + row[px]=(float) img[im].pixel[ppi].tau[0]; + else if(img[im].unit==5) + row[px]=(float) img[im].pixel[ppi].numRays; + else { + if(!silent) bail_out("Image unit number invalid"); + exit(0); + } + if (fabs(row[px]) Date: Fri, 8 Apr 2016 10:39:42 +0200 Subject: [PATCH 07/43] Added new structs. This is the first of many commits I plan for an eventual pull request which will provide LIME with a smoother raytracing capability. This added functionality requires quite a lot of new code and will also involve some cleaning up of existing code. The smoother raytracing algorithm works by linear interpolation of level populations inside the Delaunay cells of the model grid. This first commit lays the groundwork for reading and using the cell information by defining several new structs to allow convenient storage and manipulation of this data. Structs added were: cell, pop2, gridInterp, gAuxType, intersectType, faceType and triangle2D. None of these are yet active so as yet they cause no change in the task function. Also added in this commit is a new user parameter traceRayAlgorithm. This will be used to invoke the new raytracing algorithm; at present it has no effect. Finally, I have commented out the function template for FastExp() in lime.h and added the variable type of the input argument of greetings_parallel() in messages.c just to avoid the warning messages that gcc issues now on compilation. This is just for present convenience, I plan to make these minor fixes the subject of a separate PR. --- src/aux.c | 1 + src/lime.h | 49 +++++++++++++++++++++++++++++++++++++++++++++++-- src/messages.c | 2 +- 3 files changed, 49 insertions(+), 3 deletions(-) diff --git a/src/aux.c b/src/aux.c index f3d9d39..a8b4a4c 100644 --- a/src/aux.c +++ b/src/aux.c @@ -40,6 +40,7 @@ parseInput(inputPars *par, image **img, molData **m){ par->sinkPoints=0; par->doPregrid=0; par->nThreads=0; + par->traceRayAlgorithm=0; /* Allocate space for output fits images */ (*img)=malloc(sizeof(image)*MAX_NSPECIES); diff --git a/src/lime.h b/src/lime.h index f1d602c..755c2f9 100644 --- a/src/lime.h +++ b/src/lime.h @@ -81,7 +81,7 @@ /* input parameters */ typedef struct { double radius,radiusSqu,minScale,minScaleSqu,tcmb,taylorCutoff; - int ncell,sinkPoints,pIntensity,nImages,nSpecies,blend; + int ncell,sinkPoints,pIntensity,nImages,nSpecies,blend,traceRayAlgorithm; char *outputfile, *binoutputfile, *inputfile; char *gridfile; char *pregrid; @@ -173,6 +173,51 @@ typedef struct { typedef struct {double x,y, *intensity, *tau;} rayData; +/* NOTE that it is assumed that vertx[i] is opposite the face that abuts with neigh[i] for all i. +*/ +struct cell { + struct grid *vertx[DIM+1]; + struct cell *neigh[DIM+1]; /* ==NULL flags an external face. */ + unsigned long id; + double centre[DIM]; +}; + +struct pop2 { + double *specNumDens, *knu, *dust; + double binv; +}; + +typedef struct{ + double x[DIM], xCmpntRay, B[3]; + struct pop2 *mol; +} gridInterp; + +struct gAuxType{ + struct pop2 *mol; +}; + +/* This struct is meant to record all relevant information about the intersection between a ray (defined by a direction unit vector 'dir' and a starting position 'r') and a face of a Delaunay cell. +*/ +typedef struct { + int fi; + /* The index (in the range {0...DIM}) of the face (and thus of the opposite vertex, i.e. the one 'missing' from the bary[] list of this face). + */ + int orientation; + /* >0 means the ray exits, <0 means it enters, ==0 means the face is parallel to ray. + */ + double bary[DIM], dist, collPar; + /* 'dist' is defined via r_int = r + dist*dir. 'collPar' is a measure of how close to any edge of the face r_int lies. + */ +} intersectType; + +typedef struct { + double r[DIM][DIM], centre[DIM];/*, norm[3], mat[1][1], det; */ +} faceType; + +typedef struct { + double xAxis[DIM], yAxis[DIM], r[3][2]; +} triangle2D; + /* Some functions */ @@ -246,7 +291,7 @@ int factorial(const int n); double taylor(const int maxOrder, const float x); void calcFastExpRange(const int maxTaylorOrder, const int maxNumBitsPerMantField, int *numMantissaFields, int *lowestExponent, int *numExponentsUsed); void calcTableEntries(const int maxTaylorOrder, const int maxNumBitsPerMantField); -inline double FastExp(const float negarg); +//inline double FastExp(const float negarg); /* Curses functions */ diff --git a/src/messages.c b/src/messages.c index cf7eb15..df05ce9 100644 --- a/src/messages.c +++ b/src/messages.c @@ -39,7 +39,7 @@ greetings(){ } void -greetings_parallel(numThreads){ +greetings_parallel(int numThreads){ #ifdef NO_NCURSES if (numThreads>1){ From 0f1725aec2ab05793de86078e117be44bbc7f8e7 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 11:50:05 +0200 Subject: [PATCH 08/43] Moved nmol from grid to populations This attribute belongs in struct populations rather than struct grid (directly). Moving it involved making changes to 18 lines of code over 7 modules. I checked afterward that the code runs still ok with the default model. I also slightly rearranged attributes x and vel in the struct template for struct grid, and added a new attribute (as yet unused) B. These changes have no further effects on the code. --- src/LTEsolution.c | 2 +- src/grid.c | 7 +------ src/lime.h | 11 +++++------ src/molinit.c | 8 ++++---- src/popsin.c | 9 +++++---- src/popsout.c | 7 +++++-- src/predefgrid.c | 2 +- src/sourcefunc.c | 8 ++++---- 8 files changed, 26 insertions(+), 28 deletions(-) diff --git a/src/LTEsolution.c b/src/LTEsolution.c index 874213d..aa1cbeb 100644 --- a/src/LTEsolution.c +++ b/src/LTEsolution.c @@ -16,7 +16,7 @@ LTE(inputPars *par, struct grid *g, molData *m){ for(ispec=0;ispecnSpecies;ispec++){ for(id=0;idpIntensity;id++){ - g[id].nmol[ispec]=g[id].abun[ispec]*g[id].dens[0]; + g[id].mol[ispec].nmol=g[id].abun[ispec]*g[id].dens[0]; z=0; for(ilev=0;ilevcollPart); (*g)[i].abun=malloc(sizeof(double)*par->nSpecies); - (*g)[i].nmol=malloc(sizeof(double)*par->nSpecies); (*g)[i].t[0]=-1; (*g)[i].t[1]=-1; } @@ -92,7 +91,7 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ int i; if( g != NULL ) { - for(i=0;i<(par->pIntensity+par->sinkPoints); i++){ + for(i=0;i<(par->ncell); i++){ if(g[i].a0 != NULL) { free(g[i].a0); @@ -129,10 +128,6 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ { free(g[i].dens); } - if(g[i].nmol != NULL) - { - free(g[i].nmol); - } if(g[i].abun != NULL) { free(g[i].abun); diff --git a/src/lime.h b/src/lime.h index 755c2f9..9ce32e4 100644 --- a/src/lime.h +++ b/src/lime.h @@ -120,16 +120,15 @@ struct rates { struct populations { - double * pops, *knu, *dust; - double dopb, binv; + double *pops, *knu, *dust; + double dopb, binv, nmol; struct rates *partner; }; /* Grid properties */ struct grid { int id; - double x[3]; - double vel[3]; + double x[DIM], vel[DIM], B[3]; /* B field only makes physical sense in 3 dimensions. */ double *a0,*a1,*a2,*a3,*a4; int numNeigh; point *dir; @@ -138,9 +137,9 @@ struct grid { int sink; int nphot; int conv; - double *dens,t[2],*nmol,*abun, dopb; + double *dens,t[2],*abun, dopb; double *ds; - struct populations* mol; + struct populations *mol; }; typedef struct { diff --git a/src/molinit.c b/src/molinit.c index 4ce8ea2..845cd88 100644 --- a/src/molinit.c +++ b/src/molinit.c @@ -241,16 +241,16 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){ for(id=0;idncell; id++){ for(ispec=0;ispecnSpecies;ispec++){ if(m[i].npart == 1 && (count[0] == 1 || count[0] == 2 || count[0] == 3)){ - g[id].nmol[ispec]=g[id].abun[ispec]*g[id].dens[0]; + g[id].mol[ispec].nmol=g[id].abun[ispec]*g[id].dens[0]; } else if(m[i].npart == 2 && (count[0] == 2 || count[0] == 3) && (count[1] == 2 || count[1] == 3)){ if(!flag){ - g[id].nmol[ispec]=g[id].abun[ispec]*(g[id].dens[0]+g[id].dens[1]); + g[id].mol[ispec].nmol=g[id].abun[ispec]*(g[id].dens[0]+g[id].dens[1]); } else { - g[id].nmol[ispec]=g[id].abun[ispec]*g[id].dens[0]; + g[id].mol[ispec].nmol=g[id].abun[ispec]*g[id].dens[0]; if(!silent) warning("Calculating molecular density with respect to first collision partner only"); } } else if(m[i].npart > 2 && !flag){ - g[id].nmol[ispec]=g[id].abun[ispec]*(g[id].dens[0]+g[id].dens[1]); + g[id].mol[ispec].nmol=g[id].abun[ispec]*(g[id].dens[0]+g[id].dens[1]); if(!silent) warning("Calculating molecular density with respect first and second collision partner"); } } diff --git a/src/popsin.c b/src/popsin.c index 5392664..4b7e279 100644 --- a/src/popsin.c +++ b/src/popsin.c @@ -5,6 +5,8 @@ * Copyright (C) 2006-2014 Christian Brinch * Copyright (C) 2015 The LIME development team * +***TODO: + - Change the definition of the file format so that nmol is now written with the other mol[] scalars. */ #include "lime.h" @@ -74,7 +76,6 @@ popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){ (*g)[i].a4 = NULL; (*g)[i].dens = NULL; (*g)[i].abun = NULL; - (*g)[i].nmol = NULL; (*g)[i].dir = NULL; (*g)[i].neigh = NULL; (*g)[i].w = NULL; @@ -83,10 +84,10 @@ popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){ fread(&(*g)[i].x, sizeof (*g)[i].x, 1, fp); fread(&(*g)[i].vel, sizeof (*g)[i].vel, 1, fp); fread(&(*g)[i].sink, sizeof (*g)[i].sink, 1, fp); - (*g)[i].nmol=malloc(par->nSpecies*sizeof(double)); - for(j=0;jnSpecies;j++) fread(&(*g)[i].nmol[j], sizeof(double), 1, fp); - fread(&(*g)[i].dopb, sizeof (*g)[i].dopb, 1, fp); (*g)[i].mol=malloc(par->nSpecies*sizeof(struct populations)); + for(j=0;jnSpecies;j++) + fread(&(*g)[i].mol[j].nmol, sizeof(double), 1, fp); + fread(&(*g)[i].dopb, sizeof (*g)[i].dopb, 1, fp); for(j=0;jnSpecies;j++){ (*g)[i].mol[j].pops=malloc(sizeof(double)*(*m)[j].nlev); for(k=0;k<(*m)[j].nlev;k++) fread(&(*g)[i].mol[j].pops[k], sizeof(double), 1, fp); diff --git a/src/popsout.c b/src/popsout.c index 93ff3e3..c41639b 100644 --- a/src/popsout.c +++ b/src/popsout.c @@ -5,6 +5,8 @@ * Copyright (C) 2006-2014 Christian Brinch * Copyright (C) 2015 The LIME development team * +***TODO: + - Change the definition of the file format so that nmol is now read with the other mol[] scalars. */ #include "lime.h" @@ -25,7 +27,7 @@ popsout(inputPars *par, struct grid *g, molData *m){ for(j=0;jncell-par->sinkPoints;j++){ dens=0.; for(l=0;lcollPart;l++) dens+=g[j].dens[l]; - fprintf(fp,"%e %e %e %e %e %e %d ", g[j].x[0], g[j].x[1], g[j].x[2], dens, g[j].t[0], g[j].nmol[0]/dens, g[j].conv); + fprintf(fp,"%e %e %e %e %e %e %d ", g[j].x[0], g[j].x[1], g[j].x[2], dens, g[j].t[0], g[j].mol[0].nmol/dens, g[j].conv); for(k=0;knSpecies,1, fp); + for(j=0;jnSpecies;j++) + fwrite(&g[i].mol[j].nmol, sizeof(double), 1, fp); fwrite(&g[i].dopb, sizeof g[i].dopb, 1, fp); for(j=0;jnSpecies;j++){ fwrite(g[i].mol[j].pops, sizeof(double)*m[j].nlev, 1, fp); diff --git a/src/predefgrid.c b/src/predefgrid.c index 4aed807..de14811 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -40,7 +40,7 @@ predefinedGrid(inputPars *par, struct grid *g){ 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].mol[0].nmol=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); diff --git a/src/sourcefunc.c b/src/sourcefunc.c index 8377ced..58a6223 100644 --- a/src/sourcefunc.c +++ b/src/sourcefunc.c @@ -18,7 +18,7 @@ sourceFunc(double *snu, double *dtau, double ds, molData *m,double vfac,struct g jnu = g[pos].mol[ispec].dust[iline]*g[pos].mol[ispec].knu[iline]; /* Line part: j_nu = v*consts*1/b*rho*n_i*A_ij */ - if(doline) jnu += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].nmol[ispec]*g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].aeinst[iline]; + if(doline) jnu += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].mol[ispec].nmol*g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].aeinst[iline]; @@ -28,7 +28,7 @@ sourceFunc(double *snu, double *dtau, double ds, molData *m,double vfac,struct g /* Line part: alpha_nu = v*const*1/b*rho*(n_j*B_ij-n_i*B_ji) */ - if(doline) alpha += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].nmol[ispec]*(g[pos].mol[ispec].pops[m[ispec].lal[iline]]*m[ispec].beinstl[iline] + if(doline) alpha += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].mol[ispec].nmol*(g[pos].mol[ispec].pops[m[ispec].lal[iline]]*m[ispec].beinstl[iline] -g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].beinstu[iline]); @@ -47,10 +47,10 @@ void sourceFunc_line(double *jnu, double *alpha, molData *m,double vfac,struct grid *g,int pos,int ispec, int iline){ /* Line part: j_nu = v*consts*1/b*rho*n_i*A_ij */ - *jnu += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].nmol[ispec]*g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].aeinst[iline]; + *jnu += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].mol[ispec].nmol*g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].aeinst[iline]; /* Line part: alpha_nu = v*const*1/b*rho*(n_j*B_ij-n_i*B_ji) */ - *alpha += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].nmol[ispec]*(g[pos].mol[ispec].pops[m[ispec].lal[iline]]*m[ispec].beinstl[iline] + *alpha += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].mol[ispec].nmol*(g[pos].mol[ispec].pops[m[ispec].lal[iline]]*m[ispec].beinstl[iline] -g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].beinstu[iline]); return; From 0a45fd759db0588b3a8a8163b94b8fab6122bee0 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 12:00:27 +0200 Subject: [PATCH 09/43] qhull to delaunay This is the first of several commits which make changes to the function qhull(). The only one of these which introduces new functionality will be an option to read and return the information about the Delaunay cells. However I plan to use this necessity to reform the whole function, starting with the name. 'qhull' is a bad name because it is not immediately explanatory and because it can be confused with the library of that name. --- src/grid.c | 4 ++-- src/lime.h | 2 +- src/popsin.c | 2 +- src/predefgrid.c | 2 +- src/smooth.c | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/grid.c b/src/grid.c index 00440ec..81a7965 100644 --- a/src/grid.c +++ b/src/grid.c @@ -146,7 +146,7 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ } void -qhull(inputPars *par, struct grid *g){ +delaunay(inputPars *par, struct grid *g){ int i,j,k,id; char flags[255]; boolT ismalloc = False; @@ -608,7 +608,7 @@ buildGrid(inputPars *par, struct grid *g){ } /* end grid allocation */ - qhull(par, g); + delaunay(par, g); distCalc(par, g); smooth(par,g); diff --git a/src/lime.h b/src/lime.h index 9ce32e4..86bb732 100644 --- a/src/lime.h +++ b/src/lime.h @@ -260,7 +260,7 @@ void lineCount(int,molData *,int **, int **, int *); void LTE(inputPars *, struct grid *, molData *); void molinit(molData *, inputPars *, struct grid *,int); void openSocket(inputPars *par, int); -void qhull(inputPars *, struct grid *); +void delaunay(inputPars *, struct grid *); void photon(int, struct grid *, molData *, int, const gsl_rng *,inputPars *,blend *,gridPointData *,double *); void parseInput(inputPars *, image **, molData **); double planckfunc(int, double, molData *, int); diff --git a/src/popsin.c b/src/popsin.c index 4b7e279..752cc25 100644 --- a/src/popsin.c +++ b/src/popsin.c @@ -105,7 +105,7 @@ popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){ } fclose(fp); - qhull(par, *g); + delaunay(par, *g); distCalc(par, *g); getVelosplines(par,*g); *popsdone=1; diff --git a/src/predefgrid.c b/src/predefgrid.c index de14811..deec8b0 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -69,7 +69,7 @@ predefinedGrid(inputPars *par, struct grid *g){ } fclose(fp); - qhull(par,g); + delaunay(par,g); distCalc(par,g); // getArea(par,g, ran); // getMass(par,g, ran); diff --git a/src/smooth.c b/src/smooth.c index 1e98310..c84acd8 100644 --- a/src/smooth.c +++ b/src/smooth.c @@ -68,7 +68,7 @@ smooth(inputPars *par, struct grid *g){ } } - qhull(par, g); + delaunay(par, g); distCalc(par, g); if(!silent) progressbar((double)(sg+1)/(double)smooth, 5); } From cfc7ebb82ac30d4ae53f2d7b33e2c6bd70224a99 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 12:14:40 +0200 Subject: [PATCH 10/43] New arguments for delaunay() I replaced the input argument par with two, namely numDims and numPoints. This helps to make the function more generic. An additional variable ppi of type matching numPoints was added to loop over the latter quantity. Iteration variables of types other than the upper limit variable can cause strange errors! --- src/grid.c | 28 +++++++++++++--------------- src/lime.h | 2 +- src/popsin.c | 2 +- src/predefgrid.c | 2 +- src/smooth.c | 2 +- 5 files changed, 17 insertions(+), 19 deletions(-) diff --git a/src/grid.c b/src/grid.c index 81a7965..5fe10a3 100644 --- a/src/grid.c +++ b/src/grid.c @@ -146,7 +146,7 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ } void -delaunay(inputPars *par, struct grid *g){ +delaunay(const int numDims, struct grid *g, const unsigned long numPoints){ int i,j,k,id; char flags[255]; boolT ismalloc = False; @@ -155,17 +155,17 @@ delaunay(inputPars *par, struct grid *g){ coordT *pt_array; int simplex[DIM+1]; int curlong, totlong; + unsigned long ppi; - pt_array=malloc(DIM*sizeof(coordT)*par->ncell); - - for(i=0;incell;i++) { - for(j=0;jncell, pt_array, ismalloc, flags, NULL, NULL)) { + if (!qh_new_qhull(DIM, (int)numPoints, pt_array, ismalloc, flags, NULL, NULL)) { /* Identify points */ FORALLvertices { id=qh_pointid(vertex->point); @@ -203,15 +203,13 @@ delaunay(inputPars *par, struct grid *g){ exit(1); } - for(i=0;incell;i++){ + for(ppi=0;ppincell); distCalc(par, g); smooth(par,g); diff --git a/src/lime.h b/src/lime.h index 86bb732..632cf76 100644 --- a/src/lime.h +++ b/src/lime.h @@ -260,7 +260,7 @@ void lineCount(int,molData *,int **, int **, int *); void LTE(inputPars *, struct grid *, molData *); void molinit(molData *, inputPars *, struct grid *,int); void openSocket(inputPars *par, int); -void delaunay(inputPars *, struct grid *); +void delaunay(const int, struct grid *, const unsigned long); void photon(int, struct grid *, molData *, int, const gsl_rng *,inputPars *,blend *,gridPointData *,double *); void parseInput(inputPars *, image **, molData **); double planckfunc(int, double, molData *, int); diff --git a/src/popsin.c b/src/popsin.c index 752cc25..90e209d 100644 --- a/src/popsin.c +++ b/src/popsin.c @@ -105,7 +105,7 @@ popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){ } fclose(fp); - delaunay(par, *g); + delaunay(DIM, *g, (unsigned long)par->ncell); distCalc(par, *g); getVelosplines(par,*g); *popsdone=1; diff --git a/src/predefgrid.c b/src/predefgrid.c index deec8b0..1d30b8d 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -69,7 +69,7 @@ predefinedGrid(inputPars *par, struct grid *g){ } fclose(fp); - delaunay(par,g); + delaunay(DIM, g, (unsigned long)par->ncell); distCalc(par,g); // getArea(par,g, ran); // getMass(par,g, ran); diff --git a/src/smooth.c b/src/smooth.c index c84acd8..4b1383f 100644 --- a/src/smooth.c +++ b/src/smooth.c @@ -68,7 +68,7 @@ smooth(inputPars *par, struct grid *g){ } } - delaunay(par, g); + delaunay(DIM, g, (unsigned long)par->ncell); distCalc(par, g); if(!silent) progressbar((double)(sg+1)/(double)smooth, 5); } From 06529af54bfdab90ac54c39bdfa49beae8e1cf5f Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 12:31:10 +0200 Subject: [PATCH 11/43] Rearranged if block in delaunay() The block starting at line 168 in grid.c is rearranged, the if test being inverted such that the code exits if qh_new_qhull() returns a non-zero status. This involves unindenting all the code (~32 lines) which was in the block before. No other functional changes have been made to these lines in the present commit. --- src/grid.c | 56 +++++++++++++++++++++++++++--------------------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/src/grid.c b/src/grid.c index 5fe10a3..336fd12 100644 --- a/src/grid.c +++ b/src/grid.c @@ -165,42 +165,42 @@ delaunay(const int numDims, struct grid *g, const unsigned long numPoints){ } sprintf(flags,"qhull d Qbb"); - if (!qh_new_qhull(DIM, (int)numPoints, pt_array, ismalloc, flags, NULL, NULL)) { - /* Identify points */ - FORALLvertices { - id=qh_pointid(vertex->point); - g[id].numNeigh=qh_setsize(vertex->neighbors); - if( g[id].neigh != NULL ) - { - free( g[id].neigh ); - } - g[id].neigh=malloc(sizeof(struct grid *)*g[id].numNeigh); - for(k=0;kpoint); + g[id].numNeigh=qh_setsize(vertex->neighbors); + if( g[id].neigh != NULL ) + { + free( g[id].neigh ); } + g[id].neigh=malloc(sizeof(struct grid *)*g[id].numNeigh); + for(k=0;kupperdelaunay) { - j=0; - FOREACHvertex_ (facet->vertices) simplex[j++]=qh_pointid(vertex->point); - for(i=0;iid != g[simplex[j]].id) { - k++; - } - g[simplex[i]].neigh[k]=&g[simplex[j]]; + /* Identify neighbors */ + FORALLfacets { + if (!facet->upperdelaunay) { + j=0; + FOREACHvertex_ (facet->vertices) simplex[j++]=qh_pointid(vertex->point); + for(i=0;iid != g[simplex[j]].id) { + k++; } + g[simplex[i]].neigh[k]=&g[simplex[j]]; } } } } - } else { - if(!silent) bail_out("Qhull failed to triangulate"); - exit(1); } for(ppi=0;ppi Date: Fri, 8 Apr 2016 12:37:31 +0200 Subject: [PATCH 12/43] Renamed simplex in delaunay() The array 'simplex' is renamed to 'pointIdsThisFacet' for added clarity. Additional variables idI and idJ have been added. All three of these have been made type unsigned long, and variable 'id' is also now of the same type. --- src/grid.c | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/grid.c b/src/grid.c index 336fd12..b0f4d6b 100644 --- a/src/grid.c +++ b/src/grid.c @@ -147,15 +147,14 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ void delaunay(const int numDims, struct grid *g, const unsigned long numPoints){ - int i,j,k,id; + int i,j,k; char flags[255]; boolT ismalloc = False; facetT *facet; vertexT *vertex,**vertexp; coordT *pt_array; - int simplex[DIM+1]; int curlong, totlong; - unsigned long ppi; + unsigned long ppi,id,pointIdsThisFacet[numDims+1],idI,idJ; pt_array=malloc(sizeof(coordT)*numDims*numPoints); for(ppi=0;ppiupperdelaunay) { j=0; - FOREACHvertex_ (facet->vertices) simplex[j++]=qh_pointid(vertex->point); - for(i=0;ivertices) pointIdsThisFacet[j++]=(unsigned long)qh_pointid(vertex->point); + + for(i=0;iid != g[simplex[j]].id) { + /* Cycle through all the non-NULL links of g[idI], storing the link if it is new. + */ + k=0; + while(g[idI].neigh[k] != NULL && g[idI].neigh[k]->id != g[idJ].id) k++; - } - g[simplex[i]].neigh[k]=&g[simplex[j]]; + g[idI].neigh[k]=&g[idJ]; } } } From 9cceb29b2e7d46ef04f66dd2d421a5d23fee506e Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 12:46:11 +0200 Subject: [PATCH 13/43] Added the rest of delaunay arguments Three more arguments, to mediate the return of the cell information. These are not yet used, but the changed interface required some additions to every function which calls delaunay(). --- src/grid.c | 8 ++++++-- src/lime.h | 2 +- src/popsin.c | 6 +++++- src/predefgrid.c | 6 +++++- src/smooth.c | 7 +++++-- 5 files changed, 22 insertions(+), 7 deletions(-) diff --git a/src/grid.c b/src/grid.c index b0f4d6b..abfb386 100644 --- a/src/grid.c +++ b/src/grid.c @@ -146,7 +146,8 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ } void -delaunay(const int numDims, struct grid *g, const unsigned long numPoints){ +delaunay(const int numDims, struct grid *g, const unsigned long numPoints\ + , const _Bool getCells, struct cell **dc, unsigned long *numCells){ int i,j,k; char flags[255]; boolT ismalloc = False; @@ -518,6 +519,8 @@ buildGrid(inputPars *par, struct grid *g){ double temp; int k=0,i; /* counters */ int flag; + struct cell *dc=NULL; /* Not used at present. */ + unsigned long numCells; gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2); /* Random number generator */ #ifdef TEST @@ -609,7 +612,7 @@ buildGrid(inputPars *par, struct grid *g){ } /* end grid allocation */ - delaunay(DIM, g, (unsigned long)par->ncell); + delaunay(DIM, g, (unsigned long)par->ncell, 0, &dc, &numCells); distCalc(par, g); smooth(par,g); @@ -624,6 +627,7 @@ buildGrid(inputPars *par, struct grid *g){ // getMass(par,g, ran); getVelosplines(par,g); dumpGrid(par,g); + if(dc!=NULL) free(dc); gsl_rng_free(ran); if(!silent) done(5); diff --git a/src/lime.h b/src/lime.h index 632cf76..e775d03 100644 --- a/src/lime.h +++ b/src/lime.h @@ -260,7 +260,7 @@ void lineCount(int,molData *,int **, int **, int *); void LTE(inputPars *, struct grid *, molData *); void molinit(molData *, inputPars *, struct grid *,int); void openSocket(inputPars *par, int); -void delaunay(const int, struct grid *, const unsigned long); +void delaunay(const int, struct grid *, const unsigned long, const _Bool, struct cell **, unsigned long *); void photon(int, struct grid *, molData *, int, const gsl_rng *,inputPars *,blend *,gridPointData *,double *); void parseInput(inputPars *, image **, molData **); double planckfunc(int, double, molData *, int); diff --git a/src/popsin.c b/src/popsin.c index 90e209d..2be88d4 100644 --- a/src/popsin.c +++ b/src/popsin.c @@ -17,6 +17,8 @@ popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){ FILE *fp; int i,j,k; double dummy; + struct cell *dc=NULL; /* Not used at present. */ + unsigned long numCells; if((fp=fopen(par->restart, "rb"))==NULL){ if(!silent) bail_out("Error reading binary output populations file!"); @@ -105,9 +107,11 @@ popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){ } fclose(fp); - delaunay(DIM, *g, (unsigned long)par->ncell); + delaunay(DIM, *g, (unsigned long)par->ncell, 0, &dc, &numCells); distCalc(par, *g); getVelosplines(par,*g); *popsdone=1; + + if(dc!=NULL) free(dc); } diff --git a/src/predefgrid.c b/src/predefgrid.c index 1d30b8d..f278891 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -14,6 +14,9 @@ predefinedGrid(inputPars *par, struct grid *g){ FILE *fp; int i; double x,y,z,scale; + struct cell *dc=NULL; /* Not used at present. */ + unsigned long numCells; + gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2); #ifdef TEST gsl_rng_set(ran,6611304); @@ -69,11 +72,12 @@ predefinedGrid(inputPars *par, struct grid *g){ } fclose(fp); - delaunay(DIM, g, (unsigned long)par->ncell); + delaunay(DIM, g, (unsigned long)par->ncell, 0, &dc, &numCells); distCalc(par,g); // getArea(par,g, ran); // getMass(par,g, ran); getVelosplines_lin(par,g); if(par->gridfile) write_VTK_unstructured_Points(par, g); gsl_rng_free(ran); + if(dc!=NULL) free(dc); } diff --git a/src/smooth.c b/src/smooth.c index 4b1383f..9197985 100644 --- a/src/smooth.c +++ b/src/smooth.c @@ -20,7 +20,9 @@ smooth(inputPars *par, struct grid *g){ int smooth=20; /* Amount of grid smoothing */ double move[3]; /* Auxillary array for smoothing the grid */ double dist; /* Distance to a neighbor */ - + struct cell *dc=NULL; /* Not used at present. */ + unsigned long numCells; + for(sg=0;sgncell && !g[i].sink;i++){ mindist=1e30; @@ -68,9 +70,10 @@ smooth(inputPars *par, struct grid *g){ } } - delaunay(DIM, g, (unsigned long)par->ncell); + delaunay(DIM, g, (unsigned long)par->ncell, 0, &dc, &numCells); distCalc(par, g); if(!silent) progressbar((double)(sg+1)/(double)smooth, 5); + if(dc!=NULL) free(dc); } } From f0a5ba53ddbc9da1503a63c705fa6210d9339d87 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 13:17:17 +0200 Subject: [PATCH 14/43] Added cell return code to delaunay() This is not yet called anywhere however. --- src/grid.c | 56 ++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 2 deletions(-) diff --git a/src/grid.c b/src/grid.c index abfb386..b5c3959 100644 --- a/src/grid.c +++ b/src/grid.c @@ -151,11 +151,12 @@ delaunay(const int numDims, struct grid *g, const unsigned long numPoints\ int i,j,k; char flags[255]; boolT ismalloc = False; - facetT *facet; + facetT *facet, *neighbor, **neighborp; vertexT *vertex,**vertexp; coordT *pt_array; int curlong, totlong; - unsigned long ppi,id,pointIdsThisFacet[numDims+1],idI,idJ; + unsigned long ppi,id,pointIdsThisFacet[numDims+1],idI,idJ,fi,ffi; + _Bool neighbourNotFound; pt_array=malloc(sizeof(coordT)*numDims*numPoints); for(ppi=0;ppiupperdelaunay) { j=0; @@ -204,6 +206,7 @@ delaunay(const int numDims, struct grid *g, const unsigned long numPoints\ } } } + (*numCells)++; } } @@ -215,6 +218,55 @@ delaunay(const int numDims, struct grid *g, const unsigned long numPoints\ } g[ppi].numNeigh=j; } + + if(getCells){ + (*dc) = malloc(sizeof(**dc)*(*numCells)); + fi = 0; + FORALLfacets { + if (!facet->upperdelaunay) { + (*dc)[fi].id = (unsigned long)facet->id; /* Do NOT expect this to be equal to fi. */ + fi++; + } + } + + fi = 0; + FORALLfacets { + if (!facet->upperdelaunay) { + i = 0; + FOREACHneighbor_(facet) { + if(neighbor->upperdelaunay){ + (*dc)[fi].neigh[i] = NULL; + }else{ + /* Have to find the member of *dc with the same id as neighbour.*/ + ffi = 0; + neighbourNotFound=1; + while(ffi<(*numCells) && neighbourNotFound){ + if((*dc)[ffi].id==(unsigned long)neighbor->id){ + (*dc)[fi].neigh[i] = &(*dc)[ffi]; + neighbourNotFound = 0; + } + ffi++; + } + if(ffi>=(*numCells) && neighbourNotFound){ + if(!silent) bail_out("Something weird going on."); + exit(1); + } + } + i++; + } + + i = 0; + FOREACHvertex_( facet->vertices ) { + id = (unsigned long)qh_pointid(vertex->point); + (*dc)[fi].vertx[i] = &g[id]; + i++; + } + + fi++; + } + } + } + qh_freeqhull(!qh_ALL); qh_memfreeshort (&curlong, &totlong); free(pt_array); From 8b8a364455aa2e69a79475e01479c31ea93ed96a Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 13:19:24 +0200 Subject: [PATCH 15/43] Changed g to gp in delaunay() This is a purely cosmetic change, but a useful one, because searching for g returns too many false positives. I plan to do this elsewhere too. --- src/grid.c | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/grid.c b/src/grid.c index b5c3959..ac602ff 100644 --- a/src/grid.c +++ b/src/grid.c @@ -146,7 +146,7 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ } void -delaunay(const int numDims, struct grid *g, const unsigned long numPoints\ +delaunay(const int numDims, struct grid *gp, const unsigned long numPoints\ , const _Bool getCells, struct cell **dc, unsigned long *numCells){ int i,j,k; char flags[255]; @@ -161,7 +161,7 @@ delaunay(const int numDims, struct grid *g, const unsigned long numPoints\ pt_array=malloc(sizeof(coordT)*numDims*numPoints); for(ppi=0;ppipoint); - g[id].numNeigh=qh_setsize(vertex->neighbors); - if( g[id].neigh != NULL ) + gp[id].numNeigh=qh_setsize(vertex->neighbors); + if( gp[id].neigh != NULL ) { - free( g[id].neigh ); + free( gp[id].neigh ); } - g[id].neigh=malloc(sizeof(struct grid *)*g[id].numNeigh); - for(k=0;kid != g[idJ].id) + while(gp[idI].neigh[k] != NULL && gp[idI].neigh[k]->id != gp[idJ].id) k++; - g[idI].neigh[k]=&g[idJ]; + gp[idI].neigh[k]=&gp[idJ]; } } } @@ -212,11 +212,11 @@ delaunay(const int numDims, struct grid *g, const unsigned long numPoints\ for(ppi=0;ppivertices ) { id = (unsigned long)qh_pointid(vertex->point); - (*dc)[fi].vertx[i] = &g[id]; + (*dc)[fi].vertx[i] = &gp[id]; i++; } From bf9cb4d93a79e399f444d3a4ca0272f2e7ab79cd Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 13:27:41 +0200 Subject: [PATCH 16/43] Edited comments in delaunay() I also rearranged the variable declarations a little. The only substantive change is the explicit cast of the return value of qh_pointid() to unsigned long. This is the last of the delaunay-focussed commits. --- src/grid.c | 50 ++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 38 insertions(+), 12 deletions(-) diff --git a/src/grid.c b/src/grid.c index ac602ff..2f15c8e 100644 --- a/src/grid.c +++ b/src/grid.c @@ -145,17 +145,38 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ } } -void -delaunay(const int numDims, struct grid *gp, const unsigned long numPoints\ +/*....................................................................*/ +void delaunay(const int numDims, struct grid *gp, const unsigned long numPoints\ , const _Bool getCells, struct cell **dc, unsigned long *numCells){ + /* +The principal purpose of this function is to perform a Delaunay triangulation for the set of points defined by the input argument 'g'. This is achieved via routines in the 3rd-party package qhull. + +A note about qhull nomenclature: a vertex is what you think it is - i.e., a point; but a facet means in this context a Delaunay triangle (in 2D) or tetrahedron (in 3D). + +Required elements of structs: + struct grid *gp: + .id + .x + +Elements of structs are set as follows: + struct grid *gp: + .numNeigh + .neigh (this is malloc'd too large and at present not realloc'd.) + + cellType *dc (if getCells>0): + .id + .neigh + .vertx + */ + + coordT *pt_array=NULL; + unsigned long ppi,id,pointIdsThisFacet[numDims+1],idI,idJ,fi,ffi; int i,j,k; char flags[255]; boolT ismalloc = False; - facetT *facet, *neighbor, **neighborp; vertexT *vertex,**vertexp; - coordT *pt_array; + facetT *facet, *neighbor, **neighborp; int curlong, totlong; - unsigned long ppi,id,pointIdsThisFacet[numDims+1],idI,idJ,fi,ffi; _Bool neighbourNotFound; pt_array=malloc(sizeof(coordT)*numDims*numPoints); @@ -173,22 +194,27 @@ delaunay(const int numDims, struct grid *gp, const unsigned long numPoints\ /* Identify points */ FORALLvertices { - id=qh_pointid(vertex->point); + id=(unsigned long)qh_pointid(vertex->point); + /* Note this is NOT the same value as vertex->id. Only the id gained via the call to qh_pointid() is the same as the index of the point in the input list. */ + gp[id].numNeigh=qh_setsize(vertex->neighbors); - if( gp[id].neigh != NULL ) - { - free( gp[id].neigh ); - } + /* Note that vertex->neighbors refers to facets abutting the vertex, not other vertices. In general there seem to be more facets surrounding a point than vertices (in fact there seem to be exactly 2x as many). In any case, mallocing to N_facets gives extra room. */ + + if(gp[id].neigh!=NULL) + free( gp[id].neigh ); gp[id].neigh=malloc(sizeof(struct grid *)*gp[id].numNeigh); for(k=0;kupperdelaunay) { + /* Store the point IDs in a list for convenience. These ID values are conveniently ordered such that qh_pointid() returns ppi for gp[ppi]. + */ j=0; FOREACHvertex_ (facet->vertices) pointIdsThisFacet[j++]=(unsigned long)qh_pointid(vertex->point); From ffdf517979038693526c8bab8ff79a80b4b9e2c1 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 13:39:18 +0200 Subject: [PATCH 17/43] velocityspline to calcLineAmp(etc). To assist the interpolative raytracing algorithm I plan to introduce a function calcLineAmpInterp() to return the spectral line amplitude relative to line centre, a quantity which in LIME is usually signified by the variable named 'vfac'. This function is therefore analogous to velocityspline2() in raytrace.c and the similarly-named functions in photon.c. This planned new addition is a good excuse to change the names and interfaces of all three existing 'velocityspline' functions to be more similar. In the present commit I have changed the names as follows: velocityspline to calcLineAmpSpline velocityspline_lin to calcLineAmpLinear velocityspline2 to calcLineAmpSample I also made some of the arguments 'const'. This has no effect on the functioning, it is just more conservative practice. The final change is the addition of a calcLineAmpLinear() template to lime.h --- src/lime.h | 5 +++-- src/photon.c | 24 ++++++++++++++---------- src/raytrace.c | 7 ++++--- 3 files changed, 21 insertions(+), 15 deletions(-) diff --git a/src/lime.h b/src/lime.h index e775d03..85f615f 100644 --- a/src/lime.h +++ b/src/lime.h @@ -281,8 +281,9 @@ void stateq(int, struct grid *, molData *, int, inputPars *,gridPointData *,d void statistics(int, molData *, struct grid *, int *, double *, double *, int *); void stokesangles(double, double, double, double, double *); void traceray(rayData, int, int, inputPars *, 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*); +void calcLineAmpLinear(struct grid*, const int, const int, const double, const double, double*); +void calcLineAmpSample(double*, double*, const double, const double, const double, double*); +void calcLineAmpSpline(struct grid*, const int, const int, const double, const double, double*); double veloproject(double *, double *); void writefits(int, inputPars *, molData *, image *); void write_VTK_unstructured_Points(inputPars *, struct grid *); diff --git a/src/photon.c b/src/photon.c index 4533870..76f18e0 100644 --- a/src/photon.c +++ b/src/photon.c @@ -49,8 +49,10 @@ sortangles(double *inidir, int id, struct grid *g, const gsl_rng *ran) { -void -velocityspline(struct grid *g, int id, int k, double binv, double deltav, double *vfac){ +/*....................................................................*/ +void calcLineAmpSpline(struct grid *g, const int id, const int k\ + , const double binv, const double deltav, double *vfac){ + int nspline,ispline,naver,iaver; double v1,v2,s1,s2,sd,v,vfacsub,d; @@ -82,8 +84,10 @@ velocityspline(struct grid *g, int id, int k, double binv, double deltav, double } -void -velocityspline_lin(struct grid *g, int id, int k, double binv, double deltav, double *vfac){ +/*....................................................................*/ +void calcLineAmpLinear(struct grid *g, const int id, const int k\ + , const double binv, const double deltav, double *vfac){ + int nspline,ispline,naver,iaver; double v1,v2,s1,s2,sd,v,vfacsub,d; @@ -217,8 +221,8 @@ 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]); - else velocityspline_lin(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); + if(!par->doPregrid) calcLineAmpSpline(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); + else calcLineAmpLinear(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); mp[l].vfac[iphot]=vfac[0]; } for(l=0;l<3;l++) x[l]=g[here].x[l]+(g[here].dir[dir].xn[l] * g[id].ds[dir]/2.); @@ -228,8 +232,8 @@ 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]); - else velocityspline_lin(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); + if(!par->doPregrid) calcLineAmpSpline(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); + else calcLineAmpLinear(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); } for(iline=0;ilinedoPregrid) 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); + if(!par->doPregrid) calcLineAmpSpline(g,here,dir,g[id].mol[counta[jline]].binv,deltav-matrix[jline].deltav,&vblend); + else calcLineAmpLinear(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; if(dtau < -30) dtau = -30; diff --git a/src/raytrace.c b/src/raytrace.c index 8f28320..ecd2788 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -10,8 +10,9 @@ #include "lime.h" -void -velocityspline2(double x[3], double dx[3], double ds, double binv, double deltav, double *vfac){ +/*....................................................................*/ +void calcLineAmpSample(double x[3], double dx[3], const double ds\ + , const double binv, const double deltav, double *vfac){ /* The bulk velocity of the model material can vary significantly with position, thus so can the value of the line-shape function at a given frequency and direction. The present function calculates 'vfac', an approximate average of the line-shape function along a path of length ds in the direction of the line of sight. */ @@ -144,7 +145,7 @@ Note that the algorithm employed here is similar to that employed in the functio /* Line centre occurs when deltav = the recession velocity of the radiating material. Explanation of the signs of the 2nd and 3rd terms on the RHS: (i) A bulk source velocity (which is defined as >0 for the receding direction) should be added to the material velocity field; this is equivalent to subtracting it from deltav, as here. (ii) A positive value of lineRedShift means the line is red-shifted wrt to the frequency specified for the image. The effect is the same as if the line and image frequencies were the same, but the bulk recession velocity were higher. lineRedShift should thus be added to the recession velocity, which is equivalent to subtracting it from deltav, as here. */ /* Calculate an approximate average line-shape function at deltav within the Voronoi cell. */ - if(!par->pregrid) velocityspline2(x,dx,ds,g[posn].mol[molI].binv,deltav,&vfac); + if(!par->pregrid) calcLineAmpSample(x,dx,ds,g[posn].mol[molI].binv,deltav,&vfac); else vfac=gaussline(deltav+veloproject(dx,g[posn].vel),g[posn].mol[molI].binv); /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. */ From 9e3006c0a4599ea270ac7759ed1cdc8c2d40fc7f Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 15:05:41 +0200 Subject: [PATCH 18/43] New module raythrucells.c This contains the code for finding the cells traversed by a ray on the one hand and interpolating within Delaunay cells on the other. None of this code is yet used, but it compiles at this stage. --- Makefile | 4 +- src/lime.h | 9 + src/raythrucells.c | 529 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 540 insertions(+), 2 deletions(-) create mode 100644 src/raythrucells.c diff --git a/Makefile b/Makefile index 10b4ee0..9fc6f9a 100644 --- a/Makefile +++ b/Makefile @@ -56,7 +56,7 @@ SRCS = src/aux.c src/messages.c src/grid.c src/LTEsolution.c \ src/raytrace.c src/smooth.c src/sourcefunc.c \ src/stateq.c src/statistics.c src/magfieldfit.c \ src/stokesangles.c src/writefits.c src/weights.c \ - src/velospline.c src/getclosest.c \ + src/velospline.c src/getclosest.c src/raythrucells.c \ src/tcpsocket.c src/defaults.c src/fastexp.c MODELS = model.c OBJS = src/aux.o src/messages.o src/grid.o src/LTEsolution.o \ @@ -65,7 +65,7 @@ OBJS = src/aux.o src/messages.o src/grid.o src/LTEsolution.o \ src/ratranInput.o src/smooth.o src/sourcefunc.o \ src/stateq.o src/statistics.o src/magfieldfit.o \ src/stokesangles.o src/writefits.o src/weights.o \ - src/velospline.o src/getclosest.o \ + src/velospline.o src/getclosest.o src/raythrucells.o \ src/tcpsocket.o src/defaults.o src/fastexp.o MODELO = src/model.o diff --git a/src/lime.h b/src/lime.h index 85f615f..f125918 100644 --- a/src/lime.h +++ b/src/lime.h @@ -218,6 +218,15 @@ typedef struct { } triangle2D; +int followRayThroughDelCells(double*, double*, struct grid*, struct cell*, const unsigned long, const double, intersectType*, unsigned long**, intersectType**, int*); +int buildRayCellChain(double*, double*, struct grid*, struct cell*, _Bool**, unsigned long, int, int, int, const double, unsigned long**, intersectType**, int*); +int getNewEntryFaceI(const unsigned long, const struct cell); +faceType extractFace(struct grid*, struct cell*, const unsigned long, const int); +void intersectLineTriangle(double*, double*, faceType, intersectType*); +triangle2D calcTriangle2D(faceType face); +void doBaryInterp(const intersectType, struct grid*, struct gAuxType*, double*, unsigned long*, molData*, const int, gridInterp*); +void doSegmentInterp(gridInterp*, const int, molData*, const int, const double, const int); + /* Some functions */ void density(double,double,double,double *); diff --git a/src/raythrucells.c b/src/raythrucells.c new file mode 100644 index 0000000..2855e6f --- /dev/null +++ b/src/raythrucells.c @@ -0,0 +1,529 @@ +/* + * raythrucells.c + * This file is part of LIME, the versatile line modeling engine + * + * Copyright (C) 2006-2014 Christian Brinch + * Copyright (C) 2015 The LIME development team + * + */ + +#include "lime.h" + +/*....................................................................*/ +int followRayThroughDelCells(double x[DIM], double dir[DIM], struct grid *gp\ + , struct cell *dc, const unsigned long numCells, const double epsilon\ + , intersectType *entryIntcpt, unsigned long **chainOfCellIds\ + , intersectType **cellExitIntcpts, int *lenChainPtrs){ + /* +The present function follows a ray through a connected, convex set of Delaunay cells and returns information about the chain of cells it passes through. If the ray is found to pass through 1 or more cells, the function returns 0, indicating success; if not, it returns a non-zero value. The chain description consists of three pieces of information: (i) intercept information for the entry face of the first cell encountered; (ii) the IDs of the cells in the chain; (iii) intercept information for the exit face of the ith cell. + */ + + const int numFaces=DIM+1, maxNumEntryFaces=100; + int numEntryFaces, fi, entryFis[maxNumEntryFaces], i, status; + faceType face; + unsigned long dci, entryDcis[maxNumEntryFaces];//, trialDci; + intersectType intcpt, entryIntcpts[maxNumEntryFaces]; + _Bool *cellVisited=NULL; + + /* Choose a set of starting faces by testing all the 'external' faces of cells which have some. */ + numEntryFaces = 0; + for(dci=0;dci0.0){ + if(numEntryFaces>maxNumEntryFaces){ + if(!silent) bail_out("Too many entry faces."); + exit(1); + } + + entryDcis[ numEntryFaces] = dci; + entryFis[ numEntryFaces] = fi; + entryIntcpts[numEntryFaces] = intcpt; + numEntryFaces++; + } + } + } + } + } + + if(numEntryFaces<=0) + return 1; + + *lenChainPtrs = 1024; /* This can be increased within followCellChain(). */ + *chainOfCellIds = malloc(sizeof(**chainOfCellIds) *(*lenChainPtrs)); + *cellExitIntcpts = malloc(sizeof(**cellExitIntcpts)*(*lenChainPtrs)); + cellVisited = malloc(sizeof(*cellVisited)*numCells); + for(dci=0;dci0){ + status = buildRayCellChain(x, dir, gp, dc, &cellVisited, entryDcis[i], entryFis[i], 0, 0, epsilon, chainOfCellIds, cellExitIntcpts, lenChainPtrs); + i++; + } + + if(status==0) + *entryIntcpt = entryIntcpts[i-1]; + /* Note that the order of the bary coords, and the value of fi, are with reference to the vertx list of the _entered_ cell. This can't of course be any other way, because this ray enters this first cell from the exterior of the model, where there are no cells. For all the intersectType objects in the list cellExitIntcpts, the bary coords etc are with reference to the exited cell. */ + + free(cellVisited); + + return status; +} + +/*....................................................................*/ +int buildRayCellChain(double x[DIM], double dir[DIM], struct grid *gp\ + , struct cell *dc, _Bool **cellVisited, unsigned long dci, int entryFaceI\ + , int levelI, int nCellsInChain, const double epsilon\ + , unsigned long **chainOfCellIds, intersectType **cellExitIntcpts, int *lenChainPtrs){ + /* +This function is designed to follow a ray (defined by a starting locus 'x' and a direction vector 'dir') through a convex connected set of Delaunay cells. The function returns an integer status value directly, and two lists (plus their common length) via the argument interface: chainOfCellIds and cellExitIntcpts. Taken together, these lists define a chain of Delaunay cells traversed by the ray. + +The task of determining which cells are traversed by the ray is simple in principle, but complications arise in computational practice due to the finite precision of floating-point calculations. Where the ray crosses a cell face near to one of its edges, numerical calculation of the 'impact parameter' may return an answer which is erroneous either way: i.e., a ray which just misses a face may be reported as hitting it, and vice versa. To deal with this, a distinction is made between impacts which are (i) 'good', that is far from any face edge; (ii) 'bad', i.e. clearly missing the face; and (iii) 'marginal', i.e. closer to the edge than some preset cutoff which is represented in the argument list by the number 'epsilon'. Note that a tally is kept of those cells which have already been tested for the current ray, and any face which abuts a neighbouring cell which has been visited already will be flagged as 'bad'. + +The function therefore looks at all the exit faces of the cell and sorts them into these three categories. How it proceeds then depends on the relative numbers of each, as described below. + + - If there is more than 1 good exit face, an exception is generated. This is a sign that epsilon has been chosen with too small a value. + + - If there is just 1 good exit face, the marginal ones are ignored. The current cell data are appended to the chain and the cell abutting the exit face becomes the new working cell. + + - If there are no good exit faces, we look at the marginal ones. If there is only 1 of these, we loop as before. If there are more than 1, the function is called recursively for each cell on the far side of an exit face. + +Thus there are two alternate modes of operation within the function: a straightforward loop along the cells in a single chain, which will continue so long as there is only a single exit face of type either 'good' or 'marginal'; and a recursive launch of the function at a fork in the chain into each of its several possible branches. + +The function terminates under the following conditions: + - It detects that the edge of the model is reached (returns success). + - There are no exit faces (returns failure). + - There are no good exit faces and either + * all of the recursive calls to marginal faces have been unsuccessful (returns failure), or + * one of these has been successful (returns success). + +At a successful termination, therefore, details of all the cells to the edge of the model are correctly stored in chainOfCellIds and cellExitIntcpts, and the number of these cells is returned in lenChainPtrs. + +***** Note that it is assumed here that a mis-indentification of the actual cell traversed by a ray will not ultimately matter to the external calling routine. This is only reasonable if whatever function or property is being sampled by the ray does not vary in a step function across cell boundaries. ***** + */ + + const int numFaces=DIM+1; + _Bool followingSingleChain; + const int bufferSize=1024; + int numGoodExits, numMarginalExits, fi, goodExitFis[numFaces], marginalExitFis[numFaces], exitFi, i, status, newEntryFaceI; + faceType face; + intersectType intcpt[numFaces], dummy_intcpt; + + followingSingleChain = 1; /* default */ + do{ /* Follow the chain through 'good' cells, i.e. ones for which entry and exit are nicely distant from face edges. (Here we also follow marginal exits if there are no good ones, and only 1 marginal one.) */ + (*cellVisited)[dci] = 1; + + /* Store the current cell ID (we leave storing the exit face for later, when we know what it is). (If there is not enough room in chainOfCellIds and cellExitIntcpts, realloc them to new value of lenChainPtrs.) */ + if(nCellsInChain>=(*lenChainPtrs)){ + *lenChainPtrs += bufferSize; + *chainOfCellIds = realloc(*chainOfCellIds, sizeof(**chainOfCellIds) *(*lenChainPtrs)); + *cellExitIntcpts = realloc(*cellExitIntcpts, sizeof(**cellExitIntcpts)*(*lenChainPtrs)); + } + (*chainOfCellIds)[nCellsInChain] = dci; + + /* calculate num good and bad exits */ + numGoodExits = 0; + numMarginalExits = 0; + for(fi=0;fiid])){ + /* Store points for this face: */ + face = extractFace(gp, dc, dci, fi); + + /* Now calculate the intercept: */ + intersectLineTriangle(x, dir, face, &intcpt[fi]); + intcpt[fi].fi = fi; /* Ultimately we need this so we can relate the bary coords for the face back to the Delaunay cell. */ + + if(intcpt[fi].orientation>0){ /* it is an exit face. */ + if(intcpt[fi].collPar-epsilon>0.0){ + goodExitFis[numGoodExits] = fi; + numGoodExits++; + }else if (intcpt[fi].collPar+epsilon>0.0){ + marginalExitFis[numMarginalExits] = fi; + numMarginalExits++; + } + } + } + } + + if(numGoodExits>1){ + if(!silent) bail_out("Some sort of bug: more than 1 firm candidate found for ray exit from cell."); + exit(1); + + }else if(numGoodExits==1 || numMarginalExits==1){ + if(numGoodExits==1) + exitFi = goodExitFis[0]; + else + exitFi = marginalExitFis[0]; + + /* Store the exit face details: */ + (*cellExitIntcpts)[nCellsInChain] = intcpt[exitFi]; + + nCellsInChain++; + + if(dc[dci].neigh[exitFi]==NULL){ /* Signals that we have reached the edge of the model. */ + /* Realloc the ptrs to their final sizes: */ + *chainOfCellIds = realloc(*chainOfCellIds, sizeof(**chainOfCellIds) *nCellsInChain); + *cellExitIntcpts = realloc(*cellExitIntcpts, sizeof(**cellExitIntcpts)*nCellsInChain); + *lenChainPtrs = nCellsInChain; + + return 0; + } + + entryFaceI = getNewEntryFaceI(dci, *(dc[dci].neigh[exitFi])); + dci = dc[dci].neigh[exitFi]->id; + + } else + followingSingleChain = 0; + } while(followingSingleChain); + + /* Now we have run out of good (or at least single) exit-face options, let's try the marginal ones. */ + + if(numMarginalExits<1) + return 1; /* Unsuccessful end of this chain. */ + + /* If we have got to this point, we must have numMarginalExits>1; thus we have a fork in the chain, and must explore each branch. We recurse here because a recursive scheme is the best way to do that. */ + + i = 0; + status = 1; /* default */ + while(i0){ + exitFi = marginalExitFis[i]; + (*cellExitIntcpts)[nCellsInChain] = intcpt[exitFi]; + + if(dc[dci].neigh[exitFi]==NULL){ + status = 0; + }else{ + newEntryFaceI = getNewEntryFaceI(dci, *(dc[dci].neigh[exitFi])); + + /* Now we dive into the branch: */ + status = buildRayCellChain(x, dir, gp, dc, cellVisited, dc[dci].neigh[exitFi]->id\ + , newEntryFaceI, levelI+1, nCellsInChain+1, epsilon\ + , chainOfCellIds, cellExitIntcpts, lenChainPtrs); + } + i++; + } + + return status; +} + +/*....................................................................*/ +int getNewEntryFaceI(const unsigned long dci, const struct cell newCell){ + /* Finds the index of the old cell in the face list of the new cell. */ + + const int numFaces=DIM+1; + _Bool matchFound = 0; + int ffi = 0, newEntryFaceI; + + while(ffiid==dci){ + matchFound = 1; + newEntryFaceI = ffi; + } + ffi++; + } + + /* Sanity check: */ + if(!matchFound){ + bail_out("Cannot find old cell ID in new cell data."); + exit(1); + } + + return newEntryFaceI; +} + +/*....................................................................*/ +faceType extractFace(struct grid *gp, struct cell *dc, const unsigned long dci, const int fi){ + /* Given a simplex dc[dci] and the face index (in the range {0...DIM}) fi, this returns the desired information about that face. Note that the ordering of the elements of face.r[] is the same as the ordering of the vertices of the simplex, dc[dci].vertx[]; just the vertex fi is omitted. + +Note that the element 'centre' of the faceType struct is mean to contain the spatial coordinates of the centre of the simplex, not of the face. This is designed to facilitate orientation of the face and thus to help determine whether rays which cross it are entering or exiting the simplex. + */ + + const int numFaces=DIM+1; + int vi, vvi, di; + struct grid point; + faceType face; + + unsigned long gi; + + vvi = 0; + for(vi=0;viid; + point = gp[gi]; + for(di=0;difi = -1; + + /* First calculate a normal vector to the triangle from the cross product. Since we don't know a priori whether the vertices of the face are listed CW or ACW seen from inside the cell, we will work this out by dotting the norm vector with a vector from the centre of the cell to vertex 0. + */ + for(di=0;di0.0){ /* it is an exit face. */ + intcpt->orientation = 1; + }else if(normDotDx<0.0){ /* it is an entry face. */ + intcpt->orientation = -1; + }else{ /* normDotDx==0.0, i.e. line and plane are parallel. */ + intcpt->orientation = 0; + intcpt->bary[0] = 0.0; + intcpt->bary[1] = 0.0; + intcpt->bary[2] = 0.0; + intcpt->dist = 0.0; + intcpt->collPar = 0.0; + + return; + } + + /* If we've got to here, we can be sure the line and plane are not parallel, and that we therefore expect meaningful results. */ + + numerator = norm[0]*(face.r[0][0] - x[0])\ + + norm[1]*(face.r[0][1] - x[1])\ + + norm[2]*(face.r[0][2] - x[2]); + + intcpt->dist = numerator/normDotDx; + + /* In order to calculate the barycentric coordinates, we need to set up a 2D coordinate basis in the plane of the triangle. I'll define the X axis as 10^, where for short I use the generic notation ji_ to mean the vector (tri[j]-tri[i]), and ji^ to mean the unit vector in the same direction. The Y axis is the unit vector in the direction (20_ - (20_.10^)*10^). + */ + face2D = calcTriangle2D(face); + + /* Now we want to express the intersection point in these coordinates: + */ + px2D[0] = (x[0] + intcpt->dist*dir[0] - face.r[0][0])*face2D.xAxis[0]\ + + (x[1] + intcpt->dist*dir[1] - face.r[0][1])*face2D.xAxis[1]\ + + (x[2] + intcpt->dist*dir[2] - face.r[0][2])*face2D.xAxis[2]; + px2D[1] = (x[0] + intcpt->dist*dir[0] - face.r[0][0])*face2D.yAxis[0]\ + + (x[1] + intcpt->dist*dir[1] - face.r[0][1])*face2D.yAxis[1]\ + + (x[2] + intcpt->dist*dir[2] - face.r[0][2])*face2D.yAxis[2]; + + /* +The barycentric coordinates (L0,L1,L2) are given by + + (tri2D[0][0]-tri2D[2][0] tri2D[1][0]-tri2D[2][0]) (L0) (px2D[0]-tri2D[2][0]) + ( )*( ) = ( ), + (tri2D[0][1]-tri2D[2][1] tri2D[1][1]-tri2D[2][1]) (L1) (px2D[1]-tri2D[2][1]) + +with L2 = 1 - L0 - L1. + */ + mat[0][0] = face2D.r[0][0]-face2D.r[2][0]; + mat[0][1] = face2D.r[1][0]-face2D.r[2][0]; + mat[1][0] = face2D.r[0][1]-face2D.r[2][1]; + mat[1][1] = face2D.r[1][1]-face2D.r[2][1]; + vec[0] = px2D[0]-face2D.r[2][0]; + vec[1] = px2D[1]-face2D.r[2][1]; + det = mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0]; + /*** We're assuming that the triangle is not pathological, i.e that det!=0. */ + intcpt->bary[0] = ( mat[1][1]*vec[0] - mat[0][1]*vec[1])/det; + intcpt->bary[1] = (-mat[1][0]*vec[0] + mat[0][0]*vec[1])/det; + intcpt->bary[2] = 1.0 - intcpt->bary[0] - intcpt->bary[1]; + + di = 0; + if(intcpt->bary[di] < 0.5) + intcpt->collPar = intcpt->bary[di]; + else + intcpt->collPar = 1.0 - intcpt->bary[di]; + + for(di=1;dibary[di] < 0.5){ + if(intcpt->bary[di] < intcpt->collPar) + intcpt->collPar = intcpt->bary[di]; + }else{ /* intcpt->bary[di]>=0.5 */ + if(1.0 - intcpt->bary[di] < intcpt->collPar) + intcpt->collPar = 1.0 - intcpt->bary[di]; + } + } + + return; +} + +/*....................................................................*/ +triangle2D calcTriangle2D(faceType face){ +/**** all this could be precalculated (norm and mat too). */ + + triangle2D face2D; + const int numDims = 3; + int i; + double lengthX, lengthY, dotResult; + + face2D.xAxis[0] = face.r[1][0]-face.r[0][0]; + face2D.xAxis[1] = face.r[1][1]-face.r[0][1]; + face2D.xAxis[2] = face.r[1][2]-face.r[0][2]; + lengthX = sqrt(face2D.xAxis[0]*face2D.xAxis[0]\ + + face2D.xAxis[1]*face2D.xAxis[1]\ + + face2D.xAxis[2]*face2D.xAxis[2]); + face2D.xAxis[0] /= lengthX; + face2D.xAxis[1] /= lengthX; + face2D.xAxis[2] /= lengthX; + + dotResult = face2D.xAxis[0]*(face.r[2][0]-face.r[0][0])\ + + face2D.xAxis[1]*(face.r[2][1]-face.r[0][1])\ + + face2D.xAxis[2]*(face.r[2][2]-face.r[0][2]); + + face2D.yAxis[0] = (face.r[2][0]-face.r[0][0]) - dotResult*face2D.xAxis[0]; + face2D.yAxis[1] = (face.r[2][1]-face.r[0][1]) - dotResult*face2D.xAxis[1]; + face2D.yAxis[2] = (face.r[2][2]-face.r[0][2]) - dotResult*face2D.xAxis[2]; + lengthY = sqrt(face2D.yAxis[0]*face2D.yAxis[0]\ + + face2D.yAxis[1]*face2D.yAxis[1]\ + + face2D.yAxis[2]*face2D.yAxis[2]); + face2D.yAxis[0] /= lengthY; + face2D.yAxis[1] /= lengthY; + face2D.yAxis[2] /= lengthY; + + /* The expressions for the triangle vertices in this basis are simple: + */ + face2D.r[0][0] = 0.0; + face2D.r[0][1] = 0.0; + face2D.r[1][0] = lengthX; + face2D.r[1][1] = 0.0; + face2D.r[2][0] = dotResult; + face2D.r[2][1] = lengthY; + + return face2D; +} + +/*....................................................................*/ +void doBaryInterp(const intersectType intcpt, struct grid *gp\ + , struct gAuxType *gAux, double xCmpntsRay[3], unsigned long gis[3]\ + , molData *md, const int numMols, gridInterp *gip){ + + int di, molI, levelI, lineI; + + (*gip).xCmpntRay = intcpt.bary[0]*xCmpntsRay[0]\ + + intcpt.bary[1]*xCmpntsRay[1]\ + + intcpt.bary[2]*xCmpntsRay[2]; + for(di=0;di Date: Fri, 8 Apr 2016 15:14:30 +0200 Subject: [PATCH 19/43] Added freePop2 and freeGAux --- src/grid.c | 34 ++++++++++++++++++++++++++++++++++ src/lime.h | 2 ++ 2 files changed, 36 insertions(+) diff --git a/src/grid.c b/src/grid.c index 2f15c8e..a4362d6 100644 --- a/src/grid.c +++ b/src/grid.c @@ -86,6 +86,40 @@ freePopulation(const inputPars *par, const molData* m, struct populations* pop ) free(pop); } } + +/*....................................................................*/ +void +freePop2(const int numSpecies, struct pop2 *mol){ + int i; + + if(mol !=NULL){ + for(i=0;i Date: Fri, 8 Apr 2016 15:43:18 +0200 Subject: [PATCH 20/43] Changed stokesangles This function is called by sourceFunc_pol() when the user parameter par->polarization is set. In the old version of stokesangles(), the user-supplied function magfield() was called for an arbitrary position in the model. In the planned interpolation raytracing scheme however it is possible to interpolate the B field between the vertices of any Delaunay tetrahedron. It is also desirable to remove the calls to user-supplied functions from any sections of code traversed in the case that either par.doPregrid or par.restart are set (although this is not yet achieved within the alterative raytrace code, as will be seen when this is added). For these reasons I re-wrote stokesangles() to accept a precalculated value of the B field as an input argument. These values are now also calculated via calls to magfield() at the grid cell loci within grid.c. At some point the storage of grid data needs to be improved, i.e. to include both velocity and B-field values. Until that time the par.doPregrid option will not work, not even after the present commit, even though it moves things in the right direction. --- src/grid.c | 11 +++++++++++ src/lime.h | 2 +- src/sourcefunc.c | 2 +- src/stokesangles.c | 12 +++++------- 4 files changed, 18 insertions(+), 9 deletions(-) diff --git a/src/grid.c b/src/grid.c index a4362d6..b1bed64 100644 --- a/src/grid.c +++ b/src/grid.c @@ -735,6 +735,17 @@ buildGrid(inputPars *par, struct grid *g){ abundance( g[i].x[0],g[i].x[1],g[i].x[2], g[i].abun); } + if(par->polarization){ + for(i=0;ipIntensity;i++) + magfield(g[i].x[0],g[i].x[1],g[i].x[2], g[i].B); + }else{ + for(i=0;ipIntensity;i++){ + g[i].B[0]=0.0; + g[i].B[1]=0.0; + g[i].B[2]=0.0; + } + } + // getArea(par,g, ran); // getMass(par,g, ran); getVelosplines(par,g); diff --git a/src/lime.h b/src/lime.h index fbe22ae..a864393 100644 --- a/src/lime.h +++ b/src/lime.h @@ -290,7 +290,7 @@ void sourceFunc_cont(double *,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 *); void statistics(int, molData *, struct grid *, int *, double *, double *, int *); -void stokesangles(double, double, double, double, double *); +void stokesangles(const double B[3], const double, double *); void traceray(rayData, int, int, inputPars *, struct grid *, molData *, image *, int, int *, int *, double); void calcLineAmpLinear(struct grid*, const int, const int, const double, const double, double*); void calcLineAmpSample(double*, double*, const double, const double, const double, double*); diff --git a/src/sourcefunc.c b/src/sourcefunc.c index 58a6223..b1d5f8c 100644 --- a/src/sourcefunc.c +++ b/src/sourcefunc.c @@ -75,7 +75,7 @@ sourceFunc_pol(double *snu, double *dtau, double ds, molData *m,double vfac,stru double dSigma, dSigma2, dI, dQ, dU, alpha; double angle[3]; - stokesangles(g[pos].x[0],g[pos].x[1],g[pos].x[2],incl,angle); + stokesangles(g[pos].B, incl, angle); /* Emission */ /* Continuum part: j_nu = rho_dust * kappa_nu */ diff --git a/src/stokesangles.c b/src/stokesangles.c index b211964..fb910ed 100644 --- a/src/stokesangles.c +++ b/src/stokesangles.c @@ -10,14 +10,12 @@ #include "lime.h" void -stokesangles(double x, double y, double z, double incl, double *angle){ +stokesangles(const double B[3], const double incl, double *trigFuncs){ - double B[3],Bp[3]; + double Bp[3]; double sPsi, cPsi; double cGam,cosIncl,sinIncl; - magfield(x,y,z,B); - cosIncl=cos(incl); sinIncl=sin(incl); Bp[0] = B[0]; @@ -42,8 +40,8 @@ stokesangles(double x, double y, double z, double incl, double *angle){ else { cGam = 0.; } - angle[0] = cPsi; //cosinus of Psi - angle[1] = sPsi; //sinus of Psi - angle[2] = cGam; //cosinus of Gamma + trigFuncs[0] = cPsi; //cosinus of Psi + trigFuncs[1] = sPsi; //sinus of Psi + trigFuncs[2] = cGam; //cosinus of Gamma } From 52b66b20f2f2799af4898abc7d451eb43e06b909 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 16:20:30 +0200 Subject: [PATCH 21/43] Changes to sourceFunc_pol() The changes are as follows: . Changes to the arguments, to make it easier to call this routine from the planned interpolated raytrace; the change from 'struct grid *g' to 'struct pop2 gm' also allows much more compact and clear expressions. . Variable name 'angles' changed to more descriptive 'trigFuncs'. . All elements of snu[3] are now explicitly defaulted to zero. The interface changes have necessitated changes in the call to sourceFunc_pol() within traceray() and the addition of a new argument gAux (whose values are set within raytrace()) to the traceray() call. I removed raytrace_1.4() at this point. This obsolete chunk of code was supposed to have already been removed in a previous PR. --- src/lime.h | 4 +- src/raytrace.c | 274 +++++------------------------------------------ src/sourcefunc.c | 43 ++++---- 3 files changed, 54 insertions(+), 267 deletions(-) diff --git a/src/lime.h b/src/lime.h index a864393..ba4f2b5 100644 --- a/src/lime.h +++ b/src/lime.h @@ -287,11 +287,11 @@ int sortangles(double *, int, struct grid *, const gsl_rng *); void sourceFunc(double *, double *, double, molData *,double,struct grid *,int,int, int,int); void sourceFunc_line(double *,double *,molData *, double, struct grid *, int, int,int); void sourceFunc_cont(double *,double *, struct grid *, int, int,int); -void sourceFunc_pol(double *, double *, double, molData *, double, struct grid *, int, int, int, double); +void sourceFunc_pol(const double, const double*, const molData, const struct pop2, const int, const double, double*, double*); void stateq(int, struct grid *, molData *, int, inputPars *,gridPointData *,double *); void statistics(int, molData *, struct grid *, int *, double *, double *, int *); void stokesangles(const double B[3], const double, double *); -void traceray(rayData, int, int, inputPars *, struct grid *, molData *, image *, int, int *, int *, double); +void traceray(rayData, int, int, inputPars*, struct grid*, struct gAuxType*, molData*, image*, int, int*, int*, double); void calcLineAmpLinear(struct grid*, const int, const int, const double, const double, double*); void calcLineAmpSample(double*, double*, const double, const double, const double, double*); void calcLineAmpSpline(struct grid*, const int, const int, const double, const double, double*); diff --git a/src/raytrace.c b/src/raytrace.c index ecd2788..8a28bdf 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -69,7 +69,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, inputPars *par, struct grid *g, struct gAuxType *gAux, 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. @@ -115,7 +115,7 @@ Note that the algorithm employed here is similar to that employed in the functio line_plane_intersect(g,&ds,posn,&nposn,dx,x,cutoff); /* Returns a new ds equal to the distance to the next Voronoi face, and nposn, the ID of the grid cell that abuts that face. */ if(par->polarization){ for(ichan=0;ichanncell); + for(ppi=0;ppincell;ppi++){ + gAux[ppi].mol = malloc(sizeof(*(gAux[ppi].mol))*par->nSpecies); + for(molI=0;molInSpecies;molI++){ + gAux[ppi].mol[molI].specNumDens = malloc(sizeof(*(gAux[ppi].mol[molI].specNumDens))*m[molI].nlev); + gAux[ppi].mol[molI].dust = malloc(sizeof(*(gAux[ppi].mol[molI].dust)) *m[molI].nline); + gAux[ppi].mol[molI].knu = malloc(sizeof(*(gAux[ppi].mol[molI].knu)) *m[molI].nline); + + for(li=0;linSpecies, m, &counta, &countb, &nlinetot); if(img[im].doline==0) nlinetot=1; @@ -268,7 +292,7 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ ray.x = size*(gsl_rng_uniform(threadRans[threadI])+px%img[im].pxls)-size*img[im].pxls/2.; ray.y = size*(gsl_rng_uniform(threadRans[threadI])+px/img[im].pxls)-size*img[im].pxls/2.; - traceray(ray, tmptrans, im, par, g, m, img, nlinetot, counta, countb, cutoff); + traceray(ray, tmptrans, im, par, g, gAux, m, img, nlinetot, counta, countb, cutoff); #pragma omp critical { @@ -298,245 +322,3 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ gsl_rng_free(ran); } - -void -raytrace_1_4(int im, inputPars *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. -I've adapted it slightly so it makes use of the function traceray(), -which was modified from the function tracerays() in v1.4. This algorithm -is not currently used, but may be useful as an option; that's why I have -kept it. - */ - - int *counta, *countb,nlinetot; - int ichan,i,px,iline,tmptrans,count; - double size,xp,yp,minfreq,absDeltaFreq; - double cutoff; - - gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2); /* Random number generator */ -#ifdef TEST - gsl_rng_set(ran,178490); -#else - gsl_rng_set(ran,time(0)); -#endif - rayData *rays; - - int sg,n; - double cx,cy; - - double x1,x2,x3,y1,y2,y3,z1,z2,z3,xt[3],yt[3],di,p,d1,d2,d3,temp1; - int zt[3]; - int c; - - char flags[255]; - boolT ismalloc = False; - facetT *facet, *neighbor, **neighborp;; - vertexT *vertex,**vertexp; - coordT *pt_array; - - int id; - coordT point[3]; - boolT isoutside; - realT bestdist; - - size=img[im].distance*img[im].imgres; - - /* Determine whether there are blended lines or not */ - lineCount(par->nSpecies, m, &counta, &countb, &nlinetot); - if(img[im].doline==0) nlinetot=1; - - /* Fix the image parameters */ - if(img[im].freq < 0) img[im].freq=m[0].freq[img[im].trans]; - if(img[im].nchan == 0 && img[im].bandwidth>0){ - img[im].nchan=(int) (img[im].bandwidth/(img[im].velres/CLIGHT*img[im].freq)); - } else if (img[im].velres<0 && img[im].bandwidth>0){ - img[im].velres = img[im].bandwidth*CLIGHT/img[im].freq/img[im].nchan; - } else img[im].bandwidth = img[im].nchan*img[im].velres/CLIGHT * img[im].freq; - - if(img[im].trans<0){ - iline=0; - minfreq=fabs(img[im].freq-m[0].freq[iline]); - tmptrans=iline; - for(iline=1;ilinepIntensity)); - - for(i=0;ipIntensity;i++){ - rays[i].x=g[i].x[0]; - rays[i].y=g[i].x[1]; - rays[i].tau=malloc(sizeof(double) * img[im].nchan); - rays[i].intensity=malloc(sizeof(double) * img[im].nchan); - for(ichan=0;ichanpIntensity); - - for(i=0;ipIntensity;i++) { - pt_array[i*2+0]=rays[i].x; - pt_array[i*2+1]=rays[i].y; - } - - sprintf(flags,"qhull v s Qbb T0"); - if (!qh_new_qhull(2, par->pIntensity, pt_array, ismalloc, flags, NULL, NULL)) { - - qh_setvoronoi_all(); - - FORALLvertices { - i=qh_pointid(vertex->point); - - cx=0.; - cy=0.; - n=0; - FOREACHneighbor_(vertex) { - if (!neighbor->upperdelaunay) n++; - } - if(n>0){ - - - } else { - if(!silent) bail_out("Qhull error"); - exit(0); - } - - FOREACHneighbor_(vertex) { - if (!neighbor->upperdelaunay) { - cx+=neighbor->center[0]; - cy+=neighbor->center[1]; - } - } - - rays[i].x = rays[i].x - (rays[i].x-cx/ (double) n)*0.1; - rays[i].y = rays[i].y - (rays[i].y-cy/ (double) n)*0.1; - } - } else { - printf("qhull error\n"); - } - - qh_freeqhull(!qh_ALL); - free(pt_array); - } - - cutoff = par->minScale*1.0e-7; - - /* Main loop through rays */ - count=0; - for(px=0;pxpIntensity;px++){ - traceray(rays[px], tmptrans, im, par, g, m, img, nlinetot, counta, countb, cutoff); - ++count; - if(!silent) progressbar((double)(count)/(double)(par->pIntensity-1), 13); - } - - /* Remap rays onto pixel grid */ - pt_array=malloc(2*sizeof(coordT)*par->pIntensity); - - for(i=0;ipIntensity;i++) { - pt_array[i*2+0]=rays[i].x; - pt_array[i*2+1]=rays[i].y; - } - -/* This allocation belongs to "Shepard's method" below - d=malloc(sizeof(double)*par->pIntensity); -*/ - size=img[im].distance*img[im].imgres; - - sprintf(flags,"qhull d Qbb"); - if (!qh_new_qhull(2, par->pIntensity, pt_array, ismalloc, flags, NULL, NULL)) { - for(px=0;pxpIntensity;i++){ - // d[i]=1./pow(sqrt(pow(xp-rays[i].x,2)+ pow(yp-rays[i].y,2)),8.); - temp1 = (xp-rays[i].x)*(xp-rays[i].x)+ (yp-rays[i].y)*(yp-rays[i].y) - d[i]=1./(temp1*temp1*temp1*temp1); - img[im].pixel[px].intense[ichan] += rays[i].intensity[ichan]*d[i]; -**** how to handle img[im].pixel[px].tau[ichan]? - di+=d[i]; - } - img[im].pixel[px].intense[ichan] /= di; - } -*/ - - - point[0]=xp; - point[1]=yp; - point[2]=0.; - - qh_setdelaunay (3, 1, point); - facet= qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside); - - c=0; - FOREACHvertex_( facet->vertices ) { - id=qh_pointid(vertex->point); - xt[c]=rays[id].x; yt[c]=rays[id].y; zt[c]=id; - c++; - } - - - x1=xt[0];x2=xt[1];x3=xt[2]; - y1=yt[0];y2=yt[1];y3=yt[2]; - - for(ichan=0;ichanpIntensity;i++){ - free(rays[i].tau); - free(rays[i].intensity); - } - free(rays); - free(counta); - free(countb); -} - - diff --git a/src/sourcefunc.c b/src/sourcefunc.c index b1d5f8c..127bbb9 100644 --- a/src/sourcefunc.c +++ b/src/sourcefunc.c @@ -70,34 +70,39 @@ sourceFunc_cont(double *jnu, double *alpha,struct grid *g,int pos,int ispec, int return; } +/*....................................................................*/ void -sourceFunc_pol(double *snu, double *dtau, double ds, molData *m,double vfac,struct grid *g,int pos,int ispec, int iline,double incl){ +sourceFunc_pol(const double ds, const double B[3], const molData md, const struct pop2 gm\ + , const int lineI, const double incl, double *snu, double *dtau){ + double dSigma, dSigma2, dI, dQ, dU, alpha; - double angle[3]; - - stokesangles(g[pos].B, incl, angle); - + double trigFuncs[3]; + + stokesangles(B, incl, trigFuncs); + /* Emission */ /* Continuum part: j_nu = rho_dust * kappa_nu */ - dSigma = g[pos].mol[ispec].dust[iline]*g[pos].mol[ispec].knu[iline]; - dSigma2 = maxp*g[pos].mol[ispec].dust[iline]*g[pos].mol[ispec].knu[iline]*(0.5*angle[2]*angle[2]-1./3.); - dI = dSigma - dSigma2; - dQ = maxp*g[pos].mol[ispec].dust[iline]*g[pos].mol[ispec].knu[iline]*(2.*angle[0]*angle[0]-1.)*angle[2]*angle[2]; - dU = maxp*g[pos].mol[ispec].dust[iline]*g[pos].mol[ispec].knu[iline]*(2.*angle[0]*angle[1]*angle[2]*angle[2]); - + dSigma = gm.dust[lineI]*gm.knu[lineI]; + dSigma2 = maxp*gm.dust[lineI]*gm.knu[lineI]*(0.5*trigFuncs[2]*trigFuncs[2]-1./3.); + dI = dSigma - dSigma2; + dQ = maxp*gm.dust[lineI]*gm.knu[lineI]*(2.*trigFuncs[0]*trigFuncs[0]-1.)*trigFuncs[2]*trigFuncs[2]; + dU = maxp*gm.dust[lineI]*gm.knu[lineI]*(2.*trigFuncs[0]*trigFuncs[1]*trigFuncs[2]*trigFuncs[2]); + /* Absorption */ /* Continuum part: Dust opacity */ - alpha = g[pos].mol[ispec].knu[iline]; - + alpha = gm.knu[lineI]; + /* Calculate source function and tau */ - *snu=0.; + snu[0] = 0.0; + snu[1] = 0.0; + snu[2] = 0.0; *dtau=0.; if(fabs(alpha)>0.){ - snu[0]=(dI/alpha)*m[ispec].norminv; - snu[1]=-(dQ/alpha)*m[ispec].norminv; - snu[2]=-(dU/alpha)*m[ispec].norminv; - - *dtau= alpha*ds; + snu[0]= (dI/alpha)*md.norminv; + snu[1]=-(dQ/alpha)*md.norminv; + snu[2]=-(dU/alpha)*md.norminv; + + *dtau = alpha*ds; } return; } From 7cccb38be6c6a9fbb87b48b637e8de03e79ae4df Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 16:38:56 +0200 Subject: [PATCH 22/43] Changed sourcefunc_*() interfaces. I just changed sourceFunc_line() and sourceFunc_cont() to bring them more into line with sourceFunc_pol(). This necessitated changes in photon(), getjbar() and traceray(). In the last of these, the call to sourceFunc_cont() is slightly rearranged, with new variables contMolI and contLineI being set before the start of the loop. --- src/lime.h | 4 ++-- src/photon.c | 10 +++++----- src/raytrace.c | 17 ++++++++++++----- src/sourcefunc.c | 30 +++++++++++++++++------------- 4 files changed, 36 insertions(+), 25 deletions(-) diff --git a/src/lime.h b/src/lime.h index ba4f2b5..fa6bc67 100644 --- a/src/lime.h +++ b/src/lime.h @@ -285,8 +285,8 @@ void report(int, inputPars *, struct grid *); void smooth(inputPars *, 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_line(double *,double *,molData *, double, struct grid *, int, int,int); -void sourceFunc_cont(double *,double *, struct grid *, int, int,int); +void sourceFunc_line(const molData, const double, const struct populations, const int, double*, double*); +void sourceFunc_cont(const struct populations, const int, double*, double*); void sourceFunc_pol(const double, const double*, const molData, const struct pop2, const int, const double, double*, double*); void stateq(int, struct grid *, molData *, int, inputPars *,gridPointData *,double *); void statistics(int, molData *, struct grid *, int *, double *, double *, int *); diff --git a/src/photon.c b/src/photon.c index 76f18e0..9b65b9a 100644 --- a/src/photon.c +++ b/src/photon.c @@ -240,8 +240,8 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar jnu=0.; alpha=0.; - sourceFunc_line(&jnu,&alpha,m,vfac[counta[iline]],g,here,counta[iline],countb[iline]); - sourceFunc_cont(&jnu,&alpha,g,here,counta[iline],countb[iline]); + sourceFunc_line(m[counta[iline]],vfac[counta[iline]],g[here].mol[counta[iline]],countb[iline],&jnu,&alpha); + sourceFunc_cont(g[here].mol[counta[iline]],countb[iline],&jnu,&alpha); dtau=alpha*ds; if(dtau < -30) dtau = -30; @@ -265,7 +265,7 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar if(matrix[jline].line1 == jline || matrix[jline].line2 == jline){ if(!par->doPregrid) calcLineAmpSpline(g,here,dir,g[id].mol[counta[jline]].binv,deltav-matrix[jline].deltav,&vblend); else calcLineAmpLinear(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]); + sourceFunc_line(m[counta[jline]],vblend,g[here].mol[counta[jline]],countb[jline],&jnu,&alpha); dtau=alpha*ds; if(dtau < -30) dtau = -30; calcSourceFn(dtau, par, &remnantSnu, &expDTau); @@ -318,8 +318,8 @@ getjbar(int posn, molData *m, struct grid *g, inputPars *par, gridPointData *mp, jnu=0.; alpha=0.; - sourceFunc_line(&jnu,&alpha,m,mp[0].vfac[iphot],g,posn,counta[iline],countb[iline]); - sourceFunc_cont(&jnu,&alpha,g,posn,counta[iline],countb[iline]); + sourceFunc_line(m[counta[iline]],mp[0].vfac[iphot],g[posn].mol[counta[iline]],countb[iline],&jnu,&alpha); + sourceFunc_cont(g[posn].mol[counta[iline]],countb[iline],&jnu,&alpha); tau=alpha*halfFirstDs[iphot]; calcSourceFn(tau, par, &remnantSnu, &expTau); remnantSnu *= jnu*m[0].norminv*halfFirstDs[iphot]; diff --git a/src/raytrace.c b/src/raytrace.c index 8a28bdf..1daeb31 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -75,7 +75,7 @@ For a given image pixel position, this function evaluates the intensity of the t Note that the algorithm employed here is similar to that employed in the function photon() which calculates the average radiant flux impinging on a grid cell: namely the notional photon is started at the side of the model near the observer and 'propagated' in the receding direction until it 'reaches' the far side. This is rather non-physical in conception but it makes the calculation easier. */ - int ichan,posn,nposn,i,iline,molI,lineI; + int ichan,posn,nposn,i,iline,molI,lineI,contMolI,contLineI; double vfac=0.,x[3],dx[3],vThisChan; double deltav,ds,dist2,ndist2,xp,yp,zp,col,lineRedShift,jnu,alpha,remnantSnu,dtau,expDTau,snu_pol[3]; @@ -96,6 +96,15 @@ Note that the algorithm employed here is similar to that employed in the functio dx[i]= img[im].rotMat[i][2]; /* This points away from the observer. */ } + contMolI = 0; /****** Always?? */ + + if(img[im].doline && img[im].trans > -1) + contLineI = img[im].trans; + else if(img[im].doline && img[im].trans == -1) + contLineI = tmptrans; + else + contLineI = 0; + /* Find the grid point nearest to the starting x. */ i=0; dist2=(x[0]-g[i].x[0])*(x[0]-g[i].x[0]) + (x[1]-g[i].x[1])*(x[1]-g[i].x[1]) + (x[2]-g[i].x[2])*(x[2]-g[i].x[2]); @@ -149,13 +158,11 @@ Note that the algorithm employed here is similar to that employed in the functio else vfac=gaussline(deltav+veloproject(dx,g[posn].vel),g[posn].mol[molI].binv); /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. */ - sourceFunc_line(&jnu,&alpha,m,vfac,g,posn,molI,lineI); + sourceFunc_line(m[molI], vfac, g[posn].mol[molI], lineI, &jnu, &alpha); } } - if(img[im].doline && img[im].trans > -1) sourceFunc_cont(&jnu,&alpha,g,posn,0,img[im].trans); - else if(img[im].doline && img[im].trans == -1) sourceFunc_cont(&jnu,&alpha,g,posn,0,tmptrans); - else sourceFunc_cont(&jnu,&alpha,g,posn,0,0); + sourceFunc_cont(g[posn].mol[contMolI], contLineI, &jnu, &alpha); dtau=alpha*ds; calcSourceFn(dtau, par, &remnantSnu, &expDTau); remnantSnu *= jnu*m[0].norminv*ds; diff --git a/src/sourcefunc.c b/src/sourcefunc.c index 127bbb9..f0ca662 100644 --- a/src/sourcefunc.c +++ b/src/sourcefunc.c @@ -43,30 +43,34 @@ sourceFunc(double *snu, double *dtau, double ds, molData *m,double vfac,struct g } +/*....................................................................*/ void -sourceFunc_line(double *jnu, double *alpha, molData *m,double vfac,struct grid *g,int pos,int ispec, int iline){ - +sourceFunc_line(const molData md, const double vfac, const struct populations gm\ + , const int lineI, double *jnu, double *alpha){ + /* Line part: j_nu = v*consts*1/b*rho*n_i*A_ij */ - *jnu += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].mol[ispec].nmol*g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].aeinst[iline]; - + *jnu += vfac*HPIP*gm.binv*gm.nmol*gm.pops[md.lau[lineI]]*md.aeinst[lineI]; + /* Line part: alpha_nu = v*const*1/b*rho*(n_j*B_ij-n_i*B_ji) */ - *alpha += vfac*HPIP*g[pos].mol[ispec].binv*g[pos].mol[ispec].nmol*(g[pos].mol[ispec].pops[m[ispec].lal[iline]]*m[ispec].beinstl[iline] - -g[pos].mol[ispec].pops[m[ispec].lau[iline]]*m[ispec].beinstu[iline]); - + *alpha += vfac*HPIP*gm.binv*gm.nmol*(gm.pops[md.lal[lineI]]*md.beinstl[lineI] + -gm.pops[md.lau[lineI]]*md.beinstu[lineI]); + return; } +/*....................................................................*/ void -sourceFunc_cont(double *jnu, double *alpha,struct grid *g,int pos,int ispec, int iline){ - +sourceFunc_cont(const struct populations gm, const int lineI, double *jnu\ + , double *alpha){ + /* Emission */ /* Continuum part: j_nu = T_dust * kappa_nu */ - *jnu += g[pos].mol[ispec].dust[iline]*g[pos].mol[ispec].knu[iline]; - + *jnu += gm.dust[lineI]*gm.knu[lineI]; + /* Absorption */ /* Continuum part: Dust opacity */ - *alpha += g[pos].mol[ispec].knu[iline]; - + *alpha += gm.knu[lineI]; + return; } From b34bf51ceb5c2820a03947b910aabe0774756086 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 16:44:20 +0200 Subject: [PATCH 23/43] New functions sourceFunc_*_raytrace(). These duplicate the functionality of sourceFunc_cont() and sourceFunc_line() but are made necessary, or at least desirable, by the present poor way in which grid data is stored in memory structures. When this is fixed (which it will hopefully be eventually), the *_raytrace versions can be merged into the previous ones. --- src/lime.h | 2 ++ src/sourcefunc.c | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/src/lime.h b/src/lime.h index fa6bc67..9e423f6 100644 --- a/src/lime.h +++ b/src/lime.h @@ -287,6 +287,8 @@ int sortangles(double *, int, struct grid *, const gsl_rng *); void sourceFunc(double *, double *, double, molData *,double,struct grid *,int,int, int,int); void sourceFunc_line(const molData, const double, const struct populations, const int, double*, double*); void sourceFunc_cont(const struct populations, const int, double*, double*); +void sourceFunc_line_raytrace(const molData, const double, const struct pop2, const int, double*, double*); +void sourceFunc_cont_raytrace(const struct pop2, const int, double*, double*); void sourceFunc_pol(const double, const double*, const molData, const struct pop2, const int, const double, double*, double*); void stateq(int, struct grid *, molData *, int, inputPars *,gridPointData *,double *); void statistics(int, molData *, struct grid *, int *, double *, double *, int *); diff --git a/src/sourcefunc.c b/src/sourcefunc.c index f0ca662..9f73cc6 100644 --- a/src/sourcefunc.c +++ b/src/sourcefunc.c @@ -5,6 +5,8 @@ * Copyright (C) 2006-2014 Christian Brinch * Copyright (C) 2015 The LIME development team * +TODO: + - Merge sourceFunc_*_raytrace() and sourceFunc_*() after changes to the way grid data is stored makes this possible. */ #include "lime.h" @@ -58,6 +60,21 @@ sourceFunc_line(const molData md, const double vfac, const struct populations gm return; } +/*....................................................................*/ +void +sourceFunc_line_raytrace(const molData md, const double vfac\ + , const struct pop2 gm, const int lineI, double *jnu, double *alpha){ + + /* Line part: j_nu = v*consts*1/b*rho*n_i*A_ij */ + *jnu += vfac*HPIP*gm.specNumDens[md.lau[lineI]]*md.aeinst[lineI]; + + /* Line part: alpha_nu = v*const*1/b*rho*(n_j*B_ij-n_i*B_ji) */ + *alpha += vfac*HPIP*(gm.specNumDens[md.lal[lineI]]*md.beinstl[lineI] + -gm.specNumDens[md.lau[lineI]]*md.beinstu[lineI]); + + return; +} + /*....................................................................*/ void sourceFunc_cont(const struct populations gm, const int lineI, double *jnu\ @@ -74,6 +91,22 @@ sourceFunc_cont(const struct populations gm, const int lineI, double *jnu\ return; } +/*....................................................................*/ +void +sourceFunc_cont_raytrace(const struct pop2 gm, const int lineI, double *jnu\ + , double *alpha){ + + /* Emission */ + /* Continuum part: j_nu = T_dust * kappa_nu */ + *jnu += gm.dust[lineI]*gm.knu[lineI]; + + /* Absorption */ + /* Continuum part: Dust opacity */ + *alpha += gm.knu[lineI]; + + return; +} + /*....................................................................*/ void sourceFunc_pol(const double ds, const double B[3], const molData md, const struct pop2 gm\ From 7b5bf5d2b88f6bd572c392deb319db084acfb487 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 16:51:37 +0200 Subject: [PATCH 24/43] Rearranged if-block in traceray(). There is a test near the start of the function to see whether the ray will intersect the model sphere at all; if so, the following large block of ~110 statements (containing the entire functionality of the function) is entered. It is obviously much less spendthrift of indents to simply exit the function if the ray does _not_ pass through the model sphere. The code has been rearranged to implement this change. All the statements which were formerly within the block have been unindented by 2 spaces; there are no further changes in the present commit. --- src/raytrace.c | 175 +++++++++++++++++++++++++------------------------ 1 file changed, 89 insertions(+), 86 deletions(-) diff --git a/src/raytrace.c b/src/raytrace.c index 1daeb31..9fe7834 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -87,111 +87,114 @@ Note that the algorithm employed here is similar to that employed in the functio xp=ray.x; yp=ray.y; - if((xp*xp+yp*yp)/par->radiusSqu <= 1 ) { - zp=-sqrt(par->radiusSqu-(xp*xp+yp*yp)); /* There are two points of intersection between the line of sight and the spherical model surface; this is the Z coordinate (in the unrotated frame) of the one nearer to the observer. */ + /* The model is circular in projection. We only follow the ray if it will intersect the model. + */ + if((xp*xp+yp*yp)>par->radiusSqu) + return; - /* Rotate the line of sight as desired. */ - for(i=0;i<3;i++){ - x[i]=xp*img[im].rotMat[i][0] + yp*img[im].rotMat[i][1] + zp*img[im].rotMat[i][2]; - dx[i]= img[im].rotMat[i][2]; /* This points away from the observer. */ - } + zp=-sqrt(par->radiusSqu-(xp*xp+yp*yp)); /* There are two points of intersection between the line of sight and the spherical model surface; this is the Z coordinate (in the unrotated frame) of the one nearer to the observer. */ - contMolI = 0; /****** Always?? */ - - if(img[im].doline && img[im].trans > -1) - contLineI = img[im].trans; - else if(img[im].doline && img[im].trans == -1) - contLineI = tmptrans; - else - contLineI = 0; - - /* Find the grid point nearest to the starting x. */ - i=0; - dist2=(x[0]-g[i].x[0])*(x[0]-g[i].x[0]) + (x[1]-g[i].x[1])*(x[1]-g[i].x[1]) + (x[2]-g[i].x[2])*(x[2]-g[i].x[2]); - posn=i; - for(i=1;incell;i++){ - ndist2=(x[0]-g[i].x[0])*(x[0]-g[i].x[0]) + (x[1]-g[i].x[1])*(x[1]-g[i].x[1]) + (x[2]-g[i].x[2])*(x[2]-g[i].x[2]); - if(ndist2 -1) + contLineI = img[im].trans; + else if(img[im].doline && img[im].trans == -1) + contLineI = tmptrans; + else + contLineI = 0; + + /* Find the grid point nearest to the starting x. */ + i=0; + dist2=(x[0]-g[i].x[0])*(x[0]-g[i].x[0]) + (x[1]-g[i].x[1])*(x[1]-g[i].x[1]) + (x[2]-g[i].x[2])*(x[2]-g[i].x[2]); + posn=i; + for(i=1;incell;i++){ + ndist2=(x[0]-g[i].x[0])*(x[0]-g[i].x[0]) + (x[1]-g[i].x[1])*(x[1]-g[i].x[1]) + (x[2]-g[i].x[2])*(x[2]-g[i].x[2]); + if(ndist2polarization){ - for(ichan=0;ichanpolarization){ + for(ichan=0;ichan img[im].freq-img[im].bandwidth/2. - && m[molI].freq[lineI] < img[im].freq+img[im].bandwidth/2.){ - /* Calculate the red shift of the transition wrt to the frequency specified for the image. */ - if(img[im].trans > -1){ - lineRedShift=(m[molI].freq[img[im].trans]-m[molI].freq[lineI])/m[molI].freq[img[im].trans]*CLIGHT; - } else { - lineRedShift=(img[im].freq-m[molI].freq[lineI])/img[im].freq*CLIGHT; - } - - vThisChan=(ichan-(img[im].nchan-1)/2.)*img[im].velres; /* Consistent with the WCS definition in writefits(). */ - deltav = vThisChan - img[im].source_vel - lineRedShift; - /* Line centre occurs when deltav = the recession velocity of the radiating material. Explanation of the signs of the 2nd and 3rd terms on the RHS: (i) A bulk source velocity (which is defined as >0 for the receding direction) should be added to the material velocity field; this is equivalent to subtracting it from deltav, as here. (ii) A positive value of lineRedShift means the line is red-shifted wrt to the frequency specified for the image. The effect is the same as if the line and image frequencies were the same, but the bulk recession velocity were higher. lineRedShift should thus be added to the recession velocity, which is equivalent to subtracting it from deltav, as here. */ - - /* Calculate an approximate average line-shape function at deltav within the Voronoi cell. */ - if(!par->pregrid) calcLineAmpSample(x,dx,ds,g[posn].mol[molI].binv,deltav,&vfac); - else vfac=gaussline(deltav+veloproject(dx,g[posn].vel),g[posn].mol[molI].binv); - - /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. */ - sourceFunc_line(m[molI], vfac, g[posn].mol[molI], lineI, &jnu, &alpha); + ray.tau[ichan]+=dtau; + } + } else { + for(ichan=0;ichan img[im].freq-img[im].bandwidth/2. + && m[molI].freq[lineI] < img[im].freq+img[im].bandwidth/2.){ + /* Calculate the red shift of the transition wrt to the frequency specified for the image. */ + if(img[im].trans > -1){ + lineRedShift=(m[molI].freq[img[im].trans]-m[molI].freq[lineI])/m[molI].freq[img[im].trans]*CLIGHT; + } else { + lineRedShift=(img[im].freq-m[molI].freq[lineI])/img[im].freq*CLIGHT; } + + vThisChan=(ichan-(img[im].nchan-1)/2.)*img[im].velres; /* Consistent with the WCS definition in writefits(). */ + deltav = vThisChan - img[im].source_vel - lineRedShift; + /* Line centre occurs when deltav = the recession velocity of the radiating material. Explanation of the signs of the 2nd and 3rd terms on the RHS: (i) A bulk source velocity (which is defined as >0 for the receding direction) should be added to the material velocity field; this is equivalent to subtracting it from deltav, as here. (ii) A positive value of lineRedShift means the line is red-shifted wrt to the frequency specified for the image. The effect is the same as if the line and image frequencies were the same, but the bulk recession velocity were higher. lineRedShift should thus be added to the recession velocity, which is equivalent to subtracting it from deltav, as here. */ + + /* Calculate an approximate average line-shape function at deltav within the Voronoi cell. */ + if(!par->pregrid) calcLineAmpSample(x,dx,ds,g[posn].mol[molI].binv,deltav,&vfac); + else vfac=gaussline(deltav+veloproject(dx,g[posn].vel),g[posn].mol[molI].binv); + + /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. */ + sourceFunc_line(m[molI], vfac, g[posn].mol[molI], lineI, &jnu, &alpha); } + } - sourceFunc_cont(g[posn].mol[contMolI], contLineI, &jnu, &alpha); - dtau=alpha*ds; - calcSourceFn(dtau, par, &remnantSnu, &expDTau); - remnantSnu *= jnu*m[0].norminv*ds; + sourceFunc_cont(g[posn].mol[contMolI], contLineI, &jnu, &alpha); + dtau=alpha*ds; + calcSourceFn(dtau, par, &remnantSnu, &expDTau); + remnantSnu *= jnu*m[0].norminv*ds; #ifdef FASTEXP - ray.intensity[ichan]+=FastExp(ray.tau[ichan])*remnantSnu; + ray.intensity[ichan]+=FastExp(ray.tau[ichan])*remnantSnu; #else - ray.intensity[ichan]+= exp(-ray.tau[ichan])*remnantSnu; + ray.intensity[ichan]+= exp(-ray.tau[ichan])*remnantSnu; #endif - ray.tau[ichan]+=dtau; - } + ray.tau[ichan]+=dtau; } + } - /* Move the working point to the edge of the next Voronoi cell. */ - for(i=0;i<3;i++) x[i]+=ds*dx[i]; - col+=ds; - posn=nposn; - } while(col < 2.0*fabs(zp)); + /* Move the working point to the edge of the next Voronoi cell. */ + for(i=0;i<3;i++) x[i]+=ds*dx[i]; + col+=ds; + posn=nposn; + } while(col < 2.0*fabs(zp)); - /* Add or subtract cmb. */ + /* Add or subtract cmb. */ #ifdef FASTEXP - for(ichan=0;ichan Date: Fri, 8 Apr 2016 17:00:33 +0200 Subject: [PATCH 25/43] Changed g to gp and m to md in traceray(). For better searchability. --- src/raytrace.c | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/raytrace.c b/src/raytrace.c index 9fe7834..babb449 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -69,7 +69,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, struct gAuxType *gAux, molData *m, image *img, int nlinetot, int *counta, int *countb, double cutoff){ +traceray(rayData ray, int tmptrans, int im, inputPars *par, struct grid *gp, struct gAuxType *gAux, molData *md, 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. @@ -111,10 +111,10 @@ Note that the algorithm employed here is similar to that employed in the functio /* Find the grid point nearest to the starting x. */ i=0; - dist2=(x[0]-g[i].x[0])*(x[0]-g[i].x[0]) + (x[1]-g[i].x[1])*(x[1]-g[i].x[1]) + (x[2]-g[i].x[2])*(x[2]-g[i].x[2]); + dist2=(x[0]-gp[i].x[0])*(x[0]-gp[i].x[0]) + (x[1]-gp[i].x[1])*(x[1]-gp[i].x[1]) + (x[2]-gp[i].x[2])*(x[2]-gp[i].x[2]); posn=i; for(i=1;incell;i++){ - ndist2=(x[0]-g[i].x[0])*(x[0]-g[i].x[0]) + (x[1]-g[i].x[1])*(x[1]-g[i].x[1]) + (x[2]-g[i].x[2])*(x[2]-g[i].x[2]); + ndist2=(x[0]-gp[i].x[0])*(x[0]-gp[i].x[0]) + (x[1]-gp[i].x[1])*(x[1]-gp[i].x[1]) + (x[2]-gp[i].x[2])*(x[2]-gp[i].x[2]); if(ndist2polarization){ for(ichan=0;ichan img[im].freq-img[im].bandwidth/2. - && m[molI].freq[lineI] < img[im].freq+img[im].bandwidth/2.){ + if(img[im].doline && md[molI].freq[lineI] > img[im].freq-img[im].bandwidth/2. + && md[molI].freq[lineI] < img[im].freq+img[im].bandwidth/2.){ /* Calculate the red shift of the transition wrt to the frequency specified for the image. */ if(img[im].trans > -1){ - lineRedShift=(m[molI].freq[img[im].trans]-m[molI].freq[lineI])/m[molI].freq[img[im].trans]*CLIGHT; + lineRedShift=(md[molI].freq[img[im].trans]-md[molI].freq[lineI])/md[molI].freq[img[im].trans]*CLIGHT; } else { - lineRedShift=(img[im].freq-m[molI].freq[lineI])/img[im].freq*CLIGHT; + lineRedShift=(img[im].freq-md[molI].freq[lineI])/img[im].freq*CLIGHT; } vThisChan=(ichan-(img[im].nchan-1)/2.)*img[im].velres; /* Consistent with the WCS definition in writefits(). */ @@ -158,18 +158,18 @@ Note that the algorithm employed here is similar to that employed in the functio /* Line centre occurs when deltav = the recession velocity of the radiating material. Explanation of the signs of the 2nd and 3rd terms on the RHS: (i) A bulk source velocity (which is defined as >0 for the receding direction) should be added to the material velocity field; this is equivalent to subtracting it from deltav, as here. (ii) A positive value of lineRedShift means the line is red-shifted wrt to the frequency specified for the image. The effect is the same as if the line and image frequencies were the same, but the bulk recession velocity were higher. lineRedShift should thus be added to the recession velocity, which is equivalent to subtracting it from deltav, as here. */ /* Calculate an approximate average line-shape function at deltav within the Voronoi cell. */ - if(!par->pregrid) calcLineAmpSample(x,dx,ds,g[posn].mol[molI].binv,deltav,&vfac); - else vfac=gaussline(deltav+veloproject(dx,g[posn].vel),g[posn].mol[molI].binv); + if(!par->pregrid) calcLineAmpSample(x,dx,ds,gp[posn].mol[molI].binv,deltav,&vfac); + else vfac=gaussline(deltav+veloproject(dx,gp[posn].vel),gp[posn].mol[molI].binv); /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. */ - sourceFunc_line(m[molI], vfac, g[posn].mol[molI], lineI, &jnu, &alpha); + sourceFunc_line(md[molI], vfac, gp[posn].mol[molI], lineI, &jnu, &alpha); } } - sourceFunc_cont(g[posn].mol[contMolI], contLineI, &jnu, &alpha); + sourceFunc_cont(gp[posn].mol[contMolI], contLineI, &jnu, &alpha); dtau=alpha*ds; calcSourceFn(dtau, par, &remnantSnu, &expDTau); - remnantSnu *= jnu*m[0].norminv*ds; + remnantSnu *= jnu*md[0].norminv*ds; #ifdef FASTEXP ray.intensity[ichan]+=FastExp(ray.tau[ichan])*remnantSnu; #else @@ -188,11 +188,11 @@ Note that the algorithm employed here is similar to that employed in the functio /* Add or subtract cmb. */ #ifdef FASTEXP for(ichan=0;ichan Date: Fri, 8 Apr 2016 17:09:16 +0200 Subject: [PATCH 26/43] Cosmetic changes to traceray(). --- src/raytrace.c | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/raytrace.c b/src/raytrace.c index babb449..f90bba5 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -5,6 +5,7 @@ * Copyright (C) 2006-2014 Christian Brinch * Copyright (C) 2015 The LIME development team * +TODO: sort out snu_pol in traceray(). */ #include "lime.h" @@ -75,7 +76,7 @@ For a given image pixel position, this function evaluates the intensity of the t Note that the algorithm employed here is similar to that employed in the function photon() which calculates the average radiant flux impinging on a grid cell: namely the notional photon is started at the side of the model near the observer and 'propagated' in the receding direction until it 'reaches' the far side. This is rather non-physical in conception but it makes the calculation easier. */ - int ichan,posn,nposn,i,iline,molI,lineI,contMolI,contLineI; + int ichan,di,i,posn,nposn,polMolI,polLineI,contMolI,contLineI,iline,molI,lineI; double vfac=0.,x[3],dx[3],vThisChan; double deltav,ds,dist2,ndist2,xp,yp,zp,col,lineRedShift,jnu,alpha,remnantSnu,dtau,expDTau,snu_pol[3]; @@ -95,9 +96,9 @@ Note that the algorithm employed here is similar to that employed in the functio zp=-sqrt(par->radiusSqu-(xp*xp+yp*yp)); /* There are two points of intersection between the line of sight and the spherical model surface; this is the Z coordinate (in the unrotated frame) of the one nearer to the observer. */ /* Rotate the line of sight as desired. */ - for(i=0;i<3;i++){ - x[i]=xp*img[im].rotMat[i][0] + yp*img[im].rotMat[i][1] + zp*img[im].rotMat[i][2]; - dx[i]= img[im].rotMat[i][2]; /* This points away from the observer. */ + for(di=0;dipolarization){ + polMolI = 0; /****** Always?? */ + polLineI = 0; /****** Always?? */ for(ichan=0;ichan img[im].freq-img[im].bandwidth/2. && md[molI].freq[lineI] < img[im].freq+img[im].bandwidth/2.){ - /* Calculate the red shift of the transition wrt to the frequency specified for the image. */ + /* Calculate the red shift of the transition wrt to the frequency specified for the image. + */ if(img[im].trans > -1){ lineRedShift=(md[molI].freq[img[im].trans]-md[molI].freq[lineI])/md[molI].freq[img[im].trans]*CLIGHT; } else { @@ -162,11 +167,11 @@ Note that the algorithm employed here is similar to that employed in the functio else vfac=gaussline(deltav+veloproject(dx,gp[posn].vel),gp[posn].mol[molI].binv); /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. */ - sourceFunc_line(md[molI], vfac, gp[posn].mol[molI], lineI, &jnu, &alpha); + sourceFunc_line_raytrace(md[molI],vfac,gAux[posn].mol[molI],lineI,&jnu,&alpha); } } - sourceFunc_cont(gp[posn].mol[contMolI], contLineI, &jnu, &alpha); + sourceFunc_cont_raytrace(gAux[posn].mol[contMolI], contLineI, &jnu, &alpha); dtau=alpha*ds; calcSourceFn(dtau, par, &remnantSnu, &expDTau); remnantSnu *= jnu*md[0].norminv*ds; From 42ee18e8aef00bbf5eb9eeff80db8446e50e4871 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 8 Apr 2016 17:13:55 +0200 Subject: [PATCH 27/43] More rearrangements in traceray(). The substantive changes are: . Where par->polarization and FASTEXP was defined, one of the two necessary exponentials was still calculated by the native exp() function. This has been changed now such that Fastexp() is called for both. . sourceFunc_cont_raytrace() is no longer (unnecessarily) called once per channel, but only once; the results being stored and added to jnu and alpha for each new channel. --- src/raytrace.c | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/raytrace.c b/src/raytrace.c index f90bba5..6fc50d4 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -77,7 +77,7 @@ For a given image pixel position, this function evaluates the intensity of the t Note that the algorithm employed here is similar to that employed in the function photon() which calculates the average radiant flux impinging on a grid cell: namely the notional photon is started at the side of the model near the observer and 'propagated' in the receding direction until it 'reaches' the far side. This is rather non-physical in conception but it makes the calculation easier. */ int ichan,di,i,posn,nposn,polMolI,polLineI,contMolI,contLineI,iline,molI,lineI; - double vfac=0.,x[3],dx[3],vThisChan; + double vfac=0.,x[3],dx[3],vThisChan,contJnu,contAlpha; double deltav,ds,dist2,ndist2,xp,yp,zp,col,lineRedShift,jnu,alpha,remnantSnu,dtau,expDTau,snu_pol[3]; for(ichan=0;ichan Date: Fri, 8 Apr 2016 17:24:10 +0200 Subject: [PATCH 28/43] Added calcLineAmpInterp() and traceray_smooth(). I also changed the interface of traceray(). This was inadvertent and should have been in a separate commit. --- src/lime.h | 4 +- src/raytrace.c | 237 +++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 235 insertions(+), 6 deletions(-) diff --git a/src/lime.h b/src/lime.h index 9e423f6..3bb465f 100644 --- a/src/lime.h +++ b/src/lime.h @@ -293,7 +293,9 @@ void sourceFunc_pol(const double, const double*, const molData, const struct pop void stateq(int, struct grid *, molData *, int, inputPars *,gridPointData *,double *); void statistics(int, molData *, struct grid *, int *, double *, double *, int *); void stokesangles(const double B[3], const double, double *); -void traceray(rayData, int, int, inputPars*, struct grid*, struct gAuxType*, molData*, image*, int, int*, int*, double); +void traceray(rayData, inputPars*, const int, image*, const int, struct grid*, struct gAuxType*, molData*, const int, int*, int*, const double); +void traceray_smooth(rayData, inputPars*, const int, image*, const int, struct grid*, struct gAuxType*, molData*, const int, int*, int*, struct cell*, const unsigned long, const double, gridInterp gips[3], const int, const double); +void calcLineAmpInterp(const double, const double, const double, double*); void calcLineAmpLinear(struct grid*, const int, const int, const double, const double, double*); void calcLineAmpSample(double*, double*, const double, const double, const double, double*); void calcLineAmpSpline(struct grid*, const int, const int, const double, const double, double*); diff --git a/src/raytrace.c b/src/raytrace.c index 6fc50d4..2698755 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -38,7 +38,27 @@ The bulk velocity of the model material can vary significantly with position, th return; } +/*....................................................................*/ +void calcLineAmpInterp(const double velCmpntRay, const double binv\ + , const double deltav, double *vfac){ + /* +The bulk velocity of the model material can vary significantly with position, thus so can the value of the line-shape function at a given frequency and direction. The present function calculates 'vfac', an approximate average of the line-shape function along a path of length ds in the direction of the line of sight. + */ + double v,val; + v = deltav - velCmpntRay; /* velCmpntRay is the component of the local bulk velocity in the direction of the ray, whereas deltav is the recession velocity of the channel we are interested in (corrected for bulk source velocity and line displacement from the nominal frequency). Remember also that, since the ray points away from the observer, positive values of the projected velocity also represent recessions. Line centre occurs when v==0, i.e. when deltav==velCmpntRay. That is the reason for the subtraction here. */ + val = fabs(v)*binv; + if(val <= 2500.){ +#ifdef FASTEXP + *vfac = FastExp(val*val); +#else + *vfac = exp(-(val*val)); +#endif + }else + *vfac = 0.0; +} + +/*....................................................................*/ void line_plane_intersect(struct grid *g, double *ds, int posn, int *nposn, double *dx, double *x, double cutoff){ /* @@ -68,17 +88,21 @@ This function returns ds as the (always positive-valued) distance between the pr if(*nposn==-1) *nposn=posn; } - +/*....................................................................*/ void -traceray(rayData ray, int tmptrans, int im, inputPars *par, struct grid *gp, struct gAuxType *gAux, molData *md, image *img, int nlinetot, int *counta, int *countb, double cutoff){ +traceray(rayData ray, inputPars *par, const int tmptrans, image *img\ + , const int im, struct grid *gp, struct gAuxType *gAux, molData *md\ + , const int nlinetot, int *counta, int *countb, const 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. Note that the algorithm employed here is similar to that employed in the function photon() which calculates the average radiant flux impinging on a grid cell: namely the notional photon is started at the side of the model near the observer and 'propagated' in the receding direction until it 'reaches' the far side. This is rather non-physical in conception but it makes the calculation easier. */ + int ichan,di,i,posn,nposn,polMolI,polLineI,contMolI,contLineI,iline,molI,lineI; - double vfac=0.,x[3],dx[3],vThisChan,contJnu,contAlpha; - double deltav,ds,dist2,ndist2,xp,yp,zp,col,lineRedShift,jnu,alpha,remnantSnu,dtau,expDTau,snu_pol[3]; + double xp,yp,zp,x[DIM],dx[DIM],dist2,ndist2,col,ds,snu_pol[3],dtau; + double contJnu,contAlpha,jnu,alpha,lineRedShift,vThisChan,deltav,vfac=0.; + double remnantSnu,expDTau; for(ichan=0;ichanpar->radiusSqu) + return; + + zp=-sqrt(par->radiusSqu-(xp*xp+yp*yp)); /* There are two points of intersection between the line of sight and the spherical model surface; this is the Z coordinate (in the unrotated frame) of the one nearer to the observer. */ + + /* Rotate the line of sight as desired. */ + for(di=0;di -1) + contLineI = img[im].trans; + else if(img[im].doline && img[im].trans == -1) + contLineI = tmptrans; + else + contLineI = 0; + /* Find the chain of cells the ray passes through. + */ + status = followRayThroughDelCells(x, dir, gp, dc, numCells, epsilon\ + , &entryIntcptFirstCell, &chainOfCellIds, &cellExitIntcpts, &lenChainPtrs);//, 0); + + if(status!=0){ + free(chainOfCellIds); + free(cellExitIntcpts); + return; + } + + entryI = 0; + exitI = 1; + dci = chainOfCellIds[0]; + + /* Obtain the indices of the grid points on the vertices of the entry face. + */ + vvi = 0; + for(vi=0;viid; + } + } + + /* Calculate, for each of the 3 vertices of the entry face, the displacement components in the direction of 'dir'. *** NOTE *** that if all the rays are parallel, we could precalculate these for all the vertices. + */ + for(vi=0;vinSpecies, &gips[entryI]); + + for(ci=0;ciid; + } + } + + /* Calculate, for each of the 3 vertices of the exit face, the displacement components in the direction of 'dir'. *** NOTE *** that if all the rays are parallel, we could precalculate these for all the vertices. + */ + for(vi=0;vinSpecies, &gips[exitI]); + + /* At this point we have interpolated all the values of interest to both the entry and exit points of the cell. Now we break the path between entry and exit into several segments and calculate all these values at the midpoint of each segment. + +At the moment I will fix the number of segments, but it might possibly be faster to rather have a fixed segment length (in barycentric coordinates) and vary the number depending on how many of these lengths fit in the path between entry and exit. + */ + ds = (gips[exitI].xCmpntRay - gips[entryI].xCmpntRay)*oneOnNumSegments; + + for(si=0;sinSpecies, oneOnNumSegments, si); + + if(par->polarization){ //************************ WTF is with this snu_pol????? + polMolI = 0; /****** Always?? */ + polLineI = 0; /****** Always?? */ + for(ichan=0;ichan img[im].freq-img[im].bandwidth*0.5 + && md[molI].freq[lineI] < img[im].freq+img[im].bandwidth*0.5){ + /* Calculate the red shift of the transition wrt to the frequency specified for the image. + */ + if(img[im].trans > -1){ + lineRedShift=(md[molI].freq[img[im].trans]-md[molI].freq[lineI])/md[molI].freq[img[im].trans]*CLIGHT; + } else { + lineRedShift=(img[im].freq-md[molI].freq[lineI])/img[im].freq*CLIGHT; + } + + vThisChan=(ichan-(img[im].nchan-1)*0.5)*img[im].velres; /* Consistent with the WCS definition in writefits(). */ + deltav = vThisChan - img[im].source_vel - lineRedShift; + /* Line centre occurs when deltav = the recession velocity of the radiating material. Explanation of the signs of the 2nd and 3rd terms on the RHS: (i) A bulk source velocity (which is defined as >0 for the receding direction) should be added to the material velocity field; this is equivalent to subtracting it from deltav, as here. (ii) A positive value of lineRedShift means the line is red-shifted wrt to the frequency specified for the image. The effect is the same as if the line and image frequencies were the same, but the bulk recession velocity were higher. lineRedShift should thus be added to the recession velocity, which is equivalent to subtracting it from deltav, as here. */ + + /* It appears to be necessary to sample the velocity function in the following way rather than interpolating it from the vertices of the Delaunay cell in the same way as with all the other quantities of interest. Velocity varies too much across the cells, and in a nonlinear way, for linear interpolation to yield a totally satisfactory result. + */ + velocity(gips[2].x[0], gips[2].x[1], gips[2].x[2], vel); + velCmpntRay = veloproject(dir, vel); + calcLineAmpInterp(velCmpntRay, gips[2].mol[molI].binv, deltav, &vfac); + + /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. + */ + sourceFunc_line_raytrace(md[molI], vfac, gips[2].mol[molI], lineI, &jnu, &alpha); + } + } /* end loop over all lines. */ + + dtau = alpha*ds; + calcSourceFn(dtau, par, &remnantSnu, &expDTau); + remnantSnu *= jnu*md[0].norminv*ds; +#ifdef FASTEXP + brightnessIncrement = FastExp(ray.tau[ichan])*remnantSnu; +#else + brightnessIncrement = exp(-ray.tau[ichan])*remnantSnu; +#endif + ray.intensity[ichan] += brightnessIncrement; + ray.tau[ichan] += dtau; + } /* End loop over channels. */ + } /* End if(par->polarization). */ + } /* End loop over segments within cell. */ + + entryI = exitI; + exitI = 1 - exitI; + } /* End loop over cells in the chain traversed by the ray. */ + + /* Add or subtract cmb. */ +#ifdef FASTEXP + for(ichan=0;ichan Date: Fri, 8 Apr 2016 17:37:25 +0200 Subject: [PATCH 29/43] Sink point values set In the previous raytracing algorithm the values of the sink point atributes are not accessed. However they are accessed in the interpolation routine. Because of this it is necessary to be a bit more careful about setting them. Thus I have set their: . nmol to zero in LTE(); . pops to L.T.E with the CMB in molinit(); . vel to the user-specified values (instead of zero) in getVelosplines(). --- src/LTEsolution.c | 3 +++ src/molinit.c | 14 +++++++++++++- src/smooth.c | 2 +- src/velospline.c | 4 +++- 4 files changed, 20 insertions(+), 3 deletions(-) diff --git a/src/LTEsolution.c b/src/LTEsolution.c index aa1cbeb..6285458 100644 --- a/src/LTEsolution.c +++ b/src/LTEsolution.c @@ -25,6 +25,9 @@ LTE(inputPars *par, struct grid *g, molData *m){ g[id].mol[ispec].pops[ilev]=m[ispec].gstat[ilev]*exp(-100*CLIGHT*HPLANCK*m[ispec].eterm[ilev]/(KBOLTZ*g[id].t[0]))/z; } } + for(id=par->pIntensity;idncell;id++){ + g[id].mol[ispec].nmol=0.0; + } } if(par->outputfile) popsout(par,g,m); } diff --git a/src/molinit.c b/src/molinit.c index 845cd88..8132a6b 100644 --- a/src/molinit.c +++ b/src/molinit.c @@ -105,7 +105,7 @@ void molinit(molData *m, inputPars *par, struct grid *g,int i){ int id, ilev, iline, itrans, ispec, itemp, *ntemp, tnint=-1, idummy, ipart, *count,flag=0; char *collpartname[] = {"H2","p-H2","o-H2","electrons","H","He","H+"}; /* definition from LAMDA */ - double fac, uprate, downrate=0, dummy, amass; + double fac, uprate, downrate=0, dummy, amass, sum; struct data { double *colld, *temp; } *part; char string[200], specref[90], partstr[90]; @@ -303,6 +303,18 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){ for(ilev=0;ilevpIntensity;idncell;id++){ + sum = 0.0; + for(ilev=0;ilevtcmb)); + sum += g[id].mol[i].pops[ilev]; + } + for(ilev=0;ilevncell, 0, &dc, &numCells); + delaunay(DIM, g, (unsigned long)par->ncell, 0, &dc, &numCells); distCalc(par, g); if(!silent) progressbar((double)(sg+1)/(double)smooth, 5); if(dc!=NULL) free(dc); diff --git a/src/velospline.c b/src/velospline.c index 3e063bd..1449eaa 100644 --- a/src/velospline.c +++ b/src/velospline.c @@ -53,12 +53,14 @@ getVelosplines(inputPars *par, struct grid *g){ } for(i=par->pIntensity;incell;i++){ + /* Set velocity values also for sink points (otherwise Delaunay ray-tracing has problems) */ + velocity(g[i].x[0],g[i].x[1],g[i].x[2],g[i].vel); + g[i].a0=malloc(g[i].numNeigh*sizeof(double)); g[i].a1=malloc(g[i].numNeigh*sizeof(double)); g[i].a2=malloc(g[i].numNeigh*sizeof(double)); g[i].a3=malloc(g[i].numNeigh*sizeof(double)); g[i].a4=malloc(g[i].numNeigh*sizeof(double)); - for(j=0;j<3;j++) g[i].vel[j]=0.; for(j=0;j Date: Fri, 8 Apr 2016 17:57:54 +0200 Subject: [PATCH 30/43] Added interpolation machinery to raytrace(). With this commit, the new algorithm is essentially complete. I'll make 1 more with some minor cosmetic changes. --- src/raytrace.c | 78 ++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 70 insertions(+), 8 deletions(-) diff --git a/src/raytrace.c b/src/raytrace.c index 2698755..35fe8de 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -437,11 +437,19 @@ At the moment I will fix the number of segments, but it might possibly be faster /*....................................................................*/ void raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ - int *counta, *countb,nlinetot,aa,ppi,molI,li; - int ichan,px,iline,tmptrans,i,threadI,nRaysDone; - double size,minfreq,absDeltaFreq,totalNumPixelsMinus1=(double)(img[im].pxls*img[im].pxls-1); - double cutoff; + const int numFaces=1+DIM, numInterpPoints=3, numSegments=5; + const double oneOnNumSegments = 1.0/(double)numSegments, oneOnNFaces=1.0/(double)numFaces; + const double epsilon = 1.0e-6; // Needs thinking about. Double precision is much smaller than this. + const double oneOnNAlias = 1.0/(double)par->antialias; + const double oneOnTotalNumPixelsMinus1=1.0/(double)(img[im].pxls*img[im].pxls-1); + int *counta, *countb,nlinetot,aa,ii,ppi; + int ichan,px,iline,tmptrans,i,threadI,nRaysDone,molI,di,vi,li; + double size,minfreq,absDeltaFreq; + double cutoff,sum,progress; const gsl_rng_type *ranNumGenType = gsl_rng_ranlxs2; + struct cell *dc=NULL; + unsigned long numCells, dci; + static double lastProgress = 0.0; struct gAuxType *gAux=NULL; /* This will hold some precalculated values for the grid points. */ gsl_rng *ran = gsl_rng_alloc(ranNumGenType); /* Random number generator */ @@ -461,6 +469,28 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ size=img[im].distance*img[im].imgres; + if(par->traceRayAlgorithm==1){ + delaunay(DIM, g, (unsigned long)par->ncell, 1, &dc, &numCells); + + /* We need to process the list of cells a bit further - calculate their centres, and reset the id values to be the same as the index of the cell in the list. (This last because we are going to construct other lists to indicate which cells have been visited etc.) + */ + for(dci=0;dcix[di]; + } + dc[dci].centre[di] = sum*oneOnNFaces; + } + + dc[dci].id = dci; + } + + }else if(par->traceRayAlgorithm!=0){ + if(!silent) bail_out("Unrecognized value of par.traceRayAlgorithm"); + exit(1); + } + /* Precalculate binv*nmol*pops for all grid points. */ gAux = malloc(sizeof(*gAux)*par->ncell); @@ -526,9 +556,24 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ /* Declaration of thread-private pointers. */ rayData ray; + gridInterp gips[numInterpPoints]; + ray.intensity=malloc(sizeof(double) * img[im].nchan); ray.tau=malloc(sizeof(double) * img[im].nchan); + if(par->traceRayAlgorithm==1){ + /* Allocate memory for the interpolation points: + */ + for(ii=0;iinSpecies); + for(molI=0;molInSpecies;molI++){ + gips[ii].mol[molI].specNumDens = malloc(sizeof(*(gips[ii].mol[molI].specNumDens))*m[molI].nlev); + gips[ii].mol[molI].dust = malloc(sizeof(*(gips[ii].mol[molI].dust)) *m[molI].nline); + gips[ii].mol[molI].knu = malloc(sizeof(*(gips[ii].mol[molI].knu)) *m[molI].nline); + } + } + } + #pragma omp for /* Main loop through pixel grid. */ for(px=0;px<(img[im].pxls*img[im].pxls);px++){ @@ -539,27 +584,44 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ ray.x = size*(gsl_rng_uniform(threadRans[threadI])+px%img[im].pxls)-size*img[im].pxls/2.; ray.y = size*(gsl_rng_uniform(threadRans[threadI])+px/img[im].pxls)-size*img[im].pxls/2.; - traceray(ray, par, tmptrans, img, im, g, gAux, m, nlinetot, counta, countb, cutoff); + if(par->traceRayAlgorithm==0){ + traceray(ray, par, tmptrans, img, im, g, gAux, m, nlinetot, counta, countb, cutoff); + }else if(par->traceRayAlgorithm==1) + traceray_smooth(ray, par, tmptrans, img, im, g, gAux, m, nlinetot, counta, countb\ + , dc, numCells, epsilon, gips, numSegments, oneOnNumSegments); #pragma omp critical { for(ichan=0;ichanantialias; - img[im].pixel[px].tau[ichan] += ray.tau[ichan]/(double) par->antialias; + img[im].pixel[px].intense[ichan] += ray.intensity[ichan]*oneOnNAlias; + img[im].pixel[px].tau[ichan] += ray.tau[ ichan]*oneOnNAlias; } } } if (threadI == 0){ /* i.e., is master thread */ - if(!silent) progressbar((double)(nRaysDone)/totalNumPixelsMinus1, 13); + if(!silent) { + progress = ((double)nRaysDone)*oneOnTotalNumPixelsMinus1; + if(progress-lastProgress>0.002){ + lastProgress = progress; + progressbar(progress, 13); + } + } } } + if(par->traceRayAlgorithm==1){ + free(dc); + for(ii=0;iinSpecies, gips[ii].mol); + } + } free(ray.tau); free(ray.intensity); } /* End of parallel block. */ img[im].trans=tmptrans; + freeGAux((unsigned long)par->ncell, par->nSpecies, gAux); free(counta); free(countb); for (i=0;inThreads;i++){ From f8b24efd5b6e527d61291c996ed8a1d2dec361f7 Mon Sep 17 00:00:00 2001 From: ims Date: Mon, 11 Apr 2016 12:11:19 +0200 Subject: [PATCH 31/43] Minor cosmetic changes. Here I have: . Improved some variable names in raytrace(). . Ordered all the function templates alphabetically (or at least, more so) in lime.h. . Removed some tab indents in predefgrid.c which were bugging me. --- src/lime.h | 61 ++++++++++++++++++++++++------------------------ src/predefgrid.c | 16 ++++++------- src/raytrace.c | 22 ++++++++--------- 3 files changed, 49 insertions(+), 50 deletions(-) diff --git a/src/lime.h b/src/lime.h index 3bb465f..83436cf 100644 --- a/src/lime.h +++ b/src/lime.h @@ -218,19 +218,8 @@ typedef struct { } triangle2D; -int followRayThroughDelCells(double*, double*, struct grid*, struct cell*, const unsigned long, const double, intersectType*, unsigned long**, intersectType**, int*); -int buildRayCellChain(double*, double*, struct grid*, struct cell*, _Bool**, unsigned long, int, int, int, const double, unsigned long**, intersectType**, int*); -int getNewEntryFaceI(const unsigned long, const struct cell); -faceType extractFace(struct grid*, struct cell*, const unsigned long, const int); -void intersectLineTriangle(double*, double*, faceType, intersectType*); -triangle2D calcTriangle2D(faceType face); -void doBaryInterp(const intersectType, struct grid*, struct gAuxType*, double*, unsigned long*, molData*, const int, gridInterp*); -void doSegmentInterp(gridInterp*, const int, molData*, const int, const double, const int); -void freePop2(const int, struct pop2*); -void freeGAux(const unsigned long, const int, struct gAuxType*); - -/* Some functions */ +/* User-specifiable functions */ void density(double,double,double,double *); void temperature(double,double,double,double *); void abundance(double,double,double,double *); @@ -243,26 +232,45 @@ void gasIIdust(double,double,double,double *); void binpopsout(inputPars *, struct grid *, molData *); void buildGrid(inputPars *, struct grid *); -void calcSourceFn(double dTau, const inputPars *par, double *remnantSnu, double *expDTau); -void continuumSetup(int, image *, molData *, inputPars *, struct grid *); -void distCalc(inputPars *, struct grid *); +int buildRayCellChain(double*, double*, struct grid*, struct cell*, _Bool**, unsigned long, int, int, int, const double, unsigned long**, intersectType**, int*); +void calcFastExpRange(const int, const int, int*, int*, int*); +void calcLineAmpInterp(const double, const double, const double, double*); +void calcLineAmpLinear(struct grid*, const int, const int, const double, const double, double*); +void calcLineAmpSample(double*, double*, const double, const double, const double, double*); +void calcLineAmpSpline(struct grid*, const int, const int, const double, const double, double*); +void calcSourceFn(double, const inputPars*, double*, double*); +void calcTableEntries(const int, const int); +triangle2D calcTriangle2D(faceType); +void continuumSetup(int, image*, molData*, inputPars*, struct grid*); +void delaunay(const int, struct grid*, const unsigned long, const _Bool, struct cell**, unsigned long*); +void distCalc(inputPars*, struct grid*); +void doBaryInterp(const intersectType, struct grid*, struct gAuxType*, double*, unsigned long*, molData*, const int, gridInterp*); +void doSegmentInterp(gridInterp*, const int, molData*, const int, const double, const int); +faceType extractFace(struct grid*, struct cell*, const unsigned long, const int); +int factorial(const int n); +//inline double FastExp(const float negarg); 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 freeInput(inputPars *, image*, molData* m ); -void freeGrid(const inputPars * par, const molData* m, struct grid * g); -void freePopulation(const inputPars * par, const molData* m, struct populations * pop); +int followRayThroughDelCells(double*, double*, struct grid*, struct cell*, const unsigned long, const double, intersectType*, unsigned long**, intersectType**, int*); +void freeInput(inputPars*, image*, molData* m ); +void freeGAux(const unsigned long, const int, struct gAuxType*); +void freeGrid(const inputPars*, const molData*, struct grid*); +void freePopulation(const inputPars*, const molData*, struct populations*); +void freePop2(const int, struct pop2*); double gaussline(double, double); void getArea(inputPars *, 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 getmatrix(int, gsl_matrix *, molData *, struct grid *, int, gridPointData *); -void getclosest(double, double, double, long *, long *, double *, double *, double *); +int getNewEntryFaceI(const unsigned long, const struct cell); void getVelosplines(inputPars *, struct grid *); void getVelosplines_lin(inputPars *, struct grid *); void gridAlloc(inputPars *, struct grid **); +void input(inputPars *, image *); +void intersectLineTriangle(double*, double*, faceType, intersectType*); +float invSqrt(float); void kappa(molData *, struct grid *, inputPars *,int); void levelPops(molData *, inputPars *, struct grid *, int *); void line_plane_intersect(struct grid *, double *, int , int *, double *, double *, double); @@ -271,7 +279,6 @@ void lineCount(int,molData *,int **, int **, int *); void LTE(inputPars *, struct grid *, molData *); void molinit(molData *, inputPars *, struct grid *,int); void openSocket(inputPars *par, int); -void delaunay(const int, struct grid *, const unsigned long, const _Bool, struct cell **, unsigned long *); void photon(int, struct grid *, molData *, int, const gsl_rng *,inputPars *,blend *,gridPointData *,double *); void parseInput(inputPars *, image **, molData **); double planckfunc(int, double, molData *, int); @@ -293,20 +300,12 @@ void sourceFunc_pol(const double, const double*, const molData, const struct pop void stateq(int, struct grid *, molData *, int, inputPars *,gridPointData *,double *); void statistics(int, molData *, struct grid *, int *, double *, double *, int *); void stokesangles(const double B[3], const double, double *); +double taylor(const int maxOrder, const float x); void traceray(rayData, inputPars*, const int, image*, const int, struct grid*, struct gAuxType*, molData*, const int, int*, int*, const double); void traceray_smooth(rayData, inputPars*, const int, image*, const int, struct grid*, struct gAuxType*, molData*, const int, int*, int*, struct cell*, const unsigned long, const double, gridInterp gips[3], const int, const double); -void calcLineAmpInterp(const double, const double, const double, double*); -void calcLineAmpLinear(struct grid*, const int, const int, const double, const double, double*); -void calcLineAmpSample(double*, double*, const double, const double, const double, double*); -void calcLineAmpSpline(struct grid*, const int, const int, const double, const double, double*); double veloproject(double *, double *); void writefits(int, inputPars *, molData *, image *); void write_VTK_unstructured_Points(inputPars *, struct grid *); -int factorial(const int n); -double taylor(const int maxOrder, const float x); -void calcFastExpRange(const int maxTaylorOrder, const int maxNumBitsPerMantField, int *numMantissaFields, int *lowestExponent, int *numExponentsUsed); -void calcTableEntries(const int maxTaylorOrder, const int maxNumBitsPerMantField); -//inline double FastExp(const float negarg); /* Curses functions */ diff --git a/src/predefgrid.c b/src/predefgrid.c index f278891..fe85204 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -42,14 +42,14 @@ predefinedGrid(inputPars *par, struct grid *g){ g[i].sink=0; - g[i].t[1]=g[i].t[0]; - g[i].mol[0].nmol=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); + g[i].t[1]=g[i].t[0]; + g[i].mol[0].nmol=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); } for(i=par->pIntensity;incell;i++){ diff --git a/src/raytrace.c b/src/raytrace.c index 35fe8de..7e1e84b 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -442,8 +442,8 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ const double epsilon = 1.0e-6; // Needs thinking about. Double precision is much smaller than this. const double oneOnNAlias = 1.0/(double)par->antialias; const double oneOnTotalNumPixelsMinus1=1.0/(double)(img[im].pxls*img[im].pxls-1); - int *counta, *countb,nlinetot,aa,ii,ppi; - int ichan,px,iline,tmptrans,i,threadI,nRaysDone,molI,di,vi,li; + int *allLineMolIs,*allLineLineIs,nlinetot,aa,ii,ppi; + int ichan,px,iline,tmptrans,i,threadI,nPixelsDone,molI,di,vi,li; double size,minfreq,absDeltaFreq; double cutoff,sum,progress; const gsl_rng_type *ranNumGenType = gsl_rng_ranlxs2; @@ -515,7 +515,7 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ } /* Determine whether there are blended lines or not. */ - lineCount(par->nSpecies, m, &counta, &countb, &nlinetot); + lineCount(par->nSpecies, m, &allLineMolIs, &allLineLineIs, &nlinetot); if(img[im].doline==0) nlinetot=1; /* Fix the image parameters. */ @@ -548,9 +548,9 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ } } - nRaysDone=0; + nPixelsDone=0; omp_set_dynamic(0); - #pragma omp parallel private(px,aa,threadI) num_threads(par->nThreads) + #pragma omp parallel private(px,aa,threadI,ii) num_threads(par->nThreads) { threadI = omp_get_thread_num(); @@ -578,16 +578,16 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ /* Main loop through pixel grid. */ for(px=0;px<(img[im].pxls*img[im].pxls);px++){ #pragma omp atomic - ++nRaysDone; + ++nPixelsDone; for(aa=0;aaantialias;aa++){ ray.x = size*(gsl_rng_uniform(threadRans[threadI])+px%img[im].pxls)-size*img[im].pxls/2.; ray.y = size*(gsl_rng_uniform(threadRans[threadI])+px/img[im].pxls)-size*img[im].pxls/2.; if(par->traceRayAlgorithm==0){ - traceray(ray, par, tmptrans, img, im, g, gAux, m, nlinetot, counta, countb, cutoff); + traceray(ray, par, tmptrans, img, im, g, gAux, m, nlinetot, allLineMolIs, allLineLineIs, cutoff); }else if(par->traceRayAlgorithm==1) - traceray_smooth(ray, par, tmptrans, img, im, g, gAux, m, nlinetot, counta, countb\ + traceray_smooth(ray, par, tmptrans, img, im, g, gAux, m, nlinetot, allLineMolIs, allLineLineIs\ , dc, numCells, epsilon, gips, numSegments, oneOnNumSegments); #pragma omp critical @@ -600,7 +600,7 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ } if (threadI == 0){ /* i.e., is master thread */ if(!silent) { - progress = ((double)nRaysDone)*oneOnTotalNumPixelsMinus1; + progress = ((double)nPixelsDone)*oneOnTotalNumPixelsMinus1; if(progress-lastProgress>0.002){ lastProgress = progress; progressbar(progress, 13); @@ -622,8 +622,8 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ img[im].trans=tmptrans; freeGAux((unsigned long)par->ncell, par->nSpecies, gAux); - free(counta); - free(countb); + free(allLineMolIs); + free(allLineLineIs); for (i=0;inThreads;i++){ gsl_rng_free(threadRans[i]); } From 441ebd8000f9b9c90dffff76fbe5db97c265a53a Mon Sep 17 00:00:00 2001 From: ims Date: Wed, 20 Apr 2016 15:30:05 +0200 Subject: [PATCH 32/43] Smarter call to velocity in raytrace This improves the speed of raytracing by implementing a suggestion made by T Lunttila in the discussion of issue #95. In the present master branch, the user-specified function velocity() is called for every channel/line combination judged to have a significant line amplitude; but of course the velocity depends only on position in the model. The function is actually called N times (currently N=10) at N spatially separated points along the path of the ray from entry to exit of a particular Voronoi cell. In the present commit, the N velocity values have been precalculated just once per cell. This speeds raytracing up on my desktop machine from 25 sec to 11 sec. In addition, I have - cleaned up some divisions by int within velocityspline2(); - tightened up the function interface to velocityspline2(); - removed the raytrace_1.4 code, because (i) his is available elsewhere, (ii) I'll try to distil it into a PR soon, and (iii) it was generating errors now due to an altered traceray() interface. --- src/lime.h | 5 +- src/raytrace.c | 299 +++++++------------------------------------------ 2 files changed, 41 insertions(+), 263 deletions(-) diff --git a/src/lime.h b/src/lime.h index f1d602c..a50e13b 100644 --- a/src/lime.h +++ b/src/lime.h @@ -236,9 +236,10 @@ void sourceFunc_pol(double *, double *, double, molData *, double, struct gri void stateq(int, struct grid *, molData *, int, inputPars *,gridPointData *,double *); void statistics(int, molData *, struct grid *, int *, double *, double *, int *); void stokesangles(double, double, double, double, double *); -void traceray(rayData, int, int, inputPars *, struct grid *, molData *, image *, int, int *, int *, double); +void traceray(rayData, int, int, inputPars*, struct grid*, molData*, image*, int, int*, int*, double, const int, const double); void velocityspline(struct grid *, int, int, double, double, double*); -void velocityspline2(double *, double *, double, double, double, double*); +void velocityspline2(const double*, const double*, const double\ + , const double, double*, const int, const double, const double, double*); double veloproject(double *, double *); void writefits(int, inputPars *, molData *, image *); void write_VTK_unstructured_Points(inputPars *, struct grid *); diff --git a/src/raytrace.c b/src/raytrace.c index 8f28320..9a7a191 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -10,19 +10,20 @@ #include "lime.h" -void -velocityspline2(double x[3], double dx[3], double ds, double binv, double deltav, double *vfac){ +/*....................................................................*/ +void velocityspline2(const double x[3], const double dx[3], const double ds\ + , const double binv, double *projVels, const int nSteps\ + , const double oneOnNSteps, const double deltav, double *vfac){ /* The bulk velocity of the model material can vary significantly with position, thus so can the value of the line-shape function at a given frequency and direction. The present function calculates 'vfac', an approximate average of the line-shape function along a path of length ds in the direction of the line of sight. */ - int i,steps=10; + int i; double v,d,val,vel[3]; *vfac=0.; - for(i=0;iradiusSqu-(xp*xp+yp*yp)); /* There are two points of intersection between the line of sight and the spherical model surface; this is the Z coordinate (in the unrotated frame) of the one nearer to the observer. */ /* Rotate the line of sight as desired. */ - for(i=0;i<3;i++){ + for(i=0;ipolarization){ for(ichan=0;ichanpregrid){ + for(i=0;i img[im].freq-img[im].bandwidth/2. + if(img[im].doline\ + && m[molI].freq[lineI] > img[im].freq-img[im].bandwidth/2. && m[molI].freq[lineI] < img[im].freq+img[im].bandwidth/2.){ /* Calculate the red shift of the transition wrt to the frequency specified for the image. */ if(img[im].trans > -1){ @@ -144,8 +158,10 @@ Note that the algorithm employed here is similar to that employed in the functio /* Line centre occurs when deltav = the recession velocity of the radiating material. Explanation of the signs of the 2nd and 3rd terms on the RHS: (i) A bulk source velocity (which is defined as >0 for the receding direction) should be added to the material velocity field; this is equivalent to subtracting it from deltav, as here. (ii) A positive value of lineRedShift means the line is red-shifted wrt to the frequency specified for the image. The effect is the same as if the line and image frequencies were the same, but the bulk recession velocity were higher. lineRedShift should thus be added to the recession velocity, which is equivalent to subtracting it from deltav, as here. */ /* Calculate an approximate average line-shape function at deltav within the Voronoi cell. */ - if(!par->pregrid) velocityspline2(x,dx,ds,g[posn].mol[molI].binv,deltav,&vfac); - else vfac=gaussline(deltav+veloproject(dx,g[posn].vel),g[posn].mol[molI].binv); + if(!par->pregrid) + velocityspline2(x,dx,ds,g[posn].mol[molI].binv,projVels,nSteps,oneOnNSteps,deltav,&vfac); + else + vfac = gaussline(deltav+veloproject(dx,g[posn].vel),g[posn].mol[molI].binv); /* Increment jnu and alpha for this Voronoi cell by the amounts appropriate to the spectral line. */ sourceFunc_line(&jnu,&alpha,m,vfac,g,posn,molI,lineI); @@ -156,6 +172,7 @@ Note that the algorithm employed here is similar to that employed in the functio else if(img[im].doline && img[im].trans == -1) sourceFunc_cont(&jnu,&alpha,g,posn,0,tmptrans); else sourceFunc_cont(&jnu,&alpha,g,posn,0,0); dtau=alpha*ds; +//??? if(dtau < -30) dtau = -30; // as in photon()? calcSourceFn(dtau, par, &remnantSnu, &expDTau); remnantSnu *= jnu*m[0].norminv*ds; #ifdef FASTEXP @@ -187,13 +204,15 @@ 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){ +/*....................................................................*/ +void raytrace(int im, inputPars *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); double cutoff; const gsl_rng_type *ranNumGenType = gsl_rng_ranlxs2; + const int nStepsThruCell=10; + const double oneOnNSteps=1.0/(double)nStepsThruCell; gsl_rng *ran = gsl_rng_alloc(ranNumGenType); /* Random number generator */ #ifdef TEST @@ -267,7 +286,7 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ ray.x = size*(gsl_rng_uniform(threadRans[threadI])+px%img[im].pxls)-size*img[im].pxls/2.; ray.y = size*(gsl_rng_uniform(threadRans[threadI])+px/img[im].pxls)-size*img[im].pxls/2.; - traceray(ray, tmptrans, im, par, g, m, img, nlinetot, counta, countb, cutoff); + traceray(ray, tmptrans, im, par, g, m, img, nlinetot, counta, countb, cutoff, nStepsThruCell, oneOnNSteps); #pragma omp critical { @@ -297,245 +316,3 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ gsl_rng_free(ran); } - -void -raytrace_1_4(int im, inputPars *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. -I've adapted it slightly so it makes use of the function traceray(), -which was modified from the function tracerays() in v1.4. This algorithm -is not currently used, but may be useful as an option; that's why I have -kept it. - */ - - int *counta, *countb,nlinetot; - int ichan,i,px,iline,tmptrans,count; - double size,xp,yp,minfreq,absDeltaFreq; - double cutoff; - - gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2); /* Random number generator */ -#ifdef TEST - gsl_rng_set(ran,178490); -#else - gsl_rng_set(ran,time(0)); -#endif - rayData *rays; - - int sg,n; - double cx,cy; - - double x1,x2,x3,y1,y2,y3,z1,z2,z3,xt[3],yt[3],di,p,d1,d2,d3,temp1; - int zt[3]; - int c; - - char flags[255]; - boolT ismalloc = False; - facetT *facet, *neighbor, **neighborp;; - vertexT *vertex,**vertexp; - coordT *pt_array; - - int id; - coordT point[3]; - boolT isoutside; - realT bestdist; - - size=img[im].distance*img[im].imgres; - - /* Determine whether there are blended lines or not */ - lineCount(par->nSpecies, m, &counta, &countb, &nlinetot); - if(img[im].doline==0) nlinetot=1; - - /* Fix the image parameters */ - if(img[im].freq < 0) img[im].freq=m[0].freq[img[im].trans]; - if(img[im].nchan == 0 && img[im].bandwidth>0){ - img[im].nchan=(int) (img[im].bandwidth/(img[im].velres/CLIGHT*img[im].freq)); - } else if (img[im].velres<0 && img[im].bandwidth>0){ - img[im].velres = img[im].bandwidth*CLIGHT/img[im].freq/img[im].nchan; - } else img[im].bandwidth = img[im].nchan*img[im].velres/CLIGHT * img[im].freq; - - if(img[im].trans<0){ - iline=0; - minfreq=fabs(img[im].freq-m[0].freq[iline]); - tmptrans=iline; - for(iline=1;ilinepIntensity)); - - for(i=0;ipIntensity;i++){ - rays[i].x=g[i].x[0]; - rays[i].y=g[i].x[1]; - rays[i].tau=malloc(sizeof(double) * img[im].nchan); - rays[i].intensity=malloc(sizeof(double) * img[im].nchan); - for(ichan=0;ichanpIntensity); - - for(i=0;ipIntensity;i++) { - pt_array[i*2+0]=rays[i].x; - pt_array[i*2+1]=rays[i].y; - } - - sprintf(flags,"qhull v s Qbb T0"); - if (!qh_new_qhull(2, par->pIntensity, pt_array, ismalloc, flags, NULL, NULL)) { - - qh_setvoronoi_all(); - - FORALLvertices { - i=qh_pointid(vertex->point); - - cx=0.; - cy=0.; - n=0; - FOREACHneighbor_(vertex) { - if (!neighbor->upperdelaunay) n++; - } - if(n>0){ - - - } else { - if(!silent) bail_out("Qhull error"); - exit(0); - } - - FOREACHneighbor_(vertex) { - if (!neighbor->upperdelaunay) { - cx+=neighbor->center[0]; - cy+=neighbor->center[1]; - } - } - - rays[i].x = rays[i].x - (rays[i].x-cx/ (double) n)*0.1; - rays[i].y = rays[i].y - (rays[i].y-cy/ (double) n)*0.1; - } - } else { - printf("qhull error\n"); - } - - qh_freeqhull(!qh_ALL); - free(pt_array); - } - - cutoff = par->minScale*1.0e-7; - - /* Main loop through rays */ - count=0; - for(px=0;pxpIntensity;px++){ - traceray(rays[px], tmptrans, im, par, g, m, img, nlinetot, counta, countb, cutoff); - ++count; - if(!silent) progressbar((double)(count)/(double)(par->pIntensity-1), 13); - } - - /* Remap rays onto pixel grid */ - pt_array=malloc(2*sizeof(coordT)*par->pIntensity); - - for(i=0;ipIntensity;i++) { - pt_array[i*2+0]=rays[i].x; - pt_array[i*2+1]=rays[i].y; - } - -/* This allocation belongs to "Shepard's method" below - d=malloc(sizeof(double)*par->pIntensity); -*/ - size=img[im].distance*img[im].imgres; - - sprintf(flags,"qhull d Qbb"); - if (!qh_new_qhull(2, par->pIntensity, pt_array, ismalloc, flags, NULL, NULL)) { - for(px=0;pxpIntensity;i++){ - // d[i]=1./pow(sqrt(pow(xp-rays[i].x,2)+ pow(yp-rays[i].y,2)),8.); - temp1 = (xp-rays[i].x)*(xp-rays[i].x)+ (yp-rays[i].y)*(yp-rays[i].y) - d[i]=1./(temp1*temp1*temp1*temp1); - img[im].pixel[px].intense[ichan] += rays[i].intensity[ichan]*d[i]; -**** how to handle img[im].pixel[px].tau[ichan]? - di+=d[i]; - } - img[im].pixel[px].intense[ichan] /= di; - } -*/ - - - point[0]=xp; - point[1]=yp; - point[2]=0.; - - qh_setdelaunay (3, 1, point); - facet= qh_findbestfacet (point, qh_ALL, &bestdist, &isoutside); - - c=0; - FOREACHvertex_( facet->vertices ) { - id=qh_pointid(vertex->point); - xt[c]=rays[id].x; yt[c]=rays[id].y; zt[c]=id; - c++; - } - - - x1=xt[0];x2=xt[1];x3=xt[2]; - y1=yt[0];y2=yt[1];y3=yt[2]; - - for(ichan=0;ichanpIntensity;i++){ - free(rays[i].tau); - free(rays[i].intensity); - } - free(rays); - free(counta); - free(countb); -} - - From 106bd830bea5432da404e6d2ca53f09dce51d481 Mon Sep 17 00:00:00 2001 From: ims Date: Wed, 20 Apr 2016 16:35:39 +0200 Subject: [PATCH 33/43] Smarter call to sourceFunc_cont The call to this function has been moved outside the loop over channels in traceray(). It is only necessary to call it once per channel, then just keep the returned values. The change does not speed things up much, possibly because sourceFunc_cont() is not actually doing many FPOs. I also introduced some variables to make traceray() a leeetle more self-documenting and a leeetle less obscure. Lot of work still to go in this direction in LIME as a whole. :-/ A final change was removing the unused argument 'vfac' from sourceFunc_pol(). --- src/lime.h | 2 +- src/raytrace.c | 44 ++++++++++++++++++++++++++++++-------------- src/sourcefunc.c | 2 +- 3 files changed, 32 insertions(+), 16 deletions(-) diff --git a/src/lime.h b/src/lime.h index a50e13b..46ece5b 100644 --- a/src/lime.h +++ b/src/lime.h @@ -232,7 +232,7 @@ int sortangles(double *, int, struct grid *, const gsl_rng *); void sourceFunc(double *, double *, double, molData *,double,struct grid *,int,int, int,int); void sourceFunc_line(double *,double *,molData *, double, struct grid *, int, int,int); void sourceFunc_cont(double *,double *, struct grid *, int, int,int); -void sourceFunc_pol(double *, double *, double, molData *, double, struct grid *, int, int, int, double); +void sourceFunc_pol(double *, double *, double, molData *, struct grid *, int, int, int, double); void stateq(int, struct grid *, molData *, int, inputPars *,gridPointData *,double *); void statistics(int, molData *, struct grid *, int *, double *, double *, int *); void stokesangles(double, double, double, double, double *); diff --git a/src/raytrace.c b/src/raytrace.c index 9a7a191..88e8dbd 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -78,10 +78,10 @@ For a given image pixel position, this function evaluates the intensity of the t Note that the algorithm employed here is similar to that employed in the function photon() which calculates the average radiant flux impinging on a grid cell: namely the notional photon is started at the side of the model near the observer and 'propagated' in the receding direction until it 'reaches' the far side. This is rather non-physical in conception but it makes the calculation easier. */ - int ichan,posn,nposn,i,iline,molI,lineI; + int ichan,posn,nposn,i,iline,molI,lineI,contMolI,contLineI,polMolI,polLineI,stokesId; double vfac=0.,x[DIM],dx[DIM],vThisChan; double deltav,ds,dist2,ndist2,xp,yp,zp,col,lineRedShift,jnu,alpha,remnantSnu,dtau,expDTau,snu_pol[3]; - double projVels[nSteps],d,vel[DIM]; + double projVels[nSteps],d,vel[DIM],contJnu,contAlpha; for(ichan=0;ichan -1) + contLineI = img[im].trans; + else if(img[im].doline && img[im].trans == -1) + contLineI = tmptrans; + else + contLineI = 0; + /* Find the grid point nearest to the starting x. */ i=0; dist2=(x[0]-g[i].x[0])*(x[0]-g[i].x[0]) + (x[1]-g[i].x[1])*(x[1]-g[i].x[1]) + (x[2]-g[i].x[2])*(x[2]-g[i].x[2]); @@ -117,15 +126,18 @@ Note that the algorithm employed here is similar to that employed in the functio ds=-2.*zp-col; /* This default value is chosen to be as large as possible given the spherical model boundary. */ nposn=-1; line_plane_intersect(g,&ds,posn,&nposn,dx,x,cutoff); /* Returns a new ds equal to the distance to the next Voronoi face, and nposn, the ID of the grid cell that abuts that face. */ + if(par->polarization){ - for(ichan=0;ichanpregrid){ @@ -136,9 +148,15 @@ Note that the algorithm employed here is similar to that employed in the functio } } + /* Calculate first the continuum stuff because it is the same for all channels: + */ + contJnu = 0.0; + contAlpha = 0.0; + sourceFunc_cont(&contJnu,&contAlpha,g,posn,contMolI,contLineI); + for(ichan=0;ichan img[im].freq-img[im].bandwidth/2. && m[molI].freq[lineI] < img[im].freq+img[im].bandwidth/2.){ - /* Calculate the red shift of the transition wrt to the frequency specified for the image. */ + /* Calculate the red shift of the transition wrt to the frequency specified for the image. + */ if(img[im].trans > -1){ lineRedShift=(m[molI].freq[img[im].trans]-m[molI].freq[lineI])/m[molI].freq[img[im].trans]*CLIGHT; } else { @@ -168,9 +187,6 @@ Note that the algorithm employed here is similar to that employed in the functio } } - if(img[im].doline && img[im].trans > -1) sourceFunc_cont(&jnu,&alpha,g,posn,0,img[im].trans); - else if(img[im].doline && img[im].trans == -1) sourceFunc_cont(&jnu,&alpha,g,posn,0,tmptrans); - else sourceFunc_cont(&jnu,&alpha,g,posn,0,0); dtau=alpha*ds; //??? if(dtau < -30) dtau = -30; // as in photon()? calcSourceFn(dtau, par, &remnantSnu, &expDTau); @@ -185,7 +201,7 @@ Note that the algorithm employed here is similar to that employed in the functio } /* Move the working point to the edge of the next Voronoi cell. */ - for(i=0;i<3;i++) x[i]+=ds*dx[i]; + for(i=0;i Date: Fri, 24 Jun 2016 11:10:40 +0200 Subject: [PATCH 34/43] Raytrace loop dynamically scheduled Following a suggestion from @tlunttil, the parallel for loop in raytrace() is now scheduled dynamically. --- src/raytrace.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/raytrace.c b/src/raytrace.c index 08cb4f0..7b4717d 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -292,7 +292,7 @@ void raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ ray.intensity=malloc(sizeof(double) * img[im].nchan); ray.tau=malloc(sizeof(double) * img[im].nchan); - #pragma omp for + #pragma omp for schedule(dynamic) /* Main loop through pixel grid. */ for(px=0;px<(img[im].pxls*img[im].pxls);px++){ #pragma omp atomic From 0a0df711ad86e9b3796a77b69aa757c85c3a38b4 Mon Sep 17 00:00:00 2001 From: ims Date: Tue, 28 Jun 2016 14:40:02 +0200 Subject: [PATCH 35/43] Removed comment I had left a comment in the part of traceray() which dealt with polarized continuum, but this comment only indicated that I did not understand the code. I do now, so I have removed the comment. --- src/raytrace.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/raytrace.c b/src/raytrace.c index e0e0c4f..5e766d4 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -158,7 +158,7 @@ Note that the algorithm employed here is similar to that employed in the functio for(ichan=0;ichan Date: Wed, 29 Jun 2016 11:02:57 +0200 Subject: [PATCH 36/43] Fixed some errors There were some errors in raytrace() which I have here fixed. - Changed ran to randGen at line 213 - Progress is now measured relative to the total number of rays to be followed rather than (incorrectly) to the total number of image pixels. --- src/raytrace.c | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/raytrace.c b/src/raytrace.c index 99f5467..3ada806 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -189,12 +189,12 @@ 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){ - int *counta,*countb,nlinetot,aa,numRaysThisPixel; + int *counta,*countb,nlinetot,aa,numRaysThisPixel,numActiveRays; int ichan,gi,iline,tmptrans,i,threadI,nRaysDone; double size,minfreq,absDeltaFreq; double cutoff; unsigned int totalNumImagePixels,ppi; - double imgXiCentre,imgYiCentre,x[3],oneOnTotalNumPixelsMinus1; + double imgXiCentre,imgYiCentre,x[3],oneOnNumActiveRaysMinus1; int xi,yi,ri,j; const gsl_rng_type *ranNumGenType = gsl_rng_ranlxs2; @@ -210,12 +210,11 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ for (i=0;inThreads;i++){ threadRans[i] = gsl_rng_alloc(ranNumGenType); - gsl_rng_set(threadRans[i],(int)(gsl_rng_uniform(ran)*1e6)); + gsl_rng_set(threadRans[i],(int)(gsl_rng_uniform(randGen)*1e6)); } size=img[im].distance*img[im].imgres; totalNumImagePixels = img[im].pxls*img[im].pxls; - oneOnTotalNumPixelsMinus1 = 1.0/(double)totalNumImagePixels; imgXiCentre = img[im].pxls/2.0; imgYiCentre = img[im].pxls/2.0; @@ -288,6 +287,11 @@ In this algorithm we attempt to deal better with image pixels which cover areas exit(1); } + numActiveRays = 0; + for(ppi=0;ppiminScale*1.0e-7; nRaysDone=0; @@ -331,7 +335,7 @@ In this algorithm we attempt to deal better with image pixels which cover areas } } if (threadI == 0){ // i.e., is master thread - if(!silent) progressbar((double)(nRaysDone)*oneOnTotalNumPixelsMinus1, 13); + if(!silent) progressbar((double)(nRaysDone)*oneOnNumActiveRaysMinus1, 13); } } From 670ff7ce6983707f8276a5dbe66c77b2a18e99da Mon Sep 17 00:00:00 2001 From: ims Date: Wed, 29 Jun 2016 13:55:53 +0200 Subject: [PATCH 37/43] Added some comments and made some small mods. --- src/raytrace.c | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/raytrace.c b/src/raytrace.c index 188ae1c..34cbfbb 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -438,10 +438,9 @@ At the moment I will fix the number of segments, but it might possibly be faster /*....................................................................*/ void raytrace(int im, inputPars *par, struct grid *gp, molData *md, image *img){ - - /* -In this algorithm we attempt to deal better with image pixels which cover areas of the model where grid points are packed much more densely than the pixel spacing. In algorithm 0, such regions are sampled by only at par->antialias rays per pixel, which risks missing extreme values. The approach here is to generate approximately the same number of rays per pixel as the number of projected grid points in that pixel (with par->antialias providing a minimum floor). The set of values per pixel are then averaged to obtain the image value for that pixel. - */ + /* +This function constructs an image cube by following sets of rays (at least 1 per image pixel) through the model, solving the radiative transfer equations as appropriate for each ray. The ray locations within each pixel are chosen randomly within the pixel, but the number of rays per pixel is set equal to the number of projected model grid points falling within that pixel, down to a minimum equal to par->alias. + */ const gsl_rng_type *ranNumGenType=gsl_rng_ranlxs2; const double epsilon = 1.0e-6; /* Needs thinking about. Double precision is much smaller than this. */ @@ -558,6 +557,8 @@ In this algorithm we attempt to deal better with image pixels which cover areas img[im].pixel[ppi].numRays=0; } + /* Calculate the number of rays wanted per image pixel from the density of the projected model grid points. + */ for(gi=0;gipIntensity;gi++){ /* Apply the inverse (i.e. transpose) rotation matrix. (We use the inverse matrix here because here we want to rotate grid coordinates to the observer frame, whereas inside traceray() we rotate observer coordinates to the grid frame.) */ @@ -575,8 +576,13 @@ In this algorithm we attempt to deal better with image pixels which cover areas img[im].pixel[ppi].numRays++; } + /* Set a minimum number of rays per image pixel, and count the total number of rays. + */ numActiveRays = 0; for(ppi=0;ppiantialias) + img[im].pixel[ppi].numRays = par->antialias; + numActiveRays += img[im].pixel[ppi].numRays; oneOnNumActiveRaysMinus1 = 1.0/(double)(numActiveRays - 1); @@ -587,7 +593,7 @@ In this algorithm we attempt to deal better with image pixels which cover areas #pragma omp parallel private(ppi,molI,xi,yi) num_threads(par->nThreads) { /* Declaration of thread-private pointers. */ - int ai,ii,numRaysThisPixel,threadI = omp_get_thread_num(); + int ai,ii,threadI = omp_get_thread_num(); rayData ray; gridInterp gips[numInterpPoints]; double oneOnNRaysThisPixel; @@ -615,16 +621,12 @@ In this algorithm we attempt to deal better with image pixels which cover areas xi = (int)(ppi%(unsigned int)img[im].pxls); yi = floor(ppi/(double)img[im].pxls); - if(img[im].pixel[ppi].numRays>par->antialias) - numRaysThisPixel = img[im].pixel[ppi].numRays; - else - numRaysThisPixel = par->antialias; - oneOnNRaysThisPixel = 1.0/(double)numRaysThisPixel; + oneOnNRaysThisPixel = 1.0/(double)img[im].pixel[ppi].numRays; #pragma omp atomic ++nRaysDone; - for(ai=0;aitraceRayAlgorithm==1){ - for(ii=0;iinSpecies, gips[ii].mol); - } } free(ray.tau); free(ray.intensity); From 75c728219c53f300a3c6f73b26666637c3953286 Mon Sep 17 00:00:00 2001 From: ims Date: Wed, 29 Jun 2016 16:11:15 +0200 Subject: [PATCH 38/43] Removed some unused variables. --- src/photon.c | 4 +--- src/raythrucells.c | 9 +++------ src/raytrace.c | 3 +-- src/writefits.c | 2 +- 4 files changed, 6 insertions(+), 12 deletions(-) diff --git a/src/photon.c b/src/photon.c index 22cf38a..1c82cfb 100644 --- a/src/photon.c +++ b/src/photon.c @@ -173,7 +173,7 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar int iphot,iline,jline,here,there,firststep,dir,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; - double *tau,*expTau,x[3],inidir[3]; + double *tau,*expTau,inidir[3]; double remnantSnu; lineCount(par->nSpecies, m, &counta, &countb, &nlinetot); @@ -224,10 +224,8 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar else calcLineAmpLinear(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); mp[l].vfac[iphot]=vfac[0]; } - for(l=0;l<3;l++) x[l]=g[here].x[l]+(g[here].dir[dir].xn[l] * g[id].ds[dir]/2.); } else { ds=g[here].ds[dir]; - for(l=0;l<3;l++) x[l]=g[here].x[l]; } for(l=0;lnSpecies;l++){ diff --git a/src/raythrucells.c b/src/raythrucells.c index 2855e6f..4bd6724 100644 --- a/src/raythrucells.c +++ b/src/raythrucells.c @@ -117,7 +117,7 @@ At a successful termination, therefore, details of all the cells to the edge of const int bufferSize=1024; int numGoodExits, numMarginalExits, fi, goodExitFis[numFaces], marginalExitFis[numFaces], exitFi, i, status, newEntryFaceI; faceType face; - intersectType intcpt[numFaces], dummy_intcpt; + intersectType intcpt[numFaces]; followingSingleChain = 1; /* default */ do{ /* Follow the chain through 'good' cells, i.e. ones for which entry and exit are nicely distant from face edges. (Here we also follow marginal exits if there are no good ones, and only 1 marginal one.) */ @@ -297,8 +297,8 @@ This routine works best when the sides of the triangle are not too disparate in const int numDims=3; double norm[numDims], normDotDx, numerator, px2D[numDims-1], mat[numDims-1][numDims-1], vec[numDims-1], det; - int i, j, k, di; - double a, testSumForCW=0.0; + int j, k, di; + double testSumForCW=0.0; triangle2D face2D; intcpt->fi = -1; @@ -404,8 +404,6 @@ triangle2D calcTriangle2D(faceType face){ /**** all this could be precalculated (norm and mat too). */ triangle2D face2D; - const int numDims = 3; - int i; double lengthX, lengthY, dotResult; face2D.xAxis[0] = face.r[1][0]-face.r[0][0]; @@ -498,7 +496,6 @@ void doSegmentInterp(gridInterp gips[3], const int iA, molData *md\ const double fracA = (si + 0.5)*oneOnNumSegments, fracB = 1.0 - fracA; const int iB = 1 - iA; int di, molI, levelI, lineI; - double sum; gips[2].xCmpntRay = fracA*gips[iB].xCmpntRay + fracB*gips[iA].xCmpntRay; /* This does not seem to be used. */ diff --git a/src/raytrace.c b/src/raytrace.c index a94330d..d6d6b10 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -18,11 +18,10 @@ void calcLineAmpSample(const double x[3], const double dx[3], const double ds\ The bulk velocity of the model material can vary significantly with position, thus so can the value of the line-shape function at a given frequency and direction. The present function calculates 'vfac', an approximate average of the line-shape function along a path of length ds in the direction of the line of sight. */ int i; - double v,d,val,vel[3]; + double v,val; *vfac=0.; for(i=0;i Date: Wed, 29 Jun 2016 17:28:07 +0200 Subject: [PATCH 39/43] Added 1.4 raytrace algorithm The algorithm in 1.4 was much faster than the one in 1.5 because it solved the radiative transfer equation for a smaller number of rays, essentially corresponding to the locations of the grid points when projected onto the image plane. This has now been imported into the present merged raytracing code. I changed the 1.4 code in two respects: - The smoothing step is no longer done. - Instead of the somewhat arcane image-plane interpolation that is done between sparse rays in 1.4, a straightforward linear interpolation is performed. --- src/lime.h | 6 +- src/raytrace.c | 382 ++++++++++++++++++++++++++++++------------------- 2 files changed, 240 insertions(+), 148 deletions(-) diff --git a/src/lime.h b/src/lime.h index 6161330..469e49c 100644 --- a/src/lime.h +++ b/src/lime.h @@ -171,7 +171,10 @@ typedef struct { double deltav; } blend; -typedef struct {double x,y, *intensity, *tau;} rayData; +typedef struct { + double x,y, *intensity, *tau; + unsigned int ppi; +} rayData; /* NOTE that it is assumed that vertx[i] is opposite the face that abuts with neigh[i] for all i. */ @@ -241,6 +244,7 @@ void calcLineAmpSample(const double x[3], const double dx[3], const double, cons void calcLineAmpSpline(struct grid*, const int, const int, const double, const double, double*); void calcSourceFn(double, const inputPars*, double*, double*); void calcTableEntries(const int, const int); +void calcTriangleBaryCoords(double vertices[3][2], double, double, double barys[3]); triangle2D calcTriangle2D(faceType); void continuumSetup(int, image*, molData*, inputPars*, struct grid*); void delaunay(const int, struct grid*, const unsigned long, const _Bool, struct cell**, unsigned long*); diff --git a/src/raytrace.c b/src/raytrace.c index d6d6b10..c3dc715 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -188,6 +188,7 @@ Note that the algorithm employed here is similar to that employed in the functio for(iline=0;iline img[im].freq-img[im].bandwidth*0.5\ && md[molI].freq[lineI] < img[im].freq+img[im].bandwidth*0.5){ @@ -396,6 +397,7 @@ At the moment I will fix the number of segments, but it might possibly be faster for(iline=0;iline img[im].freq-img[im].bandwidth*0.5\ && md[molI].freq[lineI] < img[im].freq+img[im].bandwidth*0.5){ @@ -459,89 +461,39 @@ raytrace(int im, inputPars *par, struct grid *gp, molData *md, image *img){ This function constructs an image cube by following sets of rays (at least 1 per image pixel) through the model, solving the radiative transfer equations as appropriate for each ray. The ray locations within each pixel are chosen randomly within the pixel, but the number of rays per pixel is set equal to the number of projected model grid points falling within that pixel, down to a minimum equal to par->alias. */ - const gsl_rng_type *ranNumGenType=gsl_rng_ranlxs2; - const double epsilon = 1.0e-6; /* Needs thinking about. Double precision is much smaller than this. */ - const int numFaces=1+DIM, numInterpPoints=3, numSegments=5; - const double oneOnNFaces=1.0/(double)numFaces, oneOnNumSegments=1.0/(double)numSegments; + const int maxNumRaysPerPixel=20; /**** Arbitrary - could make this a global, or an argument. Set it to zero to indicate there is no maximum. */ + const double cutoff = par->minScale*1.0e-7; + const int numFaces=1+DIM, numInterpPoints=3, numSegments=5, minNumRaysForAverage=2; + const double oneOnNFaces=1.0/(double)numFaces, oneOnNumSegments = 1.0/(double)numSegments; + const double epsilon = 1.0e-6; // Needs thinking about. Double precision is much smaller than this. const int nStepsThruCell=10; const double oneOnNSteps=1.0/(double)nStepsThruCell; - int i,di,vi,gi,molI,ei,li,nlinetot,iline,tmptrans,ichan,xi,yi; - double size,imgCentreXPixels,imgCentreYPixels,sum,minfreq,absDeltaFreq,oneOnNumActiveRaysMinus1,cutoff,progress,x[2]; - unsigned int totalNumImagePixels,ppi,numActiveRays,nRaysDone; + double size,oneOnNumActiveRaysMinus1,imgCentreXPixels,imgCentreYPixels,minfreq,absDeltaFreq,x[2],sum,oneOnNumRays;//,oneOnTotalNumPixelsMinus1 + unsigned int totalNumImagePixels,ppi,numPixelsForInterp; + int gi,molI,ei,li,nlinetot,iline,tmptrans,ichan,numActiveRays,i,di,xi,yi,ri,c,id,ids[3],vi; + int *allLineMolIs,*allLineLineIs; + rayData *rays; struct cell *dc=NULL; - unsigned long numCells, dci; + unsigned long numCells,dci; struct gAuxType *gAux=NULL; /* This will hold some precalculated values for the grid points. */ - int *allLineMolIs,*allLineLineIs; - static double lastProgress = 0.0; - - gsl_rng *randGen = gsl_rng_alloc(ranNumGenType); /* Random number generator */ -#ifdef TEST - gsl_rng_set(randGen,178490); -#else - gsl_rng_set(randGen,time(0)); -#endif - - gsl_rng **threadRans; - threadRans = malloc(sizeof(gsl_rng *)*par->nThreads); - - for (i=0;inThreads;i++){ - threadRans[i] = gsl_rng_alloc(ranNumGenType); - gsl_rng_set(threadRans[i],(int)(gsl_rng_uniform(randGen)*1e6)); - } + coordT *pt_array, point[3]; + char flags[255]; + boolT ismalloc = False,isoutside; + realT bestdist; + facetT *facet;//, *neighbor, **neighborp; + vertexT *vertex,**vertexp; + int curlong, totlong; + double triangle[3][2],barys[3]; + _Bool isOutsideImage; size = img[im].distance*img[im].imgres; totalNumImagePixels = img[im].pxls*img[im].pxls; imgCentreXPixels = img[im].pxls/2.0; imgCentreYPixels = img[im].pxls/2.0; - if(par->traceRayAlgorithm==1){ - delaunay(DIM, gp, (unsigned long)par->ncell, 1, &dc, &numCells); - - /* We need to process the list of cells a bit further - calculate their centres, and reset the id values to be the same as the index of the cell in the list. (This last because we are going to construct other lists to indicate which cells have been visited etc.) - */ - for(dci=0;dcix[di]; - } - dc[dci].centre[di] = sum*oneOnNFaces; - } - - dc[dci].id = dci; - } - - }else if(par->traceRayAlgorithm!=0){ - if(!silent) bail_out("Unrecognized value of par.traceRayAlgorithm"); - exit(1); - } - - /* Precalculate binv*nmol*pops for all grid points. - */ - gAux = malloc(sizeof(*gAux)*par->ncell); - for(gi=0;gincell;gi++){ - gAux[gi].mol = malloc(sizeof(*(gAux[gi].mol))*par->nSpecies); - for(molI=0;molInSpecies;molI++){ - gAux[gi].mol[molI].specNumDens = malloc(sizeof(*(gAux[gi].mol[molI].specNumDens))*md[molI].nlev); - gAux[gi].mol[molI].dust = malloc(sizeof(*(gAux[gi].mol[molI].dust)) *md[molI].nline); - gAux[gi].mol[molI].knu = malloc(sizeof(*(gAux[gi].mol[molI].knu)) *md[molI].nline); - - for(ei=0;einSpecies, md, &allLineMolIs, &allLineLineIs, &nlinetot); + lineCount(par->nSpecies, md, &allLineMolIs, &allLineLineIs, &nlinetot); /* mallocs allLineMolIs, allLineLineIs */ if(img[im].doline==0) nlinetot=1; /* Fix the image parameters. */ @@ -567,111 +519,146 @@ This function constructs an image cube by following sets of rays (at least 1 per for(ppi=0;ppipIntensity); /* We may need to reallocate this later. */ + numActiveRays = 0; for(gi=0;gipIntensity;gi++){ /* Apply the inverse (i.e. transpose) rotation matrix. (We use the inverse matrix here because here we want to rotate grid coordinates to the observer frame, whereas inside traceray() we rotate observer coordinates to the grid frame.) */ for(i=0;i<2;i++){ x[i]=0.0; - for(di=0;di=0 && ppi=img[im].pxls || yi<0 || yi>=img[im].pxls){ + isOutsideImage = 1; + ppi = 0; /* Under these circumstances it ought never to be accessed, but it is not good to leave it without a value at all. */ + }else{ + isOutsideImage = 0; + ppi = (unsigned int)yi*(unsigned int)img[im].pxls + (unsigned int)xi; + } + + /* See if we want to keep the ray. For the time being we will include those outside the image bounds, but a cleverer algorithm would exclude some of them. Note that maxNumRaysPerPixel<1 is used to flag that there is no upper limit to the number of rays per pixel. + */ + if(isOutsideImage || maxNumRaysPerPixel<1 || img[im].pixel[ppi].numRayspIntensity) + rays = realloc(rays, sizeof(rayData)*numActiveRays); + + if(par->traceRayAlgorithm==1){ + delaunay(DIM, gp, (unsigned long)par->ncell, 1, &dc, &numCells); /* mallocs dc if getCells==T */ + + /* We need to process the list of cells a bit further - calculate their centres, and reset the id values to be the same as the index of the cell in the list. (This last because we are going to construct other lists to indicate which cells have been visited etc.) + */ + for(dci=0;dcix[di]; + } + dc[dci].centre[di] = sum*oneOnNFaces; + } + + dc[dci].id = dci; + } + + }else if(par->traceRayAlgorithm!=0){ + if(!silent) bail_out("Unrecognized value of par.traceRayAlgorithm"); + exit(1); } - /* Set a minimum number of rays per image pixel, and count the total number of rays. + /* Precalculate binv*nmol*pops for all grid points. */ - numActiveRays = 0; - for(ppi=0;ppiantialias) - img[im].pixel[ppi].numRays = par->antialias; + gAux = malloc(sizeof(*gAux)*par->ncell); + for(gi=0;gincell;gi++){ + gAux[gi].mol = malloc(sizeof(*(gAux[gi].mol))*par->nSpecies); + for(molI=0;molInSpecies;molI++){ + gAux[gi].mol[molI].specNumDens = malloc(sizeof(*(gAux[gi].mol[molI].specNumDens))*md[molI].nlev); + gAux[gi].mol[molI].dust = malloc(sizeof(*(gAux[gi].mol[molI].dust)) *md[molI].nline); + gAux[gi].mol[molI].knu = malloc(sizeof(*(gAux[gi].mol[molI].knu)) *md[molI].nline); - numActiveRays += img[im].pixel[ppi].numRays; - oneOnNumActiveRaysMinus1 = 1.0/(double)(numActiveRays - 1); + for(ei=0;eiminScale*1.0e-7; + /* This next is repetition. I do it in order to be able to use the same sourcefunc.c functions for the interpolated grid values as for the 'standard' algorithm. With a sensible arrangement of memory for the grid values, this would be unnecessary. + */ + for(li=0;linThreads) + #pragma omp parallel num_threads(par->nThreads) { - /* Declaration of thread-private pointers. */ - int ai,ii,threadI = omp_get_thread_num(); - rayData ray; + /* Declaration of thread-private pointers. + */ + int threadI = omp_get_thread_num(); + int ii, si, ri; gridInterp gips[numInterpPoints]; - double oneOnNRaysThisPixel; - - ray.intensity = malloc(sizeof(double)*img[im].nchan); - ray.tau = malloc(sizeof(double)*img[im].nchan); if(par->traceRayAlgorithm==1){ /* Allocate memory for the interpolation points: */ for(ii=0;iinSpecies); - for(molI=0;molInSpecies;molI++){ - gips[ii].mol[molI].specNumDens = malloc(sizeof(*(gips[ii].mol[molI].specNumDens))*md[molI].nlev); - gips[ii].mol[molI].dust = malloc(sizeof(*(gips[ii].mol[molI].dust)) *md[molI].nline); - gips[ii].mol[molI].knu = malloc(sizeof(*(gips[ii].mol[molI].knu)) *md[molI].nline); + for(si=0;sinSpecies;si++){ + gips[ii].mol[si].specNumDens = malloc(sizeof(*(gips[ii].mol[si].specNumDens))*md[si].nlev); + gips[ii].mol[si].dust = malloc(sizeof(*(gips[ii].mol[si].dust)) *md[si].nline); + gips[ii].mol[si].knu = malloc(sizeof(*(gips[ii].mol[si].knu)) *md[si].nline); } } } #pragma omp for schedule(dynamic) - /* Main loop through pixel grid. */ - for(ppi=0;ppitraceRayAlgorithm==0){ + traceray(rays[ri], par, tmptrans, img, im, gp, gAux, md, nlinetot\ + , allLineMolIs, allLineLineIs, cutoff, nStepsThruCell, oneOnNSteps); + }else if(par->traceRayAlgorithm==1) + traceray_smooth(rays[ri], par, tmptrans, img, im, gp, gAux, md, nlinetot\ + , allLineMolIs, allLineLineIs, dc, numCells, epsilon, gips, numSegments\ + , oneOnNumSegments, nStepsThruCell, oneOnNSteps); - for(ai=0;aitraceRayAlgorithm==0){ - traceray(ray, par, tmptrans, img, im, gp, gAux, md, nlinetot\ - , allLineMolIs, allLineLineIs, cutoff, nStepsThruCell, oneOnNSteps); - }else if(par->traceRayAlgorithm==1) - traceray_smooth(ray, par, tmptrans, img, im, gp, gAux, md, nlinetot\ - , allLineMolIs, allLineLineIs, dc, numCells, epsilon, gips, numSegments\ - , oneOnNumSegments, nStepsThruCell, oneOnNSteps); - - #pragma omp critical - { - for(ichan=0;ichan0.002){ - lastProgress = progress; - progressbar(progress, 13); - } - } + if(!silent) progressbar((double)(ri)*oneOnNumActiveRaysMinus1, 13); } } @@ -679,21 +666,122 @@ This function constructs an image cube by following sets of rays (at least 1 per for(ii=0;iinSpecies, gips[ii].mol); } - free(ray.tau); - free(ray.intensity); } /* End of parallel block. */ + /* For pixels with more than a cutoff number of rays, just average those rays into the pixel: + */ + for(ri=0;ri= minNumRaysForAverage){ + for(ichan=0;ichan= minNumRaysForAverage){ + oneOnNumRays = 1.0/(double)img[im].pixel[ppi].numRays; + for(ichan=0;ichan0){ + /* Now we enter main loop 3/3, in which we loop over image pixels, and for any we need to interpolate, we do so. But first we need to invoke qhull to get a Delaunay triangulation of the projected points. + */ + pt_array=malloc(sizeof(coordT)*2*numActiveRays); + + for(ri=0;rivertices ) { + id = qh_pointid(vertex->point); + triangle[c][0] = rays[id].x; + triangle[c][1] = rays[id].y; + ids[c]=id; + c++; + } + + calcTriangleBaryCoords(triangle, (double)point[0], (double)point[1], barys); + + /* Interpolate: */ + for(ichan=0;ichan0) */ + img[im].trans=tmptrans; freeGAux((unsigned long)par->ncell, par->nSpecies, gAux); - free(allLineMolIs); - free(allLineLineIs); if(par->traceRayAlgorithm==1) free(dc); - for (i=0;inThreads;i++){ - gsl_rng_free(threadRans[i]); + for(ri=0;ri Date: Thu, 30 Jun 2016 13:12:20 +0200 Subject: [PATCH 40/43] Updated documentation --- doc/usermanual.rst | 14 +++++++------- src/raytrace.c | 1 + 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/doc/usermanual.rst b/doc/usermanual.rst index d50dfbf..d0708d5 100644 --- a/doc/usermanual.rst +++ b/doc/usermanual.rst @@ -439,13 +439,7 @@ file. The default is blend=0 (no line blending). (integer) par->antialias (optional) -If set, LIME will anti-alias the output image. anti-alias can take the -value of any positive integer, with the value 1 (default) being no -anti-aliasing. Greater values correspond to stronger anti-aliasing. LIME -uses stochastic super-sampling anti- aliasing. This is very effective in -minimizing artifacts in the image, but it also slows down the ray-tracer -by a factor equal to the value of antialias. This parameter is the only -one that will not be ignored in case par->restart is set. +This parameter is no longer used, although it is retained for the present for purposes of backward compatibility. .. code:: c @@ -459,6 +453,12 @@ work, a magnetic field needs to be defined (see below). When polarization is switched on, LIME is identical to the DustPol code (Padovani et al., 2012). + (integer) par->traceRayAlgorithm (optional) + +This parameter specifies the algorithm used by LIME to solve the radiative-transfer equations during ray-tracing. The default value of zero invokes the algorithm used in LIME-1.5 and previous; a value of 1 invokes a new algorithm which is much more time-consuming but which produces much smoother images, free from step-artifacts. + +Note that there have been additional modifications to the raytracing algorithm which have significant effects on the output images since LIME-1.5. Image-plane interpolation is now employed in areas of the image where the grid point spacing is larger than the image pixel spacing. This leads both to a smoother image and a shorter processing time. + .. _par-nthreads: .. code:: c diff --git a/src/raytrace.c b/src/raytrace.c index c3dc715..efcf36d 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -5,6 +5,7 @@ * Copyright (C) 2006-2014 Christian Brinch * Copyright (C) 2015 The LIME development team * +TODO: - In raytrace(), look at rearranging the code to do the qhull step before choosing the rays. This would allow cells with all vertices outside the image boundaries to be excluded. If the image is much smaller than the model, this could lead to significant savings in time. The only downside might be memory useage... */ #include "lime.h" From a7c30b8bc1288ad32e9cc7b5eb3beb89dcdf9260 Mon Sep 17 00:00:00 2001 From: ims Date: Thu, 21 Jul 2016 16:23:47 +0200 Subject: [PATCH 41/43] Moved call to gasIIdust() outside of line loop. --- src/grid.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid.c b/src/grid.c index 72d0133..b5ef829 100644 --- a/src/grid.c +++ b/src/grid.c @@ -124,8 +124,8 @@ void calcGridDustOpacity(configInfo *par, molData *md, struct grid *gp){ for(di=0;dinumDensities;di++) densityForDust += gp[id].dens[di]*par->dustWeights[di]; + gasIIdust(gp[id].x[0],gp[id].x[1],gp[id].x[2],>d); for(iline=0;iline Date: Thu, 21 Jul 2016 16:29:27 +0200 Subject: [PATCH 42/43] Removed unnecessary NULL tests in frees.c --- src/frees.c | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/src/frees.c b/src/frees.c index 9716ade..426b519 100644 --- a/src/frees.c +++ b/src/frees.c @@ -14,9 +14,8 @@ freeGAux(const unsigned long numPoints, const int numSpecies, struct gAuxType *g unsigned long ppi; if(gAux !=NULL){ - for(ppi=0;ppi Date: Fri, 22 Jul 2016 13:17:23 +0200 Subject: [PATCH 43/43] Fixed iatrogenic bug in traceray(). --- src/frees.c | 1 - src/raytrace.c | 6 ++++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/frees.c b/src/frees.c index 426b519..84c88e0 100644 --- a/src/frees.c +++ b/src/frees.c @@ -92,7 +92,6 @@ void freeMolsWithBlends(struct molWithBlends *mols, const int numMolsWithBlends) if(mols[mi].lines != NULL){ for(li=0;li