Navigation

← Previous Change
Next Change →

iso_create.cpp

Timestamp:

05/10/08 09:03:02 (2 months ago)

Author:

peter

Message:

Merging all changes from mainline upto r2033, except r1902.

Files:

branches/c08_branch/source/iso_create.cpp (modified) (30 diffs)

Legend:

: Unmodified
: Added
: Removed
: Modified
: Copied
: Moved

branches/c08_branch/source/iso_create.cpp

r1916	r2034
9	9	#include "elementnames.h"
10	10	#include "helike.h"
11		~~#include "helike_cs.h"~~
12	11	#include "helike_einsta.h"
13	12	#include "hydro_bauman.h"
…	…
16	15	#include "iso.h"
17	16	#include "lines_service.h"
	17	#include "opacity.h"
18	18	#include "phycon.h"
19	19	#include "physconst.h"
20	20	#include "secondaries.h"
21	21	#include "taulines.h"
	22	#include "thirdparty.h"
22	23
23	24	/iso_zero zero data for hydrogen and helium /
…	…
34	35	/* calculate radiative lifetime of an individual iso state */
35	36	STATIC double iso_state_lifetime( long ipISO, long nelem, long n, long l );
36
37		~~/* calculate cascade probabilities, branching ratios, and associated errors. */~~
38		~~STATIC void iso_cascade(void);~~
39	37
40	38	STATIC void iso_satellite( void );
…	…
177	175	/***************************************************************/
178	176	/* NB - the above is all we need if we are compiling recombination tables. */
	177	iso_recomb_malloc();
179	178	iso_recomb_setup( ipH_LIKE );
180	179	iso_recomb_setup( ipHE_LIKE );
	180	iso_recomb_auxiliary_free();
181	181
182	182	/* set up helium collision tables */
…	…
304	304
305	305	/* following comes out very slightly off, correct here */
306		Transitions[ipH_LIKE][ipHELIUM][3][ipH2s].WLAng = 1640.f;
307		Transitions[ipH_LIKE][ipHELIUM][3][ipH2p].WLAng = 1640.f;
	306	Transitions[ipH_LIKE][ipHELIUM][ipH3s][ipH2s].WLAng = 1640.f;
	307	Transitions[ipH_LIKE][ipHELIUM][ipH3s][ipH2p].WLAng = 1640.f;
