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

Max ion stage #502

Merged
merged 4 commits into from
Apr 6, 2019
Merged

Max ion stage #502

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
843a979
changes to set a max ion stage if it is not Z+1
Apr 2, 2019
5bd2036
fixed if statement in matrix_ion
Apr 3, 2019
dfa54d5
fixed schoolboy error in ionization state calc
Apr 3, 2019
554aefe
Changed version to 83a
Apr 6, 2019
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
2 changes: 1 addition & 1 deletion source/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ LDFLAGS+= -L$(LIB) -lm -lgsl -lgslcblas

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

VERSION = 83
VERSION = 83a

startup:
@echo $(COMPILER_PRINT_STRING) # prints out compiler information
Expand Down
10 changes: 6 additions & 4 deletions source/atomic.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,9 +85,10 @@ typedef struct elements
provides and index to the atomic data */
char name[20]; /* Element name */
int z; /* Atomic number */
int firstion; /*Index into struct ions ion[firstion] is the lowest ionization state of this ion */
int nions; /*The number of ions actually read in from the data file for this element */
double abun; /*Abundance */
int firstion; /* Index into struct ions ion[firstion] is the lowest ionization state of this ion */
int nions; /* The number of ions actually read in from the data file for this element */
double abun; /* Abundance */
int istate_max; /* highest ionization stage of element */
}
ele_dummy, *ElemPtr;

Expand All @@ -103,7 +104,8 @@ double rho2nh; /* The conversion constant from rho to nh the to
typedef struct ions
{
int z; /* defines the element for this ion */
int istate; /*1=neutral, 2 = once ionized, etc. */
int istate; /* 1=neutral, 2 = once ionized, etc. */
int nelem; /* index to elements structure */
double ip; /* ionization potential of this ion (converted from eV to ergs by get_atomic) */
double g; /* multiplicity of ground state, note that this is not totally consistent
with energy levels and this needs to be reconciled */
Expand Down
39 changes: 26 additions & 13 deletions source/get_atomicdata.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,7 @@ get_atomic_data (masterfile)
double gstmin, gstmax; //The range of temperatures for which all ions have GS RR rates
double gsqrdtemp, gfftemp, s1temp, s2temp, s3temp; //Temporary storage for gaunt factors
int n_elec_yield_tot; //The number of inner shell cross sections with matching electron yield arrays
int inner_no_e_yield; //The number of inner shell cross sections with no yields
int inner_no_e_yield; //The number of inner shell cross sections with no yields
// int n_fluor_yield_tot; //The number of inner shell cross sections with matching fluorescent photon yield arrays
double I, Ea; //The ionization energy and mean electron energy for electron yields
// double energy; //The energy of inner shell fluorescent photons
Expand Down Expand Up @@ -250,13 +250,15 @@ get_atomic_data (masterfile)
ele[n].abun = (-1);
ele[n].firstion = (-1);
ele[n].nions = (-1);
ele[n].istate_max = (-1);
}

for (n = 0; n < NIONS; n++)
{
simple_line_ignore[n] = 0; // diagnostic counter for how many lines ignored
ion[n].z = (-1);
ion[n].istate = (-1);
ion[n].nelem = (-1);
ion[n].ip = (-1);
ion[n].g = (-1);
ion[n].nmax = (-1);
Expand Down Expand Up @@ -286,8 +288,8 @@ get_atomic_data (masterfile)
}

nlevels = nxphot = nphot_total = ntop_phot = nauger = ndrecomb = n_inner_tot = 0; //Added counter for DR//
n_elec_yield_tot = 0; //Counter for electron yield
// n_fluor_yield_tot = 0; and fluorescent photon yields
n_elec_yield_tot = 0; //Counter for electron yield
// n_fluor_yield_tot = 0; and fluorescent photon yields

/*This initializes the top_phot array - it is used for all ionization processes so some elements
are only used in some circumstances
Expand Down Expand Up @@ -1427,7 +1429,7 @@ described as macro-levels. */

