Skip to content

Fixes #104 #105

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
smaret merged 4 commits into from
Jul 12, 2016
Merged

Fixes #104 #105

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
c5a353b
Fixes #104
imckstewart Apr 13, 2016
d6e0477
Merge branch 'master' of https://github.com/lime-rt/lime into fix_104
imckstewart May 3, 2016
1dbb1c2
Merge branch 'master' of https://github.com/lime-rt/lime into fix_104
imckstewart Jun 24, 2016
cea608b
Merge branch 'master' of https://github.com/lime-rt/lime into fix_104
imckstewart Jul 11, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions src/lime.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,15 @@
#define N_SMOOTH_ITERS 20
#define TYPICAL_ISM_DENS 1000.0

/* Collision partner ID numbers from LAMDA */
#define CP_H2 1
#define CP_p_H2 2
#define CP_o_H2 3
#define CP_e 4
#define CP_H 5
#define CP_He 6
#define CP_Hplus 7


/* input parameters */
typedef struct {
Expand Down
83 changes: 50 additions & 33 deletions src/molinit.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,12 +103,13 @@ planckfunc(int iline, double temp, molData *m,int s){

void
molinit(molData *m, inputPars *par, struct grid *g,int i){
int id, ilev, iline, itrans, ispec, itemp, *ntemp, tnint=-1, idummy, ipart, *count,flag=0;
char *collpartname[] = {"H2","p-H2","o-H2","electrons","H","He","H+"}; /* definition from LAMDA */
int id, ilev, iline, itrans, ispec, itemp, *ntemp, tnint=-1, idummy, ipart, *collPartIDs,flag=0;
char *collpartnames[] = {"H2","p-H2","o-H2","electrons","H","He","H+"}; /* definition from LAMDA */
double fac, uprate, downrate=0, dummy, amass;
struct data { double *colld, *temp; } *part;
const int sizeI=200;

char string[200], specref[90], partstr[90];
char string[sizeI], specref[90], partstr[90];
FILE *fp;

if((fp=fopen(par->moldatfile[i], "r"))==NULL) {
Expand All @@ -117,27 +118,27 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){
}

/* Read the header of the data file */
fgets(string, 80, fp);
fgets(string, sizeI, fp);
fgets(specref, 90, fp);
fgets(string, 80, fp);
fgets(string, sizeI, fp);
fscanf(fp, "%lf\n", &amass);
fgets(string, 80, fp);
fgets(string, sizeI, fp);
fscanf(fp, "%d\n", &m[i].nlev);
fgets(string, 80, fp);
fgets(string, sizeI, fp);

m[i].eterm=malloc(sizeof(double)*m[i].nlev);
m[i].gstat=malloc(sizeof(double)*m[i].nlev);

/* Read the level energies and statistical weights */
for(ilev=0;ilev<m[i].nlev;ilev++){
fscanf(fp, "%d %lf %lf", &idummy, &m[i].eterm[ilev], &m[i].gstat[ilev]);
fgets(string, 80, fp);
fgets(string, sizeI, fp);
}

/* Read the number of transitions and allocate array space */
fgets(string, 80, fp);
fgets(string, sizeI, fp);
fscanf(fp, "%d\n", &m[i].nline);
fgets(string, 80, fp);
fgets(string, sizeI, fp);

m[i].lal = malloc(sizeof(int)*m[i].nline);
m[i].lau = malloc(sizeof(int)*m[i].nline);
Expand Down Expand Up @@ -168,18 +169,17 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){
g[id].mol[i].binv=1./g[id].mol[i].dopb;
}


/* Collision rates below here */
if(par->lte_only==0){
fgets(string, 80, fp);
fgets(string, sizeI, fp);
fscanf(fp,"%d\n", &m[i].npart);
count=malloc(sizeof(*count)*m[i].npart);
collPartIDs=malloc(sizeof(*collPartIDs)*m[i].npart);
m[i].down=malloc(sizeof(double*)*m[i].npart);
/* collision partner sanity check */

if(m[i].npart > par->collPart) flag=1;
if(m[i].npart < par->collPart){
if(!silent) bail_out("Error: Too many density profiles defined");
if(!silent) bail_out("Too many density profiles defined");
exit(1);
}

Expand All @@ -189,14 +189,28 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){
part = malloc(sizeof(struct data) * m[i].npart);

for(ipart=0;ipart<m[i].npart;ipart++){
fgets(string, 80, fp);
fscanf(fp,"%d\n", &count[ipart]);
fgets(string, 80, fp);
fgets(string, 80, fp);
fgets(string, sizeI, fp);
fscanf(fp,"%d\n", &collPartIDs[ipart]);

/* We want to test if the comment after the coll partner ID number is longer than the buffer size. To do this, we write a character - any character, as long as it is not \0 - to the last element of the buffer:
Copy link
Copy Markdown
Contributor

@smaret smaret Jul 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the comment after the collisional partner ID number used elsewhere in the code?

Copy link
Copy Markdown
Contributor Author

@allegroLeiden allegroLeiden Jul 7, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No. The intent of this elaborate check is purely to avoid a cock-up if the comment line exceeds the buffer length. It is essentially an error trap.

*/
string[sizeof(string)-1] = 'x';
if(fgets(string, sizeI, fp)==NULL){
if(!silent) bail_out("Read of collision-partner comment line failed.");
exit(1);
} else{
if(string[sizeof(string)-1]=='\0' && string[sizeof(string)-2]!='\n'){
/* The presence now of a final \0 means the comment string was either just long enough for the buffer, or too long; the absence of \n in the 2nd-last place means it was too long.
*/
if(!silent) bail_out("Collision-partner comment line is too long.");
exit(1);
}
}
fgets(string, sizeI, fp);
fscanf(fp,"%d\n", &m[i].ntrans[ipart]);
fgets(string, 80, fp);
fgets(string, sizeI, fp);
fscanf(fp,"%d\n", &ntemp[ipart]);
fgets(string, 80, fp);
fgets(string, sizeI, fp);

part[ipart].temp=malloc(sizeof(double)*ntemp[ipart]);

Expand All @@ -210,7 +224,7 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){
}

fscanf(fp,"\n");
fgets(string, 80, fp);
fgets(string, sizeI, fp);

part[ipart].colld=malloc(sizeof(double)*m[i].ntrans[ipart]*ntemp[ipart]);

Expand All @@ -236,10 +250,10 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){
fclose(fp);

/* Print out collision partner information */
strcpy(partstr, collpartname[count[0]-1]);
strcpy(partstr, collpartnames[collPartIDs[0]-1]);
for(ipart=1;ipart<m[i].npart;ipart++){
strcat( partstr, ", ");
strcat( partstr, collpartname[count[ipart]-1]);
strcat( partstr, collpartnames[collPartIDs[ipart]-1]);
}
if(!silent) {
collpartmesg(specref, m[i].npart);
Expand All @@ -250,9 +264,12 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){
/* Calculate molecular density */
for(id=0;id<par->ncell; id++){
for(ispec=0;ispec<par->nSpecies;ispec++){
if(m[i].npart == 1 && (count[0] == 1 || count[0] == 2 || count[0] == 3)){
if(m[i].npart == 1\
&& (collPartIDs[0]==CP_H2 || collPartIDs[0]==CP_p_H2 || collPartIDs[0]==CP_o_H2)){
g[id].nmol[ispec]=g[id].abun[ispec]*g[id].dens[0];
} else if(m[i].npart == 2 && (count[0] == 2 || count[0] == 3) && (count[1] == 2 || count[1] == 3)){
} else if(m[i].npart == 2\
&& (collPartIDs[0]==CP_p_H2 || collPartIDs[0]==CP_o_H2)\
&& (collPartIDs[1]==CP_p_H2 || collPartIDs[1]==CP_o_H2)){
if(!flag){
g[id].nmol[ispec]=g[id].abun[ispec]*(g[id].dens[0]+g[id].dens[1]);
} else {
Expand Down Expand Up @@ -283,13 +300,13 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){
g[id].mol[i].partner[ipart].t_binlow = tnint;
g[id].mol[i].partner[ipart].interp_coeff = fac;

} else if(g[id].t[0]<=part[ipart].temp[0]) {
g[id].mol[i].partner[ipart].t_binlow=0;
g[id].mol[i].partner[ipart].interp_coeff=0.0;
} else {
g[id].mol[i].partner[ipart].t_binlow=ntemp[ipart]-2;
g[id].mol[i].partner[ipart].interp_coeff=1.0;
}
} else if(g[id].t[0]<=part[ipart].temp[0]) {
g[id].mol[i].partner[ipart].t_binlow=0;
g[id].mol[i].partner[ipart].interp_coeff=0.0;
} else {
g[id].mol[i].partner[ipart].t_binlow=ntemp[ipart]-2;
g[id].mol[i].partner[ipart].interp_coeff=1.0;
}
}
}
}
Expand All @@ -299,7 +316,7 @@ molinit(molData *m, inputPars *par, struct grid *g,int i){
}
free(ntemp);
free(part);
free(count);
free(collPartIDs);
}
/* End of collision rates */

Expand Down