308	308	if( iso.n_HighestResolved_max[ipH_LIKE][ipHELIUM] >=3 )
309	309	{
310		Transitions[ipH_LIKE][ipHELIUM][4][ipH2s].WLAng = 1640.f;
311		Transitions[ipH_LIKE][ipHELIUM][4][ipH2p].WLAng = 1640.f;
312		Transitions[ipH_LIKE][ipHELIUM][5][ipH2s].WLAng = 1640.f;
313		Transitions[ipH_LIKE][ipHELIUM][5][ipH2p].WLAng = 1640.f;
	310	Transitions[ipH_LIKE][ipHELIUM][ipH3p][ipH2s].WLAng = 1640.f;
	311	Transitions[ipH_LIKE][ipHELIUM][ipH3p][ipH2p].WLAng = 1640.f;
	312	Transitions[ipH_LIKE][ipHELIUM][ipH3d][ipH2s].WLAng = 1640.f;
	313	Transitions[ipH_LIKE][ipHELIUM][ipH3d][ipH2p].WLAng = 1640.f;
314	314	}
315	315
…	…
330	330	{
331	331	StatesElem[ipISO][nelem][ipHi].lifetime = SMALLFLOAT;
332		for( ipLo=0; ipLo < ipHi; ipLo++ )
	332
	333	long ipLoStart = 0;
	334	if( opac.lgCaseB && L_(ipHi)==1 && (ipISO==ipH_LIKE \|\| S_(ipHi)==1) )
	335	ipLoStart = 1;
	336
	337	for( ipLo=ipLoStart; ipLo < ipHi; ipLo++ )
333	338	{
334	339	if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
…	…
391	396	iso_zero();
392	397
393		/* calculate cascade probabilities, branching ratios, and associated errors. */
394		iso_cascade();
	398	/* loop over iso sequences */
	399	for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
	400	{
	401	for( nelem = ipISO; nelem < LIMELM; nelem++ )
	402	{
	403	/* must always do helium even if turned off */
	404	if( nelem == ipISO \|\| dense.lgElmtOn[nelem] )
	405	{
	406	/* calculate cascade probabilities, branching ratios, and associated errors. */
	407	iso_cascade( ipISO, nelem);
	408	}
	409	}
	410	}
395	411
396	412	iso_satellite();
…	…
402	418	/***************************************/
403	419	/* fill in iso.strkar array */
404		for( ipISO=ipH_LIKE; ipISO < NISO; ++ipISO )
	420	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
405	421	{
406	422	for( nelem=ipISO; nelem < LIMELM; nelem++ )
…	…
480	496
481	497	/* main isoelectronic arrays, fill with sane values */
482		for( ipISO=ipH_LIKE; ipISO < NISO; ++ipISO )
	498	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
483	499	{
484	500	for( nelem=ipISO; nelem < LIMELM; nelem++ )
…	…
537	553	}
538	554
539
540		/Added:26th Oct 07/
541		/For Atoms & Molecules/
	555	/* database lines - DB lines */
542	556	for( i=0; i <linesAdded2; i++)
543	557	{
544		fixit(); /* why not TransitionZero here? */
545		EmLineZero(&atmolEmis[i]);
546		}
547
	558	TransitionZero( atmolEmis[i].tran );
	559	}
548	560
549	561	for( i=0; i < nCORotate; i++ )
…	…
577	589	iso.BranchRatio.reserve( NISO );
578	590	iso.CascadeProb.reserve( NISO );
	591	iso.CachedAs.reserve( NISO );
579	592
580	593	/* the hydrogen and helium like iso-sequences */
…	…
589	602	iso.BranchRatio.reserve( ipISO, LIMELM );
590	603	iso.CascadeProb.reserve( ipISO, LIMELM );
	604	iso.CachedAs.reserve( ipISO, LIMELM );
591	605
592	606	for( nelem=ipISO; nelem < LIMELM; ++nelem )
…	…
607	621	iso.DielecRecombVsTemp.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
608	622	iso.Boltzmann.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
	623
	624	iso.CachedAs.reserve( ipISO, nelem, MAX2(1, iso.nCollapsed_max[ipISO][nelem]) );
	625
	626	for( i = 0; i < iso.nCollapsed_max[ipISO][nelem]; ++i )
	627	{
	628	iso.CachedAs.reserve( ipISO, nelem, i, iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] );
	629	for( i1 = 0; i1 < iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem]; ++i1 )
	630	{
	631	/* allocate two spaces delta L +/- 1 */
	632	iso.CachedAs.reserve( ipISO, nelem, i, i1, 2 );
	633	}
	634	}
609	635
610	636	iso.QuantumNumbers2Index.reserve( ipISO, nelem, iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1 );
…	…
677	703	iso.Boltzmann.alloc();
678	704	iso.DielecRecombVsTemp.alloc();
	705	iso.CachedAs.alloc();
679	706	iso.QuantumNumbers2Index.alloc();
680	707	iso.BranchRatio.alloc();
681	708	iso.CascadeProb.alloc();
682	709
	710	iso.bnl_effective.alloc( iso.QuantumNumbers2Index.clone() );
	711
683	712	iso.DielecRecombVsTemp.zero();
684	713	iso.QuantumNumbers2Index.invalidate();
	714	iso.bnl_effective.invalidate();
