Skip to content

Commit

Permalink
Merge branch 'dev' of https://github.com/agnwinds/python into dev
Browse files Browse the repository at this point in the history
  • Loading branch information
@Higginbottom
Higginbottom committed Aug 14, 2015
2 parents c6e5f96 + 010228e commit e710fe8
Show file tree
Hide file tree
Showing 3 changed files with 95 additions and 85 deletions.
2 changes: 1 addition & 1 deletion source/atomic.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ int nions; /*The actual number of ions read from the datafile */
int nlevels; /*These are the actual number of levels which were read in */
#define NLTE_LEVELS 270 /* Maximum number of levels to treat explicitly */
int nlte_levels; /* Actual number of levels to treat explicityly */
#define NLEVELS_MACRO 150 /* Maximum number of macro atom levels. (SS, June 04) */
#define NLEVELS_MACRO 200 /* Maximum number of macro atom levels. (SS, June 04) */
int nlevels_macro; /* Actual number of macro atom levels. (SS, June 04) */
#define NLINES 200000 /* Maximum number of lines to be read */
int nlines; /* Actual number of lines that were read in */
Expand Down
167 changes: 85 additions & 82 deletions source/direct_ion.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#include "python.h"
#include "recipes.h"

#include <gsl/gsl_sf_expint.h> //We need this gsl library to evaluate the first exponential integral
#include <gsl/gsl_sf_expint.h> //We need this gsl library to evaluate the first exponential integral



Expand All @@ -17,126 +17,130 @@ int
compute_di_coeffs (T)
double T;
{
int i,n,nrec;
int i, n, nrec;
double rate, drdt, t, dt, scaled_t;
int imin,imax;
double rates[DERE_DI_PARAMS],temps[DERE_DI_PARAMS];
int imin, imax;
double rates[DERE_DI_PARAMS], temps[DERE_DI_PARAMS];
double exp_int;


for (n = 0; n < nions; n++)
{
if (ion[n].dere_di_flag == 0)
{
//printf ("NO COEFFS\n");
//printf ("NO COEFFS\n");
di_coeffs[n] = 0.0;
}

else
{
nrec=ion[n].nxderedi;
/* find the correct coefficient */

nrec = ion[n].nxderedi;
/* find the correct coefficient */


for (i = 0; i < dere_di_rate[nrec].nspline; i++)
{
rates[i] = dere_di_rate[nrec].rates[i];
temps[i] = dere_di_rate[nrec].temps[i];
}
for (i = 0; i < dere_di_rate[nrec].nspline; i++)
{
rates[i] = dere_di_rate[nrec].rates[i];
temps[i] = dere_di_rate[nrec].temps[i];
}


t=(BOLTZMANN*T)/(dere_di_rate[nrec].xi*EV2ERGS);
scaled_t=1.0-((log(2.0))/(log(2.0+t)));
t = (BOLTZMANN * T) / (dere_di_rate[nrec].xi * EV2ERGS);
scaled_t = 1.0 - ((log (2.0)) / (log (2.0 + t)));

if (scaled_t < temps[0]) //we are below the range of DI data data
{
Log_silent
("bad_gs_rr: Requested temp %e is below limit of data for ion %i(Tmin= %e)\n",
scaled_t, n, temps[0]);
// rate = rates[0];
imax=1;
imin=0;
}

else if (scaled_t >= temps[dere_di_rate[nrec].nspline - 1]) //we are above the range of GS data
{
Log_silent
("bad_gs_rr: Requested temp %e is above limit (%e) of data for ion %i\n",
scaled_t, temps[dere_di_rate[nrec].nspline - 1],n);
// rate = rates[BAD_GS_RR_PARAMS - 1];
imax=dere_di_rate[nrec].nspline - 1;
imin=dere_di_rate[nrec].nspline - 2;
//We will try to extrapolate.
if (scaled_t < temps[0]) //we are below the range of DI data data
{

Log_silent("compute_di_coeffs: Requested temp %e is below limit of data for ion %i(Tmin= %e)\n",
scaled_t, n, temps[0]);
//rate = rates[0];
imax = 1;
imin = 0;
}

else if (scaled_t >= temps[dere_di_rate[nrec].nspline - 1]) //we are above the range of GS data
{
Log_silent("compute_di_coeffs: Requested temp %e is above limit (%e) of data for ion %i\n",
scaled_t, temps[dere_di_rate[nrec].nspline - 1], n);

}
else //We must be within the range of tabulated data
{
for (i = 0; i < dere_di_rate[nrec].nspline - 1; i++)
{
if (temps[i] <= scaled_t && scaled_t < temps[i + 1]) //We have bracketed the correct temperature
//rate = rates[BAD_GS_RR_PARAMS - 1];
imax = dere_di_rate[nrec].nspline - 1;
imin = dere_di_rate[nrec].nspline - 2;
//We will try to extrapolate.
}

else //We must be within the range of tabulated data
{
imin = i;
imax = i + 1;
for (i = 0; i < dere_di_rate[nrec].nspline - 1; i++)
{
if (temps[i] <= scaled_t && scaled_t < temps[i + 1]) //We have bracketed the correct temperature
{
imin = i;
imax = i + 1;
}
}
/* NSH 140313 - changed the following lines to interpolate in log space */
}
}
/* NSH 140313 - changed the following lines to interpolate in log space */
}
//if (ion[n].z==6 && ion[n].istate==4){
//
//printf ("T=%e Interploating between %e(%e) and %e(%e)\n",scaled_t,rates[imin],temps[imin],rates[imax],temps[imax]);
//}
drdt = ((rates[imax]) - (rates[imin])) / ((temps[imax]) - (temps[imin]));
dt = ((scaled_t) - (temps[imin]));
rate = rates[imin] + drdt * dt;


//printf ("t=%e xi=%e scaled_rate=%e\n",t,dere_di_rate[nrec].xi,rate);

//if (ion[n].z==6 && ion[n].istate==4){
//
//printf ("T=%e Interploating between %e(%e) and %e(%e)\n",scaled_t,rates[imin],temps[imin],rates[imax],temps[imax]);
//}

di_coeffs[n]=pow(t,-0.5)*pow(dere_di_rate[nrec].xi,-1.5)*rate;
//printf ("Doing 1/t=%e\n",1.0/t);
drdt = ((rates[imax]) - (rates[imin])) / ((temps[imax]) - (temps[imin]));
dt = ((scaled_t) - (temps[imin]));
rate = rates[imin] + drdt * dt;

if (exp(-1.0/t) < (1.0/VERY_BIG))
{
exp_int=0.0;
}
else
{
exp_int=gsl_sf_expint_E1(1.0/t); //Evaluate the first exponential integral using gsl library.
}

di_coeffs[n]*=exp_int;
}
//printf ("t=%e xi=%e scaled_rate=%e\n",t,dere_di_rate[nrec].xi,rate);


di_coeffs[n] =
pow (t, -0.5) * pow (dere_di_rate[nrec].xi, -1.5) * rate;

//printf ("Doing 1/t=%e\n",1.0/t);

} //End of loop over ions
if (exp (-1.0 / t) < (1.0 / VERY_BIG))
{
exp_int = 0.0;
}
else
{
exp_int = gsl_sf_expint_E1 (1.0 / t); //Evaluate the first exponential integral using gsl library.
}

di_coeffs[n] *= exp_int;
}



} //End of loop over ions


