Skip to content

Commit

Permalink
Merge pull request #220 from allegroLeiden/dust_weights
Browse files Browse the repository at this point in the history 
New parameter collPartMolWeights added to help calculate dust opacity.
  • Loading branch information
@allegroLeiden
allegroLeiden authored Feb 17, 2017
2 parents bbe0039 + 051093a commit 7654865
Show file tree
Hide file tree
Showing 10 changed files with 582 additions and 483 deletions.
2 changes: 2 additions & 0 deletions Makefile.defs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@

CORESOURCES = \
src/aux.c\
src/collparts.c\
src/fastexp.c\
src/frees.c\
src/getclosest.c\
Expand Down Expand Up @@ -42,6 +43,7 @@ STDSOURCES = \

COREINCLUDES = \
src/dims.h\
src/collparts.h\
src/grid2fits.h\
src/gridio.h\
src/inpars.h\
Expand Down
201 changes: 19 additions & 182 deletions src/aux.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,20 @@
#include "gridio.h" //**** should not need this - fix up these gridio macros!


/*....................................................................*/
double planckfunc(const double freq, const double temp){
double bb=10.,wn;
if(temp<eps) bb = 0.0;
else {
wn=freq/CLIGHT;
if (HPLANCK*freq>100.*KBOLTZ*temp)
bb=2.*HPLANCK*wn*wn*freq*exp(-HPLANCK*freq/KBOLTZ/temp);
else
bb=2.*HPLANCK*wn*wn*freq/(exp(HPLANCK*freq/KBOLTZ/temp)-1);
}
return bb;
}

/*....................................................................*/
void
reportInfsAtOrigin(const int numElements, const double *values, const char *funcName){
Expand Down Expand Up @@ -201,6 +215,8 @@ The parameters visible to the user have now been strictly confined to members of
for(i=0;i<MAX_N_COLL_PART;i++) par->dustWeights[i] = inpar.dustWeights[i];
par->collPartNames = malloc(sizeof(char*)*MAX_N_COLL_PART);
for(i=0;i<MAX_N_COLL_PART;i++) copyInparStr(inpar.collPartNames[i], &(par->collPartNames[i]));
par->collPartMolWeights = malloc(sizeof(double)*MAX_N_COLL_PART);
for(i=0;i<MAX_N_COLL_PART;i++) par->collPartMolWeights[i] = inpar.collPartMolWeights[i];

/* Copy over the grid-density-maximum data and find out how many were set:
*/
Expand Down Expand Up @@ -622,188 +638,7 @@ Eventually I hope readOrBuildGrid() will be unilaterally called within LIME; if
} /* otherwise leave it at NULL - we will not be using it. */
}