for (nion = 0; nion < nions; nion++)
{
if (ion[nion].z == z && ion[nion].istate == istate)
if (ion[nion].z == z && ion[nion].istate == istate && ion[nion].macro_info != 1)
{
if (ion[nion].phot_info == -1)
{
Expand Down Expand Up @@ -1527,7 +1529,7 @@ described as macro-levels. */
}
for (nion = 0; nion < nions; nion++)
{
if (ion[nion].z == z && ion[nion].istate == istate)
if (ion[nion].z == z && ion[nion].istate == istate && ion[nion].macro_info != 1)
{
/* Then there is a match */
inner_cross[n_inner_tot].nlev = ion[nion].firstlevel; //All these are for the ground state
Expand Down Expand Up @@ -2436,7 +2438,7 @@ would like to have simple lines for macro-ions */
}
else
{
Error ("Get_atomic_data: more than one electron yield record for inner_cross %i z=%i istate=%i\n", n,z,istate);
Error ("Get_atomic_data: more than one electron yield record for inner_cross %i z=%i istate=%i\n", n, z, istate);
}
}
}
Expand Down Expand Up @@ -2622,14 +2624,14 @@ SCUPS 1.132e-01 2.708e-01 5.017e-01 8.519e-01 1.478e+00

n_elec_yield_tot = 0; //Reset this numnber, we are now going to use it to check we have yields for all inner shells
// n_fluor_yield_tot = 0; //Reset this numnber, we are now going to use it to check we have yields for all inner shells
inner_no_e_yield =0;
inner_no_e_yield = 0;

for (n = 0; n < n_inner_tot; n++)
{
if (inner_cross[n].n_elec_yield != -1)
n_elec_yield_tot++;
else
inner_no_e_yield++;
inner_no_e_yield++;
// Error_silent ("get_atomicdata: No inner electron yield data for inner cross section %i\n", n);
// if (inner_cross[n].n_fluor_yield != -1)
// n_fluor_yield_tot++;
Expand Down Expand Up @@ -2665,12 +2667,12 @@ SCUPS 1.132e-01 2.708e-01 5.017e-01 8.519e-01 1.478e+00
Error ("Ignored %d simple lines for macro-ion %d\n", simple_line_ignore[n], n);
}
/* report ignored collision strengths */
if (cstren_no_line > 0)
if (cstren_no_line > 0)
Error ("Ignored %d collision strengths with no matching line transition\n", cstren_no_line);
if (inner_no_e_yield >0)
Error ("Ignoring %d inner shell cross sections because no matching yields\n", inner_no_e_yield);
if (inner_no_e_yield > 0)
Error ("Ignoring %d inner shell cross sections because no matching yields\n", inner_no_e_yield);



/* Now begin a series of calculations with the data that has been read in in order
to prepare it for use by other programs*/
Expand All @@ -2696,9 +2698,20 @@ exit if there is an element with no ions */
n = 0;
while (ion[n].z != ele[nelem].z && n < nions)
n++; /* Find the first ion of that element for which there is data */

ele[nelem].firstion = n;
ion[n].nelem = nelem;

/* find the highest ion stage and the number of ions */
ele[nelem].istate_max = ion[n].istate;

while (ion[n].z == ele[nelem].z && n < nions)
{
if (ele[nelem].istate_max < ion[n].istate)
ele[nelem].istate_max = ion[n].istate;
ion[n].nelem = nelem;
n++;
}
ele[nelem].nions = n - ele[nelem].firstion;
if (ele[nelem].firstion == nions)
{
Expand Down
22 changes: 11 additions & 11 deletions source/matrix_ion.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,11 +121,11 @@ matrix_ion_populations (xplasma, mode)
{
newden[mm] = xplasma->density[mm] / elem_dens[ion[mm].z]; // newden is our local fractional density array
xion[mm] = mm; // xion is an array we use to track which ion is in which row of the matrix
if (ion[mm].istate != 1) // We can recombine since we are not in the first ionization stage
if (mm != ele[ion[mm].nelem].firstion) // We can recombine since we are not in the first ionization stage
{
rr_rates[mm] = total_rrate (mm, xplasma->t_e); // radiative recombination rates
}
if (ion[mm].istate != ion[mm].z + 1) // we can photoionize, since we are not in the z+1th ionization state (bare)
if (ion[mm].istate != ele[ion[mm].nelem].istate_max) // we can photoionize, since we are not in the highest ionization state
{
if (mode == NEBULARMODE_MATRIX_BB)
{
Expand Down Expand Up @@ -455,7 +455,7 @@ populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp,

for (mm = 0; mm < nions; mm++)
{
if (ion[mm].istate != ion[mm].z + 1) // we have electrons
if (ion[mm].istate != ele[ion[mm].nelem].istate_max) // we have electrons
{
rate_matrix[mm][mm] -= pi_rates[mm];
}
Expand All @@ -469,7 +469,7 @@ populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp,
for (nn = 0; nn < nions; nn++)
{

if (mm == nn + 1 && ion[nn].istate != ion[nn].z + 1 && ion[mm].z == ion[nn].z)
if (mm == nn + 1 && ion[nn].istate != ele[ion[nn].nelem].istate_max && ion[mm].z == ion[nn].z)
{
rate_matrix[mm][nn] += pi_rates[nn];
}
Expand All @@ -480,7 +480,7 @@ populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp,

for (mm = 0; mm < nions; mm++)
{
if (ion[mm].istate != ion[mm].z + 1 && ion[mm].dere_di_flag > 0) // we have electrons and a DI rate
if (ion[mm].istate != ele[ion[mm].nelem].istate_max && ion[mm].dere_di_flag > 0) // we have electrons and a DI rate
{
rate_matrix[mm][mm] -= (xne * di_coeffs[mm]);
}
Expand All @@ -492,7 +492,7 @@ populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp,
{
for (nn = 0; nn < nions; nn++)
{
if (mm == nn + 1 && ion[nn].istate != ion[nn].z + 1 && ion[mm].z == ion[nn].z && ion[nn].dere_di_flag > 0)
if (mm == nn + 1 && ion[nn].istate != ele[ion[nn].nelem].istate_max && ion[mm].z == ion[nn].z && ion[nn].dere_di_flag > 0)
{
rate_matrix[mm][nn] += (xne * di_coeffs[nn]);
}
Expand All @@ -504,7 +504,7 @@ populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp,

for (mm = 0; mm < nions; mm++)
{
if (ion[mm].istate != 1) // we have space for electrons
if (mm != ele[ion[mm].nelem].firstion) // we have space for electrons
{
rate_matrix[mm][mm] -= xne * (rr_rates[mm] + xne * qrecomb_coeffs[mm]);
}
Expand All @@ -517,7 +517,7 @@ populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp,
{
for (nn = 0; nn < nions; nn++)
{
if (mm == nn - 1 && ion[nn].istate != 1 && ion[mm].z == ion[nn].z)
if (mm == nn - 1 && nn != ele[ion[nn].nelem].firstion && ion[mm].z == ion[nn].z)
{
rate_matrix[mm][nn] += xne * (rr_rates[nn] + xne * qrecomb_coeffs[nn]);
}
Expand All @@ -528,7 +528,7 @@ populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp,

for (mm = 0; mm < nions; mm++)
{
if (ion[mm].istate != 1 && ion[mm].drflag > 0) // we have space for electrons
if (mm != ele[ion[mm].nelem].firstion && ion[mm].drflag > 0) // we have space for electrons
{
rate_matrix[mm][mm] -= (xne * dr_coeffs[mm]);
}
Expand All @@ -543,7 +543,7 @@ populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp,
{
for (nn = 0; nn < nions; nn++)
{
if (mm == nn - 1 && ion[nn].istate != 1 && ion[mm].z == ion[nn].z && ion[mm].drflag > 0)
if (mm == nn - 1 && nn != ele[ion[nn].nelem].firstion && ion[mm].z == ion[nn].z && ion[mm].drflag > 0)
{
rate_matrix[mm][nn] += (xne * dr_coeffs[nn]);
}
Expand Down Expand Up @@ -588,7 +588,7 @@ populate_ion_rate_matrix (rate_matrix, pi_rates, inner_rates, rr_rates, b_temp,
zcount = 0;
for (nn = 0; nn < nions; nn++)
{
if (ion[nn].istate == 1)
if (nn == ele[ion[nn].nelem].firstion)
{
b_temp[nn] = 1.0; //In the relative abundance schene this equals one.

Expand Down