685	715	iso.BranchRatio.invalidate();
686	716	iso.CascadeProb.invalidate();
…	…
759	789	}
760	790
761		/* malloc space for random error generation, but only if flag
762		* set true for at least one iso sequence. */
763		if( iso.lgRandErrGen[ipH_LIKE] \|\| iso.lgRandErrGen[ipHE_LIKE] )
764		{
765		/* This assert is nothing more than a reminder to adjust the above if clause
766		* if more iso sequences are added. */
767		ASSERT( NISO==2 );
768		iso.Error.reserve( NISO );
769		iso.SigmaCascadeProb.reserve( NISO );
770
771		/* the hydrogen and helium like iso-sequences */
772		for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
773		{
774		if( iso.lgRandErrGen[ipISO] )
775		{
776		iso.Error.reserve( ipISO, LIMELM );
777		iso.SigmaCascadeProb.reserve( ipISO, LIMELM );
778
779		for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
780		{
781		if( nelem == ipHELIUM \|\| dense.lgElmtOn[nelem] )
782		{
783
784		/* Here we add one slot for rates involving the continuum. */
785		iso.Error.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem]+1 );
786
787		for( i = 1; i< iso.numLevels_max[ipISO][nelem] + 1; i++ )
788		{
789		iso.Error.reserve( ipISO, nelem, i, i+1 );
790
791		for( j = 0; j<i+1; j++ )
792		iso.Error.reserve( ipISO, nelem, i, j, 3 );
793		}
794
795		/* Uncertainty in cascade probability and effective recombination,
796		* only when error generation is turned on. */
797		iso.SigmaCascadeProb.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
798		for( i=0; i < iso.numLevels_max[ipISO][nelem]; i++ )
799		/* error in cascade probabilities, for all levels */
800		iso.SigmaCascadeProb.reserve( ipISO, nelem, i, i+1 );
801		}
	791	/* malloc space for random error generation, if appropriate flags set. */
	792	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	793	{
	794	static bool lgNotSetYet = true;
	795
	796	if( iso.lgRandErrGen[ipISO] && lgNotSetYet )
	797	{
	798	iso.Error.reserve( NISO );
	799	iso.SigmaCascadeProb.reserve( NISO );
	800	lgNotSetYet = false;
	801	}
	802
	803	if( iso.lgRandErrGen[ipISO] )
	804	{
	805	ASSERT( !lgNotSetYet );
	806	iso.Error.reserve( ipISO, LIMELM );
	807	iso.SigmaCascadeProb.reserve( ipISO, LIMELM );
	808
	809	for( nelem=ipISO; nelem < LIMELM; ++nelem )
	810	{
	811	if( nelem == ipHELIUM \|\| dense.lgElmtOn[nelem] )
	812	{
	813	/* Here we add one slot for rates involving the continuum. */
	814	iso.Error.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem]+1 );
	815
	816	for( i = 1; i< iso.numLevels_max[ipISO][nelem] + 1; i++ )
	817	{
	818	iso.Error.reserve( ipISO, nelem, i, i+1 );
	819
	820	for( j = 0; j<i+1; j++ )
	821	iso.Error.reserve( ipISO, nelem, i, j, 3 );
	822	}
	823
	824	/* Uncertainty in cascade probability and effective recombination,
	825	* only when error generation is turned on. */
	826	iso.SigmaCascadeProb.reserve( ipISO, nelem, iso.numLevels_max[ipISO][nelem] );
	827	for( i=0; i < iso.numLevels_max[ipISO][nelem]; i++ )
	828	/* error in cascade probabilities, for all levels */
	829	iso.SigmaCascadeProb.reserve( ipISO, nelem, i, i+1 );
802	830	}
803	831	}
…	…
909	937	/* fill in StatesElem and iso.QuantumNumbers2Index arrays. */
910	938	/* this loop is over quantum number n */
911		for( in = 1L; in <= iso.n_HighestResolved_max[ip~~HE_LIKE~~][nelem]; ++in )
	939	for( in = 1L; in <= iso.n_HighestResolved_max[ipISO][nelem]; ++in )
912	940	{
913	941	for( il = 0L; il < in; ++il )
…	…
927	955	if( is == 1 )
928	956	{
929		StatesElem[ip~~HE_LIKE~~][nelem][i].n = in;
930		StatesElem[ip~~HE_LIKE~~][nelem][i].S = is;
931		StatesElem[ip~~HE_LIKE~~][nelem][i].l = il;
	957	StatesElem[ipISO][nelem][i].n = in;
	958	StatesElem[ipISO][nelem][i].S = is;
	959	StatesElem[ipISO][nelem][i].l = il;
932	960	/* this is not a typo, J=L for singlets. */
933		StatesElem[ip~~HE_LIKE~~][nelem][i].j = il;
934		iso.QuantumNumbers2Index[ip~~HE_LIKE~~][nelem][in][il][is] = i;
	961	StatesElem[ipISO][nelem][i].j = il;
	962	iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i;
935	963	++i;
936	964	}
…	…
941	969	do
942	970	{
943		StatesElem[ip~~HE_LIKE~~][nelem][i].n = in;
944		StatesElem[ip~~HE_LIKE~~][nelem][i].S = is;
945		StatesElem[ip~~HE_LIKE~~][nelem][i].l = il;
946		StatesElem[ip~~HE_LIKE~~][nelem][i].j = ij;
947		iso.QuantumNumbers2Index[ip~~HE_LIKE~~][nelem][in][il][is] = i;
	971	StatesElem[ipISO][nelem][i].n = in;
	972	StatesElem[ipISO][nelem][i].S = is;
	973	StatesElem[ipISO][nelem][i].l = il;
	974	StatesElem[ipISO][nelem][i].j = ij;
	975	iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i;
948	976	++i;
949	977	++ij;
…	…
953	981	else
954	982	{
955		StatesElem[ip~~HE_LIKE~~][nelem][i].n = in;
956		StatesElem[ip~~HE_LIKE~~][nelem][i].S = is;
957		StatesElem[ip~~HE_LIKE~~][nelem][i].l = il;
958		StatesElem[ip~~HE_LIKE~~][nelem][i].j = -1L;
959		iso.QuantumNumbers2Index[ip~~HE_LIKE~~][nelem][in][il][is] = i;
	983	StatesElem[ipISO][nelem][i].n = in;
	984	StatesElem[ipISO][nelem][i].S = is;
	985	StatesElem[ipISO][nelem][i].l = il;
	986	StatesElem[ipISO][nelem][i].j = -1L;
	987	iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = i;
960	988	++i;
961	989	}
…	…
965	993	if( in > 1L )
966	994	{
967		StatesElem[ip~~HE_LIKE~~][nelem][i].n = in;
968		StatesElem[ip~~HE_LIKE~~][nelem][i].S = 1L;
969		StatesElem[ip~~HE_LIKE~~][nelem][i].l = 1L;
970		StatesElem[ip~~HE_LIKE~~][nelem][i].j = 1L;
971		iso.QuantumNumbers2Index[ip~~HE_LIKE~~][nelem][in][1][1] = i;
	995	StatesElem[ipISO][nelem][i].n = in;
	996	StatesElem[ipISO][nelem][i].S = 1L;
	997	StatesElem[ipISO][nelem][i].l = 1L;
	998	StatesElem[ipISO][nelem][i].j = 1L;
	999	iso.QuantumNumbers2Index[ipISO][nelem][in][1][1] = i;
972	1000	++i;
973	1001	}
974	1002	}
975	1003	/* now do the collapsed levels */
976		in = iso.n_HighestResolved_max[ip~~HE_LIKE~~][nelem] + 1;
977		for( level = i; level< iso.numLevels_max[ip~~HE_LIKE~~][nelem]; ++level)
978		{
979		StatesElem[ip~~HE_LIKE~~][nelem][level].n = in;
980		StatesElem[ip~~HE_LIKE~~][nelem][level].S = -LONG_MAX;
981		StatesElem[ip~~HE_LIKE~~][nelem][level].l = -LONG_MAX;
982		StatesElem[ip~~HE_LIKE~~][nelem][level].j = -1;
	1004	in = iso.n_HighestResolved_max[ipISO][nelem] + 1;
	1005	for( level = i; level< iso.numLevels_max[ipISO][nelem]; ++level)
	1006	{
	1007	StatesElem[ipISO][nelem][level].n = in;
	1008	StatesElem[ipISO][nelem][level].S = -LONG_MAX;
	1009	StatesElem[ipISO][nelem][level].l = -LONG_MAX;
	1010	StatesElem[ipISO][nelem][level].j = -1;
983	1011	/* Point every l and s to same index for collapsed levels. */
984	1012	for( il = 0; il < in; ++il )
…	…
986	1014	for( is = 1; is <= 3; is += 2 )
987	1015	{
988		iso.QuantumNumbers2Index[ip~~HE_LIKE~~][nelem][in][il][is] = level;
	1016	iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] = level;
989	1017	}
990	1018	}
…	…
994	1022
995	1023	/* confirm that we did not overrun the array */
996		ASSERT( i <= iso.numLevels_max[ip~~HE_LIKE~~][nelem] );
	1024	ASSERT( i <= iso.numLevels_max[ipISO][nelem] );
997	1025
998	1026	/* confirm that n is positive and not greater than the max n. */
999		ASSERT( (in > 0) && (in < (iso.n_HighestResolved_max[ip~~HE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE~~][nelem] + 1) ) );
1000
1001		/* Verify StatesElem[ip~~HE_LIKE~~] and iso.QuantumNumbers2Index[ipISO] agree in all cases */
1002		for( in = 2L; in <= iso.n_HighestResolved_max[ip~~HE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE~~][nelem]; ++in )
	1027	ASSERT( (in > 0) && (in < (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1) ) );
	1028
	1029	/* Verify StatesElem[ipISO] and iso.QuantumNumbers2Index[ipISO] agree in all cases */
	1030	for( in = 2L; in <= iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem]; ++in )
1003	1031	{
1004	1032	for( il = 0L; il < in; ++il )
…	…
1010	1038	continue;
1011	1039
1012		ASSERT( StatesElem[ip~~HE_LIKE][nelem][ iso.QuantumNumbers2Index[ipHE_LIKE~~][nelem][in][il][is] ].n == in );
1013		if( in <= iso.n_HighestResolved_max[ip~~HE_LIKE~~][nelem] )
	1040	ASSERT( StatesElem[ipISO][nelem][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].n == in );
	1041	if( in <= iso.n_HighestResolved_max[ipISO][nelem] )
1014	1042	{
1015	1043	/* Must only check these for resolved levels...
1016		* collapsed levels in StatesElem[ip~~HE_LIKE~~] have pointers for l and s that will blow if used. */
1017		ASSERT( StatesElem[ip~~HE_LIKE][nelem][ iso.QuantumNumbers2Index[ipHE_LIKE~~][nelem][in][il][is] ].l == il );
1018		ASSERT( StatesElem[ip~~HE_LIKE][nelem][ iso.QuantumNumbers2Index[ipHE_LIKE~~][nelem][in][il][is] ].S == is );
	1044	* collapsed levels in StatesElem[ipISO] have pointers for l and s that will blow if used. */
	1045	ASSERT( StatesElem[ipISO][nelem][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].l == il );
	1046	ASSERT( StatesElem[ipISO][nelem][ iso.QuantumNumbers2Index[ipISO][nelem][in][il][is] ].S == is );
1019	1047	}
1020	1048	}
…	…
1142	1170
1143	1171	t->Hi->lifetime = iso_state_lifetime( ipISO, nelem, nHi, 1 );
1144
1145		~~if( ipISO == ipHE_LIKE )~~
1146		{
1147		~~/* iso_state_lifetime is not spin specific, must exclude helike triplet here. */~~
1148		~~t->Hi->lifetime /= 3.;~~
1149		~~/* there is also necessary to correct the helike lifetimes */~~
1150		~~t->Hi->lifetime = 1.1722 pow( (double)nelem, 0.1 );~~
1151		}
1152
1153		~~/* would probably need a new lifetime algorithm for any other iso sequences. */~~
1154		~~ASSERT( ipISO <= ipHE_LIKE );~~
1155	1172
1156	1173	t->Emis->dampXvel = (realnum)( 1.f / t->Hi->lifetime / PI4 / t->EnergyWN );
…	…
1213	1230	0.5 * eps2 - 0.025 * eps2 * eps2 ) );
1214	1231
	1232	if( ipISO == ipHE_LIKE )
	1233	{
	1234	/* iso_state_lifetime is not spin specific, must exclude helike triplet here. */
	1235	tau /= 3.;
	1236	/* this is also necessary to correct the helike lifetimes */
	1237	tau = 1.1722 pow( (double)nelem, 0.1 );
	1238	}
	1239
	1240	/* would probably need a new lifetime algorithm for any other iso sequences. */
	1241	ASSERT( ipISO <= ipHE_LIKE );