return(0);
}

return (0);
}

double
total_di (one, t_e)
WindPtr one; // Pointer to the current wind cell - we need the cell volume, this is not in the plasma structure
double t_e; //Current electron temperature of the cell

{
double x; //The returned variable
int nplasma; //The cell number in the plasma array
PlasmaPtr xplasma; //pointer to the relevant cell in the plasma structure
int n; //loop pointers
double x; //The returned variable
int nplasma; //The cell number in the plasma array
PlasmaPtr xplasma; //pointer to the relevant cell in the plasma structure
int n; //loop pointers


nplasma = one->nplasma; //Get the correct plasma cell related to this wind cell
nplasma = one->nplasma; //Get the correct plasma cell related to this wind cell
xplasma = &plasmamain[nplasma]; //copy the plasma structure for that cell to local variable
x = 0; //zero the luminosity
x = 0; //zero the luminosity


compute_di_coeffs (t_e); //Calculate the DR coefficients for this cell
compute_di_coeffs (t_e); //Calculate the DR coefficients for this cell


for (n = 0; n < nions; n++)
Expand All @@ -148,14 +152,13 @@ total_di (one, t_e)
else
{

x +=
xplasma->vol * xplasma->ne * xplasma->density[n] * di_coeffs[n] *
dere_di_rate[ion[n].nxderedi].xi*EV2ERGS;
//printf ("n=%i V=%e ne=%e rho=%e coeff=%e xi=%e cooling=%e\n",n, V , xplasma->ne , xplasma->density[n] , di_coeffs[n] ,
// dere_di_rate[ion[n].nxderedi].xi*EV2ERGS,x);
x += xplasma->vol * xplasma->ne * xplasma->density[n] * di_coeffs[n] *
dere_di_rate[ion[n].nxderedi].xi * EV2ERGS;

//printf ("n=%i V=%e ne=%e rho=%e coeff=%e xi=%e cooling=%e\n",n, V ,
//xplasma->ne , xplasma->density[n] , di_coeffs[n] ,
//dere_di_rate[ion[n].nxderedi].xi*EV2ERGS,x);
}
}
return (x);
}


11 changes: 9 additions & 2 deletions source/python.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -910,9 +910,16 @@ typedef struct photon
PTYPE_BL=1,
PTYPE_DISK=2,
PTYPE_WIND=3,
PTYPE_AGN=4
PTYPE_AGN=4,
PTYPE_STAR_MATOM=10,
PTYPE_BL_MATOM=11,
PTYPE_DISK_MATOM=12,
PTYPE_WIND_MATOM=13,
PTYPE_AGN_MATOM=14
} origin; /* Where this photon originated. If the photon has
scattered it's "origin" may be changed to "wind".*/
scattered it's "origin" may be changed to "wind".*/
/* note that we add 10 to origin when processed by a macro-atom
which means we need these values in the enum list */
int np; /*NSH 13/4/11 - an internal pointer to the photon number so
so we can write out details of where the photon goes */
double path; /* SWM - Photon path length */
Expand Down

0 comments on commit e710fe8

Please sign in to comment.