void checkUserDensWeights(configInfo *par){
/*
This deals with four user-settable vectors: par->collPartIds, par->nMolWeights, par->dustWeights and par->collPartNames. We have to see if these (optional) parameters were set, do some basic checks on them, and make sure they match the number of density values, which by this time should be stored in par->numDensities.
* par->collPartIds: this list acts as a link between the N density function returns (I'm using here N as shorthand for par->numDensities) and the M collision partner ID integers found in the moldatfiles. This allows us to associate density functions with the collision partner transition rates provided in the moldatfiles.
* par->collPartNames: essentially this has only cosmetic importance since it has no effect on the functioning of LIME, only on the names of the collision partners which are printed to stdout. It's main purpose is to reassure the user who has provided transition rates for a non-LAMDA collision species in their moldatfile that they are actually getting these values and not some mysterious reversion to LAMDA.
The user can specify either, none, or both of these two parameters, with the following effects:
Ids Names Effect
----------------------
0 0 LAMDA coll parts are assumed and the association between the density functions and the moldatfiles is essentially guessed.
0 1 par->collPartIds is constructed to contain integers in a sequence from 1 to N. Naturally the user should write matching collision partner ID integers in their moldatfiles.
1 0 LAMDA coll parts are assumed.
1 1 User will get what they ask for. Elements of par->collPartIds are taken (after subtracting 1) to be indices to par->collPartNames; thus in this situation par->collPartNames may have a different order to the density functions, or have more values than par->numDensities. Any extra values are ignored however, and par->collPartNames must contain a value for each par->collPartIds element.
----------------------
* par->nMolWeights: this list gives the weights to be applied to the N density values when calculating molecular densities from abundances.
* par->dustWeights: ditto for the calculation of dust opacity.
*/
int i,j,numUserSetCPIds,numUserSetNMWs,numUserSetDWs,numUserSetCPNames;
int *uniqueCPIds=NULL;
double sum;
char message[80];

/* Get the number of par->collPartIds set by the user:
*/
i = 0;
while(i<MAX_N_COLL_PART && par->collPartIds[i]>0) i++;
numUserSetCPIds = i;

if(numUserSetCPIds>0){
/* Check that they are unique.
*/
uniqueCPIds = malloc(sizeof(int)*numUserSetCPIds);
for(i=0;i<numUserSetCPIds;i++){
for(j=0;j<i;j++){
if(par->collPartIds[i]==uniqueCPIds[j]){
if(!silent) bail_out("Your list of par.collPartIds is not unique.");
exit(1);
}
}
uniqueCPIds[i] = par->collPartIds[i];
}
free(uniqueCPIds);
}

/* Get the number of par->nMolWeights set by the user:
*/
i = 0;
while(i<MAX_N_COLL_PART && par->nMolWeights[i]>=0.0) i++;
numUserSetNMWs = i;

if(numUserSetNMWs>0){
/* Check that they do not sum to zero.
*/
sum = 0.0;
for(i=0;i<numUserSetNMWs;i++){
sum += par->nMolWeights[i];
}
if(sum<=0.0){
if(!silent) bail_out("At least some of your par.nMolWeights must be non-zero!");
exit(1);
}
}

/* Get the number of par->dustWeights set by the user:
*/
i = 0;
while(i<MAX_N_COLL_PART && par->dustWeights[i]>=0.0) i++;
numUserSetDWs = i;

if(numUserSetDWs>0){
/* Check that they do not sum to zero.
*/
sum = 0.0;
for(i=0;i<numUserSetDWs;i++){
sum += par->dustWeights[i];
}
if(sum<=0.0){
if(!silent) bail_out("At least some of your par.dustWeights must be non-zero!");
exit(1);
}
}

/* Get the number of par->collPartNames set by the user:
*/
i = 0;
while(i<MAX_N_COLL_PART && par->collPartNames[i]!=NULL) i++;
numUserSetCPNames = i;

/*
Re the interaction between par->collPartIds and par->collPartNames: there are three possible scenarios as follows.
- User wishes to use only collision partners from the set in the LAMDA database, in which case there is no need to set par->collPartNames. Note that the integer ID values for the collision partners in the standard LAMDA moldatfiles begin at 1.
- User wishes to use some non-LAMDA collision partners: it is perfectly possible to replace the values in a moldatfile for any LAMDA species with values calculated for another; LIME has no internal knowledge of whether a collision partner is LAMDA or not, it just uses the transition rates you give it. The only use it makes of collision partner names is in the the initial greeting message it prints. If the user wishes to have the correct name(s) printed, this can be transmitted to LIME via the par->collPartNames parameter.
For added convenience, it is allowed in this case that the user can supply just par->collPartNames and not par->collPartIds.
* If this is done, the number of elements of par->collPartNames is required to be the same as the number of density values. A list of par->collPartIds of the same length will be constructed, the elements being {1,2,3...} etc. These ID integers must match the ones supplied by the user in their moldatfiles.
* If the user supplies both par->collPartIds and par->collPartNames, the values in par->collPartIds must (after subtraction of 1) be the indices of that collision partner in par->collPartNames, as well of course as matching the value the user has supplied in the moldatfiles.
*/
if(numUserSetCPNames>0){
if(numUserSetCPNames < par->numDensities){
if(!silent) bail_out("There must be at least 1 value of par.collPartNames for each density() return.");
exit(1);
}

if(numUserSetCPIds<=0){
for(i=0;i<par->numDensities;i++)
par->collPartIds[i] = i+1;
numUserSetCPIds = par->numDensities;
}else{
for(i=0;i<numUserSetCPIds;i++){
if(par->collPartIds[i]<1 || par->collPartIds[i]>numUserSetCPNames){
if(!silent){
sprintf(message, "par.collPartIds[%d] value %d is out of range (<1 or >%d)", i, par->collPartIds[i], numUserSetCPNames);
bail_out(message);
}
exit(1);
}
}
}
}

/* Check that we have either 0 par->collPartIds or the same number as the number of density values. If not, the par->collPartIds values the user set will be thrown away, the pointer will be reallocated, and new values will be written to it in setUpDensityAux(), taken from the values in the moldatfiles.
*/
if(numUserSetCPIds != par->numDensities){
free(par->collPartIds);
par->collPartIds = NULL;
/* Note that in the present case we will (for a line-emission image) look for the collision partners listed in the moldatfiles and set par->collPartIds from them. For that to happen, we require the number of collision partners found in the files to equal par->numDensities. */

/* numUserSetCPIds==0 is ok, this just means the user has not set the parameter at all, but for other values we should issue some warnings, because if the user sets any at all, they should set the same number as there are returns from density():
*/
if(numUserSetCPIds > 0)
if(!silent) warning("par.collPartIds will be ignored - there should be 1 for each density.");
}else{
par->collPartIds = realloc(par->collPartIds, sizeof(*(par->collPartIds))*par->numDensities);
}

/* Check if we have either 0 par->nMolWeights or the same number as the number of density values.
*/
if(numUserSetNMWs != par->numDensities){
free(par->nMolWeights);
par->nMolWeights = NULL;
/* Note that in the present case we will (for a line-emission image) look for the collision partners listed in the moldatfiles and set par->nMolWeights from them. */

/* numUserSetNMWs==0 is ok, this just means the user has not set the parameter at all, but for other values we should issue some warnings, because if the user sets any at all, they should set the same number as there are returns from density():
*/
if(numUserSetNMWs > 0)
if(!silent) warning("par->nMolWeights will be ignored - there should be 1 for each density() return.");
}else{
par->nMolWeights = realloc(par->nMolWeights, sizeof(*(par->nMolWeights))*par->numDensities);
}

/* Check if we have either 0 par->dustWeights or the same number as the number of density values. Note that the treatment of the dust weights is stricter, since we need knu for the continuum case, in which we may not have access to collision partner information from moldat files.
*/
if(numUserSetDWs != par->numDensities){
if(numUserSetDWs == 0){
/* This is ok, this just means the user has not set the parameter at all. Revert to the previous algorithm, but with a warning, because the previous algorithm is dangerous.
*/
par->dustWeights = realloc(par->dustWeights, sizeof(*(par->dustWeights))*par->numDensities);
par->dustWeights[0] = 1.0;
for(i=1;i<par->numDensities;i++)
par->dustWeights[i] = 0.0;

if(!silent) warning("User didn't set par.dustWeights. Using the first density to calculate k_nu.");

}else{
if(!silent) bail_out("There must be 1 value of par.dustWeights for each density() return.");
exit(1);
}
}else{
par->dustWeights = realloc(par->dustWeights, sizeof(*(par->dustWeights))*par->numDensities);
}
}

/*....................................................................*/
float
invSqrt(float x){
/* The magic Quake(TM) fast inverse square root algorithm */
Expand All @@ -816,6 +651,7 @@ invSqrt(float x){
return x;
}

/*....................................................................*/
void checkGridDensities(configInfo *par, struct grid *gp){
/* This checks that none of the density samples is too small. */
int i;
Expand All @@ -836,6 +672,7 @@ void checkGridDensities(configInfo *par, struct grid *gp){
}
}

/*....................................................................*/
void lineBlend(molData *m, configInfo *par, struct blendInfo *blends){
/*
This obtains information on all the lines of all the radiating species which have other lines within some cutoff velocity separation.
Expand Down
Loading

0 comments on commit 7654865

Please sign in to comment.