1215	1242	ASSERT( tau > 0. );
1216	1243
…	…
1219	1246
1220	1247	/* calculate cascade probabilities, branching ratios, and associated errors. */
1221		~~STATIC void iso_cascade( void~~ )
	1248	void iso_cascade( long ipISO, long nelem )
1222	1249	{
1223	1250	/* The sum of all A's coming out of a given n,
1224	1251	* Below we assert a monotonic trend. */
1225		~~multi_arr<double,3>~~ SumAPerN;
1226
1227		long int i, j, ip~~ISO, ipLo, ipHi, nelem~~;
	1252	double *SumAPerN;
	1253
	1254	long int i, j, ipLo, ipHi;
1228	1255
1229	1256	DEBUG_ENTRY( "iso_cascade()" );
1230	1257
1231		SumAPerN.reserve( NISO );
1232		for( ipISO = ipH_LIKE; ipISO < NISO; ipISO++ )
1233		{
1234		SumAPerN.reserve( ipISO, LIMELM );
1235		for( nelem=ipISO; nelem < LIMELM; ++nelem )
1236		{
1237		if( nelem == ipISO \|\| dense.lgElmtOn[nelem] )
1238		{
1239		SumAPerN.reserve( ipISO, nelem, iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1 );
1240		}
1241		}
1242		}
1243
1244		SumAPerN.alloc();
1245		SumAPerN.zero();
1246
1247		/* loop over iso sequences */
1248		for( ipISO = ipH_LIKE; ipISO < NISO; ipISO++ )
1249		{
1250		for( nelem = ipISO; nelem < LIMELM; nelem++ )
1251		{
1252		/* must always do helium even if turned off */
1253		if( nelem == ipISO \|\| dense.lgElmtOn[nelem] )
1254		{
1255		/* Initialize some ground state stuff, easier here than in loops. */
1256		iso.CascadeProb[ipISO][nelem][0][0] = 1.;
	1258	SumAPerN = ((double)MALLOC( (size_t)(iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1 )sizeof(double )));
	1259	memset( SumAPerN, 0, (iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] + 1 )*sizeof(double ) );
	1260
	1261	/* Initialize some ground state stuff, easier here than in loops. */
	1262	iso.CascadeProb[ipISO][nelem][0][0] = 1.;
	1263	if( iso.lgRandErrGen[ipISO] )
	1264	{
	1265	iso.SigmaAtot[ipISO][nelem][0] = 0.;
	1266	iso.SigmaCascadeProb[ipISO][nelem][0][0] = 0.;
	1267	}
	1268
	1269	/***************************************************************************/
	1270	/**** Cascade probabilities, Branching ratios, and associated errors ***/
	1271	/***************************************************************************/
	1272	for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
	1273	{
	1274	/** Cascade probabilities are as defined in Robbins 68,
	1275	* generalized here for cascade probability for any iso sequence.
	1276	* >>refer He triplets Robbins, R.R. 1968, ApJ 151, 497R
	1277	* >>refer He triplets Robbins, R.R. 1968a, ApJ 151, 511R */
	1278
	1279	/* initialize variables. */
	1280	iso.CascadeProb[ipISO][nelem][ipHi][ipHi] = 1.;
	1281	iso.CascadeProb[ipISO][nelem][ipHi][0] = 0.;
	1282	iso.BranchRatio[ipISO][nelem][ipHi][0] = 0.;
	1283
	1284	if( iso.lgRandErrGen[ipISO] )
	1285	{
	1286	iso.SigmaAtot[ipISO][nelem][ipHi] = 0.;
	1287	iso.SigmaCascadeProb[ipISO][nelem][ipHi][ipHi] = 0.;
	1288	}
	1289
	1290	long ipLoStart = 0;
	1291	if( opac.lgCaseB && L_(ipHi)==1 && (ipISO==ipH_LIKE \|\| S_(ipHi)==1) )
	1292	ipLoStart = 1;
	1293
	1294	for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
	1295	{
	1296	if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
	1297	{
	1298	iso.CascadeProb[ipISO][nelem][ipHi][ipLo] = 0.;
	1299	iso.BranchRatio[ipISO][nelem][ipHi][ipLo] = 0.;
	1300	continue;
	1301	}
	1302
	1303	iso.CascadeProb[ipISO][nelem][ipHi][ipLo] = 0.;
	1304	iso.BranchRatio[ipISO][nelem][ipHi][ipLo] =
	1305	Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul *
	1306	StatesElem[ipISO][nelem][ipHi].lifetime;
	1307
	1308	ASSERT( iso.BranchRatio[ipISO][nelem][ipHi][ipLo] <= 1.0000001 );
	1309
	1310	SumAPerN[N_(ipHi)] += Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul;
	1311
	1312	/* there are some negative energy transitions, where the order
	1313	* has changed, but these are not optically allowed, these are
	1314	* same n, different L, forbidden transitions */
	1315	ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN > 0. \|\|
	1316	Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA );
	1317
	1318	if( iso.lgRandErrGen[ipISO] )
	1319	{
	1320	ASSERT( iso.Error[ipISO][nelem][ipHi][ipLo][IPRAD] >= 0. );
	1321	/* Uncertainties in A's are added in quadrature, square root is taken below. */
	1322	iso.SigmaAtot[ipISO][nelem][ipHi] +=
	1323	pow( iso.Error[ipISO][nelem][ipHi][ipLo][IPRAD] *
	1324	(double)Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul, 2. );
	1325	}
	1326	}
	1327
	1328	if( iso.lgRandErrGen[ipISO] )
	1329	{
	1330	/* Uncertainties in A's are added in quadrature above, square root taken here. */
	1331	iso.SigmaAtot[ipISO][nelem][ipHi] = sqrt( iso.SigmaAtot[ipISO][nelem][ipHi] );
	1332	}
	1333
	1334	/* cascade probabilities */
	1335	for( ipLo=0; ipLo<ipHi; ipLo++ )
	1336	{
	1337	double SigmaA, SigmaCul = 0.;
	1338
	1339	for( i=ipLo; i<ipHi; i++ )
	1340	{
	1341	iso.CascadeProb[ipISO][nelem][ipHi][ipLo] += iso.BranchRatio[ipISO][nelem][ipHi][i] * iso.CascadeProb[ipISO][nelem][i][ipLo];
1257	1342	if( iso.lgRandErrGen[ipISO] )
1258	1343	{
1259		iso.SigmaAtot[ipISO][nelem][0] = 0.;
1260		iso.SigmaCascadeProb[ipISO][nelem][0][0] = 0.;
1261		}
1262
1263		/***************************************************************************/
1264		/**** Cascade probabilities, Branching ratios, and associated errors ***/
1265		/***************************************************************************/
1266		for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
1267		{
1268		/** Cascade probabilities are as defined in Robbins 68,
1269		* generalized here for cascade probability for any iso sequence.
1270		* >>refer He triplets Robbins, R.R. 1968, ApJ 151, 497R
1271		* >>refer He triplets Robbins, R.R. 1968a, ApJ 151, 511R */
1272
1273		/* initialize variables. */
1274		iso.CascadeProb[ipISO][nelem][ipHi][ipHi] = 1.;
1275		iso.BranchRatio[ipISO][nelem][ipHi][0] = 0.;
1276
1277		if( iso.lgRandErrGen[ipISO] )
1278		{
1279		iso.SigmaAtot[ipISO][nelem][ipHi] = 0.;
1280		iso.SigmaCascadeProb[ipISO][nelem][ipHi][ipHi] = 0.;
1281		}
1282
1283		for( ipLo=0; ipLo<ipHi; ipLo++ )
1284		{
1285		if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
	1344	if( Transitions[ipISO][nelem][ipHi][i].Emis->Aul > iso.SmallA )
	1345	{
	1346	/* Uncertainties in A's and cascade probabilities */
	1347	SigmaA = iso.Error[ipISO][nelem][ipHi][i][IPRAD] *
	1348	Transitions[ipISO][nelem][ipHi][i].Emis->Aul;
	1349	SigmaCul +=
	1350	pow(SigmaAiso.CascadeProb[ipISO][nelem][i][ipLo]StatesElem[ipISO][nelem][ipHi].lifetime, 2.) +
	1351	pow(iso.SigmaAtot[ipISO][nelem][ipHi]iso.BranchRatio[ipISO][nelem][ipHi][i]iso.CascadeProb[ipISO][nelem][i][ipLo]*StatesElem[ipISO][nelem][ipHi].lifetime, 2.) +
	1352	pow(iso.SigmaCascadeProb[ipISO][nelem][i][ipLo]*iso.BranchRatio[ipISO][nelem][ipHi][i], 2.);
	1353	}
	1354	}
	1355	}
	1356	if( iso.lgRandErrGen[ipISO] )
	1357	{
	1358	SigmaCul = sqrt(SigmaCul);
	1359	iso.SigmaCascadeProb[ipISO][nelem][ipHi][ipLo] = SigmaCul;
	1360	}
	1361	}
	1362	}
	1363
	1364	/************************************************************************/
	1365	/* Allowed decay conversion probabilities. See Robbins68b, Table 1. */
	1366	/************************************************************************/
	1367	{
	1368	enum {DEBUG_LOC=false};
	1369
	1370	if( DEBUG_LOC && (nelem == ipHELIUM) && (ipISO==ipHE_LIKE) )
	1371	{
	1372	/* To output Bm(n,l; ipLo), set ipLo, hi_l, and hi_s accordingly. */
	1373	long int hi_l,hi_s;
	1374	double Bm;
	1375
	1376	/* these must be set for following output to make sense
	1377	* as is, a dangerous bit of code - set NaN for safety */
	1378	hi_s = -100000;
	1379	hi_l = -100000;
	1380	ipLo = -100000;
	1381	/* tripS to 2^3P */
	1382	//hi_l = 0, hi_s = 3, ipLo = ipHe2p3P0;
	1383
	1384	/* tripD to 2^3P */
	1385	//hi_l = 2, hi_s = 3, ipLo = ipHe2p3P0;
	1386
	1387	/* tripP to 2^3S */
	1388	//hi_l = 1, hi_s = 3, ipLo = ipHe2s3S;
	1389
	1390	ASSERT( hi_l != StatesElem[ipISO][nelem][ipLo].l );
	1391
	1392	fprintf(ioQQQ,"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo);
	1393	fprintf(ioQQQ,"m\t2\t\t3\t\t4\t\t5\t\t6\n");
	1394
	1395	for( ipHi=ipHe2p3P2; ipHi<iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem]; ipHi++ )
	1396	{
	1397	/* Pick out excitations from metastable 2tripS to ntripP. */
	1398	if( (StatesElem[ipISO][nelem][ipHi].l == 1) && (StatesElem[ipISO][nelem][ipHi].S == 3) )
	1399	{
	1400	fprintf(ioQQQ,"\n%ld\t",StatesElem[ipISO][nelem][ipHi].n);
	1401	j = 0;
	1402	Bm = 0;
	1403	for( i = ipLo; i<=ipHi; i++)
	1404	{
	1405	if( (StatesElem[ipISO][nelem][i].l == hi_l) && (StatesElem[ipISO][nelem][i].S == hi_s) )
1286	1406	{
1287		iso.CascadeProb[ipISO][nelem][ipHi][ipLo] = 0.;
1288		iso.BranchRatio[ipISO][nelem][ipHi][ipLo] = 0.;
1289		continue;
1290		}
1291
1292		iso.CascadeProb[ipISO][nelem][ipHi][ipLo] = 0.;
1293		iso.BranchRatio[ipISO][nelem][ipHi][ipLo] =
1294		Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul *
1295		StatesElem[ipISO][nelem][ipHi].lifetime;
1296
1297		ASSERT( iso.BranchRatio[ipISO][nelem][ipHi][ipLo] <= 1.0000001 );
1298
1299		SumAPerN[ipISO][nelem][N_(ipHi)] += Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul;
1300
1301		/* there are some negative energy transitions, where the order
1302		* has changed, but these are not optically allowed, these are
1303		* same n, different L, forbidden transitions */
1304		ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN > 0. \|\|
1305		Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA );
1306
1307		if( iso.lgRandErrGen[ipISO] )
1308		{
1309		ASSERT( iso.Error[ipISO][nelem][ipHi][ipLo][IPRAD] >= 0. );
1310		/* Uncertainties in A's are added in quadrature, square root is taken below. */
1311		iso.SigmaAtot[ipISO][nelem][ipHi] +=
1312		pow( iso.Error[ipISO][nelem][ipHi][ipLo][IPRAD] *
1313		(double)Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul, 2. );
1314		}
1315		}
1316
1317		if( iso.lgRandErrGen[ipISO] )
1318		{
1319		/* Uncertainties in A's are added in quadrature above, square root taken here. */
1320		iso.SigmaAtot[ipISO][nelem][ipHi] = sqrt( iso.SigmaAtot[ipISO][nelem][ipHi] );
1321		}
1322
1323		/* cascade probabilities */
1324		for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
1325		{
1326		double SigmaA, SigmaCul = 0.;
1327
1328		for( i=ipLo; i<ipHi; i++ )
1329		{
1330		iso.CascadeProb[ipISO][nelem][ipHi][ipLo] += iso.BranchRatio[ipISO][nelem][ipHi][