From 843a979f0868c751ff9d1f04c3cdc6524e369ea5 Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Tue, 2 Apr 2019 13:59:56 +0100 Subject: [PATCH 1/4] changes to set a max ion stage if it is not Z+1 --- source/atomic.h | 10 ++++++---- source/get_atomicdata.c | 17 +++++++++++++++-- source/matrix_ion.c | 20 ++++++++++---------- 3 files changed, 31 insertions(+), 16 deletions(-) diff --git a/source/atomic.h b/source/atomic.h index fdf15115d..c9467a923 100755 --- a/source/atomic.h +++ b/source/atomic.h @@ -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; @@ -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 */ diff --git a/source/get_atomicdata.c b/source/get_atomicdata.c index 744cbbc18..436eb6240 100644 --- a/source/get_atomicdata.c +++ b/source/get_atomicdata.c @@ -250,6 +250,7 @@ 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++) @@ -257,6 +258,7 @@ get_atomic_data (masterfile) 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); @@ -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) { @@ -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 @@ -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) { diff --git a/source/matrix_ion.c b/source/matrix_ion.c index 1fe2fadd2..5841f170b 100644 --- a/source/matrix_ion.c +++ b/source/matrix_ion.c @@ -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) { @@ -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 (mm != ele[ion[mm].nelem].firstion) // we have electrons { rate_matrix[mm][mm] -= pi_rates[mm]; } @@ -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]; } @@ -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]); } @@ -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]); } @@ -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]); } @@ -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]); } @@ -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]); } @@ -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. From 5bd2036c67e5be8c2c320c3e0148764fcc7dddae Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Wed, 3 Apr 2019 09:59:59 +0100 Subject: [PATCH 2/4] fixed if statement in matrix_ion --- source/matrix_ion.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/matrix_ion.c b/source/matrix_ion.c index 5841f170b..40c7c3bb6 100644 --- a/source/matrix_ion.c +++ b/source/matrix_ion.c @@ -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[nn].nelem].firstion) // we have space for electrons { rate_matrix[mm][mm] -= xne * (rr_rates[mm] + xne * qrecomb_coeffs[mm]); } From dfa54d55430a0f5102bc0c179b4d798f93be9cc8 Mon Sep 17 00:00:00 2001 From: jhmatthews Date: Wed, 3 Apr 2019 13:24:02 +0100 Subject: [PATCH 3/4] fixed schoolboy error in ionization state calc --- source/get_atomicdata.c | 26 +++++++++++++------------- source/matrix_ion.c | 8 ++++---- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/source/get_atomicdata.c b/source/get_atomicdata.c index 436eb6240..a1a4f73ec 100644 --- a/source/get_atomicdata.c +++ b/source/get_atomicdata.c @@ -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 @@ -288,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 @@ -2438,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); } } } @@ -2624,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++; @@ -2667,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*/ @@ -2698,12 +2698,12 @@ 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; + ele[nelem].istate_max = ion[n].istate; while (ion[n].z == ele[nelem].z && n < nions) { diff --git a/source/matrix_ion.c b/source/matrix_ion.c index 40c7c3bb6..db8646dcf 100644 --- a/source/matrix_ion.c +++ b/source/matrix_ion.c @@ -121,7 +121,7 @@ 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 (mm != ele[ion[mm].nelem].firstion) // 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 } @@ -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 (mm != ele[ion[mm].nelem].firstion) // we have electrons + if (ion[mm].istate != ele[ion[mm].nelem].istate_max) // we have electrons { rate_matrix[mm][mm] -= pi_rates[mm]; } @@ -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 (mm != ele[ion[nn].nelem].firstion) // 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]); } @@ -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 (mm != ele[ion[mm].nelem].firstion && 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]); } From 554aefecf51e8aac42bc9473624c0e01089f881b Mon Sep 17 00:00:00 2001 From: Knox Long Date: Sat, 6 Apr 2019 09:39:04 -0500 Subject: [PATCH 4/4] Changed version to 83a --- source/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/Makefile b/source/Makefile index a0afc89d8..3b841452f 100755 --- a/source/Makefile +++ b/source/Makefile @@ -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