Skip to content
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

New scheme for two level atom #88

Merged
merged 13 commits into from
Jul 22, 2014
Merged

New scheme for two level atom #88

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
802f0b0
new zero emit call
Higginbottom Feb 13, 2014
ff3abc6
Merge branch 'dev' of https://github.com/agnwinds/python into dev
Higginbottom Apr 24, 2014
5f035a6
Modelled J_nu in lines
Higginbottom Apr 24, 2014
54e6953
test
Higginbottom May 6, 2014
2545d07
test
Higginbottom May 6, 2014
1c73c03
put mode in
Higginbottom May 9, 2014
42c297a
removed prints
Higginbottom May 9, 2014
62ed85c
Lines
Higginbottom May 9, 2014
d953780
quick test for w
Higginbottom Jun 13, 2014
ff7be1e
backstep
Higginbottom Jun 17, 2014
4d4f120
Final TLA version
Higginbottom Jul 21, 2014
667bac5
reset
Higginbottom Jul 21, 2014
3ff17ab
New two level atom scheme
Higginbottom Jul 22, 2014
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
4 changes: 3 additions & 1 deletion Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,9 @@ endif
LDFLAGS= -L$(LIB) -L$(LIB2) -lm -lcfitsio -lgsl -lgslcblas

#Note that version should be a single string without spaces.
VERSION = 78

VERSION = 77a_dev_nsh

CHOICE=1 // Compress plasma as much as possible
# CHOICE=0 // Keep relation between plasma and wind identical

