From 76e031a368bdd0f14ac952de0e653800b495f8cc Mon Sep 17 00:00:00 2001 From: ims Date: Thu, 18 Feb 2016 14:34:47 +0100 Subject: [PATCH 01/13] Read and write grid to fits A new module grid2fits.c has been added which contains code to do this. There is also a new function processFitsError() in messages.c. The appropriate necessary additions have been made to lime.h and the makefile. Note I also made some cosmetic changes to the Makefile, mostly evening out the tab counts in SRCS and OBJS lines. So far I have not added any machinery by which the user can access the new FITS i/o routines. The functionality of the task is unaltered. I just wanted to add the bulk code first (which has already been tested elsewhere) and make sure it all compiled before proceeding further. --- Makefile | 30 +- src/grid2fits.c | 1399 +++++++++++++++++++++++++++++++++++++++++++++++ src/lime.h | 27 + src/messages.c | 21 + 4 files changed, 1462 insertions(+), 15 deletions(-) create mode 100644 src/grid2fits.c diff --git a/Makefile b/Makefile index 10b4ee0..c5bc052 100644 --- a/Makefile +++ b/Makefile @@ -49,27 +49,27 @@ endif ## TARGET = lime.x -CC = gcc -fopenmp +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/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/weights.c \ - src/velospline.c src/getclosest.c \ - src/tcpsocket.c src/defaults.c src/fastexp.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/grid2fits.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 \ - 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/weights.o \ - src/velospline.o src/getclosest.o \ - src/tcpsocket.o src/defaults.o src/fastexp.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/weights.o \ + src/velospline.o src/getclosest.o src/grid2fits.o \ + src/tcpsocket.o src/defaults.o src/fastexp.o MODELO = src/model.o -#CCFLAGS = -O3 -falign-loops=16 -fno-strict-aliasing -DTEST -DFASTEXP -DNO_NCURSES +#CCFLAGS = -O3 -falign-loops=16 -fno-strict-aliasing -DTEST CCFLAGS = -O3 -falign-loops=16 -fno-strict-aliasing LDFLAGS = -lgsl -lgslcblas -l${QHULL} -lcfitsio -lncurses -lm diff --git a/src/grid2fits.c b/src/grid2fits.c new file mode 100644 index 0000000..6dcb526 --- /dev/null +++ b/src/grid2fits.c @@ -0,0 +1,1399 @@ +/* + * 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. + */ + +#include "lime.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. + +1) Different amounts of information in grid. +-------------------------------------------- +Describing the format of the FITS file is complicated by the fact that the amount of information stored in the grid data structure tends to increase during the run of the LIME program. One can distingish 4 stages of completion which are at least conceptually different, even if they are not cleanly separated in the present implementation of the code. These are summarized as follows. + + A) At this stage the vector of grid objects has been malloc'd and values have been generated for the following struct elements: + id + x + sink + This stage is flagged to the software by the element 'neigh' of the first grid point being set to NULL. + + B) This stage is entered after the Delaunay neighbours of each grid point have been determined. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: + numNeigh + neigh + This stage is flagged by the element 'vel' of the first grid point being set to NULL. + + C) This stage is entered after sampling the user-supplied functions for density, velocity etc. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: + vel + a0, a1 etc. + dens + t + dopb + abun + This stage is flagged by the element 'mol' of the first grid point being set to NULL. + + D) After 1 or more iterations of populating the levels, we enter stage D, in which all the grid struct elements have been malloc'd and given at least preliminary values; specifically now the element + mol + +(Other grid struct values are deducible from the ones given or are used for temporary storage only.) + +The dependency relationship tree for the stages may be diagrammed as follows: + + A -> B -> C -> D. + +2) The FITS file format for stage D +----------------------------------- +This is defined in the following list of extensions. 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. + + 1) GRID + Number of rows = number of grid points. + + Keywords: + RADIUS + COLLPARn # 1 for each nth collision partner. + + Columns: + ID V + Xj D # Cartesian components of the point location, 1 col per jth dimension. + IS_SINK L # =True iff the point lies on the edge of the model. + NUMNEIGH U + FIRST_NN V # See explanation in section 3 below. + VELj D # 1 col per jth dimension. + TURBDPLR E # Given Gaussian lineshape exp(-v^2/[B^2 + 2*k*T/m]), this is B. + DENSITYn D # 1 per nth collision partner. + TEMPKNTC E # From t[0]. + TEMPDUST E # From t[1]. + ABUNMOLm E # 1 per mth molecular species. + + 2) NN_INDICES (see explanation in section 3 below) + Number of rows = number of grid points * average number of Delaunay links per point. + + Columns: + LINK_I V + + 3) LINKS (see explanation in section 3 below) + Number of rows = number of Delaunay links. + + Columns: + GRID_I_1 V + GRID_I_2 V + ACOEFF_p D # 1 per pth order of the velocity polynomial. + + 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: + MOL_NAME + +The subsets of these FITS structures which are written under the other 3 stages are described in the header comment to the function writeGridToFits(). + +3) The NN_INDICES and LINKS extensions and the FIRST_NN column. +--------------------------------------------------------------- +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 FITS 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;j at least stage B. */ + constructLinkArrays((unsigned int)par.ncell, g, &links, &totalNumLinks\ + , &nnLinks, &firstNearNeigh, &totalNumNeigh); + } + + 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); + } + + /* GRID + */ + writeGridExtToFits(fptr, par, numDims, g, firstNearNeigh, collPartNames); /* Deals internally with stages B and C. */ + + if (g[0].neigh!=NULL){ /* => at least stage B. */ + /* NN_INDICES + */ + writeNnIndicesExtToFits(fptr, totalNumNeigh, nnLinks, links); + + /* LINKS + */ + if(g[0].vel!=NULL) stageC = 1; + writeLinksExtToFits(fptr, stageC, totalNumLinks, numACoeffs, links); + + if (g[0].mol!=NULL){ /* => stage D. */ + for(speciesI=0;speciesInumNeigh;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].g[0] = gAPtr; + (*links)[li].g[1] = gBPtr; + + li++; + } + + nnLinkIs[ni] = linkId; + ni++; + } + + pointIsDone[idA] = 1; + } + *totalNumLinks = li; + + *links = realloc(*links, sizeof(struct linkType)*(*totalNumLinks)); + + if(g[0].vel!=NULL){ /* => at least stage C. */ + for(li=0;li<*totalNumLinks;li++){ + gAPtr = (*links)[li].g[0]; + /* Find which neighbour of gAPtr corresponds to the link: */ + linkNotFound = 1; + for(jA=0;jAnumNeigh;jA++){ + idB = gAPtr->neigh[jA]->id; + if(linkNotFound && idB==(*links)[li].g[1]->id){ + (*links)[li].aCoeffs[0] = gAPtr->a0[jA]; + (*links)[li].aCoeffs[1] = gAPtr->a1[jA]; + (*links)[li].aCoeffs[2] = gAPtr->a2[jA]; + (*links)[li].aCoeffs[3] = gAPtr->a3[jA]; + (*links)[li].aCoeffs[4] = gAPtr->a4[jA]; + linkNotFound = 0; + } + } + if(linkNotFound){ + sprintf(message, "SNAFU with link %d grid point %d", li, gAPtr->id); + if(!silent) bail_out(message); + exit(1); + } + } + } + + *nnLinks = malloc(sizeof(**nnLinks)*(*totalNumNeigh)); + for(ni=0;ni<(*totalNumNeigh);ni++) + (*nnLinks)[ni] = &(*links)[nnLinkIs[ni]]; + + free(nnLinkIs); + free(pointIsDone); + /* The calling routine must free firstNearNeigh, nnLinks, links. */ +} + +/*....................................................................*/ +void defineGridExtColumns(unsigned short numKwdChars, unsigned short numDims\ + , unsigned short numCollPart, unsigned short numSpecies, char *ttype[]\ + , char *tform[], char *tunit[], int dataTypes[]){ + + /* Note that data types in all capitals are defined in fitsio.h. */ + + const int maxNumDims=9, maxNumSpecies=9, maxNumCollPart=9; + int colI=0, i; + char message[80]; + + if(numDims>maxNumDims){ + if(!silent){ + sprintf(message, "Caller asked for %d dims but colnames can only be written for %d.", numDims, maxNumDims); + bail_out(message); + } + exit(1); + } + + if(numSpecies>maxNumSpecies){ + if(!silent){ + sprintf(message, "Caller asked for %d species but colnames can only be written for %d.", numSpecies, maxNumSpecies); + bail_out(message); + } + exit(1); + } + + if(numCollPart>maxNumCollPart){ + if(!silent){ + sprintf(message, "Caller asked for %d coll. part. but colnames can only be written for %d.", numCollPart, maxNumCollPart); + bail_out(message); + } + exit(1); + } + + colI = 0; + ttype[colI] = malloc(numKwdChars); + sprintf(ttype[colI], "ID"); + tform[colI] = "V"; + tunit[colI] = "\0"; + dataTypes[colI] = TUINT; + + /* should rather have a vector column? */ + for(i=0;i at least stage B. */ + colI++; + numNeigh = malloc(sizeof(*numNeigh)*par.ncell); + for(i=0;i at least stage C. */ + velj = malloc(sizeof(*velj)*par.ncell); + for(j=0;jmaxNumCollPart){ + if(!silent){ + sprintf(message, "There seem to be %d collision partners but keywords can only be written for %d.", par.collPart, maxNumCollPart); + warning(message); + } + localNumCollPart = maxNumCollPart; + }else{ + localNumCollPart = par.collPart; + } + + for(i=0;iid; + fits_write_col(fptr, dataType, 1, firstRow, firstElem, (LONGLONG)totalNumNeigh, linkIs, &status); + processFitsError(status); + free(linkIs); +} + +/*....................................................................*/ +void writeLinksExtToFits(fitsfile *fptr, _Bool stageC, unsigned int totalNumLinks\ + , unsigned short numACoeffs, struct linkType *links){ + /* +See the comment at the beginning of the present module for a description of how the LINKS extension relates to the grid struct. + +Note that data types in all capitals are defined in fitsio.h. + */ + + unsigned int *ids=NULL; + double *aCoeffs=NULL; + int status=0, colI=0, i, n; + LONGLONG firstRow=1, firstElem=1; + int numCols = 2 + numACoeffs; + char extname[] = "LINKS"; + const int numKwdChars = 9; /* 8 characters + \0. */ + + /* Define the name, datatype, and physical units for the columns. + */ + char *ttype[numCols]; + char *tform[numCols]; + char *tunit[numCols]; + int dataTypes[numCols]; + + colI = 0; + ttype[colI] = malloc(numKwdChars); + sprintf(ttype[colI], "GRID_I_1"); + tform[colI] = "V"; + tunit[colI] = "\0"; + dataTypes[colI] = TUINT; + + colI++; + ttype[colI] = malloc(numKwdChars); + sprintf(ttype[colI], "GRID_I_2"); + tform[colI] = "V"; + tunit[colI] = "\0"; + dataTypes[colI] = TUINT; + + /* should rather have a vector column? */ + for(n=0;nid; + fits_write_col(fptr, dataTypes[colI], colI+1, firstRow, firstElem, (LONGLONG)totalNumLinks, ids, &status); + processFitsError(status); + + colI++; + for(i=0;iid; + fits_write_col(fptr, dataTypes[colI], colI+1, firstRow, firstElem, (LONGLONG)totalNumLinks, ids, &status); + processFitsError(status); + free(ids); + + if (stageC){ + aCoeffs = malloc(sizeof(*aCoeffs)*totalNumLinks); + for(n=0;n: + */ + fits_movnam_hdu(fptr, IMAGE_HDU, extname, 0, &status); + processFitsError(status); + + readPopsExtFromFits(fptr, (unsigned int)par.ncell, md, si, g); + } + } + } + + if(*popsValuesFound==0){ + for(i=0;imaxNumCollPart){ + if(!silent){ + sprintf(message, "There seem to be %d collision partners but keywords can only be written for %d.", par.collPart, maxNumCollPart); + warning(message); + } + *numCollPartRead = maxNumCollPart; + }else{ + *numCollPartRead = par.collPart; + } + + *collPartNames = malloc(sizeof(**collPartNames)*(*numCollPartRead)); + + for(i=0;i<(*numCollPartRead);i++){ + sprintf(genericKwd, "COLLPAR%d", i+1); + (*collPartNames)[i] = malloc(sizeof(char)*100); + fits_read_key(fptr, TSTRING, genericKwd, (*collPartNames)[i], NULL, &status); + processFitsError(status); + } + } +} + +/*....................................................................*/ +void readLinksExtFromFits(fitsfile *fptr, _Bool userValuesFound\ + , unsigned short numACoeffs, struct grid *g, struct linkType **links){ + /* +See the comment at the beginning of the present module for a description of how the LINKS extension relates to the grid struct. + +The present function mallocs *links. + */ + + LONGLONG totalNumLinks, firstRow=1, firstElem=1; + int status=0, colNum, anynul=0, li, n; + char colName[21]; + unsigned int *ids=NULL; + double *aCoeffs=NULL; + + /* Find out how many rows there are, then malloc the array. + */ + fits_get_num_rowsll(fptr, &totalNumLinks, &status); + processFitsError(status); + + *links = malloc(sizeof(**links)*totalNumLinks); + + /* Read GRID_I_1 column. + */ + fits_get_colnum(fptr, CASEINSEN, "GRID_I_1", &colNum, &status); + processFitsError(status); + + ids = malloc(sizeof(*ids)*totalNumLinks); + fits_read_col(fptr, TUINT, colNum, firstRow, firstElem, totalNumLinks, 0, ids, &anynul, &status); + processFitsError(status); + + for(li=0;li0){ /* Means we are in at least stage C. */ + aCoeffs = malloc(sizeof(*aCoeffs)*totalNumLinks); + for(n=0;ng[0]->id==i){ + (*g)[i].neigh[j] = linkPtr->g[1]; + }else{ + (*g)[i].neigh[j] = linkPtr->g[0]; + } + + (*g)[i].a0[j] = linkPtr->aCoeffs[0]; + (*g)[i].a1[j] = linkPtr->aCoeffs[1]; + (*g)[i].a2[j] = linkPtr->aCoeffs[2]; + (*g)[i].a3[j] = linkPtr->aCoeffs[3]; + (*g)[i].a4[j] = linkPtr->aCoeffs[4]; + } + } +} + +/*....................................................................*/ +void readPopsExtFromFits(fitsfile *fptr, unsigned int numGridPoints\ + , molData *md, unsigned short speciesI, struct grid **g){ + /* +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 (*g)[yi].mol[speciesI].pops for each grid point yi and species. + */ + + const int maxLenMolName = 8; + int bitpix, naxis, status=0, xi, yi, anynul=0; + long *naxes=NULL; + float *row=NULL; + long fpixels[2],lpixels[2]; + long inc[2] = {1,1}; + char molNameRead[maxLenMolName+1]; + unsigned int numEnergyLevels = md[speciesI].nlev; + + fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status); + processFitsError(status); + /****** error if naxes != {numEnergyLevels, numGridPoints} or bitpix!=FLOAT_IMG.*/ + free(naxes); + + row = malloc(sizeof(*row)*numEnergyLevels); + + /* Read FITS data. + */ + for(yi=0;yi Date: Fri, 19 Feb 2016 12:16:11 +0100 Subject: [PATCH 02/13] Added some write points Here I have added some places where the grid is written to file. Everything seems to work ok, and the output values look sensible, but I have not yet done any exhaustive tests of the values in the fits files. Some new elements of the struct inputPars have been necessary. The default values of these cause no fits output to be written. In fact there should be no changes in the way the software functions with the standard model.c file. I also added a new function printMessage() to messages.c and some lines in main.c which print the number of threads used when this is more than 1. The fits-read functionality remains unused. Adding the ability to read the grid from file should perhaps wait on discussion among the developers about the future of the pregrid etc machinery. --- Makefile | 1 + src/aux.c | 7 ++++++- src/grid.c | 17 +++++++++++++++-- src/grid2fits.c | 26 ++++++++++++++++---------- src/lime.h | 6 +++++- src/main.c | 16 ++++++++++++++-- src/messages.c | 25 ++++++++++++++++++++++--- 7 files changed, 79 insertions(+), 19 deletions(-) diff --git a/Makefile b/Makefile index c5bc052..699d10e 100644 --- a/Makefile +++ b/Makefile @@ -50,6 +50,7 @@ endif TARGET = lime.x CC = gcc -fopenmp +#CC = gcc -fopenmp -g 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 \ diff --git a/src/aux.c b/src/aux.c index f3d9d39..238e765 100644 --- a/src/aux.c +++ b/src/aux.c @@ -41,9 +41,14 @@ parseInput(inputPars *par, image **img, molData **m){ par->doPregrid=0; par->nThreads=0; + for(i=0;iwriteGridAtStage[i] = 0; + par->gridFitsOutSets[i] = ""; + }; + /* Allocate space for output fits images */ (*img)=malloc(sizeof(image)*MAX_NSPECIES); - par->moldatfile=malloc(sizeof(char *) * MAX_NSPECIES); + par->moldatfile=malloc(sizeof(char *)*MAX_NSPECIES); for(id=0;idmoldatfile[id]=NULL; diff --git a/src/grid.c b/src/grid.c index 872f301..5b58afb 100644 --- a/src/grid.c +++ b/src/grid.c @@ -515,13 +515,14 @@ getMass(inputPars *par, struct grid *g, const gsl_rng *ran){ void -buildGrid(inputPars *par, struct grid *g){ +buildGrid(inputPars *par, struct grid *g, molData *md){ double lograd; /* The logarithm of the model radius */ double logmin; /* Logarithm of par->minScale */ double r,theta,phi,sinPhi,x,y,z,semiradius; /* Coordinates */ double temp; int k=0,i; /* counters */ - int flag; + int flag,stageI,status=0; + char **collPartNames=NULL; /*** this is a placeholder until we start reading these. */ gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2); /* Random number generator */ #ifdef TEST @@ -617,6 +618,13 @@ buildGrid(inputPars *par, struct grid *g){ distCalc(par, g); smooth(par,g); +/* Can't do this yet because .dens is not NULL, but .a0, .a1 etc are. Have to rearrange the mallocs of g elements first. Simply set g[0].dens=NULL to kluge around this? + stageI = 1; + if(par->writeGridAtStage[stageI]) + status = writeGridToFits(par->gridFitsOutSets[stageI], *par, (unsigned short)DIM\ + , (unsigned short)NUM_VEL_COEFFS, g, md, collPartNames); +*/ + 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); @@ -631,6 +639,11 @@ buildGrid(inputPars *par, struct grid *g){ gsl_rng_free(ran); if(!silent) done(5); + + stageI = 2; + if(par->writeGridAtStage[stageI]) + status = writeGridToFits(par->gridFitsOutSets[stageI], *par, (unsigned short)DIM\ + , (unsigned short)NUM_VEL_COEFFS, g, md, collPartNames); } diff --git a/src/grid2fits.c b/src/grid2fits.c index 6dcb526..fd2bc27 100644 --- a/src/grid2fits.c +++ b/src/grid2fits.c @@ -30,7 +30,7 @@ Describing the format of the FITS file is complicated by the fact that the amoun B) This stage is entered after the Delaunay neighbours of each grid point have been determined. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: numNeigh neigh - This stage is flagged by the element 'vel' of the first grid point being set to NULL. + This stage is flagged by the element 'dens' of the first grid point being set to NULL. C) This stage is entered after sampling the user-supplied functions for density, velocity etc. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: vel @@ -154,7 +154,7 @@ The pseudo-code is as follows: # GRID_I_0 # GRID_I_1 - if (g[0].vel!=NULL){ # = at least stage C, write: + if (g[0].dens!=NULL){ # = at least stage C, write: # GRID # COLLPARn # -------- @@ -182,16 +182,24 @@ The pseudo-code is as follows: fitsfile *fptr; int status = 0; unsigned short speciesI; - char negfile[100]="! "; + char negfile[100]="! ", message[80]; struct linkType *links=NULL, **nnLinks=NULL; unsigned int totalNumLinks, totalNumNeigh, *firstNearNeigh=NULL; _Bool stageC=0; - if (g==NULL){ - warning("Cannot write grid list to file, there are no entries in it."); + if (outFileName==""){ + if(!silent) warning("Cannot write grid list to file, filename is blank."); return 1; } + if (g==NULL){ + if(!silent) warning("Cannot write grid list to file, there are no entries in it."); + return 2; + } + + sprintf(message, "Writing grid list to file %s", outFileName); + if(!silent) printMessage(message); + /* If we get to here, we have at least stage A info in the grid. */ if (g[0].neigh!=NULL){ /* => at least stage B. */ @@ -219,7 +227,7 @@ The pseudo-code is as follows: /* LINKS */ - if(g[0].vel!=NULL) stageC = 1; + if(g[0].dens!=NULL) stageC = 1; writeLinksExtToFits(fptr, stageC, totalNumLinks, numACoeffs, links); if (g[0].mol!=NULL){ /* => stage D. */ @@ -321,7 +329,7 @@ See the comment at the beginning of the present module for a description of how *links = realloc(*links, sizeof(struct linkType)*(*totalNumLinks)); - if(g[0].vel!=NULL){ /* => at least stage C. */ + if(g[0].dens!=NULL){ /* => at least stage C. */ for(li=0;li<*totalNumLinks;li++){ gAPtr = (*links)[li].g[0]; /* Find which neighbour of gAPtr corresponds to the link: */ @@ -497,7 +505,6 @@ void writeGridExtToFits(fitsfile *fptr, inputPars par, unsigned short numDims\ LONGLONG firstRow=1, firstElem=1; int numCols = 7 + numDims*2 + par.collPart + par.nSpecies; char extname[] = "GRID"; -// char *genericComment; char genericComment[80]; char genericKwd[numKwdChars], message[80]; @@ -558,7 +565,7 @@ void writeGridExtToFits(fitsfile *fptr, inputPars par, unsigned short numDims\ fits_write_col(fptr, dataTypes[colI], colI+1, firstRow, firstElem, (LONGLONG)par.ncell, firstNearNeigh, &status); processFitsError(status); - if (g[0].vel!=NULL){ /* => at least stage C. */ + if (g[0].dens!=NULL){ /* => at least stage C. */ velj = malloc(sizeof(*velj)*par.ncell); for(j=0;j1){ + sprintf(message, "Number of threads used: %d", par.nThreads); + printMessage(message); + } + if(par.doPregrid) { gridAlloc(&par,&g); @@ -50,7 +57,7 @@ int main () { else { gridAlloc(&par,&g); - buildGrid(&par,g); + buildGrid(&par,g,m); } for(i=0;i0) + { + printf("%s\n", message ); + } +#else + move(22,0); printw("*** %s\n",message); + refresh(); +#endif +} + void warning(char message[80]){ #ifdef NO_NCURSES @@ -341,12 +358,14 @@ processFitsError(int status){ char message[80]; if (status){ + if(!silent){ #ifdef NO_NCURSES - fits_report_error(stderr, status); /* print error report */ + fits_report_error(stderr, status); /* print error report */ #else - sprintf(message, "Error in cfitsio: status=%d", status); - bail_out(message); + sprintf(message, "Error in cfitsio: status=%d", status); + bail_out(message); #endif + } exit( status ); /* terminate the program, returning error status */ } return; From a0bf2a3785867ec879b460bbc21dfb5ed0a527d2 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 19 Feb 2016 15:32:45 +0100 Subject: [PATCH 03/13] Name of molecule now written to fits file. --- src/grid.c | 2 +- src/lime.h | 2 +- src/messages.c | 2 +- src/molinit.c | 8 +++++--- 4 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/grid.c b/src/grid.c index 5b58afb..3e16a42 100644 --- a/src/grid.c +++ b/src/grid.c @@ -618,7 +618,7 @@ buildGrid(inputPars *par, struct grid *g, molData *md){ distCalc(par, g); smooth(par,g); -/* Can't do this yet because .dens is not NULL, but .a0, .a1 etc are. Have to rearrange the mallocs of g elements first. Simply set g[0].dens=NULL to kluge around this? +/* Can't do this yet because .dens is not NULL, but .a0, .a1 etc are. Have to rearrange the mallocs of g elements first. Simply set g[0].dens=NULL to kluge around this? But then have to malloc it again (sigh). stageI = 1; if(par->writeGridAtStage[stageI]) status = writeGridToFits(par->gridFitsOutSets[stageI], *par, (unsigned short)DIM\ diff --git a/src/lime.h b/src/lime.h index e4a6be3..60dd462 100644 --- a/src/lime.h +++ b/src/lime.h @@ -101,7 +101,7 @@ typedef struct { int *lal,*lau,*lcl,*lcu; double *aeinst,*freq,*beinstu,*beinstl,*up,*down,*eterm,*gstat; double norm,norminv,*cmb,*local_cmb; - char *molName; + char molName[80]; } molData; /* Data concerning a single grid vertex which is passed from photon() to stateq(). This data needs to be thread-safe. */ diff --git a/src/messages.c b/src/messages.c index fa231a2..3f540d0 100644 --- a/src/messages.c +++ b/src/messages.c @@ -294,7 +294,7 @@ bail_out(char message[80]){ } void -collpartmesg(char molecule[90], int partners){//, int specnumber){ +collpartmesg(char molecule[80], int partners){//, int specnumber){ #ifdef NO_NCURSES printf(" Molecule: %.25s\n", molecule); if (partners==1) diff --git a/src/molinit.c b/src/molinit.c index 4ce8ea2..6915c1d 100644 --- a/src/molinit.c +++ b/src/molinit.c @@ -108,7 +108,7 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){ double fac, uprate, downrate=0, dummy, amass; struct data { double *colld, *temp; } *part; - char string[200], specref[90], partstr[90]; + char string[200], partstr[90]; FILE *fp; if((fp=fopen(par->moldatfile[i], "r"))==NULL) { @@ -118,13 +118,15 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){ /* Read the header of the data file */ fgets(string, 80, fp); - fgets(specref, 90, fp); + fgets(m[i].molName, 80, fp); fgets(string, 80, fp); fscanf(fp, "%lf\n", &amass); fgets(string, 80, fp); fscanf(fp, "%d\n", &m[i].nlev); fgets(string, 80, fp); + m[i].molName[strlen(m[i].molName)-1] = 0; /* Lops the carriage return from the end. */ + m[i].eterm=malloc(sizeof(double)*m[i].nlev); m[i].gstat=malloc(sizeof(double)*m[i].nlev); @@ -232,7 +234,7 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){ strcat( partstr, collpartname[count[ipart]-1]); } if(!silent) { - collpartmesg(specref, m[i].npart); + collpartmesg(m[i].molName, m[i].npart); collpartmesg2(partstr, ipart); collpartmesg3(par->collPart, flag); } From d0a50a9ad7733d6622802e6451bf846be183a799 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 17 Jun 2016 14:41:35 +0200 Subject: [PATCH 04/13] Many small changes This commit covers numerous minor issues to get them out of the way before the big one (or ones). A summary is as follows: - Includes fix for issue #131: where nmol was not being set for sink points in LTE. - C++->C style comments in aux.c:levelPops(). - conv renamed nItersDone in aux.c:levelPops(). - Replaced 'if(fileName)' with 'if(fileName!=NULL)' in a few places in aux.c:levelPops(). - Replaced 'g' by 'gp' and 'm' by 'md' in numerous places. - Number of density values returned is now counted (and used to set par->collPart) in aux.c:parseInput() rather than grid.c:gridAlloc(). - Added warning in aux.c:parseInput() if could not find moldata file. - Grid velocity now set at the same time as dens, abun etc in grid.c:readOrBuildGrid() rather than in getVelosplines(). - New header file gridio.h with a couple of macros relating to file types. - New structs molInfoType, gridInfoType to contain sizes of various fields in grid info which has been read from file. (Ultimately these should be merged into a single struct which contains all the grid data.) - messages.c:done() renamed to printDone(). - photon.c:photon(): unused variable x[] removed. - In smooth.c:smooth(): qhull() and distCalc() were called at the end of each iteration; now they are called at the beginning of each. In grid.c:buildGrid() the calls to qhull() and distCalc() which preceded the call to smooth() now succeed it. The functionality is unchanged but it allows one to read a grid file at stage 1 and have the triangulation performed by qhull() without perturbation of the grid point locations. --- src/LTEsolution.c | 6 ++-- src/aux.c | 53 +++++++++++++++++++------------ src/grid.c | 81 ++++++++++++++++++++++------------------------- src/gridio.h | 12 +++++++ src/lime.h | 27 +++++++++++----- src/messages.c | 2 +- src/photon.c | 4 +-- src/popsout.c | 12 +++---- src/predefgrid.c | 3 +- src/raytrace.c | 66 +++++++++++++++++++------------------- src/smooth.c | 39 ++++++++++++----------- src/sourcefunc.c | 31 +++++++++--------- src/velospline.c | 3 -- src/writefits.c | 2 +- 14 files changed, 184 insertions(+), 157 deletions(-) create mode 100644 src/gridio.h diff --git a/src/LTEsolution.c b/src/LTEsolution.c index 874213d..8ae879b 100644 --- a/src/LTEsolution.c +++ b/src/LTEsolution.c @@ -15,8 +15,10 @@ LTE(inputPars *par, struct grid *g, molData *m){ double z; for(ispec=0;ispecnSpecies;ispec++){ - for(id=0;idpIntensity;id++){ + for(id=0;idncell;id++){ g[id].nmol[ispec]=g[id].abun[ispec]*g[id].dens[0]; + } + for(id=0;idpIntensity;id++){ z=0; for(ilev=0;ilevoutputfile) popsout(par,g,m); + if(par->outputfile != NULL) popsout(par,g,m); } diff --git a/src/aux.c b/src/aux.c index b47306d..f72cfe1 100644 --- a/src/aux.c +++ b/src/aux.c @@ -19,6 +19,7 @@ parseInput(inputPars *par, image **img, molData **m){ int i,id; double BB[3]; double cosPhi,sinPhi,cosTheta,sinTheta; + double temp[99]; /* Set default values */ par->dust = NULL; @@ -79,6 +80,7 @@ parseInput(inputPars *par, image **img, molData **m){ /* Check if files exists */ for(id=0;idnSpecies;id++){ if((fp=fopen(par->moldatfile[id], "r"))==NULL) { + if(!silent) warning("Could not find moldata file, attempting to obtain it over the web."); openSocket(par, id); } else { @@ -109,6 +111,15 @@ parseInput(inputPars *par, image **img, molData **m){ par->minScaleSqu=par->minScale*par->minScale; if(par->pregrid!=NULL) par->doPregrid=1; + if(par->doPregrid || par->restart) par->collPart=1; + else{ + for(i=0;i<99;i++) temp[i]=-1; + density(AU,AU,AU,temp); + i=0; + par->collPart=0; + while(temp[i++]>-1) par->collPart++; + } + /* Now we need to calculate the cutoff value used in calcSourceFn(). The issue is to decide between @@ -350,22 +361,22 @@ invSqrt(float x){ } void -continuumSetup(int im, image *img, molData *m, inputPars *par, struct grid *g){ +continuumSetup(int im, image *img, molData *m, inputPars *par, struct grid *gp){ int id; img[im].trans=0; m[0].nline=1; m[0].freq= malloc(sizeof(double)); m[0].freq[0]=img[im].freq; for(id=0;idncell;id++) { - freePopulation( par, m, g[id].mol ); - g[id].mol=malloc(sizeof(struct populations)*1); - g[id].mol[0].dust = malloc(sizeof(double)*m[0].nline); - g[id].mol[0].knu = malloc(sizeof(double)*m[0].nline); - g[id].mol[0].pops = NULL; - g[id].mol[0].partner = NULL; + freePopulation( par, m, gp[id].mol ); + gp[id].mol=malloc(sizeof(struct populations)*1); + gp[id].mol[0].dust = malloc(sizeof(double)*m[0].nline); + gp[id].mol[0].knu = malloc(sizeof(double)*m[0].nline); + gp[id].mol[0].pops = NULL; + gp[id].mol[0].partner = NULL; } - if(par->outputfile) popsout(par,g,m); - kappa(m,g,par,0); + if(par->outputfile) popsout(par,gp,m); + kappa(m,gp,par,0); } void @@ -432,7 +443,7 @@ lineBlend(molData *m, inputPars *par, blend **matrix){ void levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ - int id,conv=0,iter,ilev,prog=0,ispec,c=0,n,i,threadI,nVerticesDone; + int id,iter,ilev,prog=0,ispec,c=0,n,i,threadI,nVerticesDone,nItersDone; double percent=0.,*median,result1=0,result2=0,snr,delta_pop; blend *matrix; struct statistics { double *pop, *ave, *sigma; } *stat; @@ -486,8 +497,7 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ } } - if(par->outputfile) popsout(par,g,m); - + if(par->outputfile != NULL) popsout(par,g,m); /* Initialize convergence flag */ for(id=0;idncell;id++){ @@ -495,6 +505,7 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ } if(par->lte_only==0){ + nItersDone=0; do{ if(!silent) progressbar2(0, prog++, 0, result1, result2); @@ -512,8 +523,8 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ threadI = omp_get_thread_num(); /* Declare and allocate thread-private variables */ - gridPointData *mp; // Could have declared them earlier - double *halfFirstDs; // and included them in private() I guess. + gridPointData *mp; /* Could have declared them earlier */ + double *halfFirstDs; /* and included them in private() I guess. */ mp=malloc(sizeof(gridPointData)*par->nSpecies); for (i=0;inSpecies;i++){ mp[i].phot = malloc(sizeof(double)*m[i].nline*max_phot); @@ -527,21 +538,21 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ #pragma omp atomic ++nVerticesDone; - if (threadI == 0){ // i.e., is master thread + if (threadI == 0){ /* i.e., is master thread. */ if(!silent) progressbar((double)nVerticesDone/par->pIntensity,10); } if(g[id].dens[0] > 0 && g[id].t[0] > 0){ photon(id,g,m,0,threadRans[threadI],par,matrix,mp,halfFirstDs); for(ispec=0;ispecnSpecies;ispec++) stateq(id,g,m,ispec,par,mp,halfFirstDs); } - if (threadI == 0){ // i.e., is master thread + if (threadI == 0){ /* i.e., is master thread. */ if(!silent) warning(""); } } freeGridPointData(par, mp); free(halfFirstDs); - } // end parallel block. + } /* end parallel block. */ for(id=0;idncell && !g[id].sink;id++){ snr=0; @@ -578,16 +589,16 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ } gsl_sort(median, 1, c); - if(conv>1){ + 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,g,m); - } while(conv++binoutputfile) binpopsout(par,g,m); + if(par->outputfile != NULL) popsout(par,g,m); + } while(nItersDone++binoutputfile != NULL) binpopsout(par,g,m); } for (i=0;inThreads;i++){ diff --git a/src/grid.c b/src/grid.c index 3e16a42..57d16a6 100644 --- a/src/grid.c +++ b/src/grid.c @@ -13,20 +13,10 @@ void gridAlloc(inputPars *par, struct grid **g){ int i; - double temp[99]; *g=malloc(sizeof(struct grid)*(par->pIntensity+par->sinkPoints)); memset(*g, 0., sizeof(struct grid) * (par->pIntensity+par->sinkPoints)); - if(par->doPregrid || par->restart) par->collPart=1; - else{ - for(i=0;i<99;i++) temp[i]=-1; - density(AU,AU,AU,temp); - i=0; - par->collPart=0; - while(temp[i++]>-1) par->collPart++; - } - for(i=0;i<(par->pIntensity+par->sinkPoints); i++){ (*g)[i].a0 = NULL; (*g)[i].a1 = NULL; @@ -151,7 +141,7 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ } void -qhull(inputPars *par, struct grid *g){ +qhull(inputPars *par, struct grid *gp){ int i,j,k,id; char flags[255]; boolT ismalloc = False; @@ -165,7 +155,7 @@ qhull(inputPars *par, struct grid *g){ for(i=0;incell;i++) { for(j=0;jpoint); - 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[simplex[j]].id) { + while(gp[simplex[i]].neigh[k] != NULL && gp[simplex[i]].neigh[k]->id != gp[simplex[j]].id) { k++; } - g[simplex[i]].neigh[k]=&g[simplex[j]]; + gp[simplex[i]].neigh[k]=&gp[simplex[j]]; } } } @@ -210,13 +200,13 @@ qhull(inputPars *par, struct grid *g){ for(i=0;incell;i++){ j=0; - for(k=0;kncell;i++){ - if( g[i].dir != NULL ) + if( gp[i].dir != NULL ) { - free( g[i].dir ); + free( gp[i].dir ); } - if( g[i].ds != NULL ) + if( gp[i].ds != NULL ) { - free( g[i].ds ); + free( gp[i].ds ); } - g[i].dir=malloc(sizeof(point)*g[i].numNeigh); - g[i].ds =malloc(sizeof(double)*g[i].numNeigh); - memset(g[i].dir, 0., sizeof(point) * g[i].numNeigh); - memset(g[i].ds, 0., sizeof(double) * g[i].numNeigh); - for(k=0;kx[l] - g[i].x[l]; - g[i].ds[k]=sqrt(g[i].dir[k].x[0]*g[i].dir[k].x[0]+g[i].dir[k].x[1]*g[i].dir[k].x[1]+g[i].dir[k].x[2]*g[i].dir[k].x[2]); - for(l=0;l<3;l++) g[i].dir[k].xn[l] = g[i].dir[k].x[l]/g[i].ds[k]; + gp[i].dir=malloc(sizeof(point)*gp[i].numNeigh); + gp[i].ds =malloc(sizeof(double)*gp[i].numNeigh); + memset(gp[i].dir, 0., sizeof(point) * gp[i].numNeigh); + memset(gp[i].ds, 0., sizeof(double) * gp[i].numNeigh); + for(k=0;kx[l] - gp[i].x[l]; + gp[i].ds[k]=sqrt(gp[i].dir[k].x[0]*gp[i].dir[k].x[0]+gp[i].dir[k].x[1]*gp[i].dir[k].x[1]+gp[i].dir[k].x[2]*gp[i].dir[k].x[2]); + for(l=0;l<3;l++) gp[i].dir[k].xn[l] = gp[i].dir[k].x[l]/gp[i].ds[k]; } - g[i].nphot=ininphot*g[i].numNeigh; + gp[i].nphot=ininphot*gp[i].numNeigh; } } @@ -520,7 +510,7 @@ buildGrid(inputPars *par, struct grid *g, molData *md){ double logmin; /* Logarithm of par->minScale */ double r,theta,phi,sinPhi,x,y,z,semiradius; /* Coordinates */ double temp; - int k=0,i; /* counters */ + int k=0,i,j; /* counters */ int flag,stageI,status=0; char **collPartNames=NULL; /*** this is a placeholder until we start reading these. */ @@ -591,7 +581,7 @@ buildGrid(inputPars *par, struct grid *g, molData *md){ if(!silent) progressbar((double) k/((double)par->pIntensity-1), 4); } /* end model grid point assignment */ - if(!silent) done(4); + if(!silent) printDone(4); /* Add surface sink particles */ for(i=0;isinkPoints;i++){ @@ -610,13 +600,15 @@ buildGrid(inputPars *par, struct grid *g, molData *md){ g[k].dens[0]=1e-30; g[k].t[0]=par->tcmb; g[k].t[1]=par->tcmb; - g[k++].dopb=0.; + g[k].dopb=0.; + for(j=0;jwriteGridAtStage[stageI]) diff --git a/src/gridio.h b/src/gridio.h new file mode 100644 index 0000000..ac5c038 --- /dev/null +++ b/src/gridio.h @@ -0,0 +1,12 @@ +/* + * gridio.h + * This file is part of LIME, the versatile line modeling engine + * + * Copyright (C) 2006-2014 Christian Brinch + * Copyright (C) 2015 The LIME development team + * + */ + +#define lime_fptr fitsfile +#define lime_FITS 1 + diff --git a/src/lime.h b/src/lime.h index 0879e30..ed96fd3 100644 --- a/src/lime.h +++ b/src/lime.h @@ -17,6 +17,7 @@ #include #include #include +#include "gridio.h" #ifdef OLD_QHULL #include @@ -115,17 +116,16 @@ typedef struct { /* Point coordinate */ typedef struct { - double x[3]; - double xn[3]; + double x[DIM]; + double xn[DIM]; } point; struct rates { double *up, *down; }; - struct populations { - double * pops, *knu, *dust; + double *pops, *knu, *dust; double dopb, binv; struct rates *partner; }; @@ -133,8 +133,8 @@ struct populations { /* Grid properties */ struct grid { int id; - double x[3]; - double vel[3]; + double x[DIM]; + double vel[DIM]; double *a0,*a1,*a2,*a3,*a4; int numNeigh; point *dir; @@ -148,10 +148,21 @@ struct grid { struct populations* mol; }; +struct molInfoType{ + char *molName; + int nLevels, nLines; +}; + +struct gridInfoType{ + unsigned int nInternalPoints, nSinkPoints, nLinks, nNNIndices; + unsigned short nDims, nSpecies, nDensities, nACoeffs; + struct molInfoType *mols; +}; + struct linkType { unsigned int id; struct grid *g[2]; - double aCoeffs[NUM_VEL_COEFFS]; + double *aCoeffs; }; typedef struct { @@ -284,7 +295,7 @@ double FastExp(const float negarg); void greetings(); void greetings_parallel(int); void screenInfo(); -void done(int); +void printDone(int); void progressbar(double,int); void progressbar2(int,int,double,double,double); void casaStyleProgressBar(const int,int); diff --git a/src/messages.c b/src/messages.c index b1cc5f9..6dbe387 100644 --- a/src/messages.c +++ b/src/messages.c @@ -106,7 +106,7 @@ screenInfo(){ } void -done(int line){ +printDone(int line){ #ifdef NO_NCURSES if (line == 4) printf( " Building grid: DONE \n\n"); diff --git a/src/photon.c b/src/photon.c index e0ef063..81930bf 100644 --- a/src/photon.c +++ b/src/photon.c @@ -169,7 +169,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); @@ -220,10 +220,8 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar else velocityspline_lin(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/popsout.c b/src/popsout.c index 93ff3e3..615d11d 100644 --- a/src/popsout.c +++ b/src/popsout.c @@ -10,7 +10,7 @@ #include "lime.h" void -popsout(inputPars *par, struct grid *g, molData *m){ +popsout(inputPars *par, struct grid *gp, molData *md){ FILE *fp; int j,k,l; double dens; @@ -22,13 +22,13 @@ popsout(inputPars *par, struct grid *g, molData *m){ exit(1); } fprintf(fp,"# Column definition: x, y, z, H2 density, kinetic gas temperature, molecular abundance, convergence flag, pops_0...pops_n\n"); - for(j=0;jncell-par->sinkPoints;j++){ + for(j=0;jpIntensity;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); - for(k=0;kcollPart;l++) dens+=gp[j].dens[l]; + fprintf(fp,"%e %e %e %e %e %e %d ", gp[j].x[0], gp[j].x[1], gp[j].x[2], dens, gp[j].t[0], gp[j].nmol[0]/dens, gp[j].conv); + for(k=0;ktcmb; g[i].t[1]=par->tcmb; g[i].dopb=0.; + for(j=0;j 0){ newdist=numerator/denominator; if(newdist<*ds && newdist > cutoff){ *ds=newdist; - *nposn=g[posn].neigh[i]->id; + *nposn=gp[posn].neigh[i]->id; } } } @@ -68,7 +68,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 *gp, 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. @@ -97,10 +97,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(). */ @@ -144,20 +144,20 @@ 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,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(&jnu,&alpha,m,vfac,g,posn,molI,lineI); + sourceFunc_line(&jnu,&alpha,md,vfac,gp,posn,molI,lineI); } } - 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); + if(img[im].doline && img[im].trans > -1) sourceFunc_cont(&jnu,&alpha,gp,posn,0,img[im].trans); + else if(img[im].doline && img[im].trans == -1) sourceFunc_cont(&jnu,&alpha,gp,posn,0,tmptrans); + else sourceFunc_cont(&jnu,&alpha,gp,posn,0,0); 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 @@ -176,11 +176,11 @@ Note that the algorithm employed here is similar to that employed in the functio /* Add or subtract cmb. */ #ifdef FASTEXP for(ichan=0;ichannSpecies, m, &counta, &countb, &nlinetot); + lineCount(par->nSpecies, md, &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].freq < 0) img[im].freq=md[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){ @@ -226,10 +226,10 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ if(img[im].trans<0){ iline=0; - minfreq=fabs(img[im].freq-m[0].freq[iline]); + minfreq=fabs(img[im].freq-md[0].freq[iline]); tmptrans=iline; - for(iline=1;ilinencell && !g[i].sink;i++){ + qhull(par, gp); + distCalc(par, gp); + + for(i=0;incell && !gp[i].sink;i++){ mindist=1e30; cn=-1; - for(k=0;ksink==0){ - if(g[i].ds[k]sink==0){ + if(gp[i].ds[k]radius-sqrt(g[i].x[0]*g[i].x[0] + g[i].x[1]*g[i].x[1] + g[i].x[2]*g[i].x[2])radius-sqrt(gp[i].x[0]*gp[i].x[0] + gp[i].x[1]*gp[i].x[1] + gp[i].x[2]*gp[i].x[2])-1) { for(k=0;kradiusSqu && (move[0]*move[0]+move[1]*move[1]+move[2]*move[2])>par->minScaleSqu){ - for(k=0;kpIntensity;incell;i++){ mindist=1e30; cn=-1; - for(k=0;ksink==1){ - if(g[i].ds[k]sink==1){ + if(gp[i].ds[k]id; + j=gp[i].neigh[cn]->id; for(k=0;kradius/sqrt(g[i].x[0]*g[i].x[0] + g[i].x[1]*g[i].x[1] + g[i].x[2]*g[i].x[2]); + dist=par->radius/sqrt(gp[i].x[0]*gp[i].x[0] + gp[i].x[1]*gp[i].x[1] + gp[i].x[2]*gp[i].x[2]); for(k=0;kpIntensity;incell;i++){ g[i].a0=malloc(g[i].numNeigh*sizeof(double)); g[i].a1=malloc(g[i].numNeigh*sizeof(double)); - for(j=0;j<3;j++) g[i].vel[j]=0.; for(j=0;j Date: Fri, 17 Jun 2016 15:23:52 +0200 Subject: [PATCH 05/13] Changed velocity interpolation indep. variable Velocity polynomial interpolation coefficients are now calculated with respect to an independent variable d/ds-0.5 instead of d. The new variable is symmetrical about the midpoint of the distance ds between grid points and so allows easy conversion of the coefficients from one member of a pair of connected points to the other. I took the opportunity to rename and rewrite functions as follows: - velospline.c:getVelosplines() renamed to calcInterpCoeffs() (same for _lin() versions). - photon.c:velocityspline() renamed to calcAvRelLineAmp() (same for _lin() versions). - Added photon.c:interpolate(). --- src/grid.c | 2 +- src/lime.h | 13 ++-- src/photon.c | 144 ++++++++++++++++++++++++--------------- src/popsin.c | 2 +- src/predefgrid.c | 2 +- src/velospline.c | 172 +++++++++++++++++++++++++++++------------------ 6 files changed, 205 insertions(+), 130 deletions(-) diff --git a/src/grid.c b/src/grid.c index 57d16a6..0d2a9d3 100644 --- a/src/grid.c +++ b/src/grid.c @@ -627,7 +627,7 @@ buildGrid(inputPars *par, struct grid *g, molData *md){ // getArea(par,g, ran); // getMass(par,g, ran); - getVelosplines(par,g); + calcInterpCoeffs(par,g); /* Mallocs and sets .a0, .a1 etc. */ dumpGrid(par,g); gsl_rng_free(ran); diff --git a/src/lime.h b/src/lime.h index ed96fd3..7ddd6c7 100644 --- a/src/lime.h +++ b/src/lime.h @@ -210,6 +210,10 @@ void gasIIdust(double,double,double,double *); void binpopsout(inputPars *, struct grid *, molData *); void buildGrid(inputPars *, struct grid *, molData *); +void calcAvRelLineAmp(struct grid*, int, int, double, double, double*); +void calcAvRelLineAmp_lin(struct grid*, int, int, double, double, double*); +void calcInterpCoeffs(inputPars*, struct grid*); +void calcInterpCoeffs_lin(inputPars*, struct grid*); void calcSourceFn(double dTau, const inputPars *par, double *remnantSnu, double *expDTau); void constructLinkArrays(unsigned int, struct grid *, struct linkType **\ , unsigned int *, struct linkType ***, unsigned int **, unsigned int *); @@ -221,8 +225,6 @@ void distCalc(inputPars *, struct grid *); 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); @@ -232,9 +234,10 @@ void getjbar(int, molData *, struct grid *, inputPars *,gridPointData *,doubl 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 *); -void getVelosplines(inputPars *, struct grid *); -void getVelosplines_lin(inputPars *, struct grid *); void gridAlloc(inputPars *, struct grid **); +void input(inputPars *, image *); +double interpolate(double, double, double, double, double, double); +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); @@ -273,8 +276,6 @@ 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*); double veloproject(double *, double *); void writeGridExtToFits(fitsfile *, inputPars, unsigned short, struct grid *, unsigned int *, char **); int writeGridToFits(char *, inputPars, unsigned short, unsigned short, struct grid *, molData *, char **); diff --git a/src/photon.c b/src/photon.c index 81930bf..9b3af09 100644 --- a/src/photon.c +++ b/src/photon.c @@ -47,70 +47,100 @@ sortangles(double *inidir, int id, struct grid *g, const gsl_rng *ran) { } } +inline double interpolate(double a0, double a1, double a2, double a3, double a4, double s){ + /* +Calculates a polynomial interpolation of velocity at location d along an interval of length ds + __4 + \ + v(d) ~ > a_j*s^j + /_j=0 -void -velocityspline(struct grid *g, int id, int k, double binv, double deltav, double *vfac){ - int nspline,ispline,naver,iaver; - double v1,v2,s1,s2,sd,v,vfacsub,d; +where it is assumed that + + d + s = --- - 0.5. + ds + */ + double vel; + vel = (((a4*s + a3)*s + a2)*s + a1)*s + a0; + return vel; +} + +void calcAvRelLineAmp(struct grid *gp, int id, int k, double oneOnB, double deltav, double *avRelLineAmp){ + /* +This performs a numerical integration of the amplitude of a spectral line along the link or edge between two grid points. The line shapes are Gaussian and are normalized such that the maximum equals unity. The integral may be expressed as + + ds + 1 / (-[deltav-v{x}]^2 ) + I = ---- | dx exp(-----------------) + ds / ( B ) + 0 + +where ds is the distance between the two grid points and v(x) is the component of vector velocity in the direction of the 'edge' between the grid points. v(x) is estimated via polynomial interpolation - see function 'interpolate()' for the exact formula. + */ + + int nCoarseSamples,i,nFineSamples,j; + double v1,v2,vOld,vNew,dFracOld,dFracNew,dFrac,v,fineRelLineAmp; - v1=deltav-veloproject(g[id].dir[k].xn,g[id].vel); - v2=deltav-veloproject(g[id].dir[k].xn,g[id].neigh[k]->vel); + v1 = interpolate(gp[id].a0[k], gp[id].a1[k], gp[id].a2[k], gp[id].a3[k], gp[id].a4[k], -0.5); + v2 = interpolate(gp[id].a0[k], gp[id].a1[k], gp[id].a2[k], gp[id].a3[k], gp[id].a4[k], 0.5); - nspline=(fabs(v1-v2)*binv < 1) ? 1 : (int)(fabs(v1-v2)*binv); - *vfac=0.; - s2=0; - v2=v1; + nCoarseSamples = (fabs(v1-v2)*oneOnB < 1) ? 1 : (int)(fabs(v1-v2)*oneOnB); + *avRelLineAmp = 0.0; + dFracNew = 0.0; + vNew = v1; - for(ispline=0;ispline fabs(v1-v2)*binv) ? 1 : (int)(fabs(v1-v2)*binv); - for(iaver=0;iaver fabs(vOld-vNew)*oneOnB) ? 1 : (int)(fabs(vOld-vNew)*oneOnB); + fineRelLineAmp = 0.0; + for(j=0;jvel); + v1 = gp[id].a0[k] - 0.5*gp[id].a1[k]; + v2 = gp[id].a0[k] + 0.5*gp[id].a1[k]; - nspline=(fabs(v1-v2)*binv < 1) ? 1 : (int)(fabs(v1-v2)*binv); - *vfac=0.; - s2=0; - v2=v1; + nCoarseSamples = (fabs(v1-v2)*oneOnB < 1) ? 1 : (int)(fabs(v1-v2)*oneOnB); + *avRelLineAmp = 0.0; + dFracNew = 0.0; + vNew = v1; - for(ispline=0;ispline fabs(v1-v2)*binv) ? 1 : (int)(fabs(v1-v2)*binv); - for(iaver=0;iaver fabs(vOld-vNew)*oneOnB) ? 1 : (int)(fabs(vOld-vNew)*oneOnB); + fineRelLineAmp = 0.0; + for(j=0;jnSpecies;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) + calcAvRelLineAmp_lin(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); + else + calcAvRelLineAmp( g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); mp[l].vfac[iphot]=vfac[0]; } } else { @@ -225,8 +257,10 @@ 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) + calcAvRelLineAmp_lin(g,here,dir,g[id].mol[l].binv,deltav,&vfac[l]); + else + calcAvRelLineAmp( 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) + calcAvRelLineAmp_lin(g,here,dir,g[id].mol[counta[jline]].binv,deltav-matrix[jline].deltav,&vblend); + else + calcAvRelLineAmp( 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/popsin.c b/src/popsin.c index 5392664..bdf25ad 100644 --- a/src/popsin.c +++ b/src/popsin.c @@ -106,7 +106,7 @@ popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){ qhull(par, *g); distCalc(par, *g); - getVelosplines(par,*g); + calcInterpCoeffs(par,*g); *popsdone=1; } diff --git a/src/predefgrid.c b/src/predefgrid.c index 6a0529c..e60912e 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -74,7 +74,7 @@ predefinedGrid(inputPars *par, struct grid *g){ distCalc(par,g); // getArea(par,g, ran); // getMass(par,g, ran); - getVelosplines_lin(par,g); + calcInterpCoeffs_lin(par,g); if(par->gridfile) write_VTK_unstructured_Points(par, g); gsl_rng_free(ran); } diff --git a/src/velospline.c b/src/velospline.c index 4c4bd2f..23c3316 100644 --- a/src/velospline.c +++ b/src/velospline.c @@ -9,92 +9,130 @@ #include "lime.h" -void -getVelosplines(inputPars *par, struct grid *g){ - int i,k,j,l,s; - double v[5], vel[3], x[3], d; - gsl_matrix *matrix = gsl_matrix_alloc(5,5); - gsl_vector *a = gsl_vector_alloc(5); - gsl_vector *y = gsl_vector_alloc(5); - gsl_permutation *p = gsl_permutation_alloc(5); - +void calcInterpCoeffs(inputPars *par, struct grid *gp){ + /* +The velocity v(d) (or rather the scalar component of vector velocity in the direction of the edge) at a distance d along the line (or 'edge' in triangulation jargon) between a given grid point and its neighbour is approximated in LIME by the polynomial expression + + __N-1 + \ + v(d) ~ > a_j*s^j + /_j=0 + +where + d + s = --- - 0.5. + ds + +The coefficients a_j are calculated in the present function for each grid point (which uses twice the memory needed). This is done by sampling velocity at N points along the edge with at distances evenly spread between 0 and ds, then solving the system of equations to obtain the interpolation coefficients. + +A Chebyshev interpolation would probably be preferred, but we wil leave that for future generations. + +Note that, given coefficients calculated for edge AB, the present definition allows easy conversion from the values calculated for grid point A to those for grid point B: coefficient j for B is just (-1)^j times the respective coefficient for A. + +The number of coefficients N is currently hard-wired at 5. + */ + + int i,k,di,ri,ci,dummy; + const int nCoeffs = NUM_VEL_COEFFS; + double vel[DIM],x[DIM],dFrac[nCoeffs],velComp,sToPower; + gsl_matrix *matrix = gsl_matrix_alloc(nCoeffs,nCoeffs); + gsl_vector *aCoeffs = gsl_vector_alloc(nCoeffs); + gsl_vector *velVector = gsl_vector_alloc(nCoeffs); + gsl_permutation *p = gsl_permutation_alloc(nCoeffs); + + for(ri=0;ripIntensity;i++){ - 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(k=0;kpIntensity;incell;i++){ - 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;jpIntensity;i++){ g[i].a0=malloc(g[i].numNeigh*sizeof(double)); - g[i].a1=malloc(g[i].numNeigh*sizeof(double)); + g[i].a1=malloc(g[i].numNeigh*sizeof(double)); for(k=0;kvel); - g[i].a1[k]=(v[0]-v[1])/(0-g[i].ds[k]); - g[i].a0[k]=v[0]; - } + g[i].a1[k] = v[1] - v[0]; + g[i].a0[k] = (v[0] + v[1])/2.0; + } } - + for(i=par->pIntensity;incell;i++){ g[i].a0=malloc(g[i].numNeigh*sizeof(double)); g[i].a1=malloc(g[i].numNeigh*sizeof(double)); - for(j=0;j Date: Fri, 17 Jun 2016 15:29:14 +0200 Subject: [PATCH 06/13] Fixed #136 --- src/photon.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/photon.c b/src/photon.c index 9b3af09..ff754f5 100644 --- a/src/photon.c +++ b/src/photon.c @@ -99,7 +99,7 @@ where ds is the distance between the two grid points and v(x) is the component o nFineSamples = (1 > fabs(vOld-vNew)*oneOnB) ? 1 : (int)(fabs(vOld-vNew)*oneOnB); fineRelLineAmp = 0.0; for(j=0;j fabs(vOld-vNew)*oneOnB) ? 1 : (int)(fabs(vOld-vNew)*oneOnB); fineRelLineAmp = 0.0; for(j=0;j Date: Fri, 17 Jun 2016 16:48:31 +0200 Subject: [PATCH 07/13] NITERATIONS now user-settable. I replaced the NITERATIONS macro by a user-settable parameter par->nSolveIters. Mostly this is to make it easier to iterate a previous population solution further, when this can eventually be implemented. --- src/aux.c | 12 +++++++----- src/lime.h | 4 ++-- src/messages.c | 11 +++++------ 3 files changed, 14 insertions(+), 13 deletions(-) diff --git a/src/aux.c b/src/aux.c index f72cfe1..aeaac96 100644 --- a/src/aux.c +++ b/src/aux.c @@ -41,6 +41,7 @@ parseInput(inputPars *par, image **img, molData **m){ par->sinkPoints=0; par->doPregrid=0; par->nThreads=0; + par->nSolveIters=17; for(i=0;iwriteGridAtStage[i] = 0; @@ -443,7 +444,7 @@ lineBlend(molData *m, inputPars *par, blend **matrix){ void levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ - int id,iter,ilev,prog=0,ispec,c=0,n,i,threadI,nVerticesDone,nItersDone; + int id,iter,ilev,ispec,c=0,n,i,threadI,nVerticesDone,nItersDone; double percent=0.,*median,result1=0,result2=0,snr,delta_pop; blend *matrix; struct statistics { double *pop, *ave, *sigma; } *stat; @@ -506,8 +507,8 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ if(par->lte_only==0){ nItersDone=0; - do{ - if(!silent) progressbar2(0, prog++, 0, result1, result2); + 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;idncell && !g[id].sink;id++){ for(ilev=0;ilevoutputfile != NULL) popsout(par,g,m); - } while(nItersDone++binoutputfile != NULL) binpopsout(par,g,m); } diff --git a/src/lime.h b/src/lime.h index 7ddd6c7..426587a 100644 --- a/src/lime.h +++ b/src/lime.h @@ -63,7 +63,6 @@ #define GRAV 6.67428e-11 /* Other constants */ -#define NITERATIONS 16 #define max_phot 10000 /* don't set this value higher unless you have enough memory. */ #define ininphot 9 #define minpop 1.e-6 @@ -94,6 +93,7 @@ typedef struct { char **moldatfile; _Bool writeGridAtStage[NUM_GRID_STAGES]; char *gridFitsOutSets[NUM_GRID_STAGES]; + int nSolveIters; } inputPars; /* Molecular data: shared attributes */ @@ -298,7 +298,7 @@ void greetings_parallel(int); void screenInfo(); void printDone(int); void progressbar(double,int); -void progressbar2(int,int,double,double,double); +void progressbar2(inputPars*, int,int,double,double,double); void casaStyleProgressBar(const int,int); void goodnight(int, char *); void quotemass(double); diff --git a/src/messages.c b/src/messages.c index 6dbe387..013e245 100644 --- a/src/messages.c +++ b/src/messages.c @@ -152,10 +152,10 @@ progressbar(double percent, int line){ } void -progressbar2(int flag, int prog, double percent, double minsnr, double median){ +progressbar2(inputPars *par, int flag, int prog, double percent, double minsnr, double median){ #ifdef NO_NCURSES if (flag == 0) { - printf(" Iteration %i / max %i: Starting\n", prog + 1, NITERATIONS + 1); + printf(" Iteration %i / max %i: Starting\n", prog + 1, par->nSolveIters); } 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 @@ -355,13 +355,12 @@ processFitsError(int status){ /* Print out cfitsio error messages and exit program */ /*****************************************************/ - char message[80]; - - if (status){ + 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 From f94c566a6b2add0c400e0c8b9da8dd4df2c147e6 Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 17 Jun 2016 17:08:38 +0200 Subject: [PATCH 08/13] New interface for IO of grid data. This is a substantial commit, but it is hard to see how it could have been easily broken into smaller ones. Here I have reworked the previous facility for IO to and from FITS-format files. There is another, generic IO layer now implemented in a new module gridio.c, with FITS relegated to a choie of file format. Hopefully this will make it easier if anybody else wants to write an interface to HDF5 or a DB or whatever. Some detail on the changes: - New pars par->gridInFile, par->gridOutFiles[] (user-settable) and par->writeGridAtStage[], par->dataStageI (not intended to be user settable). - New function grid.c:mallocAndSetDefaultGrid(). - grid.c:buildGrid() rewritten and renamed readOrBuildGrid(). - main.c changed to include the fix to #130. --- Makefile | 4 +- src/aux.c | 23 +- src/grid.c | 357 +++++++++++----- src/grid2fits.c | 1088 ++++++++++++++++++----------------------------- src/gridio.c | 798 ++++++++++++++++++++++++++++++++++ src/lime.h | 148 ++++--- src/main.c | 64 ++- 7 files changed, 1606 insertions(+), 876 deletions(-) create mode 100644 src/gridio.c diff --git a/Makefile b/Makefile index 699d10e..43274a7 100644 --- a/Makefile +++ b/Makefile @@ -58,7 +58,7 @@ SRCS = src/aux.c src/messages.c src/grid.c src/LTEsolution.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/grid2fits.c \ - src/tcpsocket.c src/defaults.c src/fastexp.c + src/tcpsocket.c src/defaults.c src/fastexp.c src/gridio.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 \ @@ -67,7 +67,7 @@ OBJS = src/aux.o src/messages.o src/grid.o src/LTEsolution.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/grid2fits.o \ - src/tcpsocket.o src/defaults.o src/fastexp.o + src/tcpsocket.o src/defaults.o src/fastexp.o src/gridio.o MODELO = src/model.o #CCFLAGS = -O3 -falign-loops=16 -fno-strict-aliasing -DTEST diff --git a/src/aux.c b/src/aux.c index aeaac96..8d92a39 100644 --- a/src/aux.c +++ b/src/aux.c @@ -29,6 +29,7 @@ parseInput(inputPars *par, image **img, molData **m){ par->gridfile = NULL; par->pregrid = NULL; par->restart = NULL; + par->gridInFile = NULL; par->tcmb = 2.728; par->lte_only=0; @@ -44,9 +45,10 @@ parseInput(inputPars *par, image **img, molData **m){ par->nSolveIters=17; for(i=0;iwriteGridAtStage[i] = 0; - par->gridFitsOutSets[i] = ""; + par->writeGridAtStage[i] = 0; /* The user is not expected to set this. */ + par->gridOutFiles[i] = ""; }; + par->dataStageI = 0; /* The user is not expected to set this. */ /* Allocate space for output fits images */ (*img)=malloc(sizeof(image)*MAX_NSPECIES); @@ -121,6 +123,11 @@ parseInput(inputPars *par, image **img, molData **m){ while(temp[i++]>-1) par->collPart++; } + for(i=0;igridOutFiles[i] != "") + par->writeGridAtStage[i] = 1; + }; + /* Now we need to calculate the cutoff value used in calcSourceFn(). The issue is to decide between @@ -147,6 +154,16 @@ The cutoff will be the value of abs(x) for which the error in the exact expressi } } + if(par->gridInFile!=NULL){ + if((fp=fopen(par->gridInFile, "r"))==NULL){ + if(!silent) bail_out("Error opening grid input file!"); + exit(1); + } + else { + fclose(fp); + } + } + /* Allocate pixel space and parse image information */ for(i=0;inImages;i++){ if((*img)[i].nchan == 0 && (*img)[i].velres<0 ){ @@ -603,6 +620,8 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ if(par->binoutputfile != NULL) binpopsout(par,g,m); } + par->dataStageI = 4; + for (i=0;inThreads;i++){ gsl_rng_free(threadRans[i]); } diff --git a/src/grid.c b/src/grid.c index 0d2a9d3..ba959a2 100644 --- a/src/grid.c +++ b/src/grid.c @@ -504,139 +504,264 @@ getMass(inputPars *par, struct grid *g, const gsl_rng *ran){ } -void -buildGrid(inputPars *par, struct grid *g, molData *md){ +void mallocAndSetDefaultGrid(struct grid **gp, const unsigned int numPoints){ + unsigned int i; + + *gp = malloc(sizeof(struct grid)*numPoints); + for(i=0;iminScale */ double r,theta,phi,sinPhi,x,y,z,semiradius; /* Coordinates */ double temp; int k=0,i,j; /* counters */ - int flag,stageI,status=0; - char **collPartNames=NULL; /*** this is a placeholder until we start reading these. */ + int flag; + struct gridInfoType gridInfoRead; + int status; + char **collPartNames; + int numCollPartRead; + char message[80]; + + par->dataStageI = 0; + if(par->gridInFile!=NULL){ + /* Set defaults for *gridInfoRead: + */ + gridInfoRead.nInternalPoints = 0; + gridInfoRead.nSinkPoints = 0; + gridInfoRead.nLinks = 0; + gridInfoRead.nNNIndices = 0; + gridInfoRead.nDims = 0; + gridInfoRead.nSpecies = 0; + gridInfoRead.nDensities = 0; + gridInfoRead.nACoeffs = 0; + gridInfoRead.mols = NULL; + + status = readGrid(par->gridInFile, lime_FITS, &gridInfoRead, gp\ + , &collPartNames, &numCollPartRead, &(par->dataStageI)); + + if(status){ + if(!silent){ + sprintf(message, "Read of grid file failed with status return %d", status); + bail_out(message); + } + exit(1); + } - gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2); /* Random number generator */ + /* Test gridInfoRead values against par values and overwrite the latter, with a warning, if necessary. + */ + if(gridInfoRead.nSinkPoints>0 && par->sinkPoints>0){ + if(gridInfoRead.nSinkPoints!=par->sinkPoints){ + if(!silent) warning("par->sinkPoints will be overwritten"); + par->sinkPoints = gridInfoRead.nSinkPoints; + } + if(gridInfoRead.nInternalPoints!=par->pIntensity){ + if(!silent) warning("par->pIntensity will be overwritten"); + par->pIntensity = 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 && 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->collPart>0 && gridInfoRead.nDensities!=par->collPart){ + if(!silent){ + sprintf(message, "Grid file had %d densities but you have provided %d."\ + , (int)gridInfoRead.nDensities, par->collPart); + 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? +*/ + if(par->dataStageI==4 && !silent) + warning("Sorry, LIME is not yet able to further process level populations read from file."); + } /* End of read grid file. Whether and what we subsequently calculate will depend on the value of par->dataStageI returned. */ + + if(par->dataStageI<1){ /* This can only happen if we did not read a file. */ + mallocAndSetDefaultGrid(gp, (unsigned int)par->ncell); + + gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2); /* Random number generator */ #ifdef TEST - gsl_rng_set(ran,342971); + gsl_rng_set(ran,342971); #else - gsl_rng_set(ran,time(0)); + gsl_rng_set(ran,time(0)); #endif - lograd=log10(par->radius); - logmin=log10(par->minScale); - - /* Sample pIntensity number of points */ - for(k=0;kpIntensity;k++){ - temp=gsl_rng_uniform(ran); - flag=0; - /* Pick a point and check if we like it or not */ - do{ - if(par->sampling==0){ - r=pow(10,logmin+gsl_rng_uniform(ran)*(lograd-logmin)); - theta=2.*PI*gsl_rng_uniform(ran); - phi=PI*gsl_rng_uniform(ran); - sinPhi=sin(phi); - x=r*cos(theta)*sinPhi; - y=r*sin(theta)*sinPhi; - if(DIM==3) z=r*cos(phi); - else z=0.; - } else if(par->sampling==1){ - x=(2*gsl_rng_uniform(ran)-1)*par->radius; - y=(2*gsl_rng_uniform(ran)-1)*par->radius; - if(DIM==3) z=(2*gsl_rng_uniform(ran)-1)*par->radius; - else z=0; - } else if(par->sampling==2){ - r=pow(10,logmin+gsl_rng_uniform(ran)*(lograd-logmin)); - theta=2.*PI*gsl_rng_uniform(ran); - if(DIM==3) { - z=2*gsl_rng_uniform(ran)-1.; - semiradius=r*sqrt(1.-z*z); - z*=r; + lograd=log10(par->radius); + logmin=log10(par->minScale); + + /* Sample pIntensity number of points */ + for(k=0;kpIntensity;k++){ + temp=gsl_rng_uniform(ran); + flag=0; + /* Pick a point and check if we like it or not */ + do{ + if(par->sampling==0){ + r=pow(10,logmin+gsl_rng_uniform(ran)*(lograd-logmin)); + theta=2.*PI*gsl_rng_uniform(ran); + phi=PI*gsl_rng_uniform(ran); + sinPhi=sin(phi); + x=r*cos(theta)*sinPhi; + y=r*sin(theta)*sinPhi; + if(DIM==3) z=r*cos(phi); + else z=0.; + } else if(par->sampling==1){ + x=(2*gsl_rng_uniform(ran)-1)*par->radius; + y=(2*gsl_rng_uniform(ran)-1)*par->radius; + if(DIM==3) z=(2*gsl_rng_uniform(ran)-1)*par->radius; + else z=0; + } else if(par->sampling==2){ + r=pow(10,logmin+gsl_rng_uniform(ran)*(lograd-logmin)); + theta=2.*PI*gsl_rng_uniform(ran); + if(DIM==3) { + z=2*gsl_rng_uniform(ran)-1.; + semiradius=r*sqrt(1.-z*z); + z*=r; + } else { + z=0.; + semiradius=r; + } + x=semiradius*cos(theta); + y=semiradius*sin(theta); } else { - z=0.; - semiradius=r; + if(!silent) bail_out("Don't know how to sample model"); + exit(1); } - x=semiradius*cos(theta); - y=semiradius*sin(theta); - } else { - if(!silent) bail_out("Don't know how to sample model"); - exit(1); - } - if((x*x+y*y+z*z)radiusSqu) flag=pointEvaluation(par,temp,x,y,z); - } while(!flag); - /* Now pointEvaluation has decided that we like the point */ - - /* Assign values to the k'th grid point */ - /* I don't think we actually need to do this here... */ - g[k].id=k; - g[k].x[0]=x; - g[k].x[1]=y; - if(DIM==3) g[k].x[2]=z; - else g[k].x[2]=0.; - - g[k].sink=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); - if(!silent) progressbar((double) k/((double)par->pIntensity-1), 4); + if((x*x+y*y+z*z)radiusSqu) flag=pointEvaluation(par,temp,x,y,z); + } while(!flag); + /* Now pointEvaluation has decided that we like the point */ + + /* Assign values to the k'th grid point */ + /* I don't think we actually need to do this here... */ + (*gp)[k].id=k; + (*gp)[k].x[0]=x; + (*gp)[k].x[1]=y; + if(DIM==3) (*gp)[k].x[2]=z; + else (*gp)[k].x[2]=0.; + + (*gp)[k].sink=0; + + if(!silent) progressbar((double) k/((double)par->pIntensity-1), 4); + } + /* end model grid point assignment */ + if(!silent) printDone(4); + + /* Add surface sink particles */ + for(k=par->pIntensity;kncell;k++){ + theta=gsl_rng_uniform(ran)*2*PI; + if(DIM==3) z=2*gsl_rng_uniform(ran)-1.; + else z=0.; + semiradius=sqrt(1.-z*z); + x=semiradius*cos(theta); + y=semiradius*sin(theta);; + (*gp)[k].id=k; + (*gp)[k].x[0]=par->radius*x; + (*gp)[k].x[1]=par->radius*y; + (*gp)[k].x[2]=par->radius*z; + (*gp)[k].sink=1; + } + /* end grid allocation */ + + gsl_rng_free(ran); + + smooth(par, *gp); + if(!silent) printDone(5); + + par->dataStageI = 1; + }else if(par->dataStageI==1 && par->writeGridAtStage[par->dataStageI-1]){ + sprintf(message, "You just read a grid file at data stage %d, now you want to write it again?", par->dataStageI); + if(!silent) warning(message); } - /* end model grid point assignment */ - if(!silent) printDone(4); - - /* Add surface sink particles */ - for(i=0;isinkPoints;i++){ - theta=gsl_rng_uniform(ran)*2*PI; - if(DIM==3) z=2*gsl_rng_uniform(ran)-1.; - else z=0.; - semiradius=sqrt(1.-z*z); - x=semiradius*cos(theta); - y=semiradius*sin(theta);; - g[k].id=k; - g[k].x[0]=par->radius*x; - g[k].x[1]=par->radius*y; - g[k].x[2]=par->radius*z; - g[k].sink=1; - g[k].abun[0]=0; - g[k].dens[0]=1e-30; - g[k].t[0]=par->tcmb; - g[k].t[1]=par->tcmb; - g[k].dopb=0.; - for(j=0;jdataStageI==1) /* Only happens if (i) we read no file, or (ii) we read a file at dataStageI==1. */ + writeGridIfRequired(par, *gp, NULL, lime_FITS); + + if(par->dataStageI<2){ + qhull(par, *gp); /* Mallocs and sets .neigh, sets .numNeigh */ + + par->dataStageI = 2; + }else if(par->dataStageI==2 && par->writeGridAtStage[par->dataStageI-1]){ + sprintf(message, "You just read a grid file at data stage %d, now you want to write it again?", par->dataStageI); + if(!silent) warning(message); } - /* end grid allocation */ - smooth(par,g); - qhull(par, g); - distCalc(par, g); + if(par->dataStageI==2) /* Only happens if (i) we read no file, or (ii) we read a file at dataStageI==2. */ + writeGridIfRequired(par, *gp, NULL, lime_FITS); -/* Can't do this yet because .dens is not NULL, but .a0, .a1 etc are. Have to rearrange the mallocs of g elements first. Simply set g[0].dens=NULL to kluge around this? But then have to malloc it again (sigh). - stageI = 1; - if(par->writeGridAtStage[stageI]) - status = writeGridToFits(par->gridFitsOutSets[stageI], *par, (unsigned short)DIM\ - , (unsigned short)NUM_VEL_COEFFS, g, md, collPartNames); -*/ + distCalc(par, *gp); /* Mallocs and sets .dir & .ds, sets .nphot */ - 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); - 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); - } - - // getArea(par,g, ran); - // getMass(par,g, ran); - calcInterpCoeffs(par,g); /* Mallocs and sets .a0, .a1 etc. */ - dumpGrid(par,g); - - gsl_rng_free(ran); - if(!silent) printDone(5); - - stageI = 2; - if(par->writeGridAtStage[stageI]) - status = writeGridToFits(par->gridFitsOutSets[stageI], *par, (unsigned short)DIM\ - , (unsigned short)NUM_VEL_COEFFS, g, md, collPartNames); + if(par->dataStageI<3){ + for(i=0;incell; i++){ + (*gp)[i].dens = malloc(sizeof(double)*par->collPart); + (*gp)[i].abun = malloc(sizeof(double)*par->nSpecies); + (*gp)[i].nmol = malloc(sizeof(double)*par->nSpecies); + } + + for(i=0;ipIntensity;i++){ + density( (*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2], (*gp)[i].dens); + temperature((*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2], (*gp)[i].t); + doppler( (*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2],&(*gp)[i].dopb); + abundance( (*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2], (*gp)[i].abun); + velocity( (*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2], (*gp)[i].vel); + } + + for(i=par->pIntensity;incell;i++){ + (*gp)[i].dens[0]=1e-30; + (*gp)[i].t[0]=par->tcmb; + (*gp)[i].t[1]=par->tcmb; + (*gp)[i].dopb=0.; + (*gp)[i].abun[0]=0; + for(j=0;jdataStageI = 3; + }else if(par->dataStageI==3 && par->writeGridAtStage[par->dataStageI-1]){ + sprintf(message, "You just read a grid file at data stage %d, now you want to write it again?", par->dataStageI); + if(!silent) warning(message); + } + + if(par->dataStageI==3) /* Only happens if (i) we read no file, or (ii) we read a file at dataStageI==3. */ + writeGridIfRequired(par, *gp, NULL, lime_FITS); + + dumpGrid(par,*gp); } diff --git a/src/grid2fits.c b/src/grid2fits.c index fd2bc27..b806e14 100644 --- a/src/grid2fits.c +++ b/src/grid2fits.c @@ -9,203 +9,92 @@ 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()? + - If any of the stage 3 columns (VEL etc) are not present, get the values via user-supplied functions. +*/ #include "lime.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. - -1) Different amounts of information in grid. --------------------------------------------- -Describing the format of the FITS file is complicated by the fact that the amount of information stored in the grid data structure tends to increase during the run of the LIME program. One can distingish 4 stages of completion which are at least conceptually different, even if they are not cleanly separated in the present implementation of the code. These are summarized as follows. - - A) At this stage the vector of grid objects has been malloc'd and values have been generated for the following struct elements: - id - x - sink - This stage is flagged to the software by the element 'neigh' of the first grid point being set to NULL. - - B) This stage is entered after the Delaunay neighbours of each grid point have been determined. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: - numNeigh - neigh - This stage is flagged by the element 'dens' of the first grid point being set to NULL. - - C) This stage is entered after sampling the user-supplied functions for density, velocity etc. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: - vel - a0, a1 etc. - dens - t - dopb - abun - This stage is flagged by the element 'mol' of the first grid point being set to NULL. +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 minimum stage integer associated with the presence of a particular extension, column or keyword is given on the leftmost place of each line. - D) After 1 or more iterations of populating the levels, we enter stage D, in which all the grid struct elements have been malloc'd and given at least preliminary values; specifically now the element - mol - -(Other grid struct values are deducible from the ones given or are used for temporary storage only.) - -The dependency relationship tree for the stages may be diagrammed as follows: - - A -> B -> C -> D. - -2) The FITS file format for stage D ------------------------------------ -This is defined in the following list of extensions. All extensions are binary table except where indicated. The letter in the second row for column descriptions gives the FITS data type. See eg +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. - 1) GRID +1 0) The primary HDU + Keywords: +1 STAGE I + +1 1) GRID Number of rows = number of grid points. Keywords: - RADIUS - COLLPARn # 1 for each nth collision partner. +1 RADIUS D +1 COLLPARn A # 1 for each nth collision partner. Columns: - ID V - Xj D # Cartesian components of the point location, 1 col per jth dimension. - IS_SINK L # =True iff the point lies on the edge of the model. - NUMNEIGH U - FIRST_NN V # See explanation in section 3 below. - VELj D # 1 col per jth dimension. - TURBDPLR E # Given Gaussian lineshape exp(-v^2/[B^2 + 2*k*T/m]), this is B. - DENSITYn D # 1 per nth collision partner. - TEMPKNTC E # From t[0]. - TEMPDUST E # From t[1]. - ABUNMOLm E # 1 per mth molecular species. - - 2) NN_INDICES (see explanation in section 3 below) +1 ID V +1 Xj D # Cartesian components of the point location, 1 col per jth dimension. +1 IS_SINK L # =True iff the point lies on the edge of the model. +2 NUMNEIGH U +2 FIRST_NN V # See explanation in section 3 below. +3 VELj D # 1 col per jth dimension. +3 DENSITYn D # 1 per nth collision partner. +3 TEMPKNTC E # From t[0]. +3 TEMPDUST E # From t[1]. +3 TURBDPLR E # Given Gaussian lineshape exp(-v^2/[B^2 + 2*k*T/m]), this is B. +3 ABUNMOLm E # 1 per mth molecular species. + +2 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: - LINK_I V +2 LINK_I V - 3) LINKS (see explanation in section 3 below) +2 3) LINKS (see explanation in the header to gridio.c) Number of rows = number of Delaunay links. Columns: - GRID_I_1 V - GRID_I_2 V - ACOEFF_p D # 1 per pth order of the velocity polynomial. +2 GRID_I_1 V +2 GRID_I_2 V +3 ACOEFF_p D # 1 per pth order of the velocity polynomial. - 4 etc) LEVEL_POPS_m (1 per mth molecular species) +4 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: - MOL_NAME - -The subsets of these FITS structures which are written under the other 3 stages are described in the header comment to the function writeGridToFits(). - -3) The NN_INDICES and LINKS extensions and the FIRST_NN column. ---------------------------------------------------------------- -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. +4 MOL_NAME -To meet the needs of FITS 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;j at least stage B. */ - constructLinkArrays((unsigned int)par.ncell, g, &links, &totalNumLinks\ - , &nnLinks, &firstNearNeigh, &totalNumNeigh); - } +/*....................................................................*/ +fitsfile *openFITSFileForWrite(char *outFileName, const int dataStageI){ + fitsfile *fptr=NULL; + int status=0; + char negfile[100]="! "; + short dsiShort = (short)dataStageI, *ptrToShort=&dsiShort; /* Have to do this otherwise the compiler complains that the const int is not an 'lvalue'. */ fits_create_file(&fptr, outFileName, &status); if(status!=0){ @@ -216,155 +105,26 @@ The pseudo-code is as follows: processFitsError(status); } - /* GRID - */ - writeGridExtToFits(fptr, par, numDims, g, firstNearNeigh, collPartNames); /* Deals internally with stages B and C. */ - - if (g[0].neigh!=NULL){ /* => at least stage B. */ - /* NN_INDICES - */ - writeNnIndicesExtToFits(fptr, totalNumNeigh, nnLinks, links); - - /* LINKS - */ - if(g[0].dens!=NULL) stageC = 1; - writeLinksExtToFits(fptr, stageC, totalNumLinks, numACoeffs, links); - - if (g[0].mol!=NULL){ /* => stage D. */ - for(speciesI=0;speciesInumNeigh;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].g[0] = gAPtr; - (*links)[li].g[1] = gBPtr; - - li++; - } - - nnLinkIs[ni] = linkId; - ni++; - } - - pointIsDone[idA] = 1; - } - *totalNumLinks = li; - - *links = realloc(*links, sizeof(struct linkType)*(*totalNumLinks)); - - if(g[0].dens!=NULL){ /* => at least stage C. */ - for(li=0;li<*totalNumLinks;li++){ - gAPtr = (*links)[li].g[0]; - /* Find which neighbour of gAPtr corresponds to the link: */ - linkNotFound = 1; - for(jA=0;jAnumNeigh;jA++){ - idB = gAPtr->neigh[jA]->id; - if(linkNotFound && idB==(*links)[li].g[1]->id){ - (*links)[li].aCoeffs[0] = gAPtr->a0[jA]; - (*links)[li].aCoeffs[1] = gAPtr->a1[jA]; - (*links)[li].aCoeffs[2] = gAPtr->a2[jA]; - (*links)[li].aCoeffs[3] = gAPtr->a3[jA]; - (*links)[li].aCoeffs[4] = gAPtr->a4[jA]; - linkNotFound = 0; - } - } - if(linkNotFound){ - sprintf(message, "SNAFU with link %d grid point %d", li, gAPtr->id); - if(!silent) bail_out(message); - exit(1); - } - } - } - - *nnLinks = malloc(sizeof(**nnLinks)*(*totalNumNeigh)); - for(ni=0;ni<(*totalNumNeigh);ni++) - (*nnLinks)[ni] = &(*links)[nnLinkIs[ni]]; +void closeFITSFile(fitsfile *fptr){ + int status=0; - free(nnLinkIs); - free(pointIsDone); - /* The calling routine must free firstNearNeigh, nnLinks, links. */ + fits_close_file(fptr, &status); + processFitsError(status); } /*....................................................................*/ -void defineGridExtColumns(unsigned short numKwdChars, unsigned short numDims\ - , unsigned short numCollPart, unsigned short numSpecies, char *ttype[]\ +void defineGridExtColumns(const unsigned short numKwdChars, inputPars par\ + , const unsigned short numDims, const int dataStageI, char *ttype[]\ , char *tform[], char *tunit[], int dataTypes[]){ /* Note that data types in all capitals are defined in fitsio.h. */ @@ -381,17 +141,17 @@ void defineGridExtColumns(unsigned short numKwdChars, unsigned short numDims\ exit(1); } - if(numSpecies>maxNumSpecies){ + if(par.nSpecies>maxNumSpecies){ if(!silent){ - sprintf(message, "Caller asked for %d species but colnames can only be written for %d.", numSpecies, maxNumSpecies); + sprintf(message, "Caller asked for %d species but colnames can only be written for %d.", par.nSpecies, maxNumSpecies); bail_out(message); } exit(1); } - if(numCollPart>maxNumCollPart){ + if(par.collPart>maxNumCollPart){ if(!silent){ - sprintf(message, "Caller asked for %d coll. part. but colnames can only be written for %d.", numCollPart, maxNumCollPart); + sprintf(message, "Caller asked for %d coll. part. but colnames can only be written for %d.", par.collPart, maxNumCollPart); bail_out(message); } exit(1); @@ -421,6 +181,9 @@ void defineGridExtColumns(unsigned short numKwdChars, unsigned short numDims\ tunit[colI] = "\0"; dataTypes[colI] = TLOGICAL; + if(dataStageI<2) + return; + colI++; ttype[colI] = malloc(numKwdChars); sprintf(ttype[colI], "NUMNEIGH"); @@ -435,6 +198,9 @@ void defineGridExtColumns(unsigned short numKwdChars, unsigned short numDims\ tunit[colI] = "\0"; dataTypes[colI] = TUINT; + if(dataStageI<3) + return; + /* should rather have a vector column? */ for(i=0;i at least stage B. */ + if (dataStageI>1){ colI++; numNeigh = malloc(sizeof(*numNeigh)*par.ncell); for(i=0;i at least stage C. */ + if (dataStageI>2){ velj = malloc(sizeof(*velj)*par.ncell); for(j=0;j2){ + /* Should rather have a vector column? */ + for(n=0;n2){ aCoeffs = malloc(sizeof(*aCoeffs)*totalNumLinks); for(n=0;n: - */ - fits_movnam_hdu(fptr, IMAGE_HDU, extname, 0, &status); - processFitsError(status); + return i-1; +} - readPopsExtFromFits(fptr, (unsigned int)par.ncell, md, si, g); - } - } - } +/*....................................................................*/ +int countKeywords(fitsfile *fptr, char *baseName){ + char kwdName[9]; + int i, status; + char kwdValue[80]; - if(*popsValuesFound==0){ - for(i=0;i1 then it also mallocs 'firstNearNeigh'. If dataStageI>2 then the function also mallocs the following elements of struct grid for all grid points: dens abun Otherwise these are set to NULL. @@ -949,33 +659,30 @@ Otherwise these are set to NULL. If a COLLPARn keywords are found in the GRID extension header then collPartNames is malloc'd to the number of these. */ - LONGLONG numGridCells, firstRow=1, firstElem=1; - int status=0, colNum, anynul=0, i, j, n;//, localNumCollPart; + LONGLONG numGridCells, firstRow=1, firstElem=1, i_LL; + int status=0, colNum, anynul=0, i; char colName[20]; - _Bool nnValuesFound; double modelRadius; - char *genericComment; - char genericKwd[9], message[80]; - const int maxNumCollPart = 9; - char dummyCollPartName[80]; + char genericKwd[9]; + char message[80]; unsigned int *ids=NULL; double *xj=NULL; _Bool *sink=NULL;/* Probably should be char* but this seems to work. */ - unsigned short *numNeigh=NULL; + unsigned short *numNeigh=NULL, i_us; double *velj=NULL, *densn=NULL; float *dopb=NULL, *t=NULL, *abunm=NULL; + /* Go to the GRID extension. + */ + fits_movnam_hdu(fptr, BINARY_TBL, "GRID", 0, &status); + processFitsError(status); + /* Find out how many rows there are, then malloc the array. */ fits_get_num_rowsll(fptr, &numGridCells, &status); processFitsError(status); - if(numGridCells!=par.ncell && !silent){ - sprintf(message, "The inputPars struct says there should be %d grid points, but %u were found.\n", par.ncell, (unsigned int)numGridCells); - warning(message); - } - - *g = malloc(sizeof(**g)*numGridCells); + mallocAndSetDefaultGrid(gp, (unsigned int)numGridCells); /* Read the columns. */ @@ -986,22 +693,34 @@ If a COLLPARn keywords are found in the GRID extension header then collPartNames fits_read_col(fptr, TUINT, colNum, firstRow, firstElem, numGridCells, 0, ids, &anynul, &status); processFitsError(status); - for(i=0;inDims = (unsigned short)countColumns(fptr, "X"); + + /* 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=0;inDims;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(j=0;jnSinkPoints = 0; + for(i_LL=0;i_LLnSinkPoints++; } free(sink); - /* Check for the existence of a NUMNEIGH column: - */ - nnValuesFound = 0; /* Default value. */ - fits_get_colnum(fptr, CASEINSEN, "NUMNEIGH", &colNum, &status); - if(status!=COL_NOT_FOUND){ + gridInfoRead->nInternalPoints = numGridCells - gridInfoRead->nSinkPoints; + + if(dataStageI>1){ + fits_get_colnum(fptr, CASEINSEN, "NUMNEIGH", &colNum, &status); processFitsError(status); - nnValuesFound = 1; numNeigh = malloc(sizeof(*numNeigh)*numGridCells); fits_read_col(fptr, TUSHORT, colNum, firstRow, firstElem, numGridCells, 0, numNeigh, &anynul, &status); processFitsError(status); - for(i=0;i2){ + /* Count the numbers of DENSITYn and ABUNMOLn columns: */ - *userValuesFound = 0; /* Default value. */ - i = 0; - sprintf(colName, "VEL%d", i+1); - fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); - if(status!=COL_NOT_FOUND){ + gridInfoRead->nDensities = (unsigned short)countColumns(fptr, "DENSITY"); + gridInfoRead->nSpecies = (unsigned short)countColumns(fptr, "ABUNMOL"); + + /* Dimension the appropriate elements of gp: + */ + for(i_LL=0;i_LLnDensities); + (*gp)[i_LL].abun = malloc(sizeof(double)*gridInfoRead->nSpecies); + (*gp)[i_LL].nmol = malloc(sizeof(double)*gridInfoRead->nSpecies); + for(i_us=0;i_usnSpecies;i_us++) + (*gp)[i_LL].nmol[i_us] = 0.0; + } + + /* Read the VEL columns: + */ + velj = malloc(sizeof(*velj)*numGridCells); + for(i_us=0;i_usnDims;i_us++){ + sprintf(colName, "VEL%d", (int)i_us+1); + fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); processFitsError(status); - *userValuesFound = 1; + fits_read_col(fptr, TDOUBLE, colNum, firstRow, firstElem, numGridCells, 0, velj, &anynul, &status); + processFitsError(status); - /* Dimension the appropriate elements of g: - */ - for(i=0;inDensities;i_us++){ + sprintf(colName, "DENSITY%d", (int)i_us+1); + fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); processFitsError(status); - dopb = malloc(sizeof(*dopb)*numGridCells); - fits_read_col(fptr, TFLOAT, colNum, firstRow, firstElem, numGridCells, 0, dopb, &anynul, &status); + fits_read_col(fptr, TDOUBLE, colNum, firstRow, firstElem, numGridCells, 0, densn, &anynul, &status); processFitsError(status); - for(i=0;inSpecies;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, t, &anynul, &status); + fits_read_col(fptr, TFLOAT, colNum, firstRow, firstElem, numGridCells, 0, abunm, &anynul, &status); processFitsError(status); - for(i=0;imaxNumCollPart){ - if(!silent){ - sprintf(message, "There seem to be %d collision partners but keywords can only be written for %d.", par.collPart, maxNumCollPart); - warning(message); - } - *numCollPartRead = maxNumCollPart; - }else{ - *numCollPartRead = par.collPart; - } - + *numCollPartRead = countKeywords(fptr, "COLLPAR"); + if(*numCollPartRead <= 0){ + *collPartNames = NULL; + }else{ *collPartNames = malloc(sizeof(**collPartNames)*(*numCollPartRead)); - for(i=0;i<(*numCollPartRead);i++){ sprintf(genericKwd, "COLLPAR%d", i+1); (*collPartNames)[i] = malloc(sizeof(char)*100); @@ -1199,27 +893,41 @@ If a COLLPARn keywords are found in the GRID extension header then collPartNames } /*....................................................................*/ -void readLinksExtFromFits(fitsfile *fptr, _Bool userValuesFound\ - , unsigned short numACoeffs, struct grid *g, struct linkType **links){ +void readLinksExtFromFits(fitsfile *fptr, struct gridInfoType *gridInfoRead\ + , struct grid *gp, struct linkType **links, const int dataStageI){ /* -See the comment at the beginning of the present module for a description of how the LINKS extension relates to the grid struct. +See the comment at the beginning of gridio.c for a description of how the LINKS extension relates to the grid struct. -The present function mallocs *links. +The present function mallocs the pointer *links. */ - LONGLONG totalNumLinks, firstRow=1, firstElem=1; - int status=0, colNum, anynul=0, li, n; + LONGLONG totalNumLinks, firstRow=1, firstElem=1, i_LL; + int status=0, colNum, anynul=0; char colName[21]; - unsigned int *ids=NULL; + unsigned int *ids=NULL, totalNumGridPoints, ppi; double *aCoeffs=NULL; + char message[80]; + unsigned short i_s; + + 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); + 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; *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, ppi, totalNumGridPoints); + bail_out(message); + } + exit(1); + } + (*links)[i_LL].g[0] = &gp[ppi]; } /* Read GRID_I_2 column. @@ -1241,56 +957,74 @@ The present function mallocs *links. fits_read_col(fptr, TUINT, colNum, firstRow, firstElem, totalNumLinks, 0, ids, &anynul, &status); processFitsError(status); - for(li=0;li=totalNumGridPoints){ + if(!silent){ + sprintf(message, "GRID_I_2 %dth-row value %ud is outside range [0,%ud]", (int)i_LL, ppi, totalNumGridPoints); + bail_out(message); + } + exit(1); + } + (*links)[i_LL].g[1] = &gp[ppi]; } free(ids); - if(userValuesFound>0){ /* Means we are in at least stage C. */ + if(dataStageI>2){ /* Means we are in at least stage C. */ + /* Find out how many ACOEFF_* columns there are. + */ + gridInfoRead->nACoeffs = (unsigned short)countColumns(fptr, "ACOEFF_"); + + for(i_LL=0;i_LLnACoeffs); + aCoeffs = malloc(sizeof(*aCoeffs)*totalNumLinks); - for(n=0;nnACoeffs;i_s++){ /* Read the ACOEFF_n columns. */ - sprintf(colName, "ACOEFF_%d", n+1); + sprintf(colName, "ACOEFF_%d", (int)i_s+1); fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); processFitsError(status); fits_read_col(fptr, TDOUBLE, colNum, firstRow, firstElem, totalNumLinks, 0, aCoeffs, &anynul, &status); processFitsError(status); - for(li=0;linACoeffs = 0; + for(i_LL=0;i_LLnNNIndices = (unsigned int)totalNumNeigh; *nnLinks = malloc(sizeof(**nnLinks)*totalNumNeigh); /* Read LINK_I column. @@ -1302,94 +1036,111 @@ The function mallocs *nnLinks. fits_read_col(fptr, TUINT, colNum, firstRow, firstElem, totalNumNeigh, 0, linkIs, &anynul, &status); processFitsError(status); - for(i=0;ig[0]->id==i){ - (*g)[i].neigh[j] = linkPtr->g[1]; - }else{ - (*g)[i].neigh[j] = linkPtr->g[0]; - } +_Bool checkPopsFitsExtExists(fitsfile *fptr, const unsigned short speciesI){ + const unsigned short maxNumSpecies = 9; + char message[80]; + char extname[13]; + int status=0; - (*g)[i].a0[j] = linkPtr->aCoeffs[0]; - (*g)[i].a1[j] = linkPtr->aCoeffs[1]; - (*g)[i].a2[j] = linkPtr->aCoeffs[2]; - (*g)[i].a3[j] = linkPtr->aCoeffs[3]; - (*g)[i].a4[j] = linkPtr->aCoeffs[4]; + 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); + 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, unsigned int numGridPoints\ - , molData *md, unsigned short speciesI, struct grid **g){ +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 (*g)[yi].mol[speciesI].pops for each grid point yi and species. +The function mallocs gp[yi].mol[speciesI].pops for each grid point yi and species. */ const int maxLenMolName = 8; - int bitpix, naxis, status=0, xi, yi, anynul=0; - long *naxes=NULL; + 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]; - unsigned int numEnergyLevels = md[speciesI].nlev; + char message[80]; + char extname[13]; + unsigned int numGridPoints, i_u; + + 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); - /****** error if naxes != {numEnergyLevels, numGridPoints} or bitpix!=FLOAT_IMG.*/ - free(naxes); - row = malloc(sizeof(*row)*numEnergyLevels); + 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(yi=0;yimols[speciesI].nLevels; + lpixels[1]=(int)i_u+1; fits_read_subset(fptr, TFLOAT, fpixels, lpixels, inc, 0, row, &anynul, &status); processFitsError(status); - (*g)[yi].mol[speciesI].pops = malloc(sizeof(double)*numEnergyLevels); - for(xi=0;ximols[speciesI].nLevels); + for(xi=0;ximols[speciesI].nLevels;xi++) + gp[i_u].mol[speciesI].pops[xi] = (double)row[xi]; } free(row); @@ -1397,9 +1148,8 @@ The function mallocs (*g)[yi].mol[speciesI].pops for each grid point yi and spec /* Read kwds: */ fits_read_key(fptr, TSTRING, "MOL_NAME", molNameRead, NULL, &status); + gridInfoRead->mols[speciesI].molName = molNameRead;//*****?? processFitsError(status); - /****** raise exception if it is not the same as md[speciesI].molName. */ } - diff --git a/src/gridio.c b/src/gridio.c new file mode 100644 index 0000000..30620ea --- /dev/null +++ b/src/gridio.c @@ -0,0 +1,798 @@ +/* + * gridio.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" + +/* +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 an attempt to regulate this, the following five states of completeness (encoded in values of the variable 'dataStageI') have been defined: + + dataStageI==0: No useable data in grid. + + dataStageI==1. At this stage the vector of grid objects has been malloc'd and values have been generated for the following struct elements: + id + x + sink + + dataStageI==2. This stage is entered after the Delaunay neighbours of each grid point have been determined. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: + numNeigh + neigh + + dataStageI==3. This is entered after sampling the user-supplied functions for density, velocity etc. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: + vel + a0, a1 etc. + dens + t + dopb + abun +Note that it is necessary to make this stage dependent on the previous because we need information about the near-neighbours to calculate the a0, a1 etc coefficients. + + dataStageI==4. After 1 or more iterations of populating the levels, we enter the final stage, in which all the grid struct elements have been malloc'd and given at least preliminary values; specifically now the element + mol + + *NOTE* that 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;jwriteGridAtStage[par->dataStageI-1]){ + status = writeGrid(par->gridOutFiles[par->dataStageI-1], fileFormatI\ + , *par, DIM, NUM_VEL_COEFFS, gp, md, collPartNames, par->dataStageI); + + if(status){ + sprintf(message, "writeGrid at data stage %d returned with status %d", par->dataStageI, status); + if(!silent) bail_out(message); + exit(1); + } + } + + free(collPartNames); +} + +/*....................................................................*/ +lime_fptr *openFileForRead(char *inFileName, const int fileFormatI, int *dataStageI){ + lime_fptr *fptr=NULL; + + if(fileFormatI==lime_FITS){ + fptr = openFITSFileForRead(inFileName, dataStageI); + }else{ + return NULL; + } + + return fptr; +} + +/*....................................................................*/ +lime_fptr *openFileForWrite(char *outFileName, const int fileFormatI, const int dataStageI){ + lime_fptr *fptr=NULL; + + if(fileFormatI==lime_FITS){ + fptr = openFITSFileForWrite(outFileName, dataStageI); + }else{ + return NULL; + } + + return fptr; +} + +/*....................................................................*/ +void closeAndFree(lime_fptr *fptr, const int fileFormatI\ + , unsigned int *firstNearNeigh, struct linkType **nnLinks\ + , struct linkType *links, const unsigned int totalNumLinks){ + + unsigned int li; + + closeFile(fptr, fileFormatI); + free(firstNearNeigh); + free(nnLinks); + if(links!=NULL){ + for(li=0;linumNeigh;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].g[0] = gAPtr; + (*links)[li].g[1] = gBPtr; + + li++; + } + + nnLinkIs[ni] = linkId; + ni++; + } + + pointIsDone[idA] = 1; + } + *totalNumLinks = li; + + *links = realloc(*links, sizeof(struct linkType)*(*totalNumLinks)); + + if(dataStageI>2){ + for(li=0;li<*totalNumLinks;li++) + (*links)[li].aCoeffs = malloc(sizeof(double)*NUM_VEL_COEFFS); + + for(li=0;li<*totalNumLinks;li++){ + if((*links)[li].g[0]->sink && (*links)[li].g[1]->sink){ + for(ci=0;cisink){ /* If we get to here, this ensures that g[1] is not a sink point. */ + nearI = 1; + evenCoeffSign = -1.0; + }else{ + nearI = 0; + evenCoeffSign = 1.0; + } + + gAPtr = (*links)[li].g[nearI]; + /* Find which neighbour of gAPtr corresponds to the link: */ + linkNotFound = 1; + for(jA=0;jAnumNeigh;jA++){ + idB = gAPtr->neigh[jA]->id; + if(linkNotFound && idB==(*links)[li].g[1-nearI]->id){ + (*links)[li].aCoeffs[0] = evenCoeffSign*gAPtr->a0[jA]; + (*links)[li].aCoeffs[1] = gAPtr->a1[jA]; + (*links)[li].aCoeffs[2] = evenCoeffSign*gAPtr->a2[jA]; + (*links)[li].aCoeffs[3] = gAPtr->a3[jA]; + (*links)[li].aCoeffs[4] = evenCoeffSign*gAPtr->a4[jA]; + linkNotFound = 0; + } + } + if(linkNotFound){ + sprintf(message, "SNAFU with link %d grid point %d", li, gAPtr->id); + if(!silent) bail_out(message); + exit(1); + } + } /* end if not both sink */ + } /* end loop over links */ + }else{ /* dataStageI<=2 */ + for(li=0;li<*totalNumLinks;li++) + (*links)[li].aCoeffs = NULL; + } + + *nnLinks = malloc(sizeof(**nnLinks)*(*totalNumNeigh)); + for(ni=0;ni<(*totalNumNeigh);ni++) + (*nnLinks)[ni] = &(*links)[nnLinkIs[ni]]; + + free(nnLinkIs); + free(pointIsDone); + /* The calling routine must free firstNearNeigh, nnLinks, links. */ +} + + +/*....................................................................*/ +int writeGrid(char *outFileName, const int fileFormatI, inputPars par\ + , unsigned short numDims, unsigned short numACoeffs, struct grid *g\ + , molData *md, char **collPartNames, const int dataStageI){ + + lime_fptr *fptr=NULL; + int status = 0; + unsigned short speciesI; + char message[80]; + struct linkType *links=NULL, **nnLinks=NULL; + unsigned int totalNumLinks, totalNumNeigh, *firstNearNeigh=NULL; + + if (outFileName==""){ + if(!silent) warning("Cannot write grid list to file, filename is blank."); + return 1; + } + + if (g==NULL || dataStageI<1){ + if(!silent) warning("Cannot write grid list to file, there are no entries in it."); + return 2; + } + + sprintf(message, "Writing grid-point list to file %s", outFileName); + if(!silent) printMessage(message); + + if (dataStageI>1){ + constructLinkArrays((unsigned int)par.ncell, g, &links, &totalNumLinks\ + , &nnLinks, &firstNearNeigh, &totalNumNeigh, dataStageI); + } + + fptr = openFileForWrite(outFileName, fileFormatI, dataStageI); + if(fptr==NULL){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); + return 3; + } + + status = writeGridBlock(fptr, fileFormatI, par, numDims, g, firstNearNeigh, collPartNames, dataStageI); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); + return 4; + } + + if (dataStageI>1){ + status = writeNnIndicesBlock(fptr, fileFormatI, totalNumNeigh, nnLinks, links); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); + return 5; + } + + status = writeLinksBlock(fptr, fileFormatI, totalNumLinks, numACoeffs, links, dataStageI); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); + return 6; + } + + if (dataStageI>3){ + for(speciesI=0;speciesINUM_GRID_STAGES){ + if(!silent) bail_out("Data stage out of range."); + exit(1); + } + + sprintf(message, "Grid-point data stage determined as %d", *dataStageI); + if(!silent) printMessage(message); + + /* Read the values which should be in grid for every stage. + */ + status = readGridBlock(fptr, fileFormatI, gridInfoRead, gp, &firstNearNeigh\ + , collPartNames, numCollPartRead, *dataStageI); + + if(status || gridInfoRead->nSinkPoints<=0\ + || gridInfoRead->nInternalPoints<=0 || gridInfoRead->nDims<=0){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, NULL, NULL, 0); + freeReadGrid(gp, *gridInfoRead); + + if(status) + return 1; + else if(gridInfoRead->nSinkPoints<=0) + return 2; + else if(gridInfoRead->nInternalPoints<=0) + return 3; + else if(gridInfoRead->nDims<=0) + return 4; + else{ + sprintf(message, "This indicates a programming error. Please contact the developer."); + if(!silent) bail_out(message); + exit(1); + } + } + + if((*dataStageI)>2 && (gridInfoRead->nDensities<=0 || gridInfoRead->nSpecies<=0)){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, NULL, NULL, 0); + freeReadGrid(gp, *gridInfoRead); + + if(gridInfoRead->nDensities<=0) + return 5; + else if(gridInfoRead->nSpecies<=0) //***** what if all continuum images?? + return 6; + else{ + sprintf(message, "This indicates a programming error. Please contact the developer."); + if(!silent) bail_out(message); + exit(1); + } + } + + totalNumGridPoints = gridInfoRead->nInternalPoints+gridInfoRead->nSinkPoints; + + if((*dataStageI)>1){ /* there should be blocks for the nnIndices and links. */ + status = readLinksBlock(fptr, fileFormatI, gridInfoRead, *gp, &links, *dataStageI); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, NULL, links, gridInfoRead->nLinks); + freeReadGrid(gp, *gridInfoRead); + return 7; + } + + if ((*dataStageI)>2 && gridInfoRead->nACoeffs<=0){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(gp, *gridInfoRead); + return 8; + } + + status = readNnIndicesBlock(fptr, fileFormatI, links, &nnLinks, gridInfoRead); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(gp, *gridInfoRead); + return 9; + } + + /* Convert the NN information back to the standard LIME grid struct format. + */ + loadNnIntoGrid(firstNearNeigh, nnLinks, links, *gridInfoRead, *gp, *dataStageI); + + if ((*dataStageI)>3){ /* there should be pops blocks. */ + status = getNumPopsBlocks(fptr, fileFormatI, &numBlocks); + + if(status || numBlocks<=0 || numBlocks!=gridInfoRead->nSpecies){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(gp, *gridInfoRead); + if(status) + return 10; + else if(numBlocks<=0) + return 11; + else if(numBlocks!=gridInfoRead->nSpecies) + return 12; + else{ + sprintf(message, "This indicates a programming error. Please contact the developer."); + if(!silent) bail_out(message); + exit(1); + } + } + + gridInfoRead->mols = malloc(sizeof(struct molInfoType)*gridInfoRead->nSpecies); + + for(i_u=0;i_unSpecies); + for(i_s=0;i_snSpecies;i_s++){ + (*gp)[i_u].mol[i_s].dust = NULL; + (*gp)[i_u].mol[i_s].knu = NULL; + (*gp)[i_u].mol[i_s].pops = NULL; + (*gp)[i_u].mol[i_s].partner = NULL; + } + } + + for(i_s=0;i_snSpecies;i_s++){ + gridInfoRead->mols[i_s].molName = NULL; + gridInfoRead->mols[i_s].nLevels = -1; + gridInfoRead->mols[i_s].nLines = -1; + + status = readPopsBlock(fptr, fileFormatI, i_s, *gp, gridInfoRead); + + if(status || gridInfoRead->mols[i_s].nLevels<=0){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(gp, *gridInfoRead); + if(status) + return 13; + else if(gridInfoRead->mols[i_s].nLevels<=0) + return 14; + else{ + sprintf(message, "This indicates a programming error. Please contact the developer."); + if(!silent) bail_out(message); + exit(1); + } + } + } + } + } + + /* Set unread pointers to NULL: + */ + if((*dataStageI)<4){ + for(i_u=0;i_unLinks); + return 0; +} + +/*....................................................................*/ +int readGridBlock(lime_fptr *fptr, const int fileFormatI\ + , struct gridInfoType *gridInfoRead, struct grid **gp\ + , unsigned int **firstNearNeigh, char ***collPartNames\ + , int *numCollPartRead, const int dataStageI){ + + int status=0; + + if(fileFormatI==lime_FITS){ + readGridExtFromFits(fptr, gridInfoRead, gp, firstNearNeigh\ + , collPartNames, numCollPartRead, dataStageI); + }else{ + status = 1; + } + + return status; +} + +/*....................................................................*/ +int readLinksBlock(lime_fptr *fptr, const int fileFormatI\ + , struct gridInfoType *gridInfoRead, struct grid *g\ + , struct linkType **links, const int dataStageI){ + + int status=0; + + if(fileFormatI==lime_FITS){ + readLinksExtFromFits(fptr, gridInfoRead, g, links, dataStageI); + }else{ + status = 1; + } + + return status; +} + +/*....................................................................*/ +int readNnIndicesBlock(lime_fptr *fptr, const int fileFormatI, struct linkType *links\ + , struct linkType ***nnLinks, struct gridInfoType *gridInfoRead){ + + int status=0; + + if(fileFormatI==lime_FITS){ + readNnIndicesExtFromFits(fptr, links, nnLinks, gridInfoRead); + }else{ + status = 1; + } + + return status; +} + +/*....................................................................*/ +void loadNnIntoGrid(unsigned int *firstNearNeigh, struct linkType **nnLinks\ + , struct linkType *links, struct gridInfoType gridInfoRead, struct grid *gp\ + , const int dataStageI){ + /* +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 the following extensions of struct g for each grid point: + neigh +if dataStageI>2: + a0 + a1 + a2 + a3 + a4 + +The business here with the evenCoeffSign is due to the fact that, for each grid point A, the 'a' coefficients are calculated for a polynomial with independent coordinate being the distance from that grid point towards the appropriate neighbour B. Point B will also store its own 'a' coefficients, which produce identically the same polynomial, but which in general will have different values, because in this case the coefficients are calculated appropriate to the distance from B towards A. Now that the definition of the independent coordinate has been defined as the fractional distance minus 0.5, thus symmetrical about the midpoint between A and B, the coefficients for B can be obtained from those of A by inverting the sign of the even-power coefficients. Note that we invert the even coefficients, not the odd, because not only has the sign of the independent coordinate (i.e the fractional distance along the link minus 0.5) reversed, but also the that of the function being interpolated (the component of velocity in the direction of the neighbour point). + +Further: the 'a' coefficients for the link or edge between A and B are stored in a struct of type linkType. They are stored appropriate to point g[0] of that struct. We cycle through all the grid points in the present function and, for each point/neighbour pair, we look first to see which end of the matching link element 'point' is - i.e., whether it is g[0] or g[1] of the link. If it is g[0], then the 'a' coefficients stored in link are appropriate for it; if g[1], then the even-power ones have to be inverted in sign. + */ + + const unsigned int totalNumGridPoints = gridInfoRead.nInternalPoints+gridInfoRead.nSinkPoints; + unsigned int i_u; + int j; + struct linkType *linkPtr; + char message[80]; + double evenCoeffSign; + + for(i_u=0;i_ug[0]->id==(int)i_u){ + gp[i_u].neigh[j] = linkPtr->g[1]; + }else{ + gp[i_u].neigh[j] = linkPtr->g[0]; + } + } + } + + if(dataStageI>2){ + /* Just for the time being: */ + if(gridInfoRead.nACoeffs!=5){ + if(!silent){ + sprintf(message, "There should be %d ACOEFF_n columns, but %d were read.", NUM_VEL_COEFFS, (int)gridInfoRead.nACoeffs); + bail_out(message); + } + exit(1); + } + + for(i_u=0;i_ug[0]->id==(int)i_u){ + evenCoeffSign = 1.0; + }else{ + evenCoeffSign = -1.0; + } + + gp[i_u].a0[j] = evenCoeffSign*linkPtr->aCoeffs[0]; + gp[i_u].a1[j] = linkPtr->aCoeffs[1]; + gp[i_u].a2[j] = evenCoeffSign*linkPtr->aCoeffs[2]; + gp[i_u].a3[j] = linkPtr->aCoeffs[3]; + gp[i_u].a4[j] = evenCoeffSign*linkPtr->aCoeffs[4]; + } + } + } +} + +/*....................................................................*/ +int getNumPopsBlocks(lime_fptr *fptr, const int fileFormatI, unsigned short *numBlocks){ + int status = 0; + _Bool blockFound = 1; + + *numBlocks = 0; + while(blockFound && !status){ + status = checkPopsBlockExists(fptr, fileFormatI, *numBlocks, &blockFound); + (*numBlocks)++; + } + + (*numBlocks)--; + + return status; +} + +/*....................................................................*/ +int checkPopsBlockExists(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 readPopsBlock(lime_fptr *fptr, const int fileFormatI\ + , const unsigned short speciesI, struct grid *g\ + , struct gridInfoType *gridInfoRead){ + + int status=0; + + if(fileFormatI==lime_FITS){ + readPopsExtFromFits(fptr, speciesI, g, gridInfoRead); + }else{ + status = 1; + } + + return status; +} + + diff --git a/src/lime.h b/src/lime.h index 426587a..e5c94a2 100644 --- a/src/lime.h +++ b/src/lime.h @@ -79,7 +79,6 @@ #define FAST_EXP_NUM_BITS 8 #define NUM_GRID_STAGES 4 - /* input parameters */ typedef struct { double radius,radiusSqu,minScale,minScaleSqu,tcmb,taylorCutoff; @@ -92,8 +91,9 @@ typedef struct { int sampling,collPart,lte_only,init_lte,antialias,polarization,doPregrid,nThreads; char **moldatfile; _Bool writeGridAtStage[NUM_GRID_STAGES]; - char *gridFitsOutSets[NUM_GRID_STAGES]; - int nSolveIters; + char *gridInFile,*gridOutFiles[NUM_GRID_STAGES]; + + int dataStageI,nSolveIters; } inputPars; /* Molecular data: shared attributes */ @@ -208,86 +208,106 @@ void gasIIdust(double,double,double,double *); /* More functions */ -void binpopsout(inputPars *, struct grid *, molData *); -void buildGrid(inputPars *, struct grid *, molData *); +void binpopsout(inputPars*, struct grid*, molData*); void calcAvRelLineAmp(struct grid*, int, int, double, double, double*); void calcAvRelLineAmp_lin(struct grid*, int, int, double, double, double*); void calcInterpCoeffs(inputPars*, struct grid*); void calcInterpCoeffs_lin(inputPars*, struct grid*); -void calcSourceFn(double dTau, const inputPars *par, double *remnantSnu, double *expDTau); -void constructLinkArrays(unsigned int, struct grid *, struct linkType **\ - , unsigned int *, struct linkType ***, unsigned int **, unsigned int *); +void calcSourceFn(double, const inputPars*, double*, double*); +int checkPopsBlockExists(lime_fptr*, const int, const unsigned short, _Bool*); +_Bool checkPopsFitsExtExists(lime_fptr*, const unsigned short); +void closeAndFree(lime_fptr*, const int, unsigned int*, struct linkType**, struct linkType*, const unsigned int); +void closeFile(lime_fptr*, const int); +void closeFITSFile(fitsfile*); +void constructLinkArrays(const unsigned int, struct grid*, struct linkType**, unsigned int*, struct linkType***, unsigned int**, unsigned int*, const int); void continuumSetup(int, image *, molData *, inputPars *, struct grid *); -void defineGridExtColumns(unsigned short, unsigned short\ - , unsigned short, unsigned short, char *ttype[]\ - , char *tform[], char *tunit[], int dataTypes[]); +int countColumns(fitsfile*, char*); +int countKeywords(fitsfile*, char*); +void defineGridExtColumns(const unsigned short, inputPars, const unsigned short, const int, char *ttype[], char *tform[], char *tunit[], int dataTypes[]); void distCalc(inputPars *, struct grid *); void fit_d1fi(double, double, double*); void fit_fi(double, double, double*); void fit_rr(double, double, double*); +void freeReadGrid(struct grid**, struct gridInfoType); 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); double gaussline(double, double); void getArea(inputPars *, struct grid *, const gsl_rng *); -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 *); -void gridAlloc(inputPars *, struct grid **); +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*); +int getNumPopsBlocks(lime_fptr*, const int, unsigned short*); +void gridAlloc(inputPars*, struct grid**); void input(inputPars *, image *); double interpolate(double, double, double, double, double, double); 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); -void lineBlend(molData *, inputPars *, blend **); -void lineCount(int,molData *,int **, int **, int *); -void loadNnIntoGrid(unsigned int *, struct linkType **, struct linkType *, unsigned int, struct grid **); -void LTE(inputPars *, struct grid *, molData *); -void molinit(molData *, inputPars *, struct grid *,int); -void openSocket(inputPars *par, int); -void qhull(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); -int pointEvaluation(inputPars *,double, double, double, double); -void popsin(inputPars *, struct grid **, molData **, int *); -void popsout(inputPars *, struct grid *, molData *); -void predefinedGrid(inputPars *, struct grid *); +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); +void lineBlend(molData*, inputPars*, blend**); +void lineCount(int,molData*,int**, int**, int*); +void loadNnIntoGrid(unsigned int*, struct linkType**, struct linkType*, struct gridInfoType, struct grid*, const int); +void LTE(inputPars*, struct grid*, molData*); +void mallocAndSetDefaultGrid(struct grid**, const unsigned int); +void molinit(molData*, inputPars*, struct grid*,int); +lime_fptr *openFileForRead(char*, const int, int*); +fitsfile *openFITSFileForRead(char*, int*); +lime_fptr *openFileForWrite(char*, const int, const int); +fitsfile *openFITSFileForWrite(char*, const int); +void openSocket(inputPars*, int); +void qhull(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); +int pointEvaluation(inputPars*,double, double, double, double); +void popsin(inputPars*, struct grid**, molData**, int *); +void popsout(inputPars*, struct grid*, molData*); +void predefinedGrid(inputPars*, struct grid*); void processFitsError(int); -double ratranInput(char *, char *, double, double, double); -void raytrace(int, inputPars *, struct grid *, molData *, image *); -void readGridExtFromFits(fitsfile *, inputPars, unsigned short\ - , struct grid **, unsigned int **, _Bool *, char ***, int *); -void readGridFromFits(char *, inputPars, unsigned short, unsigned short\ - , struct grid **, _Bool *, _Bool *, _Bool *, molData *, char ***, int *); -void readLinksExtFromFits(fitsfile *, _Bool, unsigned short, struct grid *, struct linkType **); -void readNnIndicesExtFromFits(fitsfile *, struct linkType *, struct linkType ***); -void readPopsExtFromFits(fitsfile *, unsigned int, molData *, unsigned short, struct grid **); -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_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 traceray(rayData, int, int, inputPars *, struct grid *, molData *, image *, int, int *, int *, double); -double veloproject(double *, double *); -void writeGridExtToFits(fitsfile *, inputPars, unsigned short, struct grid *, unsigned int *, char **); -int writeGridToFits(char *, inputPars, unsigned short, unsigned short, struct grid *, molData *, char **); -void writefits(int, inputPars *, molData *, image *); -void writeLinksExtToFits(fitsfile *, _Bool, unsigned int, unsigned short, struct linkType *); -void writeNnIndicesExtToFits(fitsfile *, unsigned int, struct linkType **, struct linkType *); -void writePopsExtToFits(fitsfile *, unsigned int, molData *, unsigned short, struct grid *); -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); +double ratranInput(char*, char*, double, double, double); +void raytrace(int, inputPars*, struct grid*, molData*, image*); +int readGrid(char*, const int, struct gridInfoType*, struct grid**, char***, int*, int*); +int readGridBlock(lime_fptr*, const int, struct gridInfoType*, struct grid**, unsigned int**, char***, int*, const int); +void readGridExtFromFits(fitsfile*, struct gridInfoType*, struct grid**, unsigned int**, char***, int*, const int); +void readOrBuildGrid(inputPars*, struct grid**); +int readLinksBlock(lime_fptr*, const int, struct gridInfoType*, struct grid*, struct linkType**, const int); +void readLinksExtFromFits(fitsfile*, struct gridInfoType*, struct grid*, struct linkType**, const int); +int readNnIndicesBlock(lime_fptr*, const int, struct linkType*, struct linkType***, struct gridInfoType*); +void readNnIndicesExtFromFits(fitsfile*, struct linkType*, struct linkType***, struct gridInfoType*); +void readOrBuildGrid(inputPars*, struct grid**); +int readPopsBlock(lime_fptr*, const int, const unsigned short, struct grid*, struct gridInfoType*); +void readPopsExtFromFits(fitsfile*, const unsigned short, struct grid*, struct gridInfoType*); +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_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 traceray(rayData, int, int, inputPars*, struct grid*, molData*, image*, int, int*, int*, double); +void velocityspline2(double*, double*, double, double, double, double*); +double veloproject(double*, double*); +int writeGrid(char*, const int, inputPars, unsigned short, unsigned short, struct grid*, molData*, char**, const int); +int writeGridBlock(lime_fptr*, const int, inputPars, unsigned short, struct grid*, unsigned int*, char**, const int); +void writeGridExtToFits(fitsfile*, inputPars, unsigned short, struct grid*, unsigned int*, char**, const int); +void writeGridIfRequired(inputPars*, struct grid*, molData*, const int); +void writefits(int, inputPars*, molData*, image*); +int writeLinksBlock(lime_fptr*, const int, const unsigned int, const unsigned short, struct linkType*, const int); +void writeLinksExtToFits(fitsfile*, const unsigned int, const unsigned short, struct linkType*, const int); +int writeNnIndicesBlock(lime_fptr*, const int, const unsigned int, struct linkType**, struct linkType*); +void writeNnIndicesExtToFits(fitsfile*, const unsigned int, struct linkType**, struct linkType*); +int writePopsBlock(lime_fptr*, const int, unsigned int, molData*, unsigned short, struct grid*); +void writePopsExtToFits(fitsfile*, const unsigned int, molData*, const unsigned short, struct grid*); +void write_VTK_unstructured_Points(inputPars*, struct grid*); +int factorial(const int); +double taylor(const int, const float); +void calcFastExpRange(const int, const int, int*, int*, int*); +void calcTableEntries(const int, const int); double FastExp(const float negarg); diff --git a/src/main.c b/src/main.c index f8bbc7d..ff5e8ed 100644 --- a/src/main.c +++ b/src/main.c @@ -21,14 +21,13 @@ double EXP_TABLE_3D[1][1][1]; #endif int main () { - int i,stageI,status=0; + int i,nLineImages; int initime=time(0); int popsdone=0; - molData* m = NULL; + molData* md = NULL; inputPars par; - struct grid* g = NULL; + struct grid* gp = NULL; image* img = NULL; - char **collPartNames=NULL; /*** this is a placeholder until we start reading these. */ char message[80]; if(!silent) greetings(); @@ -38,7 +37,7 @@ int main () { calcTableEntries(FAST_EXP_MAX_TAYLOR, FAST_EXP_NUM_BITS); #endif - parseInput(&par,&img,&m); + parseInput(&par,&img,&md); if(!silent && par.nThreads>1){ sprintf(message, "Number of threads used: %d", par.nThreads); @@ -47,39 +46,58 @@ int main () { if(par.doPregrid) { - gridAlloc(&par,&g); - predefinedGrid(&par,g); + gridAlloc(&par,&gp); + predefinedGrid(&par,gp); + par.dataStageI = 3; // Sort of. } else if(par.restart) { - popsin(&par,&g,&m,&popsdone); + popsin(&par,&gp,&md,&popsdone); + par.dataStageI = 4; // Sort of. } else { - gridAlloc(&par,&g); - buildGrid(&par,g,m); + readOrBuildGrid(&par,&gp); } + /* Make all the continuum images, and count the non-continuum images at the same time: + */ + nLineImages = 0; for(i=0;i0 && !popsdone){ // eventually, replace !popdone by (!popsdone || dataStageI<4)? *Really* eventually we want to get rid of popsdone. + levelPops(md,&par,gp,&popsdone); + par.dataStageI = 4; +/* Disable the next lines for now, since we have not tested dataStageI<4 in the 'if' of this block, because we can't use an input grid file at dataStageI==4 yet: we have to disentangle all the functionality of molinit() before we can contemplate doing that. + }else if(par.dataStageI==4 && par->nSolveIters>0 && par.writeGridAtStage[par.dataStageI-1]){ + sprintf(message, "You just read a grid file at data stage %d, now you want to write it again?", par.dataStageI); + if(!silent) warning(message); +*/ } - stageI = 3; - if(par.writeGridAtStage[stageI]) - status = writeGridToFits(par.gridFitsOutSets[stageI], par, (unsigned short)DIM\ - , (unsigned short)NUM_VEL_COEFFS, g, m, collPartNames); + if(par.dataStageI==4) + writeGridIfRequired(&par, gp, md, lime_FITS); + + /* Now make the line images. + */ + for(i=0;i Date: Wed, 13 Jul 2016 16:04:21 +0200 Subject: [PATCH 09/13] Added a couple of comments. --- src/grid.c | 6 +++--- src/gridio.c | 1 + 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/grid.c b/src/grid.c index ba959a2..6e9ee03 100644 --- a/src/grid.c +++ b/src/grid.c @@ -709,7 +709,7 @@ void readOrBuildGrid(inputPars *par, struct grid **gp){ if(!silent) warning(message); } - if(par->dataStageI==1) /* Only happens if (i) we read no file, or (ii) we read a file at dataStageI==1. */ + if(par->dataStageI==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); if(par->dataStageI<2){ @@ -721,7 +721,7 @@ void readOrBuildGrid(inputPars *par, struct grid **gp){ if(!silent) warning(message); } - if(par->dataStageI==2) /* Only happens if (i) we read no file, or (ii) we read a file at dataStageI==2. */ + if(par->dataStageI==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); distCalc(par, *gp); /* Mallocs and sets .dir & .ds, sets .nphot */ @@ -758,7 +758,7 @@ void readOrBuildGrid(inputPars *par, struct grid **gp){ if(!silent) warning(message); } - if(par->dataStageI==3) /* Only happens if (i) we read no file, or (ii) we read a file at dataStageI==3. */ + if(par->dataStageI==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); dumpGrid(par,*gp); diff --git a/src/gridio.c b/src/gridio.c index 30620ea..f30bc16 100644 --- a/src/gridio.c +++ b/src/gridio.c @@ -444,6 +444,7 @@ int readGrid(char *inFileName, const int fileFormatI\ sprintf(message, "Reading grid-point list from file %s", inFileName); if(!silent) printMessage(message); + /* Open the file and also return the data stage. */ fptr = openFileForRead(inFileName, fileFormatI, dataStageI); if(*dataStageI<1 || *dataStageI>NUM_GRID_STAGES){ From 053147b30b67d4607d714f0102dd871039a1c24a Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 29 Jul 2016 11:02:44 +0200 Subject: [PATCH 10/13] Added some convenience bit-mask functions. Also added data-stage flag macros. None of this stuff is used yet, I am just clearing the decks before the big change. --- .gitignore | 4 +++- src/aux.c | 34 ++++++++++++++++++++++++++++++++++ src/lime.h | 32 ++++++++++++++++++++++++++++++++ 3 files changed, 69 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 182c93e..091b3ea 100644 --- a/.gitignore +++ b/.gitignore @@ -7,4 +7,6 @@ # Ignore example files example/grid.vtk example/*.pop -example/*.fits +example/*.fits* +example/lime_*.x + diff --git a/src/aux.c b/src/aux.c index 8d92a39..b19e8af 100644 --- a/src/aux.c +++ b/src/aux.c @@ -636,5 +636,39 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ *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/lime.h b/src/lime.h index e5c94a2..3482b11 100644 --- a/src/lime.h +++ b/src/lime.h @@ -79,6 +79,33 @@ #define FAST_EXP_NUM_BITS 8 #define NUM_GRID_STAGES 4 +/* 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_ACOEFF 7 /* a0, a1, a2, a3, a4 */ +#define DS_bit_populations 8 /* mol */ + +#define DS_mask_x (1 << DS_bit_x) +#define DS_mask_neighbours (1 << DS_bit_neighbours) | DS_mask_x +#define DS_mask_velocity (1 << DS_bit_velocity) | DS_mask_x +#define DS_mask_density (1 << DS_bit_density) | DS_mask_x +#define DS_mask_abundance (1 << DS_bit_abundance) | DS_mask_x +#define DS_mask_turb_doppler (1 << DS_bit_turb_doppler) | DS_mask_x +#define DS_mask_temperatures (1 << DS_bit_temperatures) | DS_mask_x +#define DS_mask_ACOEFF (1 << DS_bit_ACOEFF) | DS_mask_neighbours | DS_mask_velocity + +#define DS_mask_1 DS_mask_x +#define DS_mask_2 DS_mask_neighbours +#define DS_mask_3 DS_mask_2|DS_mask_density|DS_mask_abundance|DS_mask_turb_doppler|DS_mask_temperatures|DS_mask_ACOEFF +#define DS_mask_populations (1 << DS_bit_populations) | DS_mask_3 +#define DS_mask_4 DS_mask_populations +#define DS_mask_all DS_mask_populations + /* input parameters */ typedef struct { double radius,radiusSqu,minScale,minScaleSqu,tcmb,taylorCutoff; @@ -208,6 +235,11 @@ void gasIIdust(double,double,double,double *); /* More functions */ +_Bool allBitsSet(const int flags, const int mask); +_Bool anyBitSet(const int flags, const int mask); +_Bool bitIsSet(const int flags, const int bitI); +_Bool onlyBitsSet(const int flags, const int mask); + void binpopsout(inputPars*, struct grid*, molData*); void calcAvRelLineAmp(struct grid*, int, int, double, double, double*); void calcAvRelLineAmp_lin(struct grid*, int, int, double, double, double*); From 749f89aef0cc3af8266c983825403852f94296de Mon Sep 17 00:00:00 2001 From: ims Date: Fri, 29 Jul 2016 11:10:17 +0200 Subject: [PATCH 11/13] Minor changes --- src/gridio.h | 5 +++++ src/velospline.c | 3 ++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/gridio.h b/src/gridio.h index ac5c038..944f217 100644 --- a/src/gridio.h +++ b/src/gridio.h @@ -7,6 +7,11 @@ * */ +#ifndef GRIDIO_H +#define GRIDIO_H + #define lime_fptr fitsfile #define lime_FITS 1 +#endif /* GRIDIO_H */ + diff --git a/src/velospline.c b/src/velospline.c index 23c3316..b92df7b 100644 --- a/src/velospline.c +++ b/src/velospline.c @@ -52,7 +52,8 @@ The number of coefficients N is currently hard-wired at 5. for(k=0;k Date: Fri, 29 Jul 2016 12:43:32 +0200 Subject: [PATCH 12/13] Grid i/o made more flexible There are two main changes here, although there are a lot of changed lines of code. - A generic i/o layer has been added, with choice of FITS just one of the (currently the only) possibilities for the format of the disk files. Hopefully this will make it easier for people to add interfaces to HDF5 or whatever at a later date. - The information which can legally be read from or written to the files has now been made much more flexible. When reading the file, if a particular element of the disk struct, say density, was not supplied, the task will automatically seek to fill in these missing values via the standard route. The completion state of the grid struct is now indicated by bit flags rather than a single 'stage' integer. --- src/aux.c | 4 +- src/grid.c | 242 +++++------ src/grid2fits.c | 1071 +++++++++++++++++++++++++++++----------------- src/gridio.c | 889 ++++++++++++++++++++------------------ src/lime.h | 38 +- src/main.c | 11 +- src/popsin.c | 13 + src/predefgrid.c | 36 +- 8 files changed, 1314 insertions(+), 990 deletions(-) diff --git a/src/aux.c b/src/aux.c index b19e8af..31b2336 100644 --- a/src/aux.c +++ b/src/aux.c @@ -48,7 +48,7 @@ parseInput(inputPars *par, image **img, molData **m){ par->writeGridAtStage[i] = 0; /* The user is not expected to set this. */ par->gridOutFiles[i] = ""; }; - par->dataStageI = 0; /* The user is not expected to set this. */ + par->dataFlags = 0; /* The user is not expected to set this. */ /* Allocate space for output fits images */ (*img)=malloc(sizeof(image)*MAX_NSPECIES); @@ -620,7 +620,7 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ if(par->binoutputfile != NULL) binpopsout(par,g,m); } - par->dataStageI = 4; + par->dataFlags |= DS_mask_4; for (i=0;inThreads;i++){ gsl_rng_free(threadRans[i]); diff --git a/src/grid.c b/src/grid.c index 6e9ee03..11a3955 100644 --- a/src/grid.c +++ b/src/grid.c @@ -9,33 +9,6 @@ #include "lime.h" - -void -gridAlloc(inputPars *par, struct grid **g){ - int i; - - *g=malloc(sizeof(struct grid)*(par->pIntensity+par->sinkPoints)); - memset(*g, 0., sizeof(struct grid) * (par->pIntensity+par->sinkPoints)); - - for(i=0;i<(par->pIntensity+par->sinkPoints); i++){ - (*g)[i].a0 = NULL; - (*g)[i].a1 = NULL; - (*g)[i].a2 = NULL; - (*g)[i].a3 = NULL; - (*g)[i].a4 = NULL; - (*g)[i].mol = NULL; - (*g)[i].dir = NULL; - (*g)[i].neigh = NULL; - (*g)[i].w = NULL; - (*g)[i].ds = NULL; - (*g)[i].dens=malloc(sizeof(double)*par->collPart); - (*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; - } -} - void freePopulation(const inputPars *par, const molData* m, struct populations* pop ) { if( pop !=NULL ) @@ -77,67 +50,29 @@ freePopulation(const inputPars *par, const molData* m, struct populations* pop ) free(pop); } } + void freeGrid(const inputPars *par, const molData* m ,struct grid* g){ int i; - if( g != NULL ) - { - for(i=0;i<(par->pIntensity+par->sinkPoints); i++){ - if(g[i].a0 != NULL) - { - free(g[i].a0); - } - if(g[i].a1 != NULL) - { - free(g[i].a1); - } - if(g[i].a2 != NULL) - { - free(g[i].a2); - } - if(g[i].a3 != NULL) - { - free(g[i].a3); - } - if(g[i].a4 != NULL) - { - free(g[i].a4); - } - if(g[i].dir != NULL) - { - free(g[i].dir); - } - if(g[i].neigh != NULL) - { - free(g[i].neigh); - } - if(g[i].w != NULL) - { - free(g[i].w); - } - if(g[i].dens != NULL) - { - free(g[i].dens); - } - if(g[i].nmol != NULL) - { - free(g[i].nmol); - } - if(g[i].abun != NULL) - { - free(g[i].abun); - } - if(g[i].ds != NULL) - { - free(g[i].ds); - } - if(g[i].mol != NULL) - { - freePopulation( par, m, g[i].mol ); - } - } - free(g); + if(g != NULL){ + for(i=0;incell;i++){ + free(g[i].a0); + free(g[i].a1); + free(g[i].a2); + free(g[i].a3); + free(g[i].a4); + free(g[i].dir); + free(g[i].neigh); + free(g[i].w); + free(g[i].dens); + free(g[i].nmol); + free(g[i].abun); + free(g[i].ds); + if(g[i].mol != NULL) + freePopulation( par, m, g[i].mol ); } + free(g); + } } void @@ -541,22 +476,10 @@ void readOrBuildGrid(inputPars *par, struct grid **gp){ int numCollPartRead; char message[80]; - par->dataStageI = 0; + par->dataFlags = 0; if(par->gridInFile!=NULL){ - /* Set defaults for *gridInfoRead: - */ - gridInfoRead.nInternalPoints = 0; - gridInfoRead.nSinkPoints = 0; - gridInfoRead.nLinks = 0; - gridInfoRead.nNNIndices = 0; - gridInfoRead.nDims = 0; - gridInfoRead.nSpecies = 0; - gridInfoRead.nDensities = 0; - gridInfoRead.nACoeffs = 0; - gridInfoRead.mols = NULL; - status = readGrid(par->gridInFile, lime_FITS, &gridInfoRead, gp\ - , &collPartNames, &numCollPartRead, &(par->dataStageI)); + , &collPartNames, &numCollPartRead, &(par->dataFlags)); if(status){ if(!silent){ @@ -566,16 +489,37 @@ void readOrBuildGrid(inputPars *par, struct grid **gp){ 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 are set as well: */ + if(bitIsSet(par->dataFlags, DS_bit_populations)\ + && !allBitsSet(par->dataFlags, (DS_mask_all & ~(1 << DS_bit_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(gridInfoRead.nSinkPoints!=par->sinkPoints){ + if((int)gridInfoRead.nSinkPoints!=par->sinkPoints){ if(!silent) warning("par->sinkPoints will be overwritten"); - par->sinkPoints = gridInfoRead.nSinkPoints; + par->sinkPoints = (int)gridInfoRead.nSinkPoints; } - if(gridInfoRead.nInternalPoints!=par->pIntensity){ + if((int)gridInfoRead.nInternalPoints!=par->pIntensity){ if(!silent) warning("par->pIntensity will be overwritten"); - par->pIntensity = gridInfoRead.nInternalPoints; + par->pIntensity = (int)gridInfoRead.nInternalPoints; } par->ncell = par->sinkPoints + par->pIntensity; } @@ -586,7 +530,7 @@ void readOrBuildGrid(inputPars *par, struct grid **gp){ } exit(1); } - if(gridInfoRead.nSpecies>0 && par->nSpecies>0 && gridInfoRead.nSpecies!=par->nSpecies){ + 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); @@ -595,7 +539,7 @@ void readOrBuildGrid(inputPars *par, struct grid **gp){ 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->collPart>0 && gridInfoRead.nDensities!=par->collPart){ + if(gridInfoRead.nDensities>0 && par->collPart>0 && (int)gridInfoRead.nDensities!=par->collPart){ if(!silent){ sprintf(message, "Grid file had %d densities but you have provided %d."\ , (int)gridInfoRead.nDensities, par->collPart); @@ -609,11 +553,11 @@ void readOrBuildGrid(inputPars *par, struct grid **gp){ **** 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? */ - if(par->dataStageI==4 && !silent) + if(allBitsSet(par->dataFlags, DS_mask_4) && !silent) warning("Sorry, LIME is not yet able to further process level populations read from file."); } /* End of read grid file. Whether and what we subsequently calculate will depend on the value of par->dataStageI returned. */ - if(par->dataStageI<1){ /* This can only happen if we did not read a file. */ + 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); gsl_rng *ran = gsl_rng_alloc(gsl_rng_ranlxs2); /* Random number generator */ @@ -703,62 +647,86 @@ void readOrBuildGrid(inputPars *par, struct grid **gp){ smooth(par, *gp); if(!silent) printDone(5); - par->dataStageI = 1; - }else if(par->dataStageI==1 && par->writeGridAtStage[par->dataStageI-1]){ - sprintf(message, "You just read a grid file at data stage %d, now you want to write it again?", par->dataStageI); - if(!silent) warning(message); + par->dataFlags |= DS_mask_1; } - if(par->dataStageI==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. */ + 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); - if(par->dataStageI<2){ + if(!allBitsSet(par->dataFlags, DS_mask_neighbours)){ qhull(par, *gp); /* Mallocs and sets .neigh, sets .numNeigh */ - par->dataStageI = 2; - }else if(par->dataStageI==2 && par->writeGridAtStage[par->dataStageI-1]){ - sprintf(message, "You just read a grid file at data stage %d, now you want to write it again?", par->dataStageI); - if(!silent) warning(message); + par->dataFlags |= DS_mask_2; } + 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. */ - if(par->dataStageI==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. */ + 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); - distCalc(par, *gp); /* Mallocs and sets .dir & .ds, sets .nphot */ + 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); + for(i=par->pIntensity;incell;i++){ + for(j=0;jdataStageI<3){ - for(i=0;incell; i++){ + par->dataFlags |= DS_mask_velocity; + } + + if(!allBitsSet(par->dataFlags, DS_mask_density)){ + for(i=0;incell; i++) (*gp)[i].dens = malloc(sizeof(double)*par->collPart); + 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++) + (*gp)[i].dens[0]=1e-30; + + par->dataFlags |= DS_mask_density; + } + + if(!allBitsSet(par->dataFlags, DS_mask_abundance)){ + for(i=0;incell; i++){ (*gp)[i].abun = malloc(sizeof(double)*par->nSpecies); - (*gp)[i].nmol = 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++) + (*gp)[i].abun[0]=0; - for(i=0;ipIntensity;i++){ - density( (*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2], (*gp)[i].dens); - temperature((*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2], (*gp)[i].t); - doppler( (*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2],&(*gp)[i].dopb); - abundance( (*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2], (*gp)[i].abun); - velocity( (*gp)[i].x[0],(*gp)[i].x[1],(*gp)[i].x[2], (*gp)[i].vel); - } + par->dataFlags |= DS_mask_abundance; + } + for(i=0;incell; i++) /* We don't store the nmol values so we have to do this malloc whether we read a file or not. */ + (*gp)[i].nmol = malloc(sizeof(double)*par->nSpecies); //**** mind you, it would be better to malloc them just before calculating them in molinit. + + 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); + for(i=par->pIntensity;incell;i++) + (*gp)[i].dopb=0.; + + par->dataFlags |= DS_mask_turb_doppler; + } + + 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].dens[0]=1e-30; (*gp)[i].t[0]=par->tcmb; (*gp)[i].t[1]=par->tcmb; - (*gp)[i].dopb=0.; - (*gp)[i].abun[0]=0; - for(j=0;jdataFlags |= DS_mask_temperatures; + } + + if(!allBitsSet(par->dataFlags, DS_mask_ACOEFF)){ calcInterpCoeffs(par,*gp); /* Mallocs and sets .a0, .a1 etc. */ - par->dataStageI = 3; - }else if(par->dataStageI==3 && par->writeGridAtStage[par->dataStageI-1]){ - sprintf(message, "You just read a grid file at data stage %d, now you want to write it again?", par->dataStageI); - if(!silent) warning(message); + par->dataFlags |= DS_mask_3; } - if(par->dataStageI==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. */ + if(onlyBitsSet(par->dataFlags, 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); dumpGrid(par,*gp); diff --git a/src/grid2fits.c b/src/grid2fits.c index b806e14..a573896 100644 --- a/src/grid2fits.c +++ b/src/grid2fits.c @@ -10,14 +10,15 @@ TODOS (some day): - 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()? - - If any of the stage 3 columns (VEL etc) are not present, get the values via user-supplied functions. */ #include "lime.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 minimum stage integer associated with the presence of a particular extension, column or keyword is given on the leftmost place of each line. +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 @@ -25,76 +26,56 @@ Note that all extensions are binary table except where indicated. The letter in for a key to these. -1 0) The primary HDU - Keywords: -1 STAGE I +0 0) The primary HDU -1 1) GRID +0 1) GRID Number of rows = number of grid points. Keywords: -1 RADIUS D -1 COLLPARn A # 1 for each nth collision partner. +0 RADIUS D +0 COLLPARn A # 1 for each nth collision partner. Columns: -1 ID V -1 Xj D # Cartesian components of the point location, 1 col per jth dimension. -1 IS_SINK L # =True iff the point lies on the edge of the model. -2 NUMNEIGH U -2 FIRST_NN V # See explanation in section 3 below. -3 VELj D # 1 col per jth dimension. +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. -3 TEMPKNTC E # From t[0]. -3 TEMPDUST E # From t[1]. -3 TURBDPLR E # Given Gaussian lineshape exp(-v^2/[B^2 + 2*k*T/m]), this is B. -3 ABUNMOLm E # 1 per mth molecular species. +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]. -2 2) NN_INDICES (see explanation in the header to gridio.c) +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: -2 LINK_I V +1 LINK_I V -2 3) LINKS (see explanation in the header to gridio.c) +1 3) LINKS (see explanation in the header to gridio.c) Number of rows = number of Delaunay links. Columns: -2 GRID_I_1 V -2 GRID_I_2 V -3 ACOEFF_p D # 1 per pth order of the velocity polynomial. +1 GRID_I_1 V +1 GRID_I_2 V +7 ACOEFF_p D # 1 per pth order of the velocity polynomial. -4 4 etc) LEVEL_POPS_m (1 per mth molecular species) +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: -4 MOL_NAME +8 MOL_NAME Note that at present the data in the 'partner' element of grid.mol is *NOT* being stored. */ /*....................................................................*/ -fitsfile *openFITSFileForRead(char *inFileName, int *dataStageI){ - fitsfile *fptr=NULL; - int status=0; - short tempStageI; - - fits_open_file(&fptr, inFileName, READONLY, &status); - processFitsError(status); - - fits_read_key(fptr, TSHORT, "STAGE ", &tempStageI, NULL, &status); - processFitsError(status); - - *dataStageI = (int)tempStageI; - - return fptr; -} - -/*....................................................................*/ -fitsfile *openFITSFileForWrite(char *outFileName, const int dataStageI){ +fitsfile *openFITSFileForWrite(char *outFileName){ fitsfile *fptr=NULL; int status=0; char negfile[100]="! "; - short dsiShort = (short)dataStageI, *ptrToShort=&dsiShort; /* Have to do this otherwise the compiler complains that the const int is not an 'lvalue'. */ fits_create_file(&fptr, outFileName, &status); if(status!=0){ @@ -108,9 +89,6 @@ fitsfile *openFITSFileForWrite(char *outFileName, const int dataStageI){ fits_create_img(fptr, 8, 0, NULL, &status); processFitsError(status); - fits_write_key(fptr, TSHORT, "STAGE ", ptrToShort, "Data stage", &status); - processFitsError(status); - return fptr; } @@ -122,274 +100,194 @@ void closeFITSFile(fitsfile *fptr){ processFitsError(status); } -/*....................................................................*/ -void defineGridExtColumns(const unsigned short numKwdChars, inputPars par\ - , const unsigned short numDims, const int dataStageI, char *ttype[]\ - , char *tform[], char *tunit[], int dataTypes[]){ - - /* Note that data types in all capitals are defined in fitsio.h. */ - - const int maxNumDims=9, maxNumSpecies=9, maxNumCollPart=9; - int colI=0, i; - char message[80]; - - if(numDims>maxNumDims){ - if(!silent){ - sprintf(message, "Caller asked for %d dims but colnames can only be written for %d.", numDims, maxNumDims); - bail_out(message); - } - exit(1); - } - - if(par.nSpecies>maxNumSpecies){ - if(!silent){ - sprintf(message, "Caller asked for %d species but colnames can only be written for %d.", par.nSpecies, maxNumSpecies); - bail_out(message); - } - exit(1); - } - - if(par.collPart>maxNumCollPart){ - if(!silent){ - sprintf(message, "Caller asked for %d coll. part. but colnames can only be written for %d.", par.collPart, maxNumCollPart); - bail_out(message); - } - exit(1); - } - - colI = 0; - ttype[colI] = malloc(numKwdChars); - sprintf(ttype[colI], "ID"); - tform[colI] = "V"; - tunit[colI] = "\0"; - dataTypes[colI] = TUINT; - - /* should rather have a vector column? */ - for(i=0;i1){ - colI++; + colI = allColNumbers[getColIndex(allColNames, maxNumCols, "NUMNEIGH")]; + if(colI>0){ numNeigh = malloc(sizeof(*numNeigh)*par.ncell); - for(i=0;i2){ - velj = malloc(sizeof(*velj)*par.ncell); - for(j=0;j0){ + velj = malloc(sizeof(*velj)*par.ncell); + for(j=0;j0){ + densn = malloc(sizeof(*densn)*par.ncell); + for(n=0;n0){ + abunm = malloc(sizeof(*abunm)*par.ncell); + for(m=0;m0){ + dopb = malloc(sizeof(*dopb)*par.ncell); + for(i=0;i0){ + t = malloc(sizeof(*t)*par.ncell); + for(i=0;imaxNumDims){ + if(!silent){ + sprintf(message, "Caller asked for %d dims but colnames can only be written for %d.", (int)numDims, (int)maxNumDims); + bail_out(message); + } + exit(1); + } + + if(numSpecies>maxNumSpecies){ + if(!silent){ + sprintf(message, "Caller asked for %d species but colnames can only be written for %d.", (int)numSpecies, (int)maxNumSpecies); + bail_out(message); + } + exit(1); + } + + if(numDensities>maxNumDensities){ + if(!silent){ + sprintf(message, "Caller asked for %d coll. part. but colnames can only be written for %d.", (int)numDensities, (int)maxNumDensities); + bail_out(message); + } + exit(1); + } + + *maxNumCols = 7 + numDims*2 + numSpecies + numDensities; + + *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=0;i0){ + 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); +} + +/*....................................................................*/ +int getColIndex(char **allColNames, const int maxNumCols, char *colName){ + int i=0; + int colFound=0; /* -> bool */ + + while(i; + + */ + const int numColNameChars=21; unsigned int *ids=NULL; double *aCoeffs=NULL; int status=0, colI=0, i, n; LONGLONG firstRow=1, firstElem=1; int numCols; char extname[] = "LINKS"; - const int numKwdChars = 9; /* 8 characters + \0. */ + char **tempColNames=NULL; - if(dataStageI<2){ - if(!silent) bail_out("Data stage indicates no link data!"); + if(links==NULL){ + if(!silent) bail_out("No link data!"); exit(1); } - if(dataStageI<3) + if(links[0].aCoeffs==NULL) numCols = 2; else numCols = 2 + numACoeffs; @@ -492,57 +649,61 @@ Note that data types in all capitals are defined in fitsio.h. char *tunit[numCols]; int dataTypes[numCols]; + tempColNames = malloc(sizeof(*tempColNames)*numCols); + for(i=0;i2){ + if(links[0].aCoeffs!=NULL){ /* Should rather have a vector column? */ for(n=0;nid; + ids[i] = (unsigned int)links[i].gp[0]->id; fits_write_col(fptr, dataTypes[colI], colI+1, firstRow, firstElem, (LONGLONG)totalNumLinks, ids, &status); processFitsError(status); colI++; for(i=0;iid; + ids[i] = (unsigned int)links[i].gp[1]->id; fits_write_col(fptr, dataTypes[colI], colI+1, firstRow, firstElem, (LONGLONG)totalNumLinks, ids, &status); processFitsError(status); free(ids); - if(dataStageI>2){ + if(links[0].aCoeffs!=NULL){ aCoeffs = malloc(sizeof(*aCoeffs)*totalNumLinks); for(n=0;n1 then it also mallocs 'firstNearNeigh'. If dataStageI>2 then the function also mallocs the following elements of struct grid for all grid points: - dens - abun -Otherwise these are set to NULL. +The present function mallocs 'gp' and sets defaults for all the simple or first-level struct elements. If a COLLPARn keywords are found in the GRID extension header then collPartNames is malloc'd to the number of these. */ @@ -681,12 +850,20 @@ If a COLLPARn keywords are found in the GRID extension header then collPartNames */ fits_get_num_rowsll(fptr, &numGridCells, &status); processFitsError(status); + if(numGridCells<=0){ + if(!silent) warning("No rows found in grid dataset."); + return; /* I.e. with dataFlags left unchanged. */ + } mallocAndSetDefaultGrid(gp, (unsigned int)numGridCells); /* Read the columns. */ fits_get_colnum(fptr, CASEINSEN, "ID", &colNum, &status); + if(status==COL_NOT_FOUND){ + if(!silent) warning("No ID column found in grid dataset."); + return; /* I.e. with dataFlags left unchanged. */ + } processFitsError(status); ids = malloc(sizeof(*ids)*numGridCells); @@ -699,6 +876,10 @@ If a COLLPARn keywords are found in the GRID extension header then collPartNames free(ids); gridInfoRead->nDims = (unsigned short)countColumns(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. */ @@ -726,6 +907,10 @@ If a COLLPARn keywords are found in the GRID extension header then collPartNames free(xj); fits_get_colnum(fptr, CASEINSEN, "IS_SINK", &colNum, &status); + if(status==COL_NOT_FOUND){ + if(!silent) warning("No IS_SINK column found in grid dataset."); + return; /* I.e. with dataFlags left unchanged. */ + } processFitsError(status); sink = malloc(sizeof(*sink)*numGridCells); @@ -742,8 +927,12 @@ If a COLLPARn keywords are found in the GRID extension header then collPartNames gridInfoRead->nInternalPoints = numGridCells - gridInfoRead->nSinkPoints; - if(dataStageI>1){ - fits_get_colnum(fptr, CASEINSEN, "NUMNEIGH", &colNum, &status); + /* 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); @@ -756,28 +945,26 @@ If a COLLPARn keywords are found in the GRID extension header then collPartNames free(numNeigh); fits_get_colnum(fptr, CASEINSEN, "FIRST_NN", &colNum, &status); - processFitsError(status); - - *firstNearNeigh = malloc(sizeof(**firstNearNeigh)*numGridCells); - fits_read_col(fptr, TUINT, colNum, firstRow, firstElem, numGridCells, 0, *firstNearNeigh, &anynul, &status); - processFitsError(status); - } + if(status!=COL_NOT_FOUND){ + processFitsError(status); - if(dataStageI>2){ - /* Count the numbers of DENSITYn and ABUNMOLn columns: - */ - gridInfoRead->nDensities = (unsigned short)countColumns(fptr, "DENSITY"); - gridInfoRead->nSpecies = (unsigned short)countColumns(fptr, "ABUNMOL"); + *firstNearNeigh = malloc(sizeof(**firstNearNeigh)*numGridCells); + fits_read_col(fptr, TUINT, colNum, firstRow, firstElem, numGridCells, 0, *firstNearNeigh, &anynul, &status); + processFitsError(status); - /* Dimension the appropriate elements of gp: - */ - for(i_LL=0;i_LLnDensities); - (*gp)[i_LL].abun = malloc(sizeof(double)*gridInfoRead->nSpecies); - (*gp)[i_LL].nmol = malloc(sizeof(double)*gridInfoRead->nSpecies); - for(i_us=0;i_usnSpecies;i_us++) - (*gp)[i_LL].nmol[i_us] = 0.0; + /* If we made it to here, we can set the neighbour bit of dataFlags. Note however that this bit is only on trial til we check for the LINKS extension. So it had better behave itself. + */ + (*dataFlags) |= (1 << DS_bit_neighbours); } + } + + /* See if there are any VEL columns: + */ + status = 0; + sprintf(colName, "VEL1"); + fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); + if(status!=COL_NOT_FOUND){ + processFitsError(status); /* Read the VEL columns: */ @@ -796,19 +983,15 @@ If a COLLPARn keywords are found in the GRID extension header then collPartNames } free(velj); - /* Read the TURBDPLR column: - */ - fits_get_colnum(fptr, CASEINSEN, "TURBDPLR", &colNum, &status); - processFitsError(status); - - dopb = malloc(sizeof(*dopb)*numGridCells); - fits_read_col(fptr, TFLOAT, colNum, firstRow, firstElem, numGridCells, 0, dopb, &anynul, &status); - processFitsError(status); + (*dataFlags) |= (1 << DS_bit_velocity); + } - for(i_LL=0;i_LLnDensities = (unsigned short)countColumns(fptr, "DENSITY"); + if(gridInfoRead->nDensities > 0){ + for(i_LL=0;i_LLnDensities); /* Read the DENSITY columns: */ @@ -827,52 +1010,96 @@ If a COLLPARn keywords are found in the GRID extension header then collPartNames } free(densn); - /* Read the TEMPKNTC column: + (*dataFlags) |= (1 << DS_bit_density); + } + + /* Count the numbers of ABUNMOLm columns: + */ + gridInfoRead->nSpecies = (unsigned short)countColumns(fptr, "ABUNMOL"); + if(gridInfoRead->nSpecies > 0){ + for(i_LL=0;i_LLnSpecies); + } + + /* Read the ABUNMOL columns: */ - fits_get_colnum(fptr, CASEINSEN, "TEMPKNTC", &colNum, &status); + 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_LLnSpecies;i_us++){ - sprintf(colName, "ABUNMOL%d", (int)i_us+1); - fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); + fits_get_colnum(fptr, CASEINSEN, "TEMPDUST", &colNum, &status); + if(status==COL_NOT_FOUND){ /* Set t0 back to defaults: */ + 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. @@ -922,6 +1160,13 @@ The present function mallocs the pointer *links. 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, ppi, totalNumGridPoints); - bail_out(message); + 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, ppi, totalNumGridPoints); + bail_out(message); + } + exit(1); } - exit(1); + (*links)[i_LL].gp[0] = &gp[ppi]; } - (*links)[i_LL].g[0] = &gp[ppi]; - } - /* Read GRID_I_2 column. - */ - fits_get_colnum(fptr, CASEINSEN, "GRID_I_2", &colNum, &status); - processFitsError(status); + /* 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); + 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, ppi, totalNumGridPoints); - bail_out(message); + 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, ppi, totalNumGridPoints); + bail_out(message); + } + exit(1); + } + (*links)[i_LL].gp[1] = &gp[ppi]; } - exit(1); } - (*links)[i_LL].g[1] = &gp[ppi]; + free(ids); } - free(ids); - if(dataStageI>2){ /* Means we are in at least stage C. */ - /* Find out how many ACOEFF_* columns there are. - */ - gridInfoRead->nACoeffs = (unsigned short)countColumns(fptr, "ACOEFF_"); + 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 ACOEFF_* columns there are. + */ + gridInfoRead->nACoeffs = (unsigned short)countColumns(fptr, "ACOEFF_"); + if(gridInfoRead->nACoeffs<=0){ for(i_LL=0;i_LLnACoeffs); + (*links)[i_LL].aCoeffs = NULL; + return; + } - aCoeffs = malloc(sizeof(*aCoeffs)*totalNumLinks); - for(i_s=0;i_snACoeffs;i_s++){ - /* Read the ACOEFF_n columns. - */ - sprintf(colName, "ACOEFF_%d", (int)i_s+1); - fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); - processFitsError(status); + for(i_LL=0;i_LLnACoeffs); - fits_read_col(fptr, TDOUBLE, colNum, firstRow, firstElem, totalNumLinks, 0, aCoeffs, &anynul, &status); - processFitsError(status); + aCoeffs = malloc(sizeof(*aCoeffs)*totalNumLinks); + for(i_s=0;i_snACoeffs;i_s++){ + /* Read the ACOEFF_n columns. + */ + sprintf(colName, "ACOEFF_%d", (int)i_s+1); + fits_get_colnum(fptr, CASEINSEN, colName, &colNum, &status); + processFitsError(status); - for(i_LL=0;i_LLnACoeffs = 0; for(i_LL=0;i_LLnNNIndices = (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); diff --git a/src/gridio.c b/src/gridio.c index f30bc16..eb1e72f 100644 --- a/src/gridio.c +++ b/src/gridio.c @@ -5,38 +5,42 @@ * Copyright (C) 2006-2014 Christian Brinch * Copyright (C) 2015 The LIME development team * +TODO: + - After merging with master, try to unify freeReadGrid and freeGrid. */ #include "lime.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 an attempt to regulate this, the following five states of completeness (encoded in values of the variable 'dataStageI') have been defined: - - dataStageI==0: No useable data in grid. - - dataStageI==1. At this stage the vector of grid objects has been malloc'd and values have been generated for the following struct elements: - id - x - sink - - dataStageI==2. This stage is entered after the Delaunay neighbours of each grid point have been determined. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: - numNeigh - neigh - - dataStageI==3. This is entered after sampling the user-supplied functions for density, velocity etc. The following further struct elements are expected to have been malloc'd (in the case of pointers) and given values: - vel - a0, a1 etc. - dens - t - dopb - abun -Note that it is necessary to make this stage dependent on the previous because we need information about the near-neighbours to calculate the a0, a1 etc coefficients. - - dataStageI==4. After 1 or more iterations of populating the levels, we enter the final stage, in which all the grid struct elements have been malloc'd and given at least preliminary values; specifically now the element - mol - - *NOTE* that LIME may run to completion without ever reaching stage 4 - if all the images required were continuum ones, for example. - +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. + +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: dataStageI=0 dataStageI=1 dataStageI=2 dataStageI=3 dataStageI=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_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': ----------------------------------------------------------------- @@ -59,13 +63,31 @@ void writeGridIfRequired(inputPars *par, struct grid *gp, molData *md, const int 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[par->dataStageI-1]){ - status = writeGrid(par->gridOutFiles[par->dataStageI-1], fileFormatI\ - , *par, DIM, NUM_VEL_COEFFS, gp, md, collPartNames, par->dataStageI); + if(par->writeGridAtStage[dataStageI-1]){ + status = writeGrid(par->gridOutFiles[dataStageI-1], fileFormatI\ + , *par, DIM, NUM_VEL_COEFFS, gp, md, collPartNames, par->dataFlags); if(status){ - sprintf(message, "writeGrid at data stage %d returned with status %d", par->dataStageI, status); + sprintf(message, "writeGrid at data stage %d returned with status %d", dataStageI, status); if(!silent) bail_out(message); exit(1); } @@ -75,11 +97,11 @@ void writeGridIfRequired(inputPars *par, struct grid *gp, molData *md, const int } /*....................................................................*/ -lime_fptr *openFileForRead(char *inFileName, const int fileFormatI, int *dataStageI){ +lime_fptr *openFileForWrite(char *outFileName, const int fileFormatI){ lime_fptr *fptr=NULL; if(fileFormatI==lime_FITS){ - fptr = openFITSFileForRead(inFileName, dataStageI); + fptr = openFITSFileForWrite(outFileName); }else{ return NULL; } @@ -88,79 +110,9 @@ lime_fptr *openFileForRead(char *inFileName, const int fileFormatI, int *dataSta } /*....................................................................*/ -lime_fptr *openFileForWrite(char *outFileName, const int fileFormatI, const int dataStageI){ - lime_fptr *fptr=NULL; - - if(fileFormatI==lime_FITS){ - fptr = openFITSFileForWrite(outFileName, dataStageI); - }else{ - return NULL; - } - - return fptr; -} - -/*....................................................................*/ -void closeAndFree(lime_fptr *fptr, const int fileFormatI\ - , unsigned int *firstNearNeigh, struct linkType **nnLinks\ - , struct linkType *links, const unsigned int totalNumLinks){ - - unsigned int li; - - closeFile(fptr, fileFormatI); - free(firstNearNeigh); - free(nnLinks); - if(links!=NULL){ - for(li=0;linumNeigh;jA++){ @@ -222,8 +177,8 @@ See the comment at the beginning of the present module for a description of how */ linkId = li; (*links)[li].id = linkId; - (*links)[li].g[0] = gAPtr; - (*links)[li].g[1] = gBPtr; + (*links)[li].gp[0] = gAPtr; + (*links)[li].gp[1] = gBPtr; li++; } @@ -238,19 +193,19 @@ See the comment at the beginning of the present module for a description of how *links = realloc(*links, sizeof(struct linkType)*(*totalNumLinks)); - if(dataStageI>2){ + if(allBitsSet(dataFlags, DS_mask_ACOEFF)){ for(li=0;li<*totalNumLinks;li++) (*links)[li].aCoeffs = malloc(sizeof(double)*NUM_VEL_COEFFS); for(li=0;li<*totalNumLinks;li++){ - if((*links)[li].g[0]->sink && (*links)[li].g[1]->sink){ + if((*links)[li].gp[0]->sink && (*links)[li].gp[1]->sink){ for(ci=0;cisink){ /* If we get to here, this ensures that g[1] is not a sink point. */ + if((*links)[li].gp[0]->sink){ /* If we get to here, this ensures that gp[1] is not a sink point. */ nearI = 1; evenCoeffSign = -1.0; }else{ @@ -258,12 +213,12 @@ See the comment at the beginning of the present module for a description of how evenCoeffSign = 1.0; } - gAPtr = (*links)[li].g[nearI]; + gAPtr = (*links)[li].gp[nearI]; /* Find which neighbour of gAPtr corresponds to the link: */ linkNotFound = 1; for(jA=0;jAnumNeigh;jA++){ idB = gAPtr->neigh[jA]->id; - if(linkNotFound && idB==(*links)[li].g[1-nearI]->id){ + if(linkNotFound && idB==(*links)[li].gp[1-nearI]->id){ (*links)[li].aCoeffs[0] = evenCoeffSign*gAPtr->a0[jA]; (*links)[li].aCoeffs[1] = gAPtr->a1[jA]; (*links)[li].aCoeffs[2] = evenCoeffSign*gAPtr->a2[jA]; @@ -279,7 +234,7 @@ See the comment at the beginning of the present module for a description of how } } /* end if not both sink */ } /* end loop over links */ - }else{ /* dataStageI<=2 */ + }else{ /* a coeffs not available */ for(li=0;li<*totalNumLinks;li++) (*links)[li].aCoeffs = NULL; } @@ -293,86 +248,40 @@ See the comment at the beginning of the present module for a description of how /* The calling routine must free firstNearNeigh, nnLinks, links. */ } - /*....................................................................*/ -int writeGrid(char *outFileName, const int fileFormatI, inputPars par\ - , unsigned short numDims, unsigned short numACoeffs, struct grid *g\ - , molData *md, char **collPartNames, const int dataStageI){ - - lime_fptr *fptr=NULL; - int status = 0; - unsigned short speciesI; - char message[80]; - struct linkType *links=NULL, **nnLinks=NULL; - unsigned int totalNumLinks, totalNumNeigh, *firstNearNeigh=NULL; - - if (outFileName==""){ - if(!silent) warning("Cannot write grid list to file, filename is blank."); - return 1; - } - - if (g==NULL || dataStageI<1){ - if(!silent) warning("Cannot write grid list to file, there are no entries in it."); - return 2; - } - - sprintf(message, "Writing grid-point list to file %s", outFileName); - if(!silent) printMessage(message); - - if (dataStageI>1){ - constructLinkArrays((unsigned int)par.ncell, g, &links, &totalNumLinks\ - , &nnLinks, &firstNearNeigh, &totalNumNeigh, dataStageI); - } - - fptr = openFileForWrite(outFileName, fileFormatI, dataStageI); - if(fptr==NULL){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); - return 3; - } +void closeFile(lime_fptr *fptr, const int fileFormatI){ + if(fileFormatI==lime_FITS){ + closeFITSFile(fptr); + }//**** error if not? +} - status = writeGridBlock(fptr, fileFormatI, par, numDims, g, firstNearNeigh, collPartNames, dataStageI); - if(status){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); - return 4; - } +/*....................................................................*/ +void closeAndFree(lime_fptr *fptr, const int fileFormatI\ + , unsigned int *firstNearNeigh, struct linkType **nnLinks\ + , struct linkType *links, const unsigned int totalNumLinks){ - if (dataStageI>1){ - status = writeNnIndicesBlock(fptr, fileFormatI, totalNumNeigh, nnLinks, links); - if(status){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); - return 5; - } + unsigned int li; - status = writeLinksBlock(fptr, fileFormatI, totalNumLinks, numACoeffs, links, dataStageI); - if(status){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); - return 6; - } + closeFile(fptr, fileFormatI); - if (dataStageI>3){ - for(speciesI=0;speciesINUM_GRID_STAGES){ - if(!silent) bail_out("Data stage out of range."); - exit(1); - } + nnIndices Basically stores the information in the element 'neigh'. - sprintf(message, "Grid-point data stage determined as %d", *dataStageI); - if(!silent) printMessage(message); + links Stores the elements a0, a1 etc. - /* Read the values which should be in grid for every stage. + pops Stores the element 'pops' of each element 'mol'. There will be a separate table for each radiating species. + --------------------------------------------------------------- */ - status = readGridBlock(fptr, fileFormatI, gridInfoRead, gp, &firstNearNeigh\ - , collPartNames, numCollPartRead, *dataStageI); - if(status || gridInfoRead->nSinkPoints<=0\ - || gridInfoRead->nInternalPoints<=0 || gridInfoRead->nDims<=0){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, NULL, NULL, 0); - freeReadGrid(gp, *gridInfoRead); + lime_fptr *fptr=NULL; + int status = 0; + unsigned short speciesI; + char message[80]; + struct linkType *links=NULL, **nnLinks=NULL; + unsigned int totalNumLinks=-1, totalNumNeigh=-1, *firstNearNeigh=NULL; - if(status) - return 1; - else if(gridInfoRead->nSinkPoints<=0) - return 2; - else if(gridInfoRead->nInternalPoints<=0) - return 3; - else if(gridInfoRead->nDims<=0) - return 4; - else{ - sprintf(message, "This indicates a programming error. Please contact the developer."); - if(!silent) bail_out(message); - exit(1); - } + if (outFileName==""){ + if(!silent) warning("Cannot write grid list to file, filename is blank."); + return 1; } - if((*dataStageI)>2 && (gridInfoRead->nDensities<=0 || gridInfoRead->nSpecies<=0)){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, NULL, NULL, 0); - freeReadGrid(gp, *gridInfoRead); + if (gp==NULL || !allBitsSet(dataFlags, DS_mask_x)){ + if(!silent) warning("Cannot write grid list to file, there are no entries in it."); + return 2; + } - if(gridInfoRead->nDensities<=0) - return 5; - else if(gridInfoRead->nSpecies<=0) //***** what if all continuum images?? - return 6; - else{ - sprintf(message, "This indicates a programming error. Please contact the developer."); - if(!silent) bail_out(message); - exit(1); - } + sprintf(message, "Writing grid-point list to file %s", outFileName); + if(!silent) printMessage(message); + + constructLinkArrays((unsigned int)par.ncell, gp, &links, &totalNumLinks\ + , &nnLinks, &firstNearNeigh, &totalNumNeigh, dataFlags); + + fptr = openFileForWrite(outFileName, fileFormatI); + if(fptr==NULL){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); + return 3; } - totalNumGridPoints = gridInfoRead->nInternalPoints+gridInfoRead->nSinkPoints; + status = writeGridTable(fptr, fileFormatI, par, numDims, gp, firstNearNeigh, collPartNames, dataFlags); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); + return 4; + } - if((*dataStageI)>1){ /* there should be blocks for the nnIndices and links. */ - status = readLinksBlock(fptr, fileFormatI, gridInfoRead, *gp, &links, *dataStageI); + if (links!=NULL && nnLinks!=NULL){ + status = writeNnIndicesTable(fptr, fileFormatI, totalNumNeigh, nnLinks, links); if(status){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, NULL, links, gridInfoRead->nLinks); - freeReadGrid(gp, *gridInfoRead); - return 7; + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); + return 5; } - if ((*dataStageI)>2 && gridInfoRead->nACoeffs<=0){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); - freeReadGrid(gp, *gridInfoRead); - return 8; + status = writeLinksTable(fptr, fileFormatI, totalNumLinks, numACoeffs, links); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, totalNumLinks); + return 6; } + } - status = readNnIndicesBlock(fptr, fileFormatI, links, &nnLinks, gridInfoRead); - if(status){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); - freeReadGrid(gp, *gridInfoRead); - return 9; + if (allBitsSet(dataFlags, DS_mask_populations)){ + for(speciesI=0;speciesI3){ /* there should be pops blocks. */ - status = getNumPopsBlocks(fptr, fileFormatI, &numBlocks); - if(status || numBlocks<=0 || numBlocks!=gridInfoRead->nSpecies){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); - freeReadGrid(gp, *gridInfoRead); - if(status) - return 10; - else if(numBlocks<=0) - return 11; - else if(numBlocks!=gridInfoRead->nSpecies) - return 12; - else{ - sprintf(message, "This indicates a programming error. Please contact the developer."); - if(!silent) bail_out(message); - exit(1); - } - } +/*....................................................................*/ +lime_fptr *openFileForRead(char *inFileName, const int fileFormatI){ + lime_fptr *fptr=NULL; - gridInfoRead->mols = malloc(sizeof(struct molInfoType)*gridInfoRead->nSpecies); + if(fileFormatI==lime_FITS){ + fptr = openFITSFileForRead(inFileName); + }else{ + return NULL; + } - for(i_u=0;i_unSpecies); - for(i_s=0;i_snSpecies;i_s++){ - (*gp)[i_u].mol[i_s].dust = NULL; - (*gp)[i_u].mol[i_s].knu = NULL; - (*gp)[i_u].mol[i_s].pops = NULL; - (*gp)[i_u].mol[i_s].partner = NULL; - } - } + return fptr; +} - for(i_s=0;i_snSpecies;i_s++){ - gridInfoRead->mols[i_s].molName = NULL; - gridInfoRead->mols[i_s].nLevels = -1; - gridInfoRead->mols[i_s].nLines = -1; - - status = readPopsBlock(fptr, fileFormatI, i_s, *gp, gridInfoRead); - - if(status || gridInfoRead->mols[i_s].nLevels<=0){ - closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); - freeReadGrid(gp, *gridInfoRead); - if(status) - return 13; - else if(gridInfoRead->mols[i_s].nLevels<=0) - return 14; - else{ - sprintf(message, "This indicates a programming error. Please contact the developer."); - if(!silent) bail_out(message); - exit(1); - } - } - } - } - } +/*....................................................................*/ +void freeReadGrid(const unsigned int numPoints, const unsigned short numSpecies\ + , struct grid *gp){ - /* Set unread pointers to NULL: - */ - if((*dataStageI)<4){ - for(i_u=0;i_unLinks); - return 0; } /*....................................................................*/ -int readGridBlock(lime_fptr *fptr, const int fileFormatI\ +int readGridTable(lime_fptr *fptr, const int fileFormatI\ , struct gridInfoType *gridInfoRead, struct grid **gp\ , unsigned int **firstNearNeigh, char ***collPartNames\ - , int *numCollPartRead, const int dataStageI){ + , int *numCollPartRead, int *dataFlags){ + /* +Individual routines called should set the appropriate bits of dataFlags; also malloc gp and set all its defaults. (Note there is a bespoke routine grid.c:mallocAndSetDefaultGrid() to do the latter.) + */ int status=0; if(fileFormatI==lime_FITS){ readGridExtFromFits(fptr, gridInfoRead, gp, firstNearNeigh\ - , collPartNames, numCollPartRead, dataStageI); + , collPartNames, numCollPartRead, dataFlags); }else{ status = 1; } @@ -632,14 +490,17 @@ int readGridBlock(lime_fptr *fptr, const int fileFormatI\ } /*....................................................................*/ -int readLinksBlock(lime_fptr *fptr, const int fileFormatI\ - , struct gridInfoType *gridInfoRead, struct grid *g\ - , struct linkType **links, const int dataStageI){ +int readLinksTable(lime_fptr *fptr, const int fileFormatI\ + , struct gridInfoType *gridInfoRead, struct grid *gp\ + , struct linkType **links, int *dataFlags){ + /* +Individual routines called should set the appropriate bits of dataFlags. + */ int status=0; if(fileFormatI==lime_FITS){ - readLinksExtFromFits(fptr, gridInfoRead, g, links, dataStageI); + readLinksExtFromFits(fptr, gridInfoRead, gp, links, dataFlags); }else{ status = 1; } @@ -648,13 +509,16 @@ int readLinksBlock(lime_fptr *fptr, const int fileFormatI\ } /*....................................................................*/ -int readNnIndicesBlock(lime_fptr *fptr, const int fileFormatI, struct linkType *links\ - , struct linkType ***nnLinks, struct gridInfoType *gridInfoRead){ +int readNnIndicesTable(lime_fptr *fptr, const int fileFormatI, struct linkType *links\ + , struct linkType ***nnLinks, struct gridInfoType *gridInfoRead, int *dataFlags){ + /* +Individual routines called should set the appropriate bits of dataFlags. + */ int status=0; if(fileFormatI==lime_FITS){ - readNnIndicesExtFromFits(fptr, links, nnLinks, gridInfoRead); + readNnIndicesExtFromFits(fptr, links, nnLinks, gridInfoRead, dataFlags); }else{ status = 1; } @@ -664,14 +528,39 @@ int readNnIndicesBlock(lime_fptr *fptr, const int fileFormatI, struct linkType * /*....................................................................*/ void loadNnIntoGrid(unsigned int *firstNearNeigh, struct linkType **nnLinks\ - , struct linkType *links, struct gridInfoType gridInfoRead, struct grid *gp\ - , const int dataStageI){ + , 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 extension 'neigh' of struct g for each grid point. + */ + + const unsigned int totalNumGridPoints = gridInfoRead.nInternalPoints+gridInfoRead.nSinkPoints; + unsigned int i_u; + int j; + struct linkType *linkPtr; + + for(i_u=0;i_ugp[0]->id==(int)i_u){ + gp[i_u].neigh[j] = linkPtr->gp[1]; + }else{ + gp[i_u].neigh[j] = linkPtr->gp[0]; + } + } + } +} + +/*....................................................................*/ +void loadACoeffsIntoGrid(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 the following extensions of struct g for each grid point: - neigh -if dataStageI>2: a0 a1 a2 @@ -690,89 +579,99 @@ Further: the 'a' coefficients for the link or edge between A and B are stored in char message[80]; double evenCoeffSign; - for(i_u=0;i_ug[0]->id==(int)i_u){ - gp[i_u].neigh[j] = linkPtr->g[1]; + if(linkPtr->gp[0]->id==(int)i_u){ + evenCoeffSign = 1.0; }else{ - gp[i_u].neigh[j] = linkPtr->g[0]; + evenCoeffSign = -1.0; } + + gp[i_u].a0[j] = evenCoeffSign*(linkPtr->aCoeffs)[0]; + gp[i_u].a1[j] = (linkPtr->aCoeffs)[1]; + gp[i_u].a2[j] = evenCoeffSign*(linkPtr->aCoeffs)[2]; + gp[i_u].a3[j] = (linkPtr->aCoeffs)[3]; + gp[i_u].a4[j] = evenCoeffSign*(linkPtr->aCoeffs)[4]; } } +} - if(dataStageI>2){ - /* Just for the time being: */ - if(gridInfoRead.nACoeffs!=5){ - if(!silent){ - sprintf(message, "There should be %d ACOEFF_n columns, but %d were read.", NUM_VEL_COEFFS, (int)gridInfoRead.nACoeffs); - bail_out(message); - } - exit(1); - } +/*....................................................................*/ +int checkPopsTableExists(lime_fptr *fptr, const int fileFormatI\ + , const unsigned short speciesI, _Bool *blockFound){ + int status=0; - for(i_u=0;i_ug[0]->id==(int)i_u){ - evenCoeffSign = 1.0; - }else{ - evenCoeffSign = -1.0; - } + *blockFound = 0; - gp[i_u].a0[j] = evenCoeffSign*linkPtr->aCoeffs[0]; - gp[i_u].a1[j] = linkPtr->aCoeffs[1]; - gp[i_u].a2[j] = evenCoeffSign*linkPtr->aCoeffs[2]; - gp[i_u].a3[j] = linkPtr->aCoeffs[3]; - gp[i_u].a4[j] = evenCoeffSign*linkPtr->aCoeffs[4]; - } - } + if(fileFormatI==lime_FITS){ + *blockFound = checkPopsFitsExtExists(fptr, speciesI); + }else{ + status = 1; } + + return status; } /*....................................................................*/ -int getNumPopsBlocks(lime_fptr *fptr, const int fileFormatI, unsigned short *numBlocks){ +int getNumPopsTables(lime_fptr *fptr, const int fileFormatI, unsigned short *numTables){ int status = 0; _Bool blockFound = 1; - *numBlocks = 0; + *numTables = 0; while(blockFound && !status){ - status = checkPopsBlockExists(fptr, fileFormatI, *numBlocks, &blockFound); - (*numBlocks)++; + status = checkPopsTableExists(fptr, fileFormatI, *numTables, &blockFound); + (*numTables)++; } - (*numBlocks)--; + (*numTables)--; return status; } /*....................................................................*/ -int checkPopsBlockExists(lime_fptr *fptr, const int fileFormatI\ - , const unsigned short speciesI, _Bool *blockFound){ +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_u; - *blockFound = 0; + totalNumGridPoints = gridInfoRead->nInternalPoints + gridInfoRead->nSinkPoints; + for(i_u=0;i_unInternalPoints = 0; + gridInfoRead->nSinkPoints = 0; + gridInfoRead->nLinks = 0; + gridInfoRead->nNNIndices = 0; + gridInfoRead->nDims = 0; + gridInfoRead->nSpecies = 0; + gridInfoRead->nDensities = 0; + gridInfoRead->nACoeffs = 0; + gridInfoRead->mols = NULL; + + /* Open the file and also return the data stage. */ + fptr = openFileForRead(inFileName, fileFormatI); + + /* 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); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 1; } - return status; + /* Some sanity checks: + */ + if((*dataFlags)!=0 && gridInfoRead->nSinkPoints<=0 || gridInfoRead->nInternalPoints<=0 || gridInfoRead->nDims<=0){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, 0); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + + if(gridInfoRead->nSinkPoints<=0) + return 2; + else if(gridInfoRead->nInternalPoints<=0) + return 3; + else if(gridInfoRead->nDims<=0) + return 4; + 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); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + + if(gridInfoRead->nDensities<=0) + return 5; + else if(gridInfoRead->nSpecies<=0) //***** what if all continuum images?? + return 6; + 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); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 7; + } + + /* Sanity check: + */ + if (allBitsSet(*dataFlags, DS_mask_ACOEFF) && gridInfoRead->nACoeffs<=0){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 8; + } + + status = readNnIndicesTable(fptr, fileFormatI, links, &nnLinks, gridInfoRead, dataFlags); /* Sets appropriate bits of dataFlags. */ + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 9; + } + + 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)) + loadACoeffsIntoGrid(firstNearNeigh, nnLinks, *gridInfoRead, *gp); + } + + status = getNumPopsTables(fptr, fileFormatI, &numTables); + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 10; + } + + /* Sanity check: + */ + if(numTables>0 && numTables!=gridInfoRead->nSpecies){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 11; + } + + if(numTables>0){ + (*dataFlags) |= (1 << DS_bit_populations); + + for(i_u=0;i_unSpecies); + + gridInfoRead->mols = malloc(sizeof(struct molInfoType)*gridInfoRead->nSpecies); + for(i_s=0;i_snSpecies;i_s++){ + gridInfoRead->mols[i_s].molName = NULL; + gridInfoRead->mols[i_s].nLevels = -1; + gridInfoRead->mols[i_s].nLines = -1; + + status = readPopsTable(fptr, fileFormatI, i_s, *gp, gridInfoRead); /* Sets defaults for all the fields under grid.mol. */ + if(status){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 12; + } + + /* Sanity check: + */ + if(gridInfoRead->mols[i_s].nLevels<=0){ + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + freeReadGrid(totalNumGridPoints, gridInfoRead->nSpecies, *gp); + return 13; + } + } + } + + closeAndFree(fptr, fileFormatI, firstNearNeigh, nnLinks, links, gridInfoRead->nLinks); + return 0; } + diff --git a/src/lime.h b/src/lime.h index 3482b11..21c6fb2 100644 --- a/src/lime.h +++ b/src/lime.h @@ -119,8 +119,7 @@ typedef struct { char **moldatfile; _Bool writeGridAtStage[NUM_GRID_STAGES]; char *gridInFile,*gridOutFiles[NUM_GRID_STAGES]; - - int dataStageI,nSolveIters; + int dataFlags,nSolveIters; } inputPars; /* Molecular data: shared attributes */ @@ -188,7 +187,7 @@ struct gridInfoType{ struct linkType { unsigned int id; - struct grid *g[2]; + struct grid *gp[2]; double *aCoeffs; }; @@ -246,8 +245,7 @@ void calcAvRelLineAmp_lin(struct grid*, int, int, double, double, double*); void calcInterpCoeffs(inputPars*, struct grid*); void calcInterpCoeffs_lin(inputPars*, struct grid*); void calcSourceFn(double, const inputPars*, double*, double*); -int checkPopsBlockExists(lime_fptr*, const int, const unsigned short, _Bool*); -_Bool checkPopsFitsExtExists(lime_fptr*, const unsigned short); +_Bool checkPopsFitsExtExists(fitsfile*, const unsigned short); void closeAndFree(lime_fptr*, const int, unsigned int*, struct linkType**, struct linkType*, const unsigned int); void closeFile(lime_fptr*, const int); void closeFITSFile(fitsfile*); @@ -255,23 +253,23 @@ void constructLinkArrays(const unsigned int, struct grid*, struct linkType**, un void continuumSetup(int, image *, molData *, inputPars *, struct grid *); int countColumns(fitsfile*, char*); int countKeywords(fitsfile*, char*); -void defineGridExtColumns(const unsigned short, inputPars, const unsigned short, const int, char *ttype[], char *tform[], char *tunit[], int dataTypes[]); +void defineAndLoadColumns(fitsfile*, const unsigned short, const unsigned short, const unsigned short, const int, const unsigned short, char***, int**, int*, int**); void distCalc(inputPars *, struct grid *); void fit_d1fi(double, double, double*); void fit_fi(double, double, double*); void fit_rr(double, double, double*); -void freeReadGrid(struct grid**, struct gridInfoType); +void freeReadGrid(const unsigned int, const unsigned short, struct grid*); 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); double gaussline(double, double); void getArea(inputPars *, struct grid *, const gsl_rng *); void getclosest(double, double, double, long *, long *, double *, double *, double *); +int getColIndex(char**, const int, char*); 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*); int getNumPopsBlocks(lime_fptr*, const int, unsigned short*); -void gridAlloc(inputPars*, struct grid**); void input(inputPars *, image *); double interpolate(double, double, double, double, double, double); float invSqrt(float); @@ -280,14 +278,11 @@ void levelPops(molData*, inputPars*, struct grid*, int*); void line_plane_intersect(struct grid*, double*, int , int*, double*, double*, double); void lineBlend(molData*, inputPars*, blend**); void lineCount(int,molData*,int**, int**, int*); -void loadNnIntoGrid(unsigned int*, struct linkType**, struct linkType*, struct gridInfoType, struct grid*, const int); void LTE(inputPars*, struct grid*, molData*); void mallocAndSetDefaultGrid(struct grid**, const unsigned int); void molinit(molData*, inputPars*, struct grid*,int); -lime_fptr *openFileForRead(char*, const int, int*); -fitsfile *openFITSFileForRead(char*, int*); -lime_fptr *openFileForWrite(char*, const int, const int); -fitsfile *openFITSFileForWrite(char*, const int); +fitsfile* openFITSFileForRead(char*); +fitsfile* openFITSFileForWrite(char*); void openSocket(inputPars*, int); void qhull(inputPars*, struct grid*); void photon(int, struct grid*, molData*, int, const gsl_rng*,inputPars*,blend*,gridPointData*,double*); @@ -301,15 +296,10 @@ void processFitsError(int); double ratranInput(char*, char*, double, double, double); void raytrace(int, inputPars*, struct grid*, molData*, image*); int readGrid(char*, const int, struct gridInfoType*, struct grid**, char***, int*, int*); -int readGridBlock(lime_fptr*, const int, struct gridInfoType*, struct grid**, unsigned int**, char***, int*, const int); -void readGridExtFromFits(fitsfile*, struct gridInfoType*, struct grid**, unsigned int**, char***, int*, const int); -void readOrBuildGrid(inputPars*, struct grid**); -int readLinksBlock(lime_fptr*, const int, struct gridInfoType*, struct grid*, struct linkType**, const int); -void readLinksExtFromFits(fitsfile*, struct gridInfoType*, struct grid*, struct linkType**, const int); -int readNnIndicesBlock(lime_fptr*, const int, struct linkType*, struct linkType***, struct gridInfoType*); -void readNnIndicesExtFromFits(fitsfile*, struct linkType*, struct linkType***, struct gridInfoType*); +void readGridExtFromFits(fitsfile*, struct gridInfoType*, struct grid**, unsigned int**, char***, int*, int*); +void readLinksExtFromFits(fitsfile*, struct gridInfoType*, struct grid*, struct linkType**, int*); +void readNnIndicesExtFromFits(fitsfile*, struct linkType*, struct linkType***, struct gridInfoType*, int*); void readOrBuildGrid(inputPars*, struct grid**); -int readPopsBlock(lime_fptr*, const int, const unsigned short, struct grid*, struct gridInfoType*); void readPopsExtFromFits(fitsfile*, const unsigned short, struct grid*, struct gridInfoType*); void report(int, inputPars*, struct grid*); void smooth(inputPars*, struct grid*); @@ -325,15 +315,11 @@ void traceray(rayData, int, int, inputPars*, struct grid*, molData*, image*, void velocityspline2(double*, double*, double, double, double, double*); double veloproject(double*, double*); int writeGrid(char*, const int, inputPars, unsigned short, unsigned short, struct grid*, molData*, char**, const int); -int writeGridBlock(lime_fptr*, const int, inputPars, unsigned short, struct grid*, unsigned int*, char**, const int); void writeGridExtToFits(fitsfile*, inputPars, unsigned short, struct grid*, unsigned int*, char**, const int); void writeGridIfRequired(inputPars*, struct grid*, molData*, const int); void writefits(int, inputPars*, molData*, image*); -int writeLinksBlock(lime_fptr*, const int, const unsigned int, const unsigned short, struct linkType*, const int); -void writeLinksExtToFits(fitsfile*, const unsigned int, const unsigned short, struct linkType*, const int); -int writeNnIndicesBlock(lime_fptr*, const int, const unsigned int, struct linkType**, struct linkType*); +void writeLinksExtToFits(fitsfile*, const unsigned int, const unsigned short, struct linkType*); void writeNnIndicesExtToFits(fitsfile*, const unsigned int, struct linkType**, struct linkType*); -int writePopsBlock(lime_fptr*, const int, unsigned int, molData*, unsigned short, struct grid*); void writePopsExtToFits(fitsfile*, const unsigned int, molData*, const unsigned short, struct grid*); void write_VTK_unstructured_Points(inputPars*, struct grid*); int factorial(const int); diff --git a/src/main.c b/src/main.c index ff5e8ed..a07f211 100644 --- a/src/main.c +++ b/src/main.c @@ -46,14 +46,12 @@ int main () { if(par.doPregrid) { - gridAlloc(&par,&gp); + mallocAndSetDefaultGrid(&gp, (unsigned int)par.ncell); predefinedGrid(&par,gp); - par.dataStageI = 3; // Sort of. } else if(par.restart) { popsin(&par,&gp,&md,&popsdone); - par.dataStageI = 4; // Sort of. } else { @@ -75,7 +73,8 @@ int main () { if(nLineImages>0 && !popsdone){ // eventually, replace !popdone by (!popsdone || dataStageI<4)? *Really* eventually we want to get rid of popsdone. levelPops(md,&par,gp,&popsdone); - par.dataStageI = 4; +// par.dataStageI = 4; + par.dataFlags |= (1 << DS_bit_populations); /* Disable the next lines for now, since we have not tested dataStageI<4 in the 'if' of this block, because we can't use an input grid file at dataStageI==4 yet: we have to disentangle all the functionality of molinit() before we can contemplate doing that. }else if(par.dataStageI==4 && par->nSolveIters>0 && par.writeGridAtStage[par.dataStageI-1]){ sprintf(message, "You just read a grid file at data stage %d, now you want to write it again?", par.dataStageI); @@ -83,8 +82,8 @@ int main () { */ } - if(par.dataStageI==4) - writeGridIfRequired(&par, gp, md, lime_FITS); +// if(par.dataStageI==4) + writeGridIfRequired(&par, gp, md, lime_FITS); /* Now make the line images. */ diff --git a/src/popsin.c b/src/popsin.c index bdf25ad..4f560d0 100644 --- a/src/popsin.c +++ b/src/popsin.c @@ -107,6 +107,19 @@ popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){ qhull(par, *g); distCalc(par, *g); calcInterpCoeffs(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_ACOEFF); + par->dataFlags |= (1 << DS_bit_populations); + +//**** should fill in any missing info via the appropriate function calls. + *popsdone=1; } diff --git a/src/predefgrid.c b/src/predefgrid.c index e60912e..6585e00 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -20,7 +20,14 @@ predefinedGrid(inputPars *par, struct grid *g){ #else gsl_rng_set(ran,time(0)); #endif - + + for(i=0;incell;i++){ + g[i].dens = malloc(sizeof(double)*par->collPart); + g[i].abun = malloc(sizeof(double)*par->nSpecies); +// g[i].nmol = malloc(sizeof(double)*par->nSpecies); +g[i].nmol = malloc(sizeof(double)*1); /* Only 1 element is accessed in the present bodgy state of affairs. */ + } + fp=fopen(par->pregrid,"r"); par->ncell=par->pIntensity+par->sinkPoints; @@ -39,14 +46,14 @@ 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].t[1]=g[i].t[0]; + g[i].nmol[0]=g[i].abun[0]*g[i].dens[0]; - /* This next step needs to be done, even though it looks stupid */ - g[i].dir=malloc(sizeof(point)*1); - g[i].ds =malloc(sizeof(double)*1); - g[i].neigh =malloc(sizeof(struct grid *)*1); - if(!silent) progressbar((double) i/((double)par->pIntensity-1), 4); + /* This next step needs to be done, even though it looks stupid */ + g[i].dir = malloc(sizeof(point)*1); + g[i].ds = malloc(sizeof(double)*1); + g[i].neigh = malloc(sizeof(struct grid *)*1); + if(!silent) progressbar((double) i/((double)par->pIntensity-1), 4); } for(i=par->pIntensity;incell;i++){ @@ -62,6 +69,7 @@ predefinedGrid(inputPars *par, struct grid *g){ g[i].sink=1; g[i].abun[0]=0; g[i].dens[0]=1e-30; + g[i].nmol[0]=0.0; /* Just to give it a value. */ g[i].t[0]=par->tcmb; g[i].t[1]=par->tcmb; g[i].dopb=0.; @@ -75,6 +83,18 @@ predefinedGrid(inputPars *par, struct grid *g){ // getArea(par,g, ran); // getMass(par,g, ran); calcInterpCoeffs_lin(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_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); } From 463002a51dde730f7794150cc4ba7ab61b545a76 Mon Sep 17 00:00:00 2001 From: ims Date: Mon, 1 Aug 2016 11:52:37 +0200 Subject: [PATCH 13/13] Fixed small errors. --- Makefile | 3 ++- src/aux.c | 23 +++++++++++++---------- src/frees.c | 1 - src/grid.c | 2 +- src/main.c | 1 - src/raytrace.c | 2 +- 6 files changed, 17 insertions(+), 15 deletions(-) diff --git a/Makefile b/Makefile index 143efc6..966be38 100644 --- a/Makefile +++ b/Makefile @@ -49,7 +49,8 @@ endif ## TARGET = lime.x -CC = gcc -fopenmp +#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 \ diff --git a/src/aux.c b/src/aux.c index 51543de..cd39f9b 100644 --- a/src/aux.c +++ b/src/aux.c @@ -17,7 +17,7 @@ void parseInput(inputPars inpar, configInfo *par, image **img, molData **m){ - int i,id,ispec; + int i,id; double BB[3],normBSquared,dens[MAX_N_COLL_PART]; double cosPhi,sinPhi,cosTheta,sinTheta,dummyVel[DIM]; FILE *fp; @@ -556,7 +556,7 @@ Pointers are indicated by a * before the attribute name and an arrow to the memo void levelPops(molData *m, configInfo *par, struct grid *g, int *popsdone){ - int id,conv=0,iter,ilev,prog=0,ispec,c=0,n,i,threadI,nVerticesDone,nlinetot,numCollParts; + int id,conv=0,iter,ilev,ispec,c=0,n,i,threadI,nVerticesDone,nItersDone,nlinetot,numCollParts; int *allCollPartIds=NULL; double percent=0.,*median,result1=0,result2=0,snr,delta_pop; int nextMolWithBlend; @@ -645,8 +645,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;idncell && !g[id].sink;id++){ for(ilev=0;ilevpIntensity,10); + if (threadI == 0){ /* i.e., is master thread. */ + if(!silent) progressbar(nVerticesDone/(double)par->pIntensity,10); } if(g[id].dens[0] > 0 && g[id].t[0] > 0){ photon(id,g,m,0,threadRans[threadI],par,nlinetot,blends,mp,halfFirstDs); @@ -739,9 +740,11 @@ While this is off however, other gsl_* etc calls will not exit if they encounter } free(median); - if(!silent) progressbar2(1, prog, percent, result1, result2); - if(par->outputfile) popsout(par,g,m); - } while(conv++outputfile != NULL) popsout(par,g,m); + nItersDone++; + } + gsl_set_error_handler(defaultErrorHandler); freeMolsWithBlends(blends.mols, blends.numMolsWithBlends); @@ -762,7 +765,7 @@ While this is off however, other gsl_* etc calls will not exit if they encounter par->dataFlags |= DS_mask_4; - if(par->binoutputfile) binpopsout(par,g,m); + if(par->binoutputfile != NULL) binpopsout(par,g,m); *popsdone=1; } diff --git a/src/frees.c b/src/frees.c index d7cda78..602dd99 100644 --- a/src/frees.c +++ b/src/frees.c @@ -13,7 +13,6 @@ void freeGrid(const unsigned int numPoints, const unsigned short numSpecies\ , struct grid *gp){ unsigned int i_u; - unsigned short i_s; if(gp != NULL){ for(i_u=0;i_unSpecies;i++){ for(id=0;idncell;id++){ diff --git a/src/main.c b/src/main.c index 5a7f3ad..b29f1ac 100644 --- a/src/main.c +++ b/src/main.c @@ -38,7 +38,6 @@ initParImg(inputPars *par, image **img) */ int i,id,nImages; - FILE *fp; /* Set 'impossible' default values for mandatory parameters */ par->radius = 0; diff --git a/src/raytrace.c b/src/raytrace.c index ae679ea..528de6f 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -203,7 +203,7 @@ Note that the algorithm employed here is similar to that employed in the functio void raytrace(int im, configInfo *par, struct grid *gp, molData *md, image *img){ - int aa,ichan,px,iline,cmbMolI,cmbLineI,i,threadI,nRaysDone,molI,lineI; + int aa,ichan,px,cmbMolI,cmbLineI,i,threadI,nRaysDone,molI,lineI; double size,minfreq,absDeltaFreq,totalNumPixelsMinus1=(double)(img[im].pxls*img[im].pxls-1); double cutoff; const gsl_rng_type *ranNumGenType = gsl_rng_ranlxs2;