From c4ffd8c4ae52bb94db6a1af8d558d2cd92e340cf Mon Sep 17 00:00:00 2001 From: Ian Stewart Date: Fri, 24 Jul 2015 12:15:27 +0200 Subject: [PATCH 1/6] New configuration struct. The new struct is typedefd as configInfo. It contains all the information which was previously stored in inputPars: the difference is that inputPars has now been trimmed to contain _only_ those parameters which are settable by the user in the model.c file. Any non-user-settable values should now be added to the config struct. The advantages of this are as follows. - We will usually want more config parameters than user-settable ones. The separation leaves it clearer which things the user needs to (or can) set. - The separation makes it easier to check for and exclude impossible combinations of parameters. - We can adopt new names (for clarity) for config parameters without bothering the user with a changed interface. Two small unrelated additional changes: - A new macro MAX_NIMAGES is defined in lime.h and used in the initial malloc of struct image instead of MAX_NSPECIES as before. - Fixed the order of functions in lime.h so they are all (instead of just mostly, as before) alphabetical. --- src/LTEsolution.c | 2 +- src/aux.c | 237 +++++++++++++++++++++++++++++----------------- src/grid.c | 20 ++-- src/lime.h | 109 ++++++++++++--------- src/main.c | 6 +- src/molinit.c | 4 +- src/photon.c | 8 +- src/popsin.c | 2 +- src/popsout.c | 4 +- src/predefgrid.c | 2 +- src/raytrace.c | 6 +- src/report.c | 2 +- src/smooth.c | 2 +- src/stateq.c | 2 +- src/velospline.c | 4 +- src/weights.c | 2 +- src/writefits.c | 2 +- 17 files changed, 245 insertions(+), 169 deletions(-) diff --git a/src/LTEsolution.c b/src/LTEsolution.c index 874213d..612a73b 100644 --- a/src/LTEsolution.c +++ b/src/LTEsolution.c @@ -10,7 +10,7 @@ #include "lime.h" void -LTE(inputPars *par, struct grid *g, molData *m){ +LTE(configInfo *par, struct grid *g, molData *m){ int id,ilev,ispec; double z; diff --git a/src/aux.c b/src/aux.c index f3d9d39..0751f77 100644 --- a/src/aux.c +++ b/src/aux.c @@ -14,77 +14,97 @@ void -parseInput(inputPars *par, image **img, molData **m){ +readUserInput(inputPars *par, image **img, int *nSpecies, int *nImages){ + /* +Here the user-set parameters (both in the inputPars struct and the image struct) are read. + */ + FILE *fp; int i,id; - double BB[3]; - double cosPhi,sinPhi,cosTheta,sinTheta; - /* Set default values */ + /* Set 'impossible' default values for mandatory parameters */ + par->radius =-1; + par->minScale =-1; + par->pIntensity=-1; + par->sinkPoints=-1; + + /* Set default values for optional parameters */ + par->sampling=2; + par->tcmb = 2.728; + par->moldatfile=malloc(sizeof(char *) * MAX_NSPECIES); + for(id=0;idmoldatfile[id]=NULL; + } par->dust = NULL; - par->inputfile = NULL; par->outputfile = NULL; par->binoutputfile= NULL; + par->restart = NULL; par->gridfile = NULL; par->pregrid = NULL; - par->restart = NULL; - - par->tcmb = 2.728; par->lte_only=0; - par->init_lte=0; - par->sampling=2; par->blend=0; par->antialias=1; par->polarization=0; - par->pIntensity=0; - par->sinkPoints=0; - par->doPregrid=0; par->nThreads=0; /* Allocate space for output fits images */ - (*img)=malloc(sizeof(image)*MAX_NSPECIES); - par->moldatfile=malloc(sizeof(char *) * MAX_NSPECIES); - for(id=0;idmoldatfile[id]=NULL; + (*img)=malloc(sizeof(image)*MAX_NIMAGES); + for(i=0;inImages=id; - if(par->nImages==0) { - if(!silent) bail_out("Error: no images defined"); + + /* Check that the mandatory parameters now have 'sensible' settings (i.e., that they have been set at all). Raise an exception if not. */ + if (par->radius<=0){ + if(!silent) bail_out("Error: you must define the radius parameter."); + exit(1); + } + if (par->minScale<=0){ + if(!silent) bail_out("Error: you must define the minScale parameter."); + exit(1); + } + if (par->pIntensity<=0){ + if(!silent) bail_out("Error: you must define the pIntensity parameter."); + exit(1); + } + if (par->sinkPoints<=0){ + if(!silent) bail_out("Error: you must define the sinkPoints parameter."); exit(1); } - *img=realloc(*img, sizeof(image)*par->nImages); - - id=-1; - while(par->moldatfile[++id]!=NULL); - par->nSpecies=id; - if( par->nSpecies == 0 ) - { - par->nSpecies = 1; - free(par->moldatfile); - par->moldatfile = NULL; - } - else - { - par->moldatfile=realloc(par->moldatfile, sizeof(char *)*par->nSpecies); - /* Check if files exists */ - for(id=0;idnSpecies;id++){ - if((fp=fopen(par->moldatfile[id], "r"))==NULL) { - openSocket(par, id); - } - else { - fclose(fp); - } + /* If the user has provided a list of moldatfile names, the corresponding elements of par->moldatfile will be non-NULL. Thus we can deduce the number of files (species) from the number of non-NULL elements. */ + *nSpecies=-1; + while(par->moldatfile[++(*nSpecies)]!=NULL); + if( *nSpecies == 0 ){ + *nSpecies = 1; + free(par->moldatfile); + par->moldatfile = NULL; + } else { + par->moldatfile=realloc(par->moldatfile, sizeof(char *)* (*nSpecies)); + /* Check if files exist */ + for(id=0;id<(*nSpecies);id++){ + if((fp=fopen(par->moldatfile[id], "r"))==NULL) { + openSocket(par, id); + } else { + fclose(fp); } } + } + + /* If the user has provided a list of image filenames, the corresponding elements of (*img).filename will be non-NULL. Thus we can deduce the number of images from the number of non-NULL elements. */ + *nImages=-1; + while((*img)[++(*nImages)].filename!=NULL); + if(*nImages==0) { + if(!silent) bail_out("Error: no images defined"); + exit(1); + } + *img=realloc(*img, sizeof(image)*(*nImages)); - /* Set defaults and read inputPars and img[] */ - for(i=0;inImages;i++) { + /* Set img defaults. */ + for(i=0;i<(*nImages);i++) { (*img)[i].source_vel=0.0; (*img)[i].phi=0.0; (*img)[i].nchan=0; @@ -93,16 +113,61 @@ parseInput(inputPars *par, image **img, molData **m){ (*img)[i].freq=-1.; (*img)[i].bandwidth=-1.; } + + /* Second-pass reading of the user-set parameters (this time just to read the par->moldatfile and img stuff). */ input(par,*img); - if(par->nThreads == 0){ // Hmm. Really ought to have a separate boolean parameter. +} + +void +setUpConfig(configInfo *par, image **img, molData **m){ + /* +Most of the work this does is setting values in the configInfo struct which contains (as one might guess) information relevant to the general configuration of the task. Non-user-settable values in the image struct are also set here; finally (and somewhat anomalously), the molData vector is initialized. + */ + + int i,id; + double BB[3]; + double cosPhi,sinPhi,cosTheta,sinTheta; + int nSpecies,nImages; + inputPars inpar; + + readUserInput(&inpar,img,&nSpecies,&nImages); + + /* Copy over user-set parameters to the configInfo versions. (This seems like duplicated effort but it is a good principle to separate the two structs, for several reasons, as follows. (i) We will usually want more config parameters than user-settable ones. The separation leaves it clearer which things the user needs to (or can) set. (ii) The separation allows checking and screening out of impossible combinations of parameters. (iii) We can adopt new names (for clarity) for config parameters without bothering the user with a changed interface.) */ + par->radius = inpar.radius; + par->minScale = inpar.minScale; + par->pIntensity = inpar.pIntensity; + par->sinkPoints = inpar.sinkPoints; + par->sampling = inpar.sampling; + par->tcmb = inpar.tcmb; + par->moldatfile = malloc(sizeof(char*)*nSpecies); + for(id=0;idmoldatfile[id] = inpar.moldatfile[id]; + } + par->dust = inpar.dust; + par->outputfile = inpar.outputfile; + par->binoutputfile= inpar.binoutputfile; + par->restart = inpar.restart; + par->gridfile = inpar.gridfile; + par->pregrid = inpar.pregrid; + par->lte_only = inpar.lte_only; + par->blend = inpar.blend; + par->antialias = inpar.antialias; + par->polarization = inpar.polarization; + + if(inpar.nThreads == 0){ par->nThreads = NTHREADS; + } else { + par->nThreads = inpar.nThreads; } - par->ncell=par->pIntensity+par->sinkPoints; - par->radiusSqu=par->radius*par->radius; - par->minScaleSqu=par->minScale*par->minScale; - if(par->pregrid!=NULL) par->doPregrid=1; + /* Now set the additional values in par. */ + par->ncell = inpar.pIntensity + inpar.sinkPoints; + par->radiusSqu = inpar.radius*inpar.radius; + par->minScaleSqu=inpar.minScale*inpar.minScale; + par->doPregrid = (inpar.pregrid==NULL)?0:1; + par->nSpecies = nSpecies; + par->nImages = nImages; /* Now we need to calculate the cutoff value used in calcSourceFn(). The issue is to decide between @@ -120,16 +185,6 @@ The cutoff will be the value of abs(x) for which the error in the exact expressi */ par->taylorCutoff = pow(24.*DBL_EPSILON, 0.25); - if(par->dust != NULL){ - if((fp=fopen(par->dust, "r"))==NULL){ - if(!silent) bail_out("Error opening dust opacity data 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 ){ @@ -137,12 +192,12 @@ The cutoff will be the value of abs(x) for which the error in the exact expressi /* Check for polarization */ BB[0]=0.; - magfield(par->minScale,par->minScale,par->minScale,BB); + magfield(inpar.minScale,inpar.minScale,inpar.minScale,BB); if(fabs(BB[0]) > 0.) par->polarization=1; if(par->polarization) (*img)[i].nchan=3; else (*img)[i].nchan=1; - if((*img)[i].trans>-1 || (*img)[i].bandwidth>-1. || (*img)[i].freq==0 || par->dust==NULL){ + if((*img)[i].trans>-1 || (*img)[i].bandwidth>-1. || (*img)[i].freq==0 || inpar.dust==NULL){ if(!silent) bail_out("Error: Image keywords are ambiguous"); exit(1); } @@ -150,7 +205,7 @@ The cutoff will be the value of abs(x) for which the error in the exact expressi } else if (((*img)[i].nchan>0 || (*img)[i].velres > 0)){ /* Assume line image */ par->polarization=0; - if(par->moldatfile==NULL){ + if(inpar.moldatfile==NULL){ if(!silent) bail_out("Error: No data file is specified for line image."); exit(1); } @@ -204,28 +259,32 @@ The cutoff will be the value of abs(x) for which the error in the exact expressi /* Allocate moldata array */ (*m)=malloc(sizeof(molData)*par->nSpecies); - for( i=0; inSpecies; i++ ) - { - (*m)[i].ntrans = NULL; - (*m)[i].lal = NULL; - (*m)[i].lau = NULL; - (*m)[i].lcl = NULL; - (*m)[i].lcu = NULL; - (*m)[i].aeinst = NULL; - (*m)[i].freq = NULL; - (*m)[i].beinstu = NULL; - (*m)[i].beinstl = NULL; - (*m)[i].up = NULL; - (*m)[i].down = NULL; - (*m)[i].eterm = NULL; - (*m)[i].gstat = NULL; - (*m)[i].cmb = NULL; - (*m)[i].local_cmb = NULL; - } + for( i=0; inSpecies; i++ ){ + (*m)[i].ntrans = NULL; + (*m)[i].lal = NULL; + (*m)[i].lau = NULL; + (*m)[i].lcl = NULL; + (*m)[i].lcu = NULL; + (*m)[i].aeinst = NULL; + (*m)[i].freq = NULL; + (*m)[i].beinstu = NULL; + (*m)[i].beinstl = NULL; + (*m)[i].up = NULL; + (*m)[i].down = NULL; + (*m)[i].eterm = NULL; + (*m)[i].gstat = NULL; + (*m)[i].cmb = NULL; + (*m)[i].local_cmb = NULL; + } + + if( inpar.moldatfile != NULL ){ + free(inpar.moldatfile); + } } + void -freeInput( inputPars *par, image* img, molData* mol ) +freeInput( configInfo *par, image* img, molData* mol ) { int i,id; if( mol!= 0 ) @@ -313,7 +372,7 @@ freeInput( inputPars *par, image* img, molData* mol ) } void -freeGridPointData(inputPars *par, gridPointData *mol){ +freeGridPointData(configInfo *par, gridPointData *mol){ int i; if (mol!= 0){ for (i=0;inSpecies;i++){ @@ -345,7 +404,7 @@ invSqrt(float x){ } void -continuumSetup(int im, image *img, molData *m, inputPars *par, struct grid *g){ +continuumSetup(int im, image *img, molData *m, configInfo *par, struct grid *g){ int id; img[im].trans=0; m[0].nline=1; @@ -386,7 +445,7 @@ lineCount(int n,molData *m,int **counta,int **countb,int *nlinetot){ } void -lineBlend(molData *m, inputPars *par, blend **matrix){ +lineBlend(molData *m, configInfo *par, blend **matrix){ int iline, jline, nlinetot=0,c; int *counta,*countb; @@ -426,7 +485,7 @@ lineBlend(molData *m, inputPars *par, blend **matrix){ } void -levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ +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; double percent=0.,*median,result1=0,result2=0,snr,delta_pop; blend *matrix; @@ -470,7 +529,7 @@ levelPops(molData *m, inputPars *par, struct grid *g, int *popsdone){ /* Check for blended lines */ lineBlend(m,par,&matrix); - if(par->lte_only || par->init_lte) LTE(par,g,m); + if(par->lte_only!=0) LTE(par,g,m); for(id=0;idpIntensity;id++){ stat[id].pop=malloc(sizeof(double)*m[0].nlev*5); diff --git a/src/grid.c b/src/grid.c index 872f301..bdf15a6 100644 --- a/src/grid.c +++ b/src/grid.c @@ -11,7 +11,7 @@ void -gridAlloc(inputPars *par, struct grid **g){ +gridAlloc(configInfo *par, struct grid **g){ int i; double temp[99]; @@ -47,7 +47,7 @@ gridAlloc(inputPars *par, struct grid **g){ } void -freePopulation(const inputPars *par, const molData* m, struct populations* pop ) { +freePopulation(const configInfo *par, const molData* m, struct populations* pop ) { if( pop !=NULL ) { int j,k; @@ -88,7 +88,7 @@ freePopulation(const inputPars *par, const molData* m, struct populations* pop ) } } void -freeGrid(const inputPars *par, const molData* m ,struct grid* g){ +freeGrid(const configInfo *par, const molData* m ,struct grid* g){ int i; if( g != NULL ) { @@ -151,7 +151,7 @@ freeGrid(const inputPars *par, const molData* m ,struct grid* g){ } void -qhull(inputPars *par, struct grid *g){ +qhull(configInfo *par, struct grid *g){ int i,j,k,id; char flags[255]; boolT ismalloc = False; @@ -224,7 +224,7 @@ qhull(inputPars *par, struct grid *g){ } void -distCalc(inputPars *par, struct grid *g){ +distCalc(configInfo *par, struct grid *g){ int i,k,l; for(i=0;incell;i++){ @@ -251,7 +251,7 @@ distCalc(inputPars *par, struct grid *g){ void -write_VTK_unstructured_Points(inputPars *par, struct grid *g){ +write_VTK_unstructured_Points(configInfo *par, struct grid *g){ FILE *fp; double length; int i,j,l=0; @@ -338,12 +338,12 @@ write_VTK_unstructured_Points(inputPars *par, struct grid *g){ } void -dumpGrid(inputPars *par, struct grid *g){ +dumpGrid(configInfo *par, struct grid *g){ if(par->gridfile) write_VTK_unstructured_Points(par, g); } void -getArea(inputPars *par, struct grid *g, const gsl_rng *ran){ +getArea(configInfo *par, struct grid *g, const gsl_rng *ran){ int i,j,k,b;//=-1; double *angle,best; /* double wsum; */ @@ -406,7 +406,7 @@ getArea(inputPars *par, struct grid *g, const gsl_rng *ran){ void -getMass(inputPars *par, struct grid *g, const gsl_rng *ran){ +getMass(configInfo *par, struct grid *g, const gsl_rng *ran){ double mass=0.,dist; double vol=0.,dp,dpbest,*farea,suma; int i,k,j,best=-1; @@ -515,7 +515,7 @@ getMass(inputPars *par, struct grid *g, const gsl_rng *ran){ void -buildGrid(inputPars *par, struct grid *g){ +buildGrid(configInfo *par, struct grid *g){ double lograd; /* The logarithm of the model radius */ double logmin; /* Logarithm of par->minScale */ double r,theta,phi,sinPhi,x,y,z,semiradius; /* Coordinates */ diff --git a/src/lime.h b/src/lime.h index f1d602c..8bd762a 100644 --- a/src/lime.h +++ b/src/lime.h @@ -73,6 +73,7 @@ #define fixset 1e-6 #define blendmask 1.e4 #define MAX_NSPECIES 100 +#define MAX_NIMAGES 100 #define N_RAN_PER_SEGMENT 3 #define FAST_EXP_MAX_TAYLOR 3 #define FAST_EXP_NUM_BITS 8 @@ -80,17 +81,32 @@ /* input parameters */ typedef struct { - double radius,radiusSqu,minScale,minScaleSqu,tcmb,taylorCutoff; - int ncell,sinkPoints,pIntensity,nImages,nSpecies,blend; - char *outputfile, *binoutputfile, *inputfile; + double radius,minScale,tcmb; + int sinkPoints,pIntensity,blend; + char *outputfile, *binoutputfile; +// char *inputfile; unused at present. char *gridfile; char *pregrid; char *restart; char *dust; - int sampling,collPart,lte_only,init_lte,antialias,polarization,doPregrid,nThreads; + int sampling,lte_only,antialias,polarization,nThreads; char **moldatfile; } inputPars; +typedef struct { + double radius,minScale,tcmb; + double radiusSqu,minScaleSqu,taylorCutoff; + int sinkPoints,pIntensity,blend; + int ncell,nImages,nSpecies,collPart,doPregrid; + char *outputfile, *binoutputfile; + char *gridfile; + char *pregrid; + char *restart; + char *dust; + int sampling,lte_only,antialias,polarization,nThreads; + char **moldatfile; +} configInfo; + /* Molecular data: shared attributes */ typedef struct { int nlev,nline,*ntrans,npart; @@ -175,7 +191,7 @@ typedef struct {double x,y, *intensity, *tau;} rayData; -/* Some functions */ +/* User-specifiable functions */ void density(double,double,double,double *); void temperature(double,double,double,double *); void abundance(double,double,double,double *); @@ -186,67 +202,68 @@ void gasIIdust(double,double,double,double *); /* More functions */ -void binpopsout(inputPars *, struct grid *, molData *); -void buildGrid(inputPars *, struct grid *); -void calcSourceFn(double dTau, const inputPars *par, double *remnantSnu, double *expDTau); -void continuumSetup(int, image *, molData *, inputPars *, struct grid *); -void distCalc(inputPars *, struct grid *); +void binpopsout(configInfo *, struct grid *, molData *); +void buildGrid(configInfo *, struct grid *); +void calcFastExpRange(const int, const int, int *, int *, int *); +void calcSourceFn(double, const configInfo *, double *, double *); +void calcTableEntries(const int maxTaylorOrder, const int maxNumBitsPerMantField); +void continuumSetup(int, image *, molData *, configInfo *, struct grid *); +void distCalc(configInfo *, struct grid *); +int factorial(const int n); +inline double FastExp(const float negarg); void fit_d1fi(double, double, double*); void fit_fi(double, double, double*); void fit_rr(double, double, double*); -void input(inputPars *, image *); -float invSqrt(float); -void freeInput(inputPars *, image*, molData* m ); -void freeGrid(const inputPars * par, const molData* m, struct grid * g); -void freePopulation(const inputPars * par, const molData* m, struct populations * pop); +void freeInput(configInfo *, image*, molData* m ); +void freeGrid(const configInfo * par, const molData* m, struct grid * g); +void freePopulation(const configInfo * 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 getArea(configInfo *, struct grid *, const gsl_rng *); +void getjbar(int, molData *, struct grid *, configInfo *,gridPointData *,double *); +void getMass(configInfo *, 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 kappa(molData *, struct grid *, inputPars *,int); -void levelPops(molData *, inputPars *, struct grid *, int *); +void getVelosplines(configInfo *, struct grid *); +void getVelosplines_lin(configInfo *, struct grid *); +void gridAlloc(configInfo *, struct grid **); +void input(inputPars *, image *); +float invSqrt(float); +void kappa(molData *, struct grid *, configInfo *,int); +void levelPops(molData *, configInfo *, struct grid *, int *); void line_plane_intersect(struct grid *, double *, int , int *, double *, double *, double); -void lineBlend(molData *, inputPars *, blend **); +void lineBlend(molData *, configInfo *, blend **); void lineCount(int,molData *,int **, int **, int *); -void LTE(inputPars *, struct grid *, molData *); -void molinit(molData *, inputPars *, struct grid *,int); -void openSocket(inputPars *par, int); -void qhull(inputPars *, struct grid *); -void photon(int, struct grid *, molData *, int, const gsl_rng *,inputPars *,blend *,gridPointData *,double *); -void parseInput(inputPars *, image **, molData **); +void LTE(configInfo *, struct grid *, molData *); +void molinit(molData *, configInfo *, struct grid *,int); +void openSocket(inputPars *, int); +void qhull(configInfo *, struct grid *); +void photon(int, struct grid *, molData *, int, const gsl_rng *,configInfo *,blend *,gridPointData *,double *); 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 *); +int pointEvaluation(configInfo *,double, double, double, double); +void popsin(configInfo *, struct grid **, molData **, int *); +void popsout(configInfo *, struct grid *, molData *); +void predefinedGrid(configInfo *, struct grid *); double ratranInput(char *, char *, double, double, double); -void raytrace(int, inputPars *, struct grid *, molData *, image *); -void report(int, inputPars *, struct grid *); -void smooth(inputPars *, struct grid *); +void raytrace(int, configInfo *, struct grid *, molData *, image *); +void readUserInput(inputPars *, image **, int *, int *); +void report(int, configInfo *, struct grid *); +void setUpConfig(configInfo *, image **, molData **); +void smooth(configInfo *, 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 stateq(int, struct grid *, molData *, int, configInfo *,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 taylor(const int maxOrder, const float x); +void traceray(rayData, int, int, configInfo *, 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 writefits(int, inputPars *, molData *, image *); -void write_VTK_unstructured_Points(inputPars *, struct grid *); -int factorial(const int n); -double taylor(const int maxOrder, const float x); -void calcFastExpRange(const int maxTaylorOrder, const int maxNumBitsPerMantField, int *numMantissaFields, int *lowestExponent, int *numExponentsUsed); -void calcTableEntries(const int maxTaylorOrder, const int maxNumBitsPerMantField); -inline double FastExp(const float negarg); +void writefits(int, configInfo *, molData *, image *); +void write_VTK_unstructured_Points(configInfo *, struct grid *); /* Curses functions */ diff --git a/src/main.c b/src/main.c index aa92ea9..0b55fb6 100644 --- a/src/main.c +++ b/src/main.c @@ -21,11 +21,11 @@ double EXP_TABLE_3D[1][1][1]; #endif int main () { - int i; + int i,nSpecies,nImages; int initime=time(0); int popsdone=0; molData* m = NULL; - inputPars par; + configInfo par; struct grid* g = NULL; image* img = NULL; @@ -36,7 +36,7 @@ int main () { calcTableEntries(FAST_EXP_MAX_TAYLOR, FAST_EXP_NUM_BITS); #endif - parseInput(&par,&img,&m); + setUpConfig(&par,&img,&m); if(par.doPregrid) { diff --git a/src/molinit.c b/src/molinit.c index d48607e..442152f 100644 --- a/src/molinit.c +++ b/src/molinit.c @@ -10,7 +10,7 @@ #include "lime.h" void -kappa(molData *m, struct grid *g, inputPars *par, int s){ +kappa(molData *m, struct grid *g, configInfo *par, int s){ FILE *fp; char string[80]; int i=0,k,j,iline,id; @@ -102,7 +102,7 @@ planckfunc(int iline, double temp, molData *m,int s){ } void -molinit(molData *m, inputPars *par, struct grid *g,int i){ +molinit(molData *m, configInfo *par, struct grid *g, int i){ int id, ilev, iline, itrans, ispec, itemp, *ntemp, tnint=-1, idummy, ipart, *count,flag=0; char *collpartname[] = {"H2","p-H2","o-H2","electrons","H","He","H+"}; /* definition from LAMDA */ double fac, uprate, downrate=0, dummy, amass; diff --git a/src/photon.c b/src/photon.c index 4533870..d56da63 100644 --- a/src/photon.c +++ b/src/photon.c @@ -131,14 +131,14 @@ double gaussline(double v, double oneOnSigma){ } -void calcSourceFn(double dTau, const inputPars *par, double *remnantSnu, double *expDTau){ +void calcSourceFn(double dTau, const configInfo *par, double *remnantSnu, double *expDTau){ /* The source function S is defined as j_nu/alpha, which is clearly not defined for alpha==0. However S is used in the algorithm only in the term (1-exp[-alpha*ds])*S, which is defined for all values of alpha. The present function calculates this term and returns it in the argument remnantSnu. For values of abs(alpha*ds) less than a pre- - calculated cutoff supplied in inputPars, a Taylor approximation is + calculated cutoff supplied in configInfo, a Taylor approximation is used. Note that the same cutoff condition holds for replacement of @@ -165,7 +165,7 @@ void calcSourceFn(double dTau, const inputPars *par, double *remnantSnu, double void -photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPars *par,blend *matrix, gridPointData *mp, double *halfFirstDs){ +photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,configInfo *par,blend *matrix, gridPointData *mp, double *halfFirstDs){ 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; @@ -300,7 +300,7 @@ photon(int id, struct grid *g, molData *m, int iter, const gsl_rng *ran,inputPar } void -getjbar(int posn, molData *m, struct grid *g, inputPars *par, gridPointData *mp, double *halfFirstDs){ +getjbar(int posn, molData *m, struct grid *g, configInfo *par, gridPointData *mp, double *halfFirstDs){ int iline,iphot; double tau, expTau, remnantSnu, vsum=0., jnu, alpha; int *counta, *countb,nlinetot; diff --git a/src/popsin.c b/src/popsin.c index 5392664..b8326f2 100644 --- a/src/popsin.c +++ b/src/popsin.c @@ -11,7 +11,7 @@ void -popsin(inputPars *par, struct grid **g, molData **m, int *popsdone){ +popsin(configInfo *par, struct grid **g, molData **m, int *popsdone){ FILE *fp; int i,j,k; double dummy; diff --git a/src/popsout.c b/src/popsout.c index 93ff3e3..7152564 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(configInfo *par, struct grid *g, molData *m){ FILE *fp; int j,k,l; double dens; @@ -35,7 +35,7 @@ popsout(inputPars *par, struct grid *g, molData *m){ void -binpopsout(inputPars *par, struct grid *g, molData *m){ +binpopsout(configInfo *par, struct grid *g, molData *m){ FILE *fp; int i,j; diff --git a/src/predefgrid.c b/src/predefgrid.c index 4aed807..992195f 100644 --- a/src/predefgrid.c +++ b/src/predefgrid.c @@ -10,7 +10,7 @@ #include "lime.h" void -predefinedGrid(inputPars *par, struct grid *g){ +predefinedGrid(configInfo *par, struct grid *g){ FILE *fp; int i; double x,y,z,scale; diff --git a/src/raytrace.c b/src/raytrace.c index a0ca190..62f1d48 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -62,7 +62,7 @@ line_plane_intersect(struct grid *g, double *ds, int posn, int *nposn, double *d 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, configInfo *par, struct grid *g, molData *m, image *img, int nlinetot, int *counta, int *countb, double cutoff){ int ichan,posn,nposn,i,iline; double vfac=0.,x[3],dx[3],vThisChan; double deltav,ds,dist2,ndist2,xp,yp,zp,col,shift,jnu,alpha,remnantSnu,dtau,expDTau,snu_pol[3]; @@ -169,7 +169,7 @@ traceray(rayData ray, int tmptrans, int im, inputPars *par, struct grid *g, molD void -raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ +raytrace(int im, configInfo *par, struct grid *g, molData *m, image *img){ int *counta, *countb,nlinetot,aa; int ichan,px,iline,tmptrans,i,threadI,nRaysDone; double size,minfreq,absDeltaFreq,totalNumPixelsMinus1=(double)(img[im].pxls*img[im].pxls-1); @@ -280,7 +280,7 @@ raytrace(int im, inputPars *par, struct grid *g, molData *m, image *img){ void -raytrace_1_4(int im, inputPars *par, struct grid *g, molData *m, image *img){ +raytrace_1_4(int im, configInfo *par, struct grid *g, molData *m, image *img){ /* This is an alternative raytracing algorithm which was implemented by C Brinch in version 1.4 (the original parallelized version) of LIME. diff --git a/src/report.c b/src/report.c index bd8c6e0..46211ac 100644 --- a/src/report.c +++ b/src/report.c @@ -13,7 +13,7 @@ void -report(int i, inputPars *par, struct grid *g){ +report(int i, configInfo *par, struct grid *g){ FILE *fp; int j,k,p,q; const int bins = 100; diff --git a/src/smooth.c b/src/smooth.c index 1e98310..df3911b 100644 --- a/src/smooth.c +++ b/src/smooth.c @@ -12,7 +12,7 @@ /* Based on Lloyds Algorithm (Lloyd, S. IEEE, 1982) */ void -smooth(inputPars *par, struct grid *g){ +smooth(configInfo *par, struct grid *g){ double mindist; /* Distance to closest neighbor */ int k=0,j,i; /* counters */ int sg; /* counter for smoothing the grid */ diff --git a/src/stateq.c b/src/stateq.c index 3588864..4896c30 100644 --- a/src/stateq.c +++ b/src/stateq.c @@ -13,7 +13,7 @@ void -stateq(int id, struct grid *g, molData *m, int ispec, inputPars *par, gridPointData *mp, double *halfFirstDs){ +stateq(int id, struct grid *g, molData *m, int ispec, configInfo *par, gridPointData *mp, double *halfFirstDs){ int t,s,iter; double *opop, *oopop; double diff; diff --git a/src/velospline.c b/src/velospline.c index 3e063bd..0aea5ec 100644 --- a/src/velospline.c +++ b/src/velospline.c @@ -10,7 +10,7 @@ #include "lime.h" void -getVelosplines(inputPars *par, struct grid *g){ +getVelosplines(configInfo *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); @@ -76,7 +76,7 @@ getVelosplines(inputPars *par, struct grid *g){ void -getVelosplines_lin(inputPars *par, struct grid *g){ +getVelosplines_lin(configInfo *par, struct grid *g){ int i,k,j; double v[2]; diff --git a/src/weights.c b/src/weights.c index f4de61a..c0be334 100644 --- a/src/weights.c +++ b/src/weights.c @@ -10,7 +10,7 @@ #include "lime.h" int -pointEvaluation(inputPars *par,double ran, double x, double y, double z){ +pointEvaluation(configInfo *par,double ran, double x, double y, double z){ double weight1, weight2, val[99],normalizer=0.0,totalDensity=0.0; int i; diff --git a/src/writefits.c b/src/writefits.c index 2ba9971..0c77fe8 100644 --- a/src/writefits.c +++ b/src/writefits.c @@ -10,7 +10,7 @@ #include "lime.h" void -writefits(int im, inputPars *par, molData *m, image *img){ +writefits(int im, configInfo *par, molData *m, image *img){ double bscale,bzero,epoch,lonpole,equinox,restfreq; double cdelt1,crpix1,crval1,cdelt2,crpix2,crval2; double cdelt3,crpix3,crval3,ru3; From 430f1514bedc35fa3eb44cabdb7957ceaf15f5ff Mon Sep 17 00:00:00 2001 From: ims Date: Tue, 24 Nov 2015 17:49:37 +0100 Subject: [PATCH 2/6] Removed reference to init_lte from documantation This input parameter had no present function in the code, so was removed from there in the previous commit. Now I also removed mention of it from the user manual. --- doc/usermanual.rst | 7 ------- 1 file changed, 7 deletions(-) diff --git a/doc/usermanual.rst b/doc/usermanual.rst index d50dfbf..e6c9578 100644 --- a/doc/usermanual.rst +++ b/doc/usermanual.rst @@ -417,13 +417,6 @@ abundance, velocity etc. There is no default value. If set, LIME performs an LTE calculation. Useful for quick checks. The default lte\_only=0, i.e., full non-LTE calculation. -.. code:: c - - (integer) par->init_lte (optional) - -If set, LIME use LTE approximation as initial one for subsequent non-LTE calculations. The -default init\_lte=0, i.e., the code will use constant value for level populations as initial solution. - .. code:: c (integer) par->blend (optional) From fcd7fc6a745d43a778a2a7481f6cb2ea99cde99e Mon Sep 17 00:00:00 2001 From: ims Date: Tue, 3 May 2016 15:09:55 +0200 Subject: [PATCH 3/6] Fixed mistake made in manual merging. --- src/lime.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lime.h b/src/lime.h index 8bd762a..e97656e 100644 --- a/src/lime.h +++ b/src/lime.h @@ -210,7 +210,7 @@ void calcTableEntries(const int maxTaylorOrder, const int maxNumBitsPerMantField void continuumSetup(int, image *, molData *, configInfo *, struct grid *); void distCalc(configInfo *, struct grid *); int factorial(const int n); -inline double FastExp(const float negarg); +double FastExp(const float negarg); void fit_d1fi(double, double, double*); void fit_fi(double, double, double*); void fit_rr(double, double, double*); From 0873706e209b95f65a3f20ed74b0add144ec16c0 Mon Sep 17 00:00:00 2001 From: ims Date: Thu, 23 Jun 2016 10:02:19 +0200 Subject: [PATCH 4/6] Added header file for inpars Following advice from @rschaaf I have created a separate header file inpars.h just to define the struct for the user-input parameters. I also included macros in both lime.h and inpars.h to prevent multiple inclusion of these files. The final change was to the interface of openSocket(). This had the inpars struct as one of its arguments, the other being the species ID, but it turns out that all it needs is the name of the moldatfile for that species. Thus I replaced the former two arguements by a single char * moldatfile. --- .gitignore | 6 ++++-- src/aux.c | 2 +- src/inpars.h | 18 ++++++++++++++++++ src/lime.h | 13 +++++++++---- src/tcpsocket.c | 16 ++++++++-------- 5 files changed, 40 insertions(+), 15 deletions(-) create mode 100644 src/inpars.h diff --git a/.gitignore b/.gitignore index 182c93e..d929ee0 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,8 @@ .DS_Store # Ignore example files -example/grid.vtk +example/*.vtk example/*.pop -example/*.fits +example/*.fits* +example/*.dat + diff --git a/src/aux.c b/src/aux.c index 9ba8d8d..b73b593 100644 --- a/src/aux.c +++ b/src/aux.c @@ -86,7 +86,7 @@ Here the user-set parameters (both in the inputPars struct and the image struct) /* Check if files exist */ for(id=0;id<(*nSpecies);id++){ if((fp=fopen(par->moldatfile[id], "r"))==NULL) { - openSocket(par, id); + openSocket(par->moldatfile[id]); } else { fclose(fp); } diff --git a/src/inpars.h b/src/inpars.h new file mode 100644 index 0000000..99003ea --- /dev/null +++ b/src/inpars.h @@ -0,0 +1,18 @@ +#ifndef INPARS_H +#define INPARS_H + +/* input parameters */ +typedef struct { + double radius,minScale,tcmb; + int sinkPoints,pIntensity,blend; + char *outputfile, *binoutputfile; +// char *inputfile; unused at present. + char *gridfile; + char *pregrid; + char *restart; + char *dust; + int sampling,lte_only,antialias,polarization,nThreads; + char **moldatfile; +} inputPars; + +#endif /* INPARS_H */ diff --git a/src/lime.h b/src/lime.h index e97656e..688ccc0 100644 --- a/src/lime.h +++ b/src/lime.h @@ -7,6 +7,11 @@ * */ +#ifndef LIME_H +#define LIME_H + +#include "inpars.h" + #include #include #include @@ -79,7 +84,7 @@ #define FAST_EXP_NUM_BITS 8 -/* input parameters */ +/* input parameters typedef struct { double radius,minScale,tcmb; int sinkPoints,pIntensity,blend; @@ -92,7 +97,7 @@ typedef struct { int sampling,lte_only,antialias,polarization,nThreads; char **moldatfile; } inputPars; - + */ typedef struct { double radius,minScale,tcmb; double radiusSqu,minScaleSqu,taylorCutoff; @@ -235,7 +240,7 @@ void lineBlend(molData *, configInfo *, blend **); void lineCount(int,molData *,int **, int **, int *); void LTE(configInfo *, struct grid *, molData *); void molinit(molData *, configInfo *, struct grid *,int); -void openSocket(inputPars *, int); +void openSocket(char*); void qhull(configInfo *, struct grid *); void photon(int, struct grid *, molData *, int, const gsl_rng *,configInfo *,blend *,gridPointData *,double *); double planckfunc(int, double, molData *, int); @@ -283,5 +288,5 @@ void collpartmesg(char *, int); void collpartmesg2(char *, int); void collpartmesg3(int, int); - +#endif /* LIME_H */ diff --git a/src/tcpsocket.c b/src/tcpsocket.c index b1a151e..87fa138 100644 --- a/src/tcpsocket.c +++ b/src/tcpsocket.c @@ -16,7 +16,7 @@ int close(int); void -openSocket(inputPars *par, int id){ +openSocket(char *moldatfile){ struct sockaddr_in *remote; int tmpres; int sock; @@ -30,20 +30,20 @@ openSocket(inputPars *par, int id){ FILE *fp; // Check if moldatfile contains .dat - if(strstr(par->moldatfile[id], ".dat") == NULL){ - size_t length = strlen(par->moldatfile[id]); + if(strstr(moldatfile, ".dat") == NULL){ + size_t length = strlen(moldatfile); s = (char*)malloc(sizeof(char) * (length + 5)); - strcpy(s,par->moldatfile[id]); + strcpy(s,moldatfile); strcat(s, ".dat"); s[length+5]='\0'; - par->moldatfile[id]=s; + moldatfile=s; } size_t length1 = strlen(page); - size_t length2 = strlen(par->moldatfile[id]); + size_t length2 = strlen(moldatfile); t = (char*)malloc(sizeof(char) * (length1 + length2 + 1)); strcpy(t,page); - strcat(t, par->moldatfile[id]); + strcat(t, moldatfile); page=t; @@ -88,7 +88,7 @@ openSocket(inputPars *par, int id){ } memset(buf, 0, sizeof(buf)); - if((fp=fopen(par->moldatfile[id], "w"))==NULL) { + if((fp=fopen(moldatfile, "w"))==NULL) { if(!silent) bail_out("Failed to write moldata!"); exit(1); } From 64d9710f292567683c10dc0c5b63235bc21cff49 Mon Sep 17 00:00:00 2001 From: ims Date: Mon, 20 Jun 2016 12:14:42 +0200 Subject: [PATCH 5/6] Fixes #138 --- src/raytrace.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/raytrace.c b/src/raytrace.c index bbf86ca..817d2cd 100644 --- a/src/raytrace.c +++ b/src/raytrace.c @@ -207,7 +207,7 @@ raytrace(int im, configInfo *par, struct grid *g, molData *m, image *img){ for (i=0;inThreads;i++){ threadRans[i] = gsl_rng_alloc(ranNumGenType); - gsl_rng_set(threadRans[i],(int)gsl_rng_uniform(ran)*1e6); + gsl_rng_set(threadRans[i],(int)(gsl_rng_uniform(ran)*1e6)); } size=img[im].distance*img[im].imgres; From 5b7abe043c07ab054fdd47eb85f92ced622f221c Mon Sep 17 00:00:00 2001 From: ims Date: Tue, 12 Jul 2016 18:10:22 +0200 Subject: [PATCH 6/6] junk --- src/main.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.c b/src/main.c index 54a4fe1..f88d84a 100644 --- a/src/main.c +++ b/src/main.c @@ -103,7 +103,7 @@ initParImg(inputPars *par, image **img) nImages++; if(nImages==0) { - if(!silent) bail_out("No images defined (or you haven't set at least the 1st filename)."); + if(!silent) bail_out("No images defined (or you haven't set the 1st filename)."); exit(1); }