diff --git a/.gitignore b/.gitignore index a2d083a..1792bba 100644 --- a/.gitignore +++ b/.gitignore @@ -8,9 +8,9 @@ example/*.vtk example/*.pop example/*.fits* +example/lime_*.x example/*.dat example/*_archive -example/lime_*.x example/vlog.txt # Ignore html doc files diff --git a/Makefile b/Makefile index d61daf0..e9b1170 100644 --- a/Makefile +++ b/Makefile @@ -49,26 +49,26 @@ endif ## TARGET = lime.x -#CC = gcc -fopenmp -g +#CC = gcc -fopenmp -g -Wunused -Wno-unused-result CC = gcc -fopenmp -SRCS = src/aux.c src/messages.c src/grid.c src/LTEsolution.c \ - src/main.c src/molinit.c src/photon.c src/popsin.c \ - src/popsout.c src/predefgrid.c src/ratranInput.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/tree_random.c \ - src/velospline.c src/getclosest.c src/frees.c \ - src/tcpsocket.c src/defaults.c src/fastexp.c \ +SRCS = src/aux.c src/messages.c src/grid.c src/LTEsolution.c \ + src/main.c src/molinit.c src/photon.c src/popsin.c \ + src/popsout.c src/predefgrid.c src/ratranInput.c \ + src/raytrace.c src/smooth.c src/sourcefunc.c src/frees.c \ + src/stateq.c src/statistics.c src/magfieldfit.c \ + src/stokesangles.c src/writefits.c src/tree_random.c \ + src/velospline.c src/getclosest.c src/grid2fits.c \ + src/tcpsocket.c src/defaults.c src/fastexp.c src/gridio.c \ src/raythrucells.c MODELS = model.c -OBJS = src/aux.o src/messages.o src/grid.o src/LTEsolution.o \ - src/main.o src/molinit.o src/photon.o src/popsin.o \ - src/popsout.o src/predefgrid.o src/raytrace.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/tree_random.o\ - src/velospline.o src/getclosest.o src/frees.o \ - src/tcpsocket.o src/defaults.o src/fastexp.o \ +OBJS = src/aux.o src/messages.o src/grid.o src/LTEsolution.o \ + src/main.o src/molinit.o src/photon.o src/popsin.o \ + src/popsout.o src/predefgrid.o src/raytrace.o \ + src/ratranInput.o src/smooth.o src/sourcefunc.o src/frees.o \ + src/stateq.o src/statistics.o src/magfieldfit.o \ + src/stokesangles.o src/writefits.o src/tree_random.o \ + src/velospline.o src/getclosest.o src/grid2fits.o \ + src/tcpsocket.o src/defaults.o src/fastexp.o src/gridio.o \ src/raythrucells.o MODELO = src/model.o diff --git a/example/model.c b/example/model.c index c0a359c..d9502ad 100644 --- a/example/model.c +++ b/example/model.c @@ -26,6 +26,7 @@ input(inputPars *par, image *img){ par->moldatfile[0] = "hco+@xpol.dat"; par->antialias = 4; par->sampling = 2; // log distr. for radius, directions distr. uniformly on a sphere. + par->nSolveIters = 14; par->outputfile = "populations.pop"; par->binoutputfile = "restart.pop"; @@ -81,6 +82,17 @@ input(inputPars *par, image *img){ par->nMolWeights[0] = 1.0; par->dustWeights[0] = 1.0; + /* Set one or more of the following parameters for full output of the grid-specific data at any of 4 stages during the processing. (See the header of gridio.c for information about the stages.) + par->gridOutFiles[0] = "grid_1.ds"; + par->gridOutFiles[1] = "grid_2.ds"; + par->gridOutFiles[2] = "grid_3.ds"; + par->gridOutFiles[3] = "grid_4.ds"; + */ + + /* You can also optionally read in a FITS file stored via the previous parameters, or prepared externally. See the header of grid2fits.c for information about the correct file format. LIME can cope with almost any sensible subset of the recognized columns; it will use the file values if they are present, then calculate the missing ones. + par->gridInFile = "grid_4.ds"; + */ + /* * Definitions for image #0. Add blocks with successive values of i for additional images. */ diff --git a/src/aux.c b/src/aux.c index 0fc4fe5..eff8e9c 100644 --- a/src/aux.c +++ b/src/aux.c @@ -47,14 +47,24 @@ parseInput(inputPars inpar, configInfo *par, image **img, molData **m){ par->antialias = inpar.antialias; par->polarization = inpar.polarization; par->nThreads = inpar.nThreads; + par->gridInFile = inpar.gridInFile; + par->nSolveIters = inpar.nSolveIters; par->traceRayAlgorithm = inpar.traceRayAlgorithm; + par->gridOutFiles = malloc(sizeof(char *)*NUM_GRID_STAGES); + for(i=0;igridOutFiles[i] = inpar.gridOutFiles[i]; + /* Now set the additional values in par. */ par->ncell = inpar.pIntensity + inpar.sinkPoints; par->radiusSqu = inpar.radius*inpar.radius; par->minScaleSqu=inpar.minScale*inpar.minScale; par->doPregrid = (inpar.pregrid==NULL)?0:1; + for(i=0;iwriteGridAtStage[i] = 0; + par->dataFlags = 0; + /* If the user has provided a list of moldatfile names, the corresponding elements of par->moldatfile will be non-NULL. Thus we can deduce the number of files (species) from the number of non-NULL elements. */ par->nSpecies=0; @@ -104,6 +114,11 @@ parseInput(inputPars inpar, configInfo *par, image **img, molData **m){ while(igridDensMaxValues[i]>=0) i++; par->numGridDensMaxima = i; + /* Check that the user has supplied the velocity function (needed in raytracing unless par->doPregrid). Note that the other previously mandatory functions (density, abundance, doppler and temperature) may not be necessary if the user reads in the appropriate values from a file. This is tested at the appropriate place in readOrBuildGrid(). + */ + if(!par->doPregrid || par->traceRayAlgorithm==1) + velocity(0.0,0.0,0.0, dummyVel); + /* Calculate par->numDensities. */ if(!(par->doPregrid || par->restart)){ /* These switches cause par->numDensities to be set in routines they call. */ @@ -143,10 +158,10 @@ parseInput(inputPars inpar, configInfo *par, image **img, molData **m){ while((*img)[par->nImages].filename!=NULL && par->nImagesnImages++; - /* Check that the user has supplied this function (needed unless par->doPregrid): - */ - if(!par->doPregrid || par->traceRayAlgorithm==1) - velocity(0.0,0.0,0.0, dummyVel); + for(i=0;igridOutFiles[i] != NULL) + par->writeGridAtStage[i] = 1; + }; /* Now we need to calculate the cutoff value used in calcSourceFn(). The issue is to decide between @@ -680,7 +695,7 @@ Pointers are indicated by a * before the attribute name and an arrow to the memo void levelPops(molData *md, configInfo *par, struct grid *gp, int *popsdone, double *lamtab, double *kaptab, const int nEntries){ - int id,conv=0,iter,ilev,prog=0,ispec,c=0,n,i,threadI,nVerticesDone,nlinetot; + int id,iter,ilev,ispec,c=0,n,i,threadI,nVerticesDone,nItersDone,nlinetot;//,conv=0 double percent=0.,*median,result1=0,result2=0,snr,delta_pop; int nextMolWithBlend; struct statistics { double *pop, *ave, *sigma; } *stat; @@ -749,8 +764,9 @@ This is done to allow proper handling of errors which may arise in the LU solver While this is off however, other gsl_* etc calls will not exit if they encounter a problem. We may need to pay some attention to trapping their errors. */ - do{ - if(!silent) progressbar2(0, prog++, 0, result1, result2); + nItersDone=0; + while(nItersDone < par->nSolveIters){ /* Not a 'for' loop because we will probably later want to add a convergence criterion. */ + if(!silent) progressbar2(par, 0, nItersDone, 0, result1, result2); for(id=0;idpIntensity;id++){ for(ilev=0;ilev1){ + if(nItersDone>1){ result1=median[0]; result2 =gsl_stats_median_from_sorted_data(median, 1, c); } free(median); - if(!silent) progressbar2(1, prog, percent, result1, result2); - if(par->outputfile) popsout(par,gp,md); - } while(conv++outputfile != NULL) popsout(par,gp,md); + nItersDone++; + } gsl_set_error_handler(defaultErrorHandler); freeMolsWithBlends(blends.mols, blends.numMolsWithBlends); @@ -867,8 +884,46 @@ While this is off however, other gsl_* etc calls will not exit if they encounter free(stat); } - if(par->binoutputfile) binpopsout(par,gp,md); + par->dataFlags |= (1 << DS_bit_populations); + + if(par->binoutputfile != NULL) binpopsout(par,gp,md); *popsdone=1; } +_Bool allBitsSet(const int flags, const int mask){ + /* Returns true only if all the masked bits of flags are set. */ + + if(~flags & mask) + return 0; + else + return 1; +} + +_Bool anyBitSet(const int flags, const int mask){ + /* Returns true if any of the masked bits of flags are set. */ + + if(flags & mask) + return 1; + else + return 0; +} + +_Bool bitIsSet(const int flags, const int bitI){ + /* Returns true if the designated bit of flags is set. */ + + if(flags & (1 << bitI)) + return 1; + else + return 0; +} + +_Bool onlyBitsSet(const int flags, const int mask){ + /* Returns true if flags has no bits set apart from those which are true in mask. */ + + if(flags & ~mask) + return 0; + else + return 1; +} + diff --git a/src/defaults.c b/src/defaults.c index 3b636c4..dca3be5 100644 --- a/src/defaults.c +++ b/src/defaults.c @@ -72,7 +72,6 @@ The grid points within the model are chosen randomly via the rejection method wi density(r[0],r[1],r[2],val); for (i=0;inumDensities;i++) totalDensity += val[i]; -// fracDensity = pow(totalDensity/par->gridDensGlobalMax,DENSITY_POWER); fracDensity = pow(totalDensity,DENSITY_POWER)/par->gridDensGlobalMax; return fracDensity; diff --git a/src/frees.c b/src/frees.c index f2eb193..91884cb 100644 --- a/src/frees.c +++ b/src/frees.c @@ -25,6 +25,9 @@ void freeGrid(const unsigned int numPoints, const unsigned short numSpecies\ if(gp != NULL){ for(i_u=0;i_ucollPartIds); free(par->nMolWeights); free(par->dustWeights); + free(par->gridOutFiles); } void freePopulation(const unsigned short numSpecies, struct populations *pop){ diff --git a/src/grid.c b/src/grid.c index f41f04b..a3d93a3 100644 --- a/src/grid.c +++ b/src/grid.c @@ -9,6 +9,7 @@ #include "lime.h" #include "tree_random.h" +#include "gridio.h" /*....................................................................*/ void mallocAndSetDefaultGrid(struct grid **gp, const unsigned int numPoints){ @@ -83,9 +84,10 @@ void calcGridMolSpecNumDens(configInfo *par, molData *md, struct grid *gp){ for(gi=0;gincell;gi++){ for(ispec=0;ispecnSpecies;ispec++){ - for(ei=0;eiminScale. */ double r,theta,phi,sinPhi,z,semiradius; double uniformRandom; - int i,j,di; + int j,di; unsigned int i_u; // _Bool pointIsAccepted; int pointIsAccepted; @@ -826,180 +830,394 @@ void randomsViaRejection(configInfo *par, const unsigned int desiredNumPoints, g } } +/*....................................................................*/ +void writeGridIfRequired(configInfo *par, struct grid *gp, molData *md, const int fileFormatI){ + int status = 0; + char **collPartNames=NULL; /*** this is a placeholder until we start reading these. */ + char message[80]; + int dataStageI=0; + + /* Work out the data stage: + */ + if(!allBitsSet(par->dataFlags, DS_mask_1)){ + if(!silent) warning("Trying to write at data stage 0."); + return; + } + + if( allBitsSet(par->dataFlags, DS_mask_4)){ + dataStageI = 4; + }else if(allBitsSet(par->dataFlags, DS_mask_3)){ + dataStageI = 3; + }else if(allBitsSet(par->dataFlags, DS_mask_2)){ + dataStageI = 2; + }else{ + dataStageI = 1; + } + + if(par->writeGridAtStage[dataStageI-1]){ + struct gridInfoType gridInfo; + unsigned short i_us; + struct keywordType *primaryKwds=malloc(sizeof(struct keywordType)*1); + + gridInfo.nInternalPoints = par->pIntensity; + gridInfo.nSinkPoints = par->sinkPoints; + gridInfo.nLinks = 0; /* This quantity is calculated when writing to file. */ + gridInfo.nNNIndices = 0; /* This quantity is calculated when writing to file. */ + gridInfo.nDims = DIM; + gridInfo.nSpecies = par->nSpecies; + gridInfo.nDensities = par->numDensities; + gridInfo.nLinkVels = NUM_VEL_COEFFS; + if(md==NULL) + gridInfo.mols = NULL; + else{ + gridInfo.mols = malloc(sizeof(*(gridInfo.mols))*gridInfo.nSpecies); + for(i_us=0;i_usradius; + primaryKwds[0].comment = "[m] Model radius."; + + status = writeGrid(par->gridOutFiles[dataStageI-1], fileFormatI\ + , gridInfo, primaryKwds, 1, gp, collPartNames, par->dataFlags); + + free(primaryKwds); + free(gridInfo.mols); + + if(status){ + sprintf(message, "writeGrid at data stage %d returned with status %d", dataStageI, status); + if(!silent) bail_out(message); + exit(1); + } + } + + free(collPartNames); +} + +/*....................................................................*/ void -buildGrid(configInfo *par, struct grid *g){ +readOrBuildGrid(configInfo *par, struct grid **gp){ const gsl_rng_type *ranNumGenType = gsl_rng_ranlxs2; - int i, k, di, numSubFields, levelI=0; - double theta,semiradius,z,fieldVolume,sumDensity,maxDensity,minNumDensity,nToP,sphereVolume; - double vals[MAX_N_COLL_PART]; + int i,j,k,di,si,levelI=0,status=0,numCollPartRead; + double theta,semiradius,z; double *outRandDensities=NULL; double (*outRandLocations)[DIM]=NULL; - extern double densityNormalizer, minDensity; - extern int numCollisionPartners; treeRandConstType rinc; treeRandVarType rinv; struct cell *dc=NULL; /* Not used at present. */ unsigned long numCells; + struct gridInfoType gridInfoRead; + char **collPartNames; + char message[80]; double x[DIM]; - unsigned int startI=0; treeType tree; - rinc.randGen = gsl_rng_alloc(ranNumGenType); /* Random number generator */ + par->dataFlags = 0; + if(par->gridInFile!=NULL){ + const int numDesiredKwds=1; + struct keywordType *desiredKwds=malloc(sizeof(struct keywordType)*numDesiredKwds); + + initializeKeyword(&desiredKwds[0]); + desiredKwds[0].datatype = TDOUBLE; + desiredKwds[0].keyname = "RADIUS "; + /* Currently not doing anything with the read keyword. */ + + status = readGrid(par->gridInFile, lime_FITS, &gridInfoRead, desiredKwds\ + , numDesiredKwds, gp, &collPartNames, &numCollPartRead, &(par->dataFlags)); + + if(status){ + if(!silent){ + sprintf(message, "Read of grid file failed with status return %d", status); + bail_out(message); + } + exit(1); + } + + /* Test that dataFlags obeys the rules. */ + /* No other bit may be set if DS_bit_x is not: */ + if(anyBitSet(par->dataFlags, (DS_mask_all & ~(1 << DS_bit_x))) && !bitIsSet(par->dataFlags, DS_bit_x)){ + if(!silent) bail_out("You may not read a grid file without X, ID or IS_SINK data."); + exit(1); + } + + /* DS_bit_ACOEFF may not be set if either DS_bit_neighbours or DS_bit_velocity is not: */ + if(bitIsSet(par->dataFlags, DS_bit_ACOEFF)\ + && !(bitIsSet(par->dataFlags, DS_bit_neighbours) && bitIsSet(par->dataFlags, DS_bit_velocity))){ + if(!silent) bail_out("You may not read a grid file with ACOEFF but no VEL or neighbour data."); + exit(1); + } + + /* DS_bit_populations may not be set unless all the others (except DS_bit_magfield) are set as well: */ + if(bitIsSet(par->dataFlags, DS_bit_populations)\ + && !allBitsSet(par->dataFlags & DS_mask_all_but_mag, DS_mask_populations)){ + if(!silent) bail_out("You may not read a grid file with pop data unless all other data is present."); + exit(1); + } + + /* Test gridInfoRead values against par values and overwrite the latter, with a warning, if necessary. + */ + if(gridInfoRead.nSinkPoints>0 && par->sinkPoints>0){ + if((int)gridInfoRead.nSinkPoints!=par->sinkPoints){ + if(!silent) warning("par->sinkPoints will be overwritten"); + par->sinkPoints = (int)gridInfoRead.nSinkPoints; + } + if((int)gridInfoRead.nInternalPoints!=par->pIntensity){ + if(!silent) warning("par->pIntensity will be overwritten"); + par->pIntensity = (int)gridInfoRead.nInternalPoints; + } + par->ncell = par->sinkPoints + par->pIntensity; + } + if(gridInfoRead.nDims!=DIM){ /* At present this situation is already detected and handled inside readGridExtFromFits(), but it may not be in future. The test here has no present functionality but saves trouble later if we change grid.x from an array to a pointer. */ + if(!silent){ + sprintf(message, "Grid file had %d dimensions but there should be %d.", (int)gridInfoRead.nDims, DIM); + bail_out(message); + } + exit(1); + } + if(gridInfoRead.nSpecies>0 && par->nSpecies>0 && (int)gridInfoRead.nSpecies!=par->nSpecies){ + if(!silent){ + sprintf(message, "Grid file had %d species but you have provided moldata files for %d."\ + , (int)gridInfoRead.nSpecies, par->nSpecies); + bail_out(message); + } + exit(1); +/**** should compare name to name - at some later time after we have read these from the moldata files? */ + } + if(gridInfoRead.nDensities>0 && par->numDensities>0 && (int)gridInfoRead.nDensities!=par->numDensities){ + if(!silent){ + sprintf(message, "Grid file had %d densities but you have provided %d."\ + , (int)gridInfoRead.nDensities, par->numDensities); + bail_out(message); + } + exit(1); + } + +/* +**** Ideally we should also have a test on nACoeffs. + +**** Ideally we should also have a test on the mols entries - at some later time after we have read the corresponding values from the moldata files? +*/ + } /* End of read grid file. Whether and what we subsequently calculate will depend on the value of par->dataStageI returned. */ + + if(!anyBitSet(par->dataFlags, DS_mask_x)){ /* This should only happen if we did not read a file. Generate the grid point locations. */ + mallocAndSetDefaultGrid(gp, (unsigned int)par->ncell); + + rinc.randGen = gsl_rng_alloc(ranNumGenType); /* Random number generator */ #ifdef TEST - gsl_rng_set(rinc.randGen,342971); + gsl_rng_set(rinc.randGen,342971); #else - gsl_rng_set(rinc.randGen,time(0)); + gsl_rng_set(rinc.randGen,time(0)); #endif - outRandDensities = malloc(sizeof(double )*par->pIntensity); /* Not used at present; and in fact they are not useful outside this routine, because they are not the values of the physical density at that point, just what densityFunc3D() returns, which is not necessarily the same thing. */ - outRandLocations = malloc(sizeof(*outRandLocations)*par->pIntensity); - - if(par->samplingAlgorithm==0){ - randomsViaRejection(par, (unsigned int)par->pIntensity, rinc.randGen, outRandLocations); - - } else if(par->samplingAlgorithm==1){ - rinc.par = *par; - rinc.verbosity = 0; - rinc.numInRandoms = TREE_N_RANDOMS; - rinc.maxRecursion = TREE_MAX_RECURSION; - rinc.maxNumTrials = TREE_MAX_N_TRIALS; - rinc.dither = TREE_DITHER; - rinc.maxNumTrialsDbl = (double)rinc.maxNumTrials; - rinc.doShuffle = 1; - rinc.doQuasiRandom = 1; - - for(di=0;diradius; - rinv.fieldWidth[di] = 2.0*par->radius; - } - rinv.expectedDesNumPoints = (double)par->pIntensity; - rinc.desiredNumPoints = (unsigned int)par->pIntensity; - - rinv.numHighPoints = par->numGridDensMaxima; - if(par->numGridDensMaxima>0){ - rinv.highPointLocations = malloc(sizeof(*(rinv.highPointLocations))*par->numGridDensMaxima); - rinv.highPointDensities = malloc(sizeof(double )*par->numGridDensMaxima); - for(i=0;inumGridDensMaxima;i++){ - for(di=0;digridDensMaxLoc[i][di]; + outRandDensities = malloc(sizeof(double )*par->pIntensity); /* Not used at present; and in fact they are not useful outside this routine, because they are not the values of the physical density at that point, just what densityFunc3D() returns, which is not necessarily the same thing. */ + outRandLocations = malloc(sizeof(*outRandLocations)*par->pIntensity); + + if(par->samplingAlgorithm==0){ + randomsViaRejection(par, (unsigned int)par->pIntensity, rinc.randGen, outRandLocations); + + } else if(par->samplingAlgorithm==1){ + rinc.par = *par; + rinc.verbosity = 0; + rinc.numInRandoms = TREE_N_RANDOMS; + rinc.maxRecursion = TREE_MAX_RECURSION; + rinc.maxNumTrials = TREE_MAX_N_TRIALS; + rinc.dither = TREE_DITHER; + rinc.maxNumTrialsDbl = (double)rinc.maxNumTrials; + rinc.doShuffle = 1; + rinc.doQuasiRandom = 1; + + for(di=0;diradius; + rinv.fieldWidth[di] = 2.0*par->radius; + } + rinv.expectedDesNumPoints = (double)par->pIntensity; + rinc.desiredNumPoints = (unsigned int)par->pIntensity; + + rinv.numHighPoints = par->numGridDensMaxima; + if(par->numGridDensMaxima>0){ + rinv.highPointLocations = malloc(sizeof(*(rinv.highPointLocations))*par->numGridDensMaxima); + rinv.highPointDensities = malloc(sizeof(double )*par->numGridDensMaxima); + for(i=0;inumGridDensMaxima;i++){ + for(di=0;digridDensMaxLoc[i][di]; + } + rinv.highPointDensities[i] = par->gridDensMaxValues[i]; } - rinv.highPointDensities[i] = par->gridDensMaxValues[i]; + }else{ + rinv.highPointLocations = NULL; + rinv.highPointDensities = NULL; } - }else{ - rinv.highPointLocations = NULL; - rinv.highPointDensities = NULL; + + initializeTree(&rinc, &rinv, gridDensity, &tree); + constructTheTree(&rinc, &rinv, levelI, gridDensity, &tree); + fillTheTree(&rinc, &tree, gridDensity, outRandLocations, outRandDensities); + + free(tree.leaves); + freeRinv(rinv); + free(rinc.inRandLocations); + + } else { + if(!silent) bail_out("Unrecognized sampling algorithm."); + exit(1); } - initializeTree(&rinc, &rinv, gridDensity, &tree); - constructTheTree(&rinc, &rinv, levelI, gridDensity, &tree); - fillTheTree(&rinc, &tree, gridDensity, outRandLocations, outRandDensities); + for(k=0;kpIntensity;k++){ + /* Assign values to the k'th grid point */ + (*gp)[k].id=k; + (*gp)[k].x[0]=outRandLocations[k][0]; + (*gp)[k].x[1]=outRandLocations[k][1]; + if(DIM==3) (*gp)[k].x[2]=outRandLocations[k][2]; + (*gp)[k].sink=0; + } - free(tree.leaves); - freeRinv(rinv); - free(rinc.inRandLocations); + /* end model grid point assignment */ + if(!silent) printDone(4); - } else { - if(!silent) bail_out("Unrecognized sampling algorithm."); - exit(1); - } + /* Add surface sink particles */ + for(k=par->pIntensity;kncell;k++){ + theta=gsl_rng_uniform(rinc.randGen)*2*PI; - for(k=0;kpIntensity;k++){ - /* Assign values to the k'th grid point */ - g[k].id=k; - g[k].x[0]=outRandLocations[k][0]; - g[k].x[1]=outRandLocations[k][1]; - if(DIM==3) g[k].x[2]=outRandLocations[k][2]; - g[k].sink=0; + if(DIM==3) { + z=2*gsl_rng_uniform(rinc.randGen)-1.; + semiradius=sqrt(1.-z*z); + x[2]=z; + } else { + semiradius=1.0; + } - /* This next step needs to be done, even though it looks stupid */ - g[k].dir=malloc(sizeof(point)*1); - g[k].ds =malloc(sizeof(double)*1); - g[k].neigh =malloc(sizeof(struct grid *)*1); + x[0]=semiradius*cos(theta); + x[1]=semiradius*sin(theta);; + (*gp)[k].id=k; + (*gp)[k].x[0]=par->radius*x[0]; + (*gp)[k].x[1]=par->radius*x[1]; + if(DIM==3) (*gp)[k].x[2]=par->radius*x[2]; + (*gp)[k].sink=1; + } + /* end grid allocation */ + + free(outRandLocations); + free(outRandDensities); + gsl_rng_free(rinc.randGen); + + if(par->samplingAlgorithm==0){ + smooth(par,*gp); + if(!silent) printDone(5); + } + + par->dataFlags |= DS_mask_1; } - /* end model grid point assignment */ - if(!silent) printDone(4); + if(onlyBitsSet(par->dataFlags, DS_mask_1)) /* Only happens if (i) we read no file and have constructed this data within LIME, or (ii) we read a file at dataStageI==1. */ + writeGridIfRequired(par, *gp, NULL, lime_FITS); - for(i=0;incell;i++){ - g[i].dens = malloc(sizeof(double)*par->numDensities); - g[i].abun = malloc(sizeof(double)*par->nSpecies); + if(!allBitsSet(par->dataFlags, DS_mask_neighbours)){ + delaunay(DIM, *gp, (unsigned long)par->ncell, 0, &dc, &numCells); + + par->dataFlags |= DS_mask_neighbours; } + distCalc(par, *gp); /* Mallocs and sets .dir & .ds, sets .nphot. We don't store these values so we have to calculate them whether we read a file or not. */ - /* Add surface sink particles */ - for(i=0;isinkPoints;i++){ - theta=gsl_rng_uniform(rinc.randGen)*2*PI; + if(onlyBitsSet(par->dataFlags, DS_mask_2)) /* Only happens if (i) we read no file and have constructed this data within LIME, or (ii) we read a file at dataStageI==2. */ + writeGridIfRequired(par, *gp, NULL, lime_FITS); - if(DIM==3) { - z=2*gsl_rng_uniform(rinc.randGen)-1.; - semiradius=sqrt(1.-z*z); - x[2]=z; - } else { - semiradius=1.0; + if(!allBitsSet(par->dataFlags, DS_mask_velocity)){ + for(i=0;ipIntensity;i++) + velocity((*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2],(*gp)[i].vel); + + /* Set velocity values also for sink points (otherwise Delaunay ray-tracing has problems) */ + for(i=par->pIntensity;incell;i++) + velocity((*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2],(*gp)[i].vel); + + par->dataFlags |= DS_mask_velocity; + } + + if(!allBitsSet(par->dataFlags, DS_mask_density)){ + for(i=0;incell; i++) + (*gp)[i].dens = malloc(sizeof(double)*par->numDensities); + for(i=0;ipIntensity;i++) + density((*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2],(*gp)[i].dens); + for(i=par->pIntensity;incell;i++){ + for(j=0;jnumDensities;j++) + (*gp)[i].dens[j]=1e-30;//************** what is the low but non zero value for? } - x[0]=semiradius*cos(theta); - x[1]=semiradius*sin(theta);; - g[k].id=k; - g[k].x[0]=par->radius*x[0]; - g[k].x[1]=par->radius*x[1]; - if(DIM==3) g[k].x[2]=par->radius*x[2]; - g[k].sink=1; - g[k].abun[0]=0; - g[k].dens[0]=1e-30;//************** what is the low but non zero value for? - g[k].t[0]=par->tcmb; - g[k].t[1]=par->tcmb; - g[k++].dopb_turb=0.; + par->dataFlags |= DS_mask_density; } - /* end grid allocation */ - /* Check that the user has supplied all necessary functions: - */ - density( 0.0,0.0,0.0, g[0].dens); - temperature(0.0,0.0,0.0, g[0].t); - doppler( 0.0,0.0,0.0,&g[0].dopb_turb); - abundance( 0.0,0.0,0.0, g[0].abun); - /* Note that velocity() is the only one of the 5 mandatory functions which is still needed (in raytrace) unless par->doPregrid. Therefore we test it already in parseInput(). */ + checkGridDensities(par, *gp); - delaunay(DIM, g, (unsigned long)par->ncell, 0, &dc, &numCells); - distCalc(par, g); + if(!allBitsSet(par->dataFlags, DS_mask_abundance)){ + for(i=0;incell; i++) + (*gp)[i].abun = malloc(sizeof(double)*par->nSpecies); + for(i=0;ipIntensity;i++) + abundance((*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2],(*gp)[i].abun); + for(i=par->pIntensity;incell;i++){ + for(si=0;sinSpecies;si++) + (*gp)[i].abun[si]=0; + } - if(par->samplingAlgorithm==0){ - smooth(par,g); - if(!silent) printDone(5); + par->dataFlags |= DS_mask_abundance; } - for(i=0;ipIntensity;i++){ - density( g[i].x[0],g[i].x[1],g[i].x[2], g[i].dens); - temperature(g[i].x[0],g[i].x[1],g[i].x[2], g[i].t); - doppler( g[i].x[0],g[i].x[1],g[i].x[2],&g[i].dopb_turb); - abundance( g[i].x[0],g[i].x[1],g[i].x[2], g[i].abun); - velocity( g[i].x[0],g[i].x[1],g[i].x[2], g[i].vel); + if(!allBitsSet(par->dataFlags, DS_mask_turb_doppler)){ + for(i=0;ipIntensity;i++) + doppler((*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2],&(*gp)[i].dopb_turb); + for(i=par->pIntensity;incell;i++) + (*gp)[i].dopb_turb=0.; + + par->dataFlags |= DS_mask_turb_doppler; } - /* Set velocity values also for sink points (otherwise Delaunay ray-tracing has problems) */ - for(i=par->pIntensity;incell;i++) - velocity(g[i].x[0],g[i].x[1],g[i].x[2],g[i].vel); + if(!allBitsSet(par->dataFlags, DS_mask_temperatures)){ + for(i=0;ipIntensity;i++) + temperature((*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2],(*gp)[i].t); + for(i=par->pIntensity;incell;i++){ + (*gp)[i].t[0]=par->tcmb; + (*gp)[i].t[1]=par->tcmb; + } - checkGridDensities(par, g); + par->dataFlags |= DS_mask_temperatures; + } - 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; + if(!allBitsSet(par->dataFlags, DS_mask_magfield)){ + if(par->polarization){ + for(i=0;ipIntensity;i++) + magfield((*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2],(*gp)[i].B); + + par->dataFlags |= DS_mask_magfield; + + }else{ + for(i=0;ipIntensity;i++){ + (*gp)[i].B[0]=0.0; + (*gp)[i].B[1]=0.0; + (*gp)[i].B[2]=0.0; + } + } + + for(i=par->pIntensity;incell;i++){ + (*gp)[i].B[0]=0.0; + (*gp)[i].B[1]=0.0; + (*gp)[i].B[2]=0.0; } } - getVelocities(par,g); - dumpGrid(par,g); - free(dc); + if(!allBitsSet(par->dataFlags, DS_mask_ACOEFF)){ + getVelocities(par,*gp); /* Mallocs and sets .v1, .v2, .v3 */ + + par->dataFlags |= DS_mask_ACOEFF; + } + + if(onlyBitsSet(par->dataFlags & DS_mask_all_but_mag, DS_mask_3)) /* Only happens if (i) we read no file and have constructed this data within LIME, or (ii) we read a file at dataStageI==3. */ + writeGridIfRequired(par, *gp, NULL, lime_FITS); - free(outRandLocations); - free(outRandDensities); - gsl_rng_free(rinc.randGen); + dumpGrid(par,*gp); + free(dc); } diff --git a/src/grid2fits.c b/src/grid2fits.c new file mode 100644 index 0000000..9af0a6e --- /dev/null +++ b/src/grid2fits.c @@ -0,0 +1,1617 @@ +/* + * grid2fits.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 + * +TODOS (some day): + - Change the type of g[i].id to unsigned int. + - Change the type of g[i].numNeigh to unsigned short. + - Change g[i].t, g[i].dopb and g[i].abun to float. + - Vector columns in defineGridExtColumns()? +*/ + +#include "lime.h" +#include "gridio.h" + +/* +The present module contains routines for transferring the LIME grid point data to or from a FITS format file. The purpose of the present comment block is to describe the FITS file format. Note that the amount and type of information stored depends on the 'data stage' of the grid struct, as described in the header remarks to module gridio.c. + +In the description below, the data stage bit associated with the presence of a particular extension, column or keyword is given on the leftmost place of each line. When the file is read, all the objects associated with a bit must be present for the bit to be set. + +Note that all extensions are binary table except where indicated. The letter in the second row for column descriptions gives the FITS data type. See eg + + https://heasarc.gsfc.nasa.gov/docs/software/fitsio/c/c_user/node20.html + +for a key to these. + +Where column names contain a lower-case letter, this is a placeholder for a digit as explained in the repective comment. + +0 0) The primary HDU + + Keywords: +0 RADIUS D # The model radius in metres. + +0 1) GRID + Number of rows = number of grid points. + + Keywords: +0 COLLPARn A # 1 for each nth collision partner. + + Columns: +0 ID V +0 Xj D # Cartesian components of the point location, 1 col per jth dimension. +0 IS_SINK L # =True iff the point lies on the edge of the model. +1 NUMNEIGH U +1 FIRST_NN V # See explanation in section 3 below. +2 VELj D # 1 col per jth dimension. +3 DENSITYn D # 1 per nth collision partner. +4 ABUNMOLm E # 1 per mth molecular species. +5 TURBDPLR E # Given Gaussian lineshape exp(-v^2/[B^2 + 2*k*T/m]), this is B. +6 TEMPKNTC E # From t[0]. +6 TEMPDUST E # From t[1]. + +1 2) NN_INDICES (see explanation in the header to gridio.c) + Number of rows = number of grid points * average number of Delaunay links per point. + + Columns: +1 LINK_I V + +1 3) LINKS (see explanation in the header to gridio.c) + Number of rows = number of Delaunay links. + + Columns: +1 GRID_I_1 V +1 GRID_I_2 V +7 V_p_j D # 1 per pth velocity sample per jth dimension. + +8 4 etc) LEVEL_POPS_m (1 per mth molecular species) + This is an image extension of size (number of grid cells)*(number of energy levels this species). + + Keywords: +8 MOL_NAME + +Note that at present the data in the 'partner' element of grid.mol is *NOT* being stored. +*/ + +/*....................................................................*/ +fitsfile * +openFITSFileForWrite(char *outFileName){ + fitsfile *fptr=NULL; + int status=0; + char negfile[100]="! "; + + fits_create_file(&fptr, outFileName, &status); + if(status!=0){ + if(!silent) warning("Overwriting existing fits file"); + status=0; + strcat(negfile,outFileName); + fits_create_file(&fptr, negfile, &status); + processFitsError(status); + } + + fits_create_img(fptr, 8, 0, NULL, &status); + processFitsError(status); + + return fptr; +} + +/*....................................................................*/ +void +initializeKeyword(struct keywordType *kwd){ + (*kwd).datatype = 0; + (*kwd).keyname = NULL; + (*kwd).comment = NULL; + (*kwd).intValue = 0; + (*kwd).floatValue = 0.0; + (*kwd).doubleValue = 0.0; + (*kwd).charValue = NULL; +} + +/*....................................................................*/ +void +closeFITSFile(fitsfile *fptr){ + int status=0; + + fits_close_file(fptr, &status); + processFitsError(status); +} + +/*....................................................................*/ +void +writeKeywordsToFits(lime_fptr *fptr, struct keywordType *kwds\ + , const int numKeywords){ + + int i, status; + char message[80]; + + for(i=0;imaxNumDims){ + if(!silent){ + sprintf(message, "Caller asked for %d dims but colnames can only be written for %d.", (int)gridInfo.nDims, (int)maxNumDims); + bail_out(message); + } + exit(1); + } + + if(gridInfo.nSpecies>maxNumSpecies){ + if(!silent){ + sprintf(message, "Caller asked for %d species but colnames can only be written for %d.", (int)gridInfo.nSpecies, (int)maxNumSpecies); + bail_out(message); + } + exit(1); + } + + if(gridInfo.nDensities>maxNumDensities){ + if(!silent){ + sprintf(message, "Caller asked for %d coll. part. but colnames can only be written for %d.", (int)gridInfo.nDensities, (int)maxNumDensities); + bail_out(message); + } + exit(1); + } + + *maxNumCols = 10 + gridInfo.nDims*2 + gridInfo.nSpecies + gridInfo.nDensities; + + *allColNames = malloc(sizeof(**allColNames) *(*maxNumCols)); + *allColNumbers = malloc(sizeof(**allColNumbers) *(*maxNumCols)); + tformAllCols = malloc(sizeof(*tformAllCols) *(*maxNumCols)); + tunitAllCols = malloc(sizeof(*tunitAllCols) *(*maxNumCols)); + dataTypeAllCols = malloc(sizeof(*dataTypeAllCols)*(*maxNumCols)); + + for(i=0;i<(*maxNumCols);i++){ /* Set default: */ + (*allColNames)[i] = malloc(sizeof(char)*numColNameChars); + (*allColNumbers)[i] = 0; + } + + colToWriteI = 0; + colI = 0; + + sprintf((*allColNames)[colI], "ID"); + if(bitIsSet(dataFlags, DS_bit_x)){ + colToWriteI++; + (*allColNumbers)[colI] = colToWriteI; + } + tformAllCols[colI] = "V"; + tunitAllCols[colI] = "\0"; + dataTypeAllCols[colI] = TUINT; + + /* should rather have a vector column? */ + for(i_us=0;i_us0){ + ttype[i-1] = (*allColNames)[colI]; + tform[i-1] = tformAllCols[colI]; + tunit[i-1] = tunitAllCols[colI]; + (*colDataTypes)[i-1] = dataTypeAllCols[colI]; + } + } + + /* Append a new empty binary table onto the FITS file. + */ + fits_create_tbl(fptr, BINARY_TBL, 0, colToWriteI, ttype, tform, tunit, "GRID", &status); + processFitsError(status); + + free(tformAllCols); + free(tunitAllCols); + free(dataTypeAllCols); +} + +/*....................................................................*/ +void +writeGridExtToFits(fitsfile *fptr, struct gridInfoType gridInfo\ + , struct grid *gp, unsigned int *firstNearNeigh\ + , char **collPartNames, const int dataFlags){ + /* +This writes whatever information is in the grid struct (as specified by the dataFlags) and which also has a dimensionality which is some simple multiple of the number of grid points, to a single FITS binary table extension called GRID. The function tries to be fairly forgiving of screwy situations but it will exit if the minimum information is not present (defined as allBitsSet(dataFlags, DS_mask_x), which implies that elements .id, .x and .sink should all contain valid values). + +Note that data types in all capitals are defined in fitsio.h. + */ + + const unsigned int totalNumGridPoints = gridInfo.nInternalPoints+gridInfo.nSinkPoints; + const unsigned short numKwdChars=9; /* 8 characters + \0. */ + const unsigned short numColNameChars=21; /* 20 characters + \0. */ + const unsigned short maxNumCollPart = 9; + unsigned int *ids=NULL,i_ui; + double *xj=NULL; + _Bool *sink=NULL;/* Probably should be char* but this seems to work. */ + unsigned short *numNeigh=NULL,i_us,localNumCollPart; + double *velj=NULL,*densn=NULL; + float *dopb=NULL, *t=NULL, *abunm=NULL, *bField=NULL; + int status=0, colI=0, i, di, maxNumCols; + LONGLONG firstRow=1, firstElem=1; + char genericComment[80]; + char genericKwd[numKwdChars], message[80]; + char colName[numColNameChars]; + char **allColNames=NULL; + int *allColNumbers=NULL, *colDataTypes=NULL; + + if(!allBitsSet(dataFlags, DS_mask_x)){ + if(!silent) bail_out("Data stage indicates no grid data!"); + exit(1); + } + + /* +Ok we have a bit of a tricky situation here in that the number of columns we write is going to depend on the information available in gp, as encoded in the dataFlags. We need to work out which columns we are going to write ahead of time because we need the appropriate data on ALL the columns to set up the table size before we can start to write their individual data values. I also want to avoid checking dataFlags twice in two different contexts - that is how errors arise. So I've arranged that the following routine will do all the donkey work of setting up only those columns we can write and then using that information to define the table size. The routine also returns three more vectors: + + allColNames - This is returned to remove what would otherwise be a hard-wired dependence that the column ordering was the same in the present routine as in defineAndLoadColumns(). With this vector, the present routine can search for a column name in it and then use the returned vector index to access (from the next vector) the number of the column in the (smaller) sequence of valid columns. + + allColNumbers - This contains the number of a column in the sequence (beginning at 1) of those for which gp has data. If gp contains no data for a given column name, its entry in allColNumbers will be 0. + + colDataTypes - This just contains the data types for the valid columns. + */ + defineAndLoadColumns(fptr, gridInfo\ + , dataFlags, numColNameChars, &allColNames, &allColNumbers, &maxNumCols, &colDataTypes); + + /* Write the columns: + */ + colI = allColNumbers[getColIndex(allColNames, maxNumCols, "ID")]; + if(colI<=0){ + if(!silent) bail_out("This should not occur, it is some sort of bug."); + exit(1); + } + ids = malloc(sizeof(*ids)*totalNumGridPoints); + for(i_ui=0;i_ui0){ + numNeigh = malloc(sizeof(*numNeigh)*totalNumGridPoints); + for(i_ui=0;i_ui0){ + velj = malloc(sizeof(*velj)*totalNumGridPoints); + for(i_us=0;i_us0){ + densn = malloc(sizeof(*densn)*totalNumGridPoints); + for(i_us=0;i_us0){ + abunm = malloc(sizeof(*abunm)*totalNumGridPoints); + for(i_us=0;i_us0){ + dopb = malloc(sizeof(*dopb)*totalNumGridPoints); + for(i_ui=0;i_ui0){ + t = malloc(sizeof(*t)*totalNumGridPoints); + for(i_ui=0;i_ui0){ + bField = malloc(sizeof(*bField)*totalNumGridPoints); + for(di=0;di<3;di++){//**** keep this hard-wired or rather test that gridInfo.nDims==3?? + sprintf(colName, "B_FIELD%d", di+1); + colI = allColNumbers[getColIndex(allColNames, maxNumCols, colName)]; + if(colI<=0){ + if(!silent) bail_out("This should not occur, it is some sort of bug."); + exit(1); + } + + for(i_ui=0;i_uimaxNumCollPart){ + if(!silent){ + sprintf(message, "There seem to be %d collision partners but keywords can only be written for %d.", (int)gridInfo.nDensities, (int)maxNumCollPart); + warning(message); + } + localNumCollPart = maxNumCollPart; + }else{ + localNumCollPart = gridInfo.nDensities; + } + + for(i_us=0;i_us bool */ + + while(iid; + fits_write_col(fptr, dataType, 1, firstRow, firstElem, (LONGLONG)gridInfo.nNNIndices, linkIs, &status); + processFitsError(status); + free(linkIs); +} + +/*....................................................................*/ +void +writeLinksExtToFits(fitsfile *fptr, struct gridInfoType gridInfo\ + , struct linkType *links){ + /* +See the comment at the beginning of module gridio.c for a description of how the LINKS extension relates to the grid struct. + +Notes: + - Data types in all capitals are defined in fitsio.h. + + - The business where the column names are first loaded into tempColNames, then the pointers to each one into ttype, is done because: + . fits_create_tbl() demands a 5th argument of type char*[] and will not accept char**; + + . sprintf() cannot be used to write a string to an unallocated char*, thus something like the following is illegal: + sprintf(ttype[colI], "ACOEFF_%d", n+1); + + . The following is legal, but cannot be used for all column names, because we have to generate some on the fly: + ttype[colI] = ; + + */ + const int numColNameChars=21; + unsigned int *ids=NULL,i_ui; + unsigned short i_us,j_us; + double *vels=NULL; + int status=0, colI=0, i; + LONGLONG firstRow=1, firstElem=1; + int numCols; + char extname[] = "LINKS"; + char **tempColNames=NULL; + + if(links==NULL){ + if(!silent) bail_out("No link data!"); + exit(1); + } + + if(links[0].vels==NULL) + numCols = 2; + else + numCols = 2 + gridInfo.nDims*gridInfo.nLinkVels; + + /* Define the name, datatype, and physical units for the columns. + */ + char *ttype[numCols]; + char *tform[numCols]; + char *tunit[numCols]; + int dataTypes[numCols]; + + tempColNames = malloc(sizeof(*tempColNames)*numCols); + for(i=0;inDims = (unsigned short)countColsBasePlusInt(fptr, "X"); + if(gridInfoRead->nDims<=0){ + if(!silent) warning("No X columns found in grid dataset."); + return; /* I.e. with dataFlags left unchanged. */ + } + + /* We have to do this here (as well after the call to readGrid()) because grid.x is a pre-sized array rather than a pointer we can malloc. Later this should be changed to allow us to define the sizes of all arrays in grid purely from the data in the file. + */ + if(gridInfoRead->nDims!=DIM){ + if(!silent){ + sprintf(message, "%d Xn columns read, but there should be %d.", (int)gridInfoRead->nDims, DIM); + bail_out(message); + } + exit(1); + } + + xj = malloc(sizeof(*xj)*numGridCells); + for(i_us=0;i_usnDims;i_us++){ + sprintf(colName, "X%d", (int)i_us+1); + fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); + processFitsError(status); + + fits_read_col(fptr, TDOUBLE, colNum, firstRow, firstElem, numGridCells, 0, xj, &anynul, &status); + processFitsError(status); + + for(i_LL=0;i_LLnSinkPoints = 0; + for(i_LL=0;i_LLnSinkPoints++; + } + free(sink); + + gridInfoRead->nInternalPoints = numGridCells - gridInfoRead->nSinkPoints; + + /* If we have made it this far, we can set the first bit of dataFlags. Woot! + */ + (*dataFlags) |= (1 << DS_bit_x); + + fits_get_colnum(fptr, CASEINSEN, "NUMNEIGH", &colNum, &status); + if(status!=COL_NOT_FOUND){ + processFitsError(status); + + numNeigh = malloc(sizeof(*numNeigh)*numGridCells); + fits_read_col(fptr, TUSHORT, colNum, firstRow, firstElem, numGridCells, 0, numNeigh, &anynul, &status); + processFitsError(status); + + for(i_LL=0;i_LLnDims;i_us++){ + sprintf(colName, "VEL%d", (int)i_us+1); + fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); + processFitsError(status); + + fits_read_col(fptr, TDOUBLE, colNum, firstRow, firstElem, numGridCells, 0, velj, &anynul, &status); + processFitsError(status); + + for(i_LL=0;i_LLnDensities = (unsigned short)countColsBasePlusInt(fptr, "DENSITY"); + if(gridInfoRead->nDensities > 0){ + for(i_LL=0;i_LLnDensities); + + /* Read the DENSITY columns: + */ + densn = malloc(sizeof(*densn)*numGridCells); + for(i_us=0;i_usnDensities;i_us++){ + sprintf(colName, "DENSITY%d", (int)i_us+1); + fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); + processFitsError(status); + + fits_read_col(fptr, TDOUBLE, colNum, firstRow, firstElem, numGridCells, 0, densn, &anynul, &status); + processFitsError(status); + + for(i_LL=0;i_LLnSpecies = (unsigned short)countColsBasePlusInt(fptr, "ABUNMOL"); + if(gridInfoRead->nSpecies > 0){ + for(i_LL=0;i_LLnSpecies); + } + + /* Read the ABUNMOL columns: + */ + abunm = malloc(sizeof(*abunm)*numGridCells); + for(i_us=0;i_usnSpecies;i_us++){ + sprintf(colName, "ABUNMOL%d", (int)i_us+1); + fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); + processFitsError(status); + + fits_read_col(fptr, TFLOAT, colNum, firstRow, firstElem, numGridCells, 0, abunm, &anynul, &status); + processFitsError(status); + + for(i_LL=0;i_LLbool + + if(!bitIsSet(*dataFlags, DS_bit_neighbours)) + return; + + totalNumGridPoints = gridInfoRead->nInternalPoints + gridInfoRead->nSinkPoints; /* Just for a bit more brevity. */ + + /* Go to the LINKS extension. + */ + fits_movnam_hdu(fptr, BINARY_TBL, "LINKS", 0, &status); + if(status==BAD_HDU_NUM){ + if(!silent) warning("No LINKS extension found in grid dataset."); + + /* Unset the DS_mask_neighbours bit: */ + *dataFlags &= ~(1 << DS_bit_neighbours); + return; + } + processFitsError(status); + + /* Find out how many rows there are, then malloc the array. + */ + fits_get_num_rowsll(fptr, &totalNumLinks, &status); + processFitsError(status); + + gridInfoRead->nLinks = (unsigned int)totalNumLinks; + if(gridInfoRead->nLinks<=0){ + if(!silent) warning("No rows in LINKS extension of grid dataset."); + + /* Unset the neighbours bit: */ + *dataFlags &= ~(1 << DS_bit_neighbours); + return; + } + *links = malloc(sizeof(**links)*totalNumLinks); + + for(i_LL=0;i_LL=totalNumGridPoints){ + if(!silent){ + sprintf(message, "GRID_I_1 %dth-row value %ud is outside range [0,%ud]", (int)i_LL, i_ui, totalNumGridPoints); + bail_out(message); + } + exit(1); + } + (*links)[i_LL].gis[0] = gp[i_ui].id; + } + + /* Read GRID_I_2 column. + */ + fits_get_colnum(fptr, CASEINSEN, "GRID_I_2", &colNum, &status); + if(status!=COL_NOT_FOUND){ + processFitsError(status); + colGrid2Found = 1; + + fits_read_col(fptr, TUINT, colNum, firstRow, firstElem, totalNumLinks, 0, ids, &anynul, &status); + processFitsError(status); + + for(i_LL=0;i_LL=totalNumGridPoints){ + if(!silent){ + sprintf(message, "GRID_I_2 %dth-row value %ud is outside range [0,%ud]", (int)i_LL, i_ui, totalNumGridPoints); + bail_out(message); + } + exit(1); + } + (*links)[i_LL].gis[1] = gp[i_ui].id; + } + } + free(ids); + } + + if(!(colGrid1Found && colGrid2Found)){ + if(!silent) warning("Both GRID_I columns not found in LINKS extension of grid dataset."); + free(*links); + *links = NULL; + /* Unset the neighbours bit: */ + *dataFlags &= ~(1 << DS_bit_neighbours); + return; + } + + /* Find out how many V_* columns there are. + */ + i = 0; + status = 0; + while(!status){ + sprintf(colName, "V_%d_", i+1); + if(countColsBasePlusInt(fptr, colName)!=gridInfoRead->nDims) + status = 1; + + i++; + } + gridInfoRead->nLinkVels = i - 1; + status = 0; + + if(gridInfoRead->nLinkVels<=0){ + for(i_LL=0;i_LLnLinkVels*gridInfoRead->nDims); + + vels = malloc(sizeof(*vels)*totalNumLinks); + for(i_us=0;i_usnLinkVels;i_us++){ + for(j_us=0;j_usnDims;j_us++){ + /* Read the V_n_d columns. + */ + sprintf(colName, "V_%d_%d", (int)i_us+1, (int)j_us+1); + fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); + processFitsError(status); + + fits_read_col(fptr, TDOUBLE, colNum, firstRow, firstElem, totalNumLinks, 0, vels, &anynul, &status); + processFitsError(status); + + for(i_LL=0;i_LLnDims*j_us + i_us] = vels[i_LL]; + } + } + free(vels); + + (*dataFlags) |= (1 << DS_bit_ACOEFF); + +} + +/*....................................................................*/ +void +readNnIndicesExtFromFits(fitsfile *fptr, struct linkType *links\ + , struct linkType ***nnLinks, struct gridInfoType *gridInfoRead, int *dataFlags){ + /* +See the comment at the beginning of gridio.c for a description of how the NN_INDICES extension relates to the grid struct. + +The function mallocs the pointer *nnLinks. + */ + + LONGLONG totalNumNeigh, firstRow=1, firstElem=1, i_LL; + int status=0, colNum, anynul=0; + unsigned int *linkIs=NULL; + + if(!bitIsSet(*dataFlags, DS_bit_neighbours)) + return; + + /* Go to the NN_INDICES extension. + */ + fits_movnam_hdu(fptr, BINARY_TBL, "NN_INDICES", 0, &status); + if(status==BAD_HDU_NUM){ + if(!silent) warning("No NN_INDICES extension found in grid dataset."); + + /* Unset the neighbours bit: */ + *dataFlags &= ~(1 << DS_bit_neighbours); + return; + } + processFitsError(status); + + /* Find out how many rows there are, then malloc the array. + */ + fits_get_num_rowsll(fptr, &totalNumNeigh, &status); + processFitsError(status); + + gridInfoRead->nNNIndices = (unsigned int)totalNumNeigh; + if(gridInfoRead->nNNIndices<=0){ + if(!silent) warning("No rows in NN_INDICES extension of grid dataset."); + + /* Unset the neighbours bit: */ + *dataFlags &= ~(1 << DS_bit_neighbours); + return; + } + *nnLinks = malloc(sizeof(**nnLinks)*totalNumNeigh); + + /* Read LINK_I column. + */ + fits_get_colnum(fptr, CASEINSEN, "LINK_I", &colNum, &status); + if(status==COL_NOT_FOUND){ + if(!silent) warning("LINK_I column not found in NN_INDICES extension of grid dataset."); + + /* Unset the neighbours bit: */ + *dataFlags &= ~(1 << DS_bit_neighbours); + return; + } + processFitsError(status); + + linkIs = malloc(sizeof(*linkIs)*totalNumNeigh); + fits_read_col(fptr, TUINT, colNum, firstRow, firstElem, totalNumNeigh, 0, linkIs, &anynul, &status); + processFitsError(status); + + for(i_LL=0;i_LLmaxNumSpecies){ + if(!silent){ + sprintf(message, "Species block %d is greater than the limit %d", (int)speciesI+1, (int)maxNumSpecies); + bail_out(message); + } + exit(1); + } + + sprintf(extname, "LEVEL_POPS_%d", (int)speciesI+1); + + /* Try to move to extension LEVEL_POPS_: + */ + fits_movnam_hdu(fptr, IMAGE_HDU, extname, 0, &status); + if(status==BAD_HDU_NUM) + return 0; + + /* If we reached here it means we found the extension, but now we check for other errors: + */ + processFitsError(status); + + return 1; +} + +/*....................................................................*/ +void +readPopsExtFromFits(fitsfile *fptr, const unsigned short speciesI\ + , struct grid *gp, struct gridInfoType *gridInfoRead){ + /* +See the comment at the beginning of the present module for a description of how the LEVEL_POPS_m extensions relate to the grid struct. + +The function mallocs gp[yi].mol[speciesI].pops for each grid point yi and species. + */ + + const int maxLenMolName = 8; + const unsigned short maxNumSpecies = 9; + int naxis, status=0, xi, anynul=0, bitpix; +/*The interface to fits_get_img_param() says long* for argument 5 but I get a seg fault unless I use an array. + long *naxes; */ +long naxes[2]; + float *row=NULL; + long fpixels[2],lpixels[2]; + long inc[2] = {1,1}; + char molNameRead[maxLenMolName+1]; + char message[80]; + char extname[13]; + unsigned int numGridPoints, i_ui; + + if(speciesI+1>maxNumSpecies){ + if(!silent){ + sprintf(message, "Species block %d is greater than the limit %d", (int)speciesI+1, (int)maxNumSpecies); + bail_out(message); + } + exit(1); + } + sprintf(extname, "LEVEL_POPS_%d", (int)speciesI+1); + + /* Try to move to extension LEVEL_POPS_: + */ + fits_movnam_hdu(fptr, IMAGE_HDU, extname, 0, &status); + processFitsError(status); + + fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status); + processFitsError(status); + + numGridPoints = gridInfoRead->nInternalPoints + gridInfoRead->nSinkPoints; + if((long)numGridPoints != naxes[1]){ + if(!silent){ + sprintf(message, "Expected %ld grid points but extension %s has %ld"\ + , (long)(gridInfoRead->nInternalPoints + gridInfoRead->nSinkPoints), extname, naxes[1]); + bail_out(message); + } + exit(1); + } + gridInfoRead->mols[speciesI].nLevels = (int)naxes[0]; + gridInfoRead->mols[speciesI].nLines = -1; +/* free(naxes); */ + + row = malloc(sizeof(*row)*gridInfoRead->mols[speciesI].nLevels); + + /* Read FITS data. + */ + for(i_ui=0;i_uimols[speciesI].nLevels; + lpixels[1]=(int)i_ui+1; + + fits_read_subset(fptr, TFLOAT, fpixels, lpixels, inc, 0, row, &anynul, &status); + processFitsError(status); + + gp[i_ui].mol[speciesI].pops = malloc(sizeof(double)*gridInfoRead->mols[speciesI].nLevels); + for(xi=0;ximols[speciesI].nLevels;xi++) + gp[i_ui].mol[speciesI].pops[xi] = (double)row[xi]; + + } + + free(row); + + /* Read kwds: + */ + fits_read_key(fptr, TSTRING, "MOL_NAME", molNameRead, NULL, &status); + gridInfoRead->mols[speciesI].molName = molNameRead;//*****?? + processFitsError(status); +} + + diff --git a/src/grid2fits.h b/src/grid2fits.h new file mode 100644 index 0000000..d4b6157 --- /dev/null +++ b/src/grid2fits.h @@ -0,0 +1,42 @@ +/* + * grid2fits.h + * This file is part of LIME, the versatile line modeling engine + * + * Copyright (C) 2006-2014 Christian Brinch + * Copyright (C) 2015-2016 The LIME development team + * + */ + +#ifndef GRID2FITS_H +#define GRID2FITS_H + +#define lime_fptr fitsfile + +struct keywordType{ + int datatype; /* Accepted: TSTRING, TINT, TFLOAT, TDOUBLE */ + char *keyname, *comment; + int intValue; + float floatValue; + double doubleValue; + char *charValue; +}; + + +_Bool checkPopsFitsExtExists(fitsfile*, const unsigned short); +void closeFITSFile(fitsfile*); +void initializeKeyword(struct keywordType*); +void readGridExtFromFits(fitsfile*, struct gridInfoType*, struct grid**, unsigned int**, char***, int*, int*); +void readKeywordsFromFits(lime_fptr*, struct keywordType*, const int); +void readLinksExtFromFits(fitsfile*, struct gridInfoType*, struct grid*, struct linkType**, int*); +void readNnIndicesExtFromFits(fitsfile*, struct linkType*, struct linkType***, struct gridInfoType*, int*); +void readPopsExtFromFits(fitsfile*, const unsigned short, struct grid*, struct gridInfoType*); +fitsfile *openFITSFileForRead(char*); +fitsfile *openFITSFileForWrite(char*); +void writeKeywordsToFits(lime_fptr*, struct keywordType*, const int); +void writeGridExtToFits(fitsfile*, struct gridInfoType, struct grid*, unsigned int*, char**, const int); +void writeLinksExtToFits(fitsfile*, struct gridInfoType, struct linkType*); +void writeNnIndicesExtToFits(fitsfile*, struct gridInfoType, struct linkType**);//, struct linkType*); +void writePopsExtToFits(fitsfile*, struct gridInfoType, const unsigned short, struct grid*); + +#endif /* GRID2FITS_H */ + diff --git a/src/gridio.c b/src/gridio.c new file mode 100644 index 0000000..543dc62 --- /dev/null +++ b/src/gridio.c @@ -0,0 +1,855 @@ +/* + * gridio.c + * This file is part of LIME, the versatile line modeling engine + * + * Copyright (C) 2006-2014 Christian Brinch + * Copyright (C) 2016 The LIME development team + * +TODO: + */ + +#include "lime.h" +#include "gridio.h" + +/* +This module contains generic routines for writing grid data to, and reading it from, a file on disk. The problem with doing this is that the grid struct contains different amounts of information at different times in the running of the code. In order to quantify and regulate this, a dataFlags integer is used to record the presence or absence (as indicated by the value of the appropriate bit in the mask) of particular types of information. The bits associated with certain fields of struct grid are given in the lime.h header. + +Notes: + - If several fields of the struct are listed for a given bit, a value of 1 for that bit indicates that all fields are expected to be present, and a value of 0 indicates that all fields of that bit will be ignored. + - Only a few combinations are not allowed, as follows: + * No other bit may be set if DS_bit_x is not. + * DS_bit_ACOEFF may not be set if either DS_bit_neighbours or DS_bit_velocity is not. + * DS_bit_populations may not be set unless all the others are set as well (except DS_bit_magfield may be set or not). + +Writing the grid file: +---------------------- +To make things simpler, four stages have been defined at which the user may write the grid data to file. These are described in the following table: + + Data mask bits set: stage 0 stage 1 stage 2 stage 3 stage 4 + ..................................................................................................... + DS_bit_x 0 1 1 1 1 + DS_bit_neighbours 0 0 1 1 1 + DS_bit_velocity 0 0 0 1 1 + DS_bit_density 0 0 0 1 1 + DS_bit_abundance 0 0 0 1 1 + DS_bit_turb_doppler 0 0 0 1 1 + DS_bit_temperatures 0 0 0 1 1 + DS_bit_magfield 0 0 0 x x + DS_bit_ACOEFF 0 0 0 1 1 + DS_bit_populations 0 0 0 0 1 + ..................................................................................................... + +Notes: + - dataStageI==0 has been included for completeness/robustness but the user may not write a file with nothing in it. + - LIME may run to completion without ever reaching stage 4 - if all the images required were continuum ones, for example. + +A note about the vectors 'links', 'nnLinks' and 'firstNearNeigh': +----------------------------------------------------------------- +The 'neigh' element of struct grid is a clunky way to store the information about the nearest neighbours of each grid point, because it leads to a host of separate little mallocs, thus memory for the grid is spread all over the place. There is also some duplication of quantities which refer to the link between 2 grid points rather than to the points themselves. This arrangement does not lend itself easily to storage in a file. + +To meet the needs of file storage we divide the information contained in 'neigh' between 3 new vectors 'links', 'nnLinks' and 'firstNearNeigh'. The 1st of these is a list of all Delaunay edges between grid points. The 2nd contains indices to the 1st, and is arranged such that indices to all the links connected to a single grid point occur in a connected sequence. The 3rd vector has the same number of entries as the total number of grid points. It contains the index of the first entry in the 2nd vector which corresponds to that grid point. Thus in order, for example, to iterate over all links connected to the ith grid point, we need only do something as follows: + + for(j=0;jid; + (*firstNearNeigh)[idA] = ni; + + for(jA=0;jAnumNeigh;jA++){ + /* Check to see if the NN has been done. + */ + gBPtr = gAPtr->neigh[jA]; + idB = gBPtr->id; + if(pointIsDone[idB]){ + /* The link exists; we need to find a pointer to it for the next entry of nnLinks. To do that, we need to go through the list of NN links of the point whose id is idB and identify that link which connects back to idA. + */ + linkNotFound = 1; /* default */ + for(jB=0;jBnumNeigh;jB++){ + trialIdA = gBPtr->neigh[jB]->id; + if(trialIdA==idA){ + linkNotFound = 0; + k = (*firstNearNeigh)[idB] + jB; + linkId = nnLinkIs[k]; + } + } + if(linkNotFound){ + if(!silent) bail_out("Link not found."); + exit(1); + } + + }else{ + /* The link does not yet exist; we must create it, load it into links, and find a pointer to it for the next entry of nnLinks. + */ + linkId = li; + (*links)[li].id = linkId; + (*links)[li].gis[0] = idA; + (*links)[li].gis[1] = idB; + + li++; + } + + nnLinkIs[ni] = linkId; + ni++; + } + + pointIsDone[idA] = 1; + } + *totalNumLinks = li; + + *links = realloc(*links, sizeof(struct linkType)*(*totalNumLinks)); + + if(allBitsSet(dataFlags, DS_mask_ACOEFF)){ + for(i_ui=0;i_ui<*totalNumLinks;i_ui++) + (*links)[i_ui].vels = malloc(sizeof(double)*gridInfo.nLinkVels*gridInfo.nDims); + + for(i_ui=0;i_ui<*totalNumLinks;i_ui++){ + gi0 = (*links)[i_ui].gis[0]; + gi1 = (*links)[i_ui].gis[1]; + if(gp[gi0].sink && gp[gi1].sink){ + for(i_us=0;i_usnumNeigh;jA++){ + idB = gAPtr->neigh[jA]->id; + if(linkNotFound && idB==(*links)[i_ui].gis[1-nearI]){ + if(nearI==0){ + for(j_us=0;j_usv1[gridInfo.nDims*jA + j_us]; + (*links)[i_ui].vels[gridInfo.nDims*1 + j_us] = gAPtr->v2[gridInfo.nDims*jA + j_us]; + (*links)[i_ui].vels[gridInfo.nDims*2 + j_us] = gAPtr->v3[gridInfo.nDims*jA + j_us]; + } + }else{ + for(j_us=0;j_usv3[gridInfo.nDims*jA + j_us]; + (*links)[i_ui].vels[gridInfo.nDims*1 + j_us] = gAPtr->v2[gridInfo.nDims*jA + j_us]; + (*links)[i_ui].vels[gridInfo.nDims*2 + j_us] = gAPtr->v1[gridInfo.nDims*jA + j_us]; + } + } + linkNotFound = 0; + } + } + if(linkNotFound){ + sprintf(message, "SNAFU with link %d grid point %d", (int)i_ui, gAPtr->id); + if(!silent) bail_out(message); + exit(1); + } + } /* end if not both sink */ + } /* end loop over links */ + }else{ /* a coeffs not available */ + for(i_ui=0;i_ui<*totalNumLinks;i_ui++) + (*links)[i_ui].vels = NULL; + } + + *nnLinks = malloc(sizeof(**nnLinks)*(*totalNumNeigh)); + for(i_ui=0;i_ui<(*totalNumNeigh);i_ui++) + (*nnLinks)[i_ui] = &(*links)[nnLinkIs[i_ui]]; + + free(nnLinkIs); + free(pointIsDone); + /* The calling routine must free firstNearNeigh, nnLinks, links. */ +} + +/*....................................................................*/ +void +closeFile(lime_fptr *fptr, const int fileFormatI){ + if(fileFormatI==lime_FITS) + closeFITSFile(fptr); + //**** error if not? +} + +/*....................................................................*/ +void +closeAndFree(lime_fptr *fptr, const int fileFormatI\ + , unsigned int *firstNearNeigh, struct linkType **nnLinks\ + , struct linkType *links, const unsigned int totalNumLinks){ + + unsigned int i_ui; + + closeFile(fptr, fileFormatI); + + free(firstNearNeigh); + free(nnLinks); + if(links!=NULL){ + for(i_ui=0;i_uigis[0]==i_ui){ + gp[i_ui].neigh[j] = &gp[linkPtr->gis[1]]; + }else{ + gp[i_ui].neigh[j] = &gp[linkPtr->gis[0]]; + } + } + } +} + +/*....................................................................*/ +void +loadLinkVelsIntoGrid(unsigned int *firstNearNeigh, struct linkType **nnLinks\ + , struct gridInfoType gridInfoRead, struct grid *gp){ + /* +See the comment at the beginning of the present module for a description of how the pointers 'links', 'nnLinks' and 'firstNearNeigh' relate to the grid struct. + +The function mallocs (then writes to) the following extensions of struct g for each grid point: + v1 + v2 + v3 + +The values are obtained from those stored in the .vels extension of elements of a list of edges (aka links) which are pointed to by elements of the input argument nnLinks. Storing these values in the grid points is x2 redundant, because if grid points A and B are at the ends of a given edge, the same set of values is stored in each. + +There is a subtle point about ordering to grasp. The order of the .vels entries in each edge element is prescribed to be the same as the order as stored in the vertex which is listed first in the edge's .gis attribute, i.e. .gis[0]. As we work through the list of grid points, if a given point corresponds to .gis[0] of the edge in question, obvisously we just copy over the velocities; but if it corresponds to .gis[1], we have to reverse the velocity order. + */ + + const unsigned int totalNumGridPoints = gridInfoRead.nInternalPoints+gridInfoRead.nSinkPoints; + unsigned int i_ui; + unsigned short i_us; + int j; + struct linkType *linkPtr; + char message[80]; + + /* Just for the time being: */ + if(gridInfoRead.nLinkVels!=NUM_VEL_COEFFS){ + if(!silent){ + sprintf(message, "There should be %d VEL_n columns, but %d were read.", NUM_VEL_COEFFS, (int)gridInfoRead.nLinkVels); + bail_out(message); + } + exit(1); + } + + for(i_ui=0;i_uigis[0]==i_ui){ + for(i_us=0;i_usvels)[gridInfoRead.nDims*0 + i_us]; + gp[i_ui].v2[gridInfoRead.nDims*j+i_us] = (linkPtr->vels)[gridInfoRead.nDims*1 + i_us]; + gp[i_ui].v3[gridInfoRead.nDims*j+i_us] = (linkPtr->vels)[gridInfoRead.nDims*2 + i_us]; + } + }else{ + for(i_us=0;i_usvels)[gridInfoRead.nDims*2 + i_us]; + gp[i_ui].v2[gridInfoRead.nDims*j+i_us] = (linkPtr->vels)[gridInfoRead.nDims*1 + i_us]; + gp[i_ui].v3[gridInfoRead.nDims*j+i_us] = (linkPtr->vels)[gridInfoRead.nDims*0 + i_us]; + } + } + } + } +} + +/*....................................................................*/ +int +checkPopsTableExists(lime_fptr *fptr, const int fileFormatI\ + , const unsigned short speciesI, _Bool *blockFound){ + int status=0; + + *blockFound = 0; + + if(fileFormatI==lime_FITS){ + *blockFound = checkPopsFitsExtExists(fptr, speciesI); + }else{ + status = 1; + } + + return status; +} + +/*....................................................................*/ +int +getNumPopsTables(lime_fptr *fptr, const int fileFormatI, unsigned short *numTables){ + int status = 0; + _Bool blockFound = 1; + + *numTables = 0; + while(blockFound && !status){ + status = checkPopsTableExists(fptr, fileFormatI, *numTables, &blockFound); + (*numTables)++; + } + + (*numTables)--; + + return status; +} + +/*....................................................................*/ +int +readPopsTable(lime_fptr *fptr, const int fileFormatI\ + , const unsigned short speciesI, struct grid *gp\ + , struct gridInfoType *gridInfoRead){ + + int status=0; + unsigned int totalNumGridPoints, i_ui; + + totalNumGridPoints = gridInfoRead->nInternalPoints + gridInfoRead->nSinkPoints; + for(i_ui=0;i_uinInternalPoints = 0; + gridInfoRead->nSinkPoints = 0; + gridInfoRead->nLinks = 0; + gridInfoRead->nNNIndices = 0; + gridInfoRead->nDims = 0; + gridInfoRead->nSpecies = 0; + gridInfoRead->nDensities = 0; + gridInfoRead->nLinkVels = 0; + gridInfoRead->mols = NULL; + + /* Open the file and also return the data stage. */ + fptr = openFileForRead(inFileName, fileFormatI); + + /* Read header keywords: + */ + status = readKeywords(fptr, fileFormatI, primaryKwds, numKeywords); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, 0); + return 1; + } + + /* Read the values which should be in grid for every stage. + */ + status = readGridTable(fptr, fileFormatI, gridInfoRead, gp, &firstNearNeigh\ + , collPartNames, numCollPartRead, dataFlags); /* Sets appropriate bits of dataFlags; also mallocs gp and sets all its defaults. */ + totalNumGridPoints = gridInfoRead->nSinkPoints + gridInfoRead->nInternalPoints; + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, 0); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 2; + } + + /* Some sanity checks: + */ + if((*dataFlags)!=0 && gridInfoRead->nSinkPoints<=0 || gridInfoRead->nInternalPoints<=0 || gridInfoRead->nDims<=0){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, 0); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + + if(gridInfoRead->nSinkPoints<=0) + return 3; + else if(gridInfoRead->nInternalPoints<=0) + return 4; + else if(gridInfoRead->nDims<=0) + return 5; + else{ + sprintf(message, "This indicates a programming error. Please contact the developer."); + if(!silent) bail_out(message); + exit(1); + } + } + + /* Some more sanity checks: + */ + if((allBitsSet(*dataFlags, DS_mask_density) && gridInfoRead->nDensities<=0)\ + || (allBitsSet(*dataFlags, DS_mask_abundance) && gridInfoRead->nSpecies<=0)){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, 0); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + + if(gridInfoRead->nDensities<=0) + return 6; + else if(gridInfoRead->nSpecies<=0) //***** what if all continuum images?? + return 7; + else{ + sprintf(message, "This indicates a programming error. Please contact the developer."); + if(!silent) bail_out(message); + exit(1); + } + } + + status = readLinksTable(fptr, fileFormatI, gridInfoRead, *gp, &links, dataFlags); /* Sets appropriate bits of dataFlags. */ + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 8; + } + + /* Sanity check: + */ + if (allBitsSet(*dataFlags, DS_mask_ACOEFF) && gridInfoRead->nLinkVels<=0){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 9; + } + + status = readNnIndicesTable(fptr, fileFormatI, links, &nnLinks, gridInfoRead, dataFlags); /* Sets appropriate bits of dataFlags. */ + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 10; + } + + if(allBitsSet(*dataFlags, DS_mask_neighbours)){ + /* Convert the NN information back to the standard LIME grid struct format. + */ + loadNnIntoGrid(firstNearNeigh, nnLinks, *gridInfoRead, *gp); /* mallocs extension 'neigh' of struct g for each grid point. */ + + if(allBitsSet(*dataFlags, DS_mask_ACOEFF)) + loadLinkVelsIntoGrid(firstNearNeigh, nnLinks, *gridInfoRead, *gp); + } + + status = getNumPopsTables(fptr, fileFormatI, &numTables); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 11; + } + + /* Sanity check: + */ + if(numTables>0 && numTables!=gridInfoRead->nSpecies){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 12; + } + + if(numTables>0){ + (*dataFlags) |= (1 << DS_bit_populations); + + for(i_ui=0;i_uinSpecies); + + gridInfoRead->mols = malloc(sizeof(struct molInfoType)*gridInfoRead->nSpecies); + for(i_us=0;i_usnSpecies;i_us++){ + gridInfoRead->mols[i_us].molName = NULL; + gridInfoRead->mols[i_us].nLevels = -1; + gridInfoRead->mols[i_us].nLines = -1; + + status = readPopsTable(fptr, fileFormatI, i_us, *gp, gridInfoRead); /* Sets defaults for all the fields under grid.mol. */ + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 13; + } + + /* Sanity check: + */ + if(gridInfoRead->mols[i_us].nLevels<=0){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 14; + } + } + } + + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + return 0; +} + + + diff --git a/src/gridio.h b/src/gridio.h new file mode 100644 index 0000000..c17e726 --- /dev/null +++ b/src/gridio.h @@ -0,0 +1,37 @@ +/* + * gridio.h + * This file is part of LIME, the versatile line modeling engine + * + * Copyright (C) 2006-2014 Christian Brinch + * Copyright (C) 2015-2016 The LIME development team + * + */ + +#ifndef GRIDIO_H +#define GRIDIO_H + +#define lime_FITS 1 + +struct linkType { + unsigned int id, gis[2]; + double *vels; +}; + +struct molInfoType{ + char *molName; + int nLevels, nLines; +}; + +struct gridInfoType{ + unsigned int nInternalPoints, nSinkPoints, nLinks, nNNIndices; + unsigned short nDims, nSpecies, nDensities, nLinkVels; + struct molInfoType *mols; +}; + +#include "grid2fits.h" + +int readGrid(char*, const int, struct gridInfoType*, struct keywordType*, const int, struct grid**, char***, int*, int*); +int writeGrid(char*, const int, struct gridInfoType, struct keywordType*, const int, struct grid*, char**, const int); + +#endif /* GRIDIO_H */ + diff --git a/src/inpars.h b/src/inpars.h index c78d2c5..c93e348 100644 --- a/src/inpars.h +++ b/src/inpars.h @@ -22,6 +22,8 @@ typedef struct { char *dust; int sampling,lte_only,init_lte,antialias,polarization,nThreads; char **moldatfile; + char *gridInFile,**gridOutFiles; + int nSolveIters; double (*gridDensMaxLoc)[DIM], *gridDensMaxValues; } inputPars; diff --git a/src/lime.h b/src/lime.h index 650cfd2..1e928ac 100644 --- a/src/lime.h +++ b/src/lime.h @@ -84,6 +84,7 @@ #define N_RAN_PER_SEGMENT 3 #define FAST_EXP_MAX_TAYLOR 3 #define FAST_EXP_NUM_BITS 8 +#define NUM_GRID_STAGES 4 #define MAX_N_COLL_PART 7 #define N_SMOOTH_ITERS 20 #define TYPICAL_ISM_DENS 1000.0 @@ -94,6 +95,8 @@ #define ERF_TABLE_SIZE 6145 #define BIN_WIDTH (ERF_TABLE_LIMIT/(ERF_TABLE_SIZE-1.)) #define IBIN_WIDTH 1./BIN_WIDTH +#define N_VEL_SEG_PER_HALF 1 +#define NUM_VEL_COEFFS 1+2*N_VEL_SEG_PER_HALF /* This is the number of velocity samples per edge (not including the grid vertices at each end of the edge). Currently this is elsewhere hard-wired at 3, the macro just being used in the file I/O modules. Note that we want an odd number of velocity samples per edge if we want to have the ability to do 2nd-order interpolation of velocity within Delaunay tetrahedra. */ /* Collision partner ID numbers from LAMDA */ #define CP_H2 1 @@ -104,6 +107,36 @@ #define CP_He 6 #define CP_Hplus 7 +/* Bit locations for the grid data-stage mask, that records the information which is present in the grid struct: */ +#define DS_bit_x 0 /* id, x, sink */ +#define DS_bit_neighbours 1 /* neigh, dir, ds, numNeigh */ +#define DS_bit_velocity 2 /* vel */ +#define DS_bit_density 3 /* dens */ +#define DS_bit_abundance 4 /* abun, nmol */ +#define DS_bit_turb_doppler 5 /* dopb */ +#define DS_bit_temperatures 6 /* t */ +#define DS_bit_magfield 7 /* B */ +#define DS_bit_ACOEFF 8 /* a0, a1, a2, a3, a4 */ +#define DS_bit_populations 9 /* mol */ + +#define DS_mask_x 1< #include "lime.h" +#include "gridio.h" /* Forward declaration of functions only used in this file */ @@ -24,7 +25,7 @@ double EXP_TABLE_3D[256][2][10]; calcFastExpRange(FAST_EXP_MAX_TAYLOR, FAST_EXP_NUM_BITS, &numMantissaFields, &lowestExponent, &numExponentsUsed) */ #else -double EXP_TABLE_2D[1][1]; // nominal definitions so the fastexp.c module will compile. +double EXP_TABLE_2D[1][1]; /* nominal definitions so the fastexp.c module will compile. */ double EXP_TABLE_3D[1][1][1]; #endif @@ -57,6 +58,7 @@ initParImg(inputPars *par, image **img) par->gridfile = NULL; par->pregrid = NULL; par->restart = NULL; + par->gridInFile = NULL; par->collPartIds = malloc(sizeof(int)*MAX_N_COLL_PART); for(i=0;icollPartIds[i] = 0; /* Possible values start at 1. */ @@ -81,8 +83,13 @@ initParImg(inputPars *par, image **img) par->antialias=1; par->polarization=0; par->nThreads = NTHREADS; + par->nSolveIters=17; par->traceRayAlgorithm=0; + par->gridOutFiles = malloc(sizeof(char *)*NUM_GRID_STAGES); + for(i=0;igridOutFiles[i] = NULL; + /* Allocate initial space for molecular data file names */ par->moldatfile=malloc(sizeof(char *)*MAX_NSPECIES); for(id=0;id0){ + for(i=0;i0){ molInit(&par, md); - if(!popsdone){ + if(!popsdone && !allBitsSet(par.dataFlags, DS_mask_populations)){ for(gi=0;gi0){ + for(i=0;inSolveIters); } else if (flag == 1){ if (minsnr < 1000) printf(" Statistics: Min(SNR) %3.3f \n", minsnr); @@ -167,7 +167,7 @@ progressbar2(int flag, int prog, double percent, double minsnr, double median){ else printf(" Statistics: Median(SNR) %.3e \n", median); - printf(" Iteration %i / max %i: DONE\n\n", prog, NITERATIONS + 1); + printf(" Iteration %i / max %i: DONE\n\n", prog+1, par->nSolveIters); } #else @@ -348,3 +348,25 @@ collpartmesg3(int number, int flag){ refresh(); #endif } + +void +processFitsError(int status){ + /*****************************************************/ + /* Print out cfitsio error messages and exit program */ + /*****************************************************/ + + if(status){ + if(!silent){ +#ifdef NO_NCURSES + fits_report_error(stderr, status); /* print error report */ +#else + char message[80]; + sprintf(message, "Error in cfitsio: status=%d", status); + bail_out(message); +#endif + } + exit( status ); /* terminate the program, returning error status */ + } + return; +} + diff --git a/src/photon.c b/src/photon.c index 93e0765..3ea6cc0 100644 --- a/src/photon.c +++ b/src/photon.c @@ -48,7 +48,6 @@ double gaussline(const double v, const double oneOnSigma){ #endif } - int getNextEdge(double *inidir, int id, struct grid *gp, const gsl_rng *ran){ int i,iOfLargest,iOfNextLargest,numPositive; double cosAngle,largest,nextLargest,mytest; @@ -103,7 +102,6 @@ int getNextEdge(double *inidir, int id, struct grid *gp, const gsl_rng *ran){ return iOfLargest; } - void calcLineAmpPWLin(struct grid *g, const int id, const int k\ , const int molI, const double deltav, double *inidir, double *vfac_in, double *vfac_out){ @@ -198,7 +196,8 @@ void calcSourceFn(double dTau, const configInfo *par, double *remnantSnu, double } /*....................................................................*/ -void photon(int id, struct grid *gp, molData *md, int iter, const gsl_rng *ran\ +void +photon(int id, struct grid *gp, molData *md, int iter, const gsl_rng *ran\ , configInfo *par, const int nlinetot, struct blendInfo blends\ , gridPointData *mp, double *halfFirstDs){ @@ -273,7 +272,6 @@ void photon(int id, struct grid *gp, molData *md, int iter, const gsl_rng *ran\ calcLineAmpLin(gp,here,neighI,molI,deltav[molI],inidir,&vfac_in[molI],&vfac_out[molI]); mp[molI].vfac[iphot]=vfac_out[molI]; - } /* Contribution of the local cell to emission and absorption is done in getjbar. @@ -359,6 +357,7 @@ void photon(int id, struct grid *gp, molData *md, int iter, const gsl_rng *ran\ here=there; } while(!gp[here].sink); + /* Add cmb contribution. */ iline = 0; @@ -372,7 +371,8 @@ void photon(int id, struct grid *gp, molData *md, int iter, const gsl_rng *ran\ } /*....................................................................*/ -void getjbar(int posn, molData *md, struct grid *gp, const int molI\ +void +getjbar(int posn, molData *md, struct grid *gp, const int molI\ , configInfo *par, struct blendInfo blends, int nextMolWithBlend\ , gridPointData *mp, double *halfFirstDs){ diff --git a/src/popsin.c b/src/popsin.c index 51ec410..b62d33a 100644 --- a/src/popsin.c +++ b/src/popsin.c @@ -106,6 +106,20 @@ popsin(configInfo *par, struct grid **gp, molData **md, int *popsdone){ delaunay(DIM, *gp, (unsigned long)par->ncell, 0, &dc, &numCells); distCalc(par, *gp); getVelocities(par,*gp); + + par->dataFlags |= (1 << DS_bit_x); + par->dataFlags |= (1 << DS_bit_neighbours); + par->dataFlags |= (1 << DS_bit_velocity); + par->dataFlags |= (1 << DS_bit_density); + par->dataFlags |= (1 << DS_bit_abundance); + par->dataFlags |= (1 << DS_bit_turb_doppler); + par->dataFlags |= (1 << DS_bit_temperatures); +/* par->dataFlags |= (1 << DS_bit_magfield); commented out because we are not yet reading it in popsin (and may never do so) */ + par->dataFlags |= (1 << DS_bit_ACOEFF); + par->dataFlags |= (1 << DS_bit_populations); + +//**** should fill in any missing info via the appropriate function calls. + *popsdone=1; free(dc); diff --git a/src/predefgrid.c b/src/predefgrid.c index 4e6e6cb..0812685 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -83,7 +83,22 @@ predefinedGrid(configInfo *par, struct grid *g){ delaunay(DIM, g, (unsigned long)par->ncell, 0, &dc, &numCells); distCalc(par,g); + // getArea(par,g, ran); + // getMass(par,g, ran); getVelocities_pregrid(par,g); + + par->dataFlags |= (1 << DS_bit_x); + par->dataFlags |= (1 << DS_bit_neighbours); + par->dataFlags |= (1 << DS_bit_velocity); + par->dataFlags |= (1 << DS_bit_density); + par->dataFlags |= (1 << DS_bit_abundance); + par->dataFlags |= (1 << DS_bit_turb_doppler); + par->dataFlags |= (1 << DS_bit_temperatures); + par->dataFlags |= (1 << DS_bit_magfield); + par->dataFlags |= (1 << DS_bit_ACOEFF); + +//**** should fill in any missing info via the appropriate function calls. + if(par->gridfile) write_VTK_unstructured_Points(par, g); gsl_rng_free(ran); free(dc); diff --git a/src/raytrace.c b/src/raytrace.c index 81f80a3..8e838d8 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -265,7 +265,7 @@ This version of traceray implements a new algorithm in which the population valu const int stokesIi=0; const int numFaces = DIM+1, nVertPerFace=3; int ichan,stokesId,di,status,lenChainPtrs,entryI,exitI,vi,vvi,ci; - int si, molI, lineI; + int si,molI,lineI; double xp,yp,zp,x[DIM],dir[DIM],projVelRay,vel[DIM]; double xCmpntsRay[nVertPerFace], ds, snu_pol[3], dtau, contJnu, contAlpha; double jnu, alpha, lineRedShift, vThisChan, deltav, vfac, remnantSnu, expDTau; @@ -530,7 +530,6 @@ Note that the argument 'md', and the grid element '.mol', are only accessed for int curlong, totlong; double triangle[3][2],barys[3],local_cmb,cmbFreq,circleSpacing,scale,angle; double *xySquared=NULL; - _Bool isOutsideImage; gsl_error_handler_t *defaultErrorHandler=NULL; size = img[im].distance*img[im].imgres; diff --git a/src/smooth.c b/src/smooth.c index f64687c..f7428bf 100644 --- a/src/smooth.c +++ b/src/smooth.c @@ -23,6 +23,9 @@ smooth(configInfo *par, struct grid *gp){ unsigned long numCells; for(sg=0;sgncell, 0, &dc, &numCells); + distCalc(par, gp); + for(i=0;ipIntensity;i++){ mindist=1e30; cn=-1; @@ -69,8 +72,6 @@ smooth(configInfo *par, struct grid *gp){ } } - delaunay(DIM, gp, (unsigned long)par->ncell, 0, &dc, &numCells); - distCalc(par, gp); if(!silent) progressbar((double)(sg+1)/(double)N_SMOOTH_ITERS, 5); free(dc); } diff --git a/src/stateq.c b/src/stateq.c index 64e1a80..a96afc2 100644 --- a/src/stateq.c +++ b/src/stateq.c @@ -14,8 +14,8 @@ void -stateq(int id, struct grid *gp, molData *md, const int ispec\ - , configInfo *par, struct blendInfo blends, int nextMolWithBlend, gridPointData *mp\ +stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\ + , struct blendInfo blends, int nextMolWithBlend, gridPointData *mp\ , double *halfFirstDs, _Bool *luWarningGiven){ int t,s,iter,status; diff --git a/src/tree_random.c b/src/tree_random.c index 17906d8..c08c274 100644 --- a/src/tree_random.c +++ b/src/tree_random.c @@ -71,8 +71,6 @@ See general remarks about the algorithm at the head of the present module. */ int di, ppi; - double x; - double r[DIM]; rinc->numSubFields = 1; for(di=0;dilastLeafI],axesMid[DIM],scales[DIM],sign,point[DIM]; int i,di,maxInRandI,lastInRandI,piIn,indices[DIM]; - double (*inRandLocations)[DIM]=NULL,r[DIM],density,acceptanceChance,desNumPtsDbl,stdDev,shrinkFrac; + double (*inRandLocations)[DIM]=NULL,r[DIM],density,acceptanceChance,shrinkFrac; unsigned int desiredNumsPoints[tree->lastLeafI]; - unsigned int remainingDesNum,piOutStart,piOut,lastPiOut; + unsigned int piOutStart,piOut,lastPiOut; gsl_qrng *qrSeqGen = NULL; if(rinc->doQuasiRandom) diff --git a/src/velospline.c b/src/velospline.c index 8b5bd47..feea5d7 100644 --- a/src/velospline.c +++ b/src/velospline.c @@ -10,59 +10,59 @@ #include "lime.h" void -getVelocities(configInfo *par, struct grid *g){ +getVelocities(configInfo *par, struct grid *gp){ int i,k,j,l; double vel[3], x[3]; for(i=0;ipIntensity;i++){ + gp[i].v1=malloc(3*gp[i].numNeigh*sizeof(double)); + gp[i].v2=malloc(3*gp[i].numNeigh*sizeof(double)); + gp[i].v3=malloc(3*gp[i].numNeigh*sizeof(double)); - g[i].v1=malloc(3*g[i].numNeigh*sizeof(double)); - g[i].v2=malloc(3*g[i].numNeigh*sizeof(double)); - g[i].v3=malloc(3*g[i].numNeigh*sizeof(double)); - - velocity(g[i].x[0],g[i].x[1],g[i].x[2],g[i].vel); + velocity(gp[i].x[0],gp[i].x[1],gp[i].x[2],gp[i].vel); - for(k=0;kpIntensity;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); + velocity(gp[i].x[0],gp[i].x[1],gp[i].x[2],gp[i].vel); - g[i].v1=malloc(3*g[i].numNeigh*sizeof(double)); - g[i].v2=malloc(3*g[i].numNeigh*sizeof(double)); - g[i].v3=malloc(3*g[i].numNeigh*sizeof(double)); + gp[i].v1=malloc(3*gp[i].numNeigh*sizeof(double)); + gp[i].v2=malloc(3*gp[i].numNeigh*sizeof(double)); + gp[i].v3=malloc(3*gp[i].numNeigh*sizeof(double)); - for(j=0;jpIntensity;incell;i++){ - for(j=0;j<3;j++) g[i].vel[j]=0.; + for(j=0;j<3;j++) gp[i].vel[j]=0.; } } +