Expand Down
125 changes: 60 additions & 65 deletions compton.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,79 +104,74 @@ kappa_ind_comp (xplasma, freq)
{
double x; // The opacity of the cell by the time we return it.
double sigma; /*The cross section, thompson, or KN if hnu/mec2 > 0.01 */
double J, expo; //The estimated intensity in the cell
int i;
double J; //The estimated intensity in the cell
// int i;

/* Previously, NSH had used the following formula whixch required ds and w to work
J=(4*PI*w*ds)/(C*xplasma->vol); //Calcuate the intensity NSH This works for a thin shell... Why? Dont know.
*/

J = 0.0; /* NSH 130605 to remove o3 compile error */

if (geo.ioniz_mode == 5 || geo.ioniz_mode == 7) /*If we are using power law ionization, use PL estimators */
{
if (geo.wcycle > 0) /* there is only any point in worrying if we have had at least one cycle otherwise there is no model */
{
for (i = 0; i < geo.nxfreq; i++)
{
if (geo.xfreq[i] < freq && freq <= geo.xfreq[i + 1]) //We have found the correct model band
{
if (xplasma->spec_mod_type[i] > 0) //Only bother if we have a model in this band
{
if (freq > xplasma->fmin_mod[i] && freq < xplasma->fmax_mod[i]) //The spectral model is defined for the frequency in question
{
if (xplasma->spec_mod_type[i] == SPEC_MOD_PL) //Power law model
{
J = pow(10,(xplasma->pl_log_w[i]+log10(freq)*xplasma->pl_alpha[i]));
}
else if (xplasma->spec_mod_type[i] == SPEC_MOD_EXP) //Exponential model
{
J = xplasma->exp_w[i] * exp ((-1 * H * freq) / (BOLTZMANN * xplasma->exp_temp[i]));
}
else
{
Error("kappa_ind_comp - unknown spectral model (%i) in band %i\n",xplasma->spec_mod_type[i], i);
J = 0.0; //Something has gone wrong
}
}
else
{
/* We have a spectral model, but it doesn't apply to the frequency in question.
clearly this is a slightly odd situation, where last time we didnt get a photon of this frequency,
but this time we did. Still this should only happen in very sparse cells, so induced compton is
unlikely to be important in such cells. We generate a warning, just so we can see if this is happening a lot */

Error ("kappa_ind_comp: frequency of photon (%e) is outside frequency range (%e - %e) of spectral model in cell %i band %i\n",
freq,xplasma->fmin_mod[i],xplasma->fmax_mod[i],xplasma->nplasma,i);

J = 0.0; //We
}
}
else /* There is no model in this band - this should not happen very often */
{
J = 0.0; //There is no modelin this band, so the best we can do is assume zero J
Error ("kappa_ind_comp: no model exists in cell %i band %i\n",xplasma->nplasma,i);
}



}
}
}
else //We have not completed an ionization cycle, so no chance of a model
{
J=0.0;
}
}

else /*Else, use BB estimator of J */
{
expo = (H * freq) / (BOLTZMANN * xplasma->t_r);
J = (2 * H * freq * freq * freq) / (C * C);
J *= 1 / (exp (expo) - 1);
J *= xplasma->w;
}

J = mean_intensity (xplasma, freq, 2); //Obtain a model for the mean intensity - we call this with mode=2, which means that if we have not yet completed a cycle, dont return a dilute blackbody estimate if we are in PL mode.

//if (geo.ioniz_mode == 5 || geo.ioniz_mode == 7) /*If we are using power law ionization, use PL estimators */
// {
// if (geo.wcycle > 0) /* there is only any point in worrying if we have had at least one cycle otherwise there is no model */
// {
// for (i = 0; i < geo.nxfreq; i++)
// {
// if (geo.xfreq[i] < freq && freq <= geo.xfreq[i + 1]) //We have found the correct model band
// {
// if (xplasma->spec_mod_type[i] > 0) //Only bother if we have a model in this band
// {
// if (freq > xplasma->fmin_mod[i] && freq < xplasma->fmax_mod[i]) //The spectral model is defined for the frequency in question
// {
// if (xplasma->spec_mod_type[i] == SPEC_MOD_PL) //Power law model
// {
// J = pow(10,(xplasma->pl_log_w[i]+log10(freq)*xplasma->pl_alpha[i]));
// }
// else if (xplasma->spec_mod_type[i] == SPEC_MOD_EXP) //Exponential model
// {
// J = xplasma->exp_w[i] * exp ((-1 * H * freq) / (BOLTZMANN * xplasma->exp_temp[i]));
// }
// else
// {
// Error("kappa_ind_comp - unknown spectral model (%i) in band %i\n",xplasma->spec_mod_type[i], i);
// J = 0.0; //Something has gone wrong
// }
// }
// else /*We have a spectral model, but it doesnt apply to the frequency in question. clearly this is a slightly odd situation, where last time we didnt get a photon of this frequency, but this time we did. Still this should only happen in very sparse cells, so induced compton is unlikely to be important in such cells. We generate a warning, just so we can see if this is happening a lot*/
// {
// Warning ("kappa_ind_comp: frequency of photon (%e) is outside frequency range (%e - %e) of spectral model in cell %i band %i\n",freq,xplasma->fmin_mod[i],xplasma->fmax_mod[i],xplasma->nplasma,i); //This is unlikely to happen very often, but if it does, we should probably know about it
// J = 0.0; //We
// }
// }
// else /* There is no model in this band - this should not happen very often */
// {
// J = 0.0; //THere is no modelin this band, so the best we can do is assume zero J
// Warning ("kappa_ind_comp: no model exists in cell %i band %i\n",xplasma->nplasma,i); //This is unlikely to happen very often, but if it does, we should probably know about it
// }
//
//
//
// }
// }
// }
// else //We have not completed an ionization cycle, so no chance of a model
// {
// J=0.0;
// }
//}
//
//else /*Else, use BB estimator of J */
// {
// / / expo = (H * freq) / (BOLTZMANN * xplasma->t_r);
// J = (2 * H * freq * freq * freq) / (C * C);
// J *= 1 / (exp (expo) - 1);
// J *= xplasma->w;
// }

sigma = klein_nishina (freq); //NSH 130214 - full KN formula

Expand Down
1 change: 1 addition & 0 deletions foo.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,7 @@ double den_config(PlasmaPtr xplasma, int nconf);
double pop_kappa_ff_array(void);
int update_banded_estimators(PlasmaPtr xplasma, PhotPtr p, double ds, double w_ave);
int save_photon_stats(WindPtr one, PhotPtr p, double ds);
double mean_intensity(PlasmaPtr xplasma, double freq, int mode);
/* wind_updates2d.c */
int wind_update(WindPtr (w));
int wind_rad_init(void);
Expand Down
31 changes: 20 additions & 11 deletions lines.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,6 @@ lum_lines (one, nmin, nmax)
if (dd > LDEN_MIN)
{ /* potentially dangerous step to avoid lines with no power */
two_level_atom (lin_ptr[n], xplasma, &d1, &d2);

x = foo1 = lin_ptr[n]->gu / lin_ptr[n]->gl * d1 - d2;

z = exp (-H_OVER_K * lin_ptr[n]->freq / t_e);
Expand Down Expand Up @@ -341,6 +340,9 @@ a21 (line_ptr)
07mar ksl 58c -- Tried to address problems associated with our
inconsistent treatment of level populations for two
level atoms, This is at best a bandaide.
14jul nsh 78 -- changed to allow the use of a computed model for
the mean intensity in a cell to calualate influence of radiation
on the upper state population of a two level atom.
*/

struct lines *old_line_ptr;
Expand All @@ -366,7 +368,7 @@ two_level_atom (line_ptr, xplasma, d1, d2)
double xw;
double ne, te, w, tr, dd;
int nion;

double J; //Model of the specific intensity


//Check and exit if this routine is called for a macro atom, since this should never happen
Expand Down Expand Up @@ -421,15 +423,22 @@ in the configuration structure. 01dec ksl */
c21 = ne * q;
c12 = c21 * g2_over_g1 * exp (-H_OVER_K * freq / te);

if (w < 1.e-6)
{ // Radiation is unimportant
n2_over_n1 = c12 / (c21 + a);
}
else
{ //Include effects of stimulated emission
z = w / (exp (H_OVER_K * freq / tr) - 1.);
n2_over_n1 = (c12 + g2_over_g1 * a * z) / (c21 + a * (1. + z));
}
// if (w < 1.e-6) //NSH this if block removed to simplify code - if can be reinstated if runtimes are a problem, but it will need an additional if statement to avoid missing it out if we are using a modelled specific intensity.
// { // Radiation is unimportant
// n2_over_n1 = c12 / (c21 + a);
// }
// else
// { //Include effects of stimulated emission
// z = w / (exp (H_OVER_K * freq / tr) - 1.); //original
z=(C*C)/(2.*H*freq*freq*freq); //This is the factor which relates the A coefficient to the b coefficient
// n2_over_n1 = (c12 + g2_over_g1 * a * z) / (c21 + a * (1. + z)); //original

J = mean_intensity (xplasma, freq, 1);/* we call mean intensity with mode 1 - this means we are happy to use the dilute blackbody approximation even if we havent run enough spectral cycles to have a model for J*/

n2_over_n1 = (c12 + g2_over_g1 * a * z * J) / (c21 + a*(1. + (J * z))); //this equation is equivalent to equation 4.29 in NSH's thesis with the einstein b coefficients replaced by a multiplied by suitable conversion factors from the einstein relations.


// }


*d1 = dd;
Expand Down
Loading