/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
#include "cddefines.h" 
#include "physconst.h" 
#include "lines_service.h"
#include "elementnames.h"
#include "taulines.h"
#include "path.h"
#include "trace.h"
#include "phycon.h"
#include "thermal.h"
#include "dense.h"
#include "iso.h"
#include "helike.h"
#include "helike_recom.h"
#include "helike_cs.h"
#include "helike_einsta.h"
#include "hydroeinsta.h"
#include "opacity.h"

static double EthRyd; /*kTRyd,,Xn_S59*/

static FILE *ofp;  

/* the designation of the levels, chLevel[n][string] */
static char **chLevel;

/* Ionization potentials (in wavenumber) for each ion in the iso-seq up to Z=30.		*/
/* These are exactly what you get if you take EionRYD below and multiply by RYD_INF.	*/
static double EionWN[LIMELM] =       
	/*198307.729760353*/
	/* the first, or H-like, element is not defined for he-like species */
	 {-DBL_MAX,
	 198310.6679     ,610003.839889137,1241136.72201499,2091948.45665631,3162116.52584231,
	 4452446.95015668,5962133.81875305,7692790.05069734,9645221.44709864,11814589.7994457,
	 14209766.0528639,16822685.5022862,19661412.9625169,22717883.6187518,26000162.0663204,
	 29508248.5246975,33234078.1790787,37185715.7345311,41363161.0813172,45766414.4389118,
	 50395475.4781030,55258409.0136949,60339085.8550283,65653635.1927626,71202056.8074231,
	 76976286.4328920,82984388.3352872,89194104.5722390,95726403.3055320};
	
/* Ionization potentials (in Rydbergs) for each ion in the iso-seq up to Z=30.	*/
/* These are exactly what you get if you take EionEV below and divide by EVRYD.	*/
static double EionRYD[LIMELM] =       
	/*1.807113*/
	/* the first, or H-like, element is not defined for he-like species */
	{-DBL_MAX,
	1.807387521,5.558764,11.310070,19.063237,28.815326,40.573682,54.330961,70.101861,
	87.893725,107.662464,129.488916,153.299590,179.167978,207.020588,236.930910,
	268.898946,302.851204,338.861175,376.928858,417.054255,459.237363,503.551674,
	549.850208,598.279945,648.840883,701.459535,756.209388,812.796486,872.323172};


/* experimental energies, in wavenumbers, for atomic helium */
#define NHE1LEVELS 111
/* These energies from Drake 1996, except these divided by c */
static double He1Energies[NHE1LEVELS] =
{0.0       , 159855.9734, 166277.4390, 169087.8298, 169086.8417, 169086.7652, 171134.8957,
183236.7908, 184864.8281, 185564.6657, 186101.5615, 186104.9656, 186209.3638, 190298.6619,
190940.6075, 191217.0826, 191444.4868, 191446.4547, 191451.8805, 191451.8964, 191492.7108,
193346.9900, 193663.5106, 193800.7280, 193917.1538, 193918.2888, 193921.1207, 193921.1298,
193921.6166, 193921.6209, 193942.4612, 194936.1184, 195114.8672, 195192.7542, 195260.0724,
195260.7694, 195262.4251, 195262.4307, 195262.7236, 195262.7261, 195262.7930, 195262.7947, 
195274.9074, 195868.2357, 195978.8938, 196027.3216, 196069.6730, 196070.1273, 196071.1763,
196071.1800, 196071.3686, 196071.3702, 196071.4141, 196071.4151, 196071.4283, 196071.4290,
196079.0865, 196461.3605, 196534.5628, 196566.7159, 196595.0620, 196595.3730, 196596.0785,
196596.0810, 196596.2092, 196596.2103, 196596.2404, 196596.2411, 196596.2503, 196596.2508,
196596.2541, 196596.2544, 196601.3992, 196861.9861, 196912.9014, 196935.3339, 196955.2261,
196955.4477, 196955.9445, 196955.9463, 196956.0373, 196956.0380, 196956.0595, 196956.0600,
196956.0666, 196956.0670, 196956.0693, 196956.0696, 196956.0705, 196956.0707, 196959.6917,
197145.2320, 197182.0643, 197198.3343, 197212.8252, 197212.9885, 197213.3513, 197213.3527,
197213.4194, 197213.4200, 197213.4358, 197213.4362, 197213.4411, 197213.4414, 197213.4431,
197213.4433, 197213.4440, 197213.4442, 197213.4445, 197213.4446, 197216.0885};
/* Last energy is 10^1P.	*/

/* >>chng 05 sep 06, RP - update energies of excited states */
#define NIONLEVELS 31
/* These energies are from Chianti 5, up to 5^1P 
 * -1 is flag to use quantum defect routines instead 
 * >>refer He-like	Energies	Landi et al., 2005; ApJSS, in press 
 * >>refer He-like	Energies	Dere, K. P. et al. 1997, 
 * >>refercon	 Astronomy and Astrophysics Suppl. Ser., Vol. 125, pp. 149-173 */
static double IonEnergies[LIMELM-2][NIONLEVELS] = 
{
{0.00,	476034.98,	491374.60,	494266.57,	494261.17,	494263.44,	501808.59,
	554754.45,	558777.88,	559501.16,	561243.67,	561273.62,	561752.82,
	579981.33,	581596.77,	581886.34,	582613.64,	582630.95,	582642.97,
	582644.04,	582830.11,	591184.26,	591989.55,	592134.36,	592504.32,
	592514.43,	592520.11,	592521.11,	-1.00	,	-1.00	,	592634.91},
{0.00,	956502.00,	981178.00,	983366.00,	983355.00,	983370.00,	997454.00,
	1121184.00,	1127705.00,	1128300.00,	1131383.00,	1131462.00,	1132390.00,
	1175295.00,	1178005.00,	1178174.00,	1179451.00,	1179495.00,	1179515.00,
	1179514.00,	1179830.00,	1199650.00,	-1.00	,	1201060.00,	1201702.00,
	1201800.00,	1201730.00,	1201742.00,	-1.00	,	-1.00	,	1201894.00},
{0.00,	1601545.00,	1635720.00,	1636938.00,	1636922.00,	1636975.00,	1657980.00,
	1882740.00,	1891790.00,	1892221.00,	1896710.00,	1896836.00,	1898063.00,
	1976420.00,	-1.00	,	1980291.11,	1982132.67,	1982220.00,	1982262.67,
	1982240.00,	1982762.00,	-1.00	,	-1.00	,	2020730.00,	2021700.00,
	-1.00	,	2021665.71,	2021760.00,	2021770.00,	2021770.00,	2022044.00},
{0.00,	2411262.00,	2455024.00,	2455162.74,	2455150.23,	2455286.01,	2483371.00,
	2839562.00,	2851180.00,	2851418.00,	2857309.67,	2857529.00,	2859375.00,
	2983541.00,	-1.00	,	2988359.00,	2990776.00,	2990923.00,	2990923.40,
	2990923.40,	2991710.00,	3048927.00,	-1.00	,	3051332.00,	3052589.00,
	3052656.00,	3052653.30,	3052653.30,	3052659.40,	3052659.40,	3053044.00},
{0.00,	3385890.00,	3439274.00,	3438312.46,	3438321.13,	3438612.15,	3473790.00,
	3991860.00,	-1.00	,	4006160.00,	4013460.00,	4013770.00,	4016390.00,
	4196800.00,	4202520.00,	4202620.00,	4205820.00,	4205830.00,	4205810.00,
	4205820.00,	4206810.00,	4290150.00,	4293020.00,	4293080.00,	4294570.00,
	4294670.00,	4294700.00,	4294700.00,	-1.00	,	-1.00	,	4296090.00},
{0.00,	4524640.00,	4588380.00,	4585620.76,	4585679.58,	4586231.19,	4629201.00,
	5338820.00,	5356420.00,	5355670.00,	5364422.67,	5365470.00,	5368550.00,
	5616340.00,	5623100.00,	5622600.00,	5626225.33,	5626670.00,	5626210.00,
	5626840.00,	5628100.00,	5742610.00,	-1.00	,	5745440.00,	5747509.33,
	5748230.00,	5747200.00,	5747820.00,	-1.00	,	-1.00	,	5748450.00},
{0.00,	5830040.00,	5903100.00,	5900600.00,	5900750.00,	5901700.00,	5949900.00,
	6885090.00,	6903270.00,	6902560.00,	6914073.33,	6915900.00,	6916590.00,
	7244270.00,	-1.00	,	7250390.00,	7255960.00,	7254240.00,	7256750.00,
	7257260.00,	7256370.00,	-1.00	,	-1.00	,	7410270.00,	7413940.00,
	7412290.00,	7414760.00,	7415300.00,	7414780.00,	-1.00	,	7414000.00},
{0.00,	7299940.00,	7382680.00,	7378205.53,	7378506.43,	7380050.00,	7436560.00,
	8623000.00,	8644880.00,	8644744.44,	8657128.67,	8662400.00,	8660530.00,
	9075200.00,	9084060.00,	9084141.11,	9090355.33,	-1.00	,	9089800.00,
	9094400.00,	9090630.00,	9282200.00,	9286650.00,	9286713.33,	9288500.00,
	-1.00	,	9289800.00,	9294400.00,	-1.00	,	-1.00	,	9290000.00},
{0.00,	8935337.00,	9027981.00,	9022354.10,	9022876.10,	9025284.70,	9088700.00,
	10558946.00,	10583431.00,	10583323.56,	10596783.40,	10597475.00,	10601080.00,
	11115065.00,	11124986.00,	11125102.78,	11130639.00,	11131017.00,	11131051.00,
	11131056.00,	11132393.00,	11369887.00,	11374868.00,	11374959.89,	11377767.00,
	11377984.00,	11377987.00,	11377991.00,	-1.00	,	-1.00	,	11378646.00},
{0.00,	10736136.00,	10838778.00,	10831985.83,	10832819.18,	10836391.13,	10906612.00,
	12691170.00,	12718304.00,	12718286.89,	12733392.33,	12734298.00,	12738006.00,
	13361991.00,	13372977.00,	13373168.22,	13379472.60,	13379830.00,	13379893.00,
	13379898.00,	13381265.00,	13669618.00,	13675137.00,	13675269.22,	13678467.13,
	13678680.00,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	13679363.00},
{0.00,	12703061.00,	12815760.00,	12807847.00,	12809088.00,	12814213.00,	12891081.00,
	15020463.00,	15050257.00,	15050434.00,	15067287.07,	15068371.00,	15072141.00,
	15816791.00,	15828851.00,	15829158.67,	15836125.13,	15836581.00,	-1.00	,
	-1.00	,	15838068.00,	16182216.00,	16188281.00,	16188471.33,	16192010.13,
	16192244.00,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	16192975.00},
{0.00,	14835945.00,	14958753.00,	14949756.42,	14951532.63,	14958690.57,	15042040.00,
	17546734.00,	17579166.00,	17579686.44,	17598406.93,	17599605.00,	17603422.00,
	18479389.00,	18492532.00,	18493007.56,	18500821.00,	18501245.00,	-1.00	,
	-1.00	,	18502736.00,	18907613.00,	18914246.00,	18914502.78,	18918476.30,
	18918694.00,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	18919421.00},
{0.00,	17135768.00,	17268828.00,	17258746.00,	17261164.00,	17270908.00,	17360546.00,
	20271100.00,	20306284.00,	20307209.11,	20327865.87,	20329412.00,	20332952.00,
	21350958.00,	21365192.00,	21365892.89,	21374428.00,	21375044.00,	21375302.00,
	21375302.00,	21376454.00,	21846994.00,	21854144.00,	21854552.89,	21858894.80,
	21859210.00,	21859340.00,	21859340.00,	-1.00	,	-1.00	,	21859464.00},
{0.00,	19602076.00,	19745473.00,	19734297.61,	19737518.84,	19750576.04,	19846285.00,
	23193163.00,	23231087.00,	23232596.56,	23255347.53,	23257195.00,	23260416.00,
	24431101.00,	24446439.00,	24447429.78,	24456981.30,	24457576.00,	-1.00	,
	-1.00	,	24458842.00,	24999972.00,	25007605.00,	25008238.11,	25013102.40,
	25013407.00,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	25014007.00},
{0.00,	22236180.00,	22390000.00,	22377820.00,	22381940.00,	22399100.00,	22500680.00,
	26314360.00,	26355050.00,	26357324.44,	26382328.67,	26384530.00,	26387270.00,
	27720900.00,	27738000.00,	27738966.67,	27749520.00,	27750400.00,	27757331.14,
	27757178.00,	27751600.00,	28367700.00,	28376500.00,	28376977.78,	28382286.67,
	28382800.00,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	28383400.00},
{0.00,	25038230.00,	25202480.00,	25189388.10,	25194588.99,	25216810.57,	25323950.00,
	29634850.00,	29678210.00,	29681541.11,	29713920.00,	29715070.00,	29713740.00,
	31221700.00,	31239280.00,	31241087.78,	-1.00	,	-1.00	,	-1.00	,
	-1.00	,	31254280.00,	31951370.00,	31960150.00,	31961182.22,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	31967860.00},
{0.00,	28004980.00,	28180480.00,	28165880.00,	28172670.00,	28200800.00,	28312910.00,
	33151930.00,	33198090.00,	33202622.22,	-1.00	,	-1.00	,	33237140.00,
	-1.00	,	34948420.00,	-1.00	,	-1.00	,	-1.00	,	-1.00	,
	-1.00	,	34964920.00,	35747360.00,	35756712.00,	35758032.00,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	35765192.00},
{0.00,	31142150.00,	31328450.00,	31312818.51,	31320486.08,	31356326.72,	31473810.00,
	36870940.00,	36919930.00,	36925900.00,	-1.00	,	-1.00	,	36962850.00,
	38850670.00,	38870530.00,	38873536.67,	-1.00	,	-1.00	,	-1.00	,
	-1.00	,	38888680.00,	39761380.00,	39771310.00,	39772968.89,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	39780630.00},
{0.00,	34448120.00,	34645360.00,	34628770.00,	34638550.00,	34682810.00,	34805000.00,
	40790620.00,	40842480.00,	40850158.89,	-1.00	,	-1.00	,	40889690.00,
	42983370.00,	43004390.00,	43008165.56,	-1.00	,	-1.00	,	-1.00	,
	-1.00	,	43024380.00,	43992240.00,	44002740.00,	44004803.33,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	44013010.00},
{0.00,	37923880.00,	38131640.00,	38114760.00,	38125260.00,	38180620.00,	38308340.00,
	44911910.00,	44966970.00,	44976472.22,	45016028.67,	45021140.00,	45018670.00,
	47328500.00,	47351600.00,	47355644.44,	47372326.67,	47374500.00,	-1.00	,
	-1.00	,	47373500.00,	48440800.00,	48452600.00,	48454581.11,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	48463700.00},
{0.00,	41568880.00,	41787830.00,	41770130.00,	41782100.00,	41849950.00,	41982380.00,
	49234710.00,	49292760.00,	49304574.44,	49347922.67,	49353910.00,	49349740.00,
	51886600.00,	51910900.00,	51915955.56,	51934113.33,	51936800.00,	-1.00	,
	-1.00	,	51935100.00,	53107300.00,	53119700.00,	53122217.78,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	53132000.00},
{0.00,	45384110.00,	45614410.00,	45595910.00,	45609360.00,	45691820.00,	45828830.00,
	53760100.00,	53821190.00,	53835660.00,	53883181.33,	53890160.00,	53884060.00,
	56658500.00,	56684100.00,	56690288.89,	56710306.67,	56713200.00,	-1.00	,
	-1.00	,	56710700.00,	57992700.00,	58005800.00,	58008904.44,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	58019400.00},
{0.00,	49370240.00,	49612040.00,	49592800.00,	49607700.00,	49707130.00,	49848620.00,
	58488800.00,	58553000.00,	58570522.22,	58622620.00,	58630700.00,	58622500.00,
	61644700.00,	61671800.00,	61679233.33,	61701013.33,	61704700.00,	-1.00	,
	-1.00	,	61701200.00,	63097900.00,	63111600.00,	63115327.78,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	63126600.00},
{0.00,	53527760.00,	53781230.00,	53761280.00,	53777570.00,	53896550.00,	54042490.00,
	63421700.00,	63489000.00,	63509966.67,	63567120.00,	63576500.00,	63565800.00,
	66846900.00,	66875000.00,	66884011.11,	66908120.00,	66912100.00,	-1.00	,
	-1.00	,	66907600.00,	68423800.00,	68438100.00,	68442730.00,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	68454800.00},
{0.00,	57857380.00,	58122700.00,	58102090.00,	58119680.00,	58261180.00,	58411430.00,
	68560000.00,	68630600.00,	68655444.44,	68718160.00,	68728900.00,	68715500.00,
	72266000.00,	72295500.00,	72306144.44,	72332586.67,	72337100.00,	-1.00	,
	-1.00	,	72331500.00,	73972000.00,	73987000.00,	73992318.89,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	74005400.00},
{0.00,	62358960.00,	62637200.00,	62615022.80,	62633778.00,	62800884.44,	62952670.00,
	73903340.00,	73976370.00,	74005924.44,	-1.00	,	-1.00	,	74070580.00,
	77900890.00,	77930480.00,	77943808.89,	-1.00	,	-1.00	,	-1.00	,
	-1.00	,	77970500.00,	79740940.00,	79755710.00,	79762730.00,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	79776290.00},
{0.00,	67035380.00,	67324970.00,	67303150.00,	67322950.00,	67519170.00,	67678080.00,
	79453200.00,	79530300.00,	79564500.00,	79640046.67,	79654200.00,	79634300.00,
	83754400.00,	83786500.00,	83801155.56,	83832993.33,	83838900.00,	-1.00	,
	-1.00	,	83830600.00,	85734300.00,	85750700.00,	85759175.00,	-1.00	,
	-1.00	,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	85773200.00},
{0.00,	71886300.00,	72188400.00,	72166200.00,	72186600.00,	72415600.00,	72579000.00,
	85212700.00,	85293200.00,	85332955.56,	85415826.67,	85431900.00,	85408300.00,
	89828600.00,	89862000.00,	89879022.22,	89913666.67,	89920800.00,	-1.00	,
	-1.00	,	89910900.00,	91953400.00,	91970500.00,	91979222.22,	91995073.33,
	92000500.00,	-1.00	,	-1.00	,	-1.00	,	-1.00	,	91994500.00}};

/*he_energy calculates energy of a given level.	*/
static double he_energy(double Eff_n,long int nelem, long int ipLo )
{
	double Ef;
	
	/* for atomic helium use experimental energies */
	if( nelem == ipHELIUM && ipLo < NHE1LEVELS )
	{
		Ef = EionWN[ipHELIUM] - He1Energies[ipLo];
	} 
	/* >>chng 05 sep 06, RP - update energies of excited states */
	else if( nelem > ipHELIUM && nelem <= ipZINC &&
		ipLo < NIONLEVELS && IonEnergies[nelem-2][ipLo] > 0. )
	{
		Ef = EionWN[nelem] - IonEnergies[nelem-2][ipLo];
	}
	else
	/* Calculate energy difference for all other trans.	*/
	{
		Ef = RYD_INF * (nelem/Eff_n) * (nelem/Eff_n);  
	}

	ASSERT(Ef > 0.);
	return Ef;
}

/*defect - calculate quantum defect for a given level and nuclear charge. */
static double defect(long int nelem, long int ipLo)  
{  
	/* The quantum defect, and parameters a,b, and c	*/
	double qd,a,b,c;	

	double HeDefectAsymptotes[2][10] = {
		{1.40005E-01,-1.20673E-02,2.08056E-03,4.21484E-04,1.14868E-04,
			4.08648E-05,1.73548E-05,8.33891E-06,4.39680E-06,2.42075E-06},
		{2.97063E-01,6.81567E-02,2.82381E-03,4.27703E-04,1.17319E-04,
			4.25254E-05,1.85549E-05,9.24641E-06,5.30882E-06,3.02877E-06}
	};

	/* Parameters for fits to quantum defects for	*/
	/* P triplet and S orbitals. The dimensions are	*/
	/*	first: l									*/
	/*	second: n									*/
	/* 	third: parameters a,b,and c.				*/
	double param[3][4][3]=  
	{ 
		{{0.6451941,0.3119437,-1.2722842},		/* ^3S	*/
			{0.7664874,0.3455675,-1.3976462}, 
			{0.8247101,0.3603131,-1.4520500}, 
			{0.8878402,0.3714450,-1.4995732}}, 

		{{1.4203514,0.5311096,-2.6728087},		/* ^1S	*/
			{1.5733513,0.5997339,-2.9253834}, 
			{1.4531025,0.5924751,-2.8662756}, 
			{1.6038999,0.6342552,-3.0298071}}, 

		{{-2.2323488,0.0890840,-0.5166053},		/* ^3P	*/
			{-2.0463691,0.1222081,-0.6672983}, 
			{-1.9904104,0.1328918,-0.7150879}, 
			{-1.9500974,0.1452111,-0.7649031}} 
	};   

	/* Because they cannot be fit to a funtion of the same form as the other orbitals,
	 * the P singlets are fit to a different function, with these parameters and dimensions */
	/*	first: n									*/
	/*	second: parameters a and b.					*/
	double P1[4][2]=
	{
		{-56.65245,-3.661923},
		{-52.03411,-4.941075},
		{-50.43744,-5.525750},
		{-49.45137,-5.908615}
	};

	long int n = iso.quant_desig[ipHE_LIKE][nelem][ipLo].n;
	long int lqn = iso.quant_desig[ipHE_LIKE][nelem][ipLo].l;
	long int s = iso.quant_desig[ipHE_LIKE][nelem][ipLo].s;
	
	ASSERT(n >= 1L);
	ASSERT(n > lqn);
	/* Only Helium and up, and only those turned on.	*/
	ASSERT((nelem >= ipHELIUM) && (nelem < LIMELM));       

	if ( n > iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
	{
		/* collapsed levels are assumed to have zero quantum defect. */
		qd = 0.;
	}
	else if ( nelem == ipHELIUM )
	{
		if( ipLo<NHE1LEVELS && n<=iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
		{
			qd = n-sqrt(0.999862926*RYD_INF/(EionWN[ipHELIUM] - He1Energies[ipLo]));
		} 
		else if ( lqn<=9 )
		{
			qd = HeDefectAsymptotes[s][lqn];
		}
		else if ( s == 0 )
		{
			qd = 0.0497*pow((double)lqn, -4.4303);
		}
		else
		{
			qd = 0.0656*pow((double)lqn, -4.5606);
		}
	}
	else if ( ipLo == ipHe1s1S )
	{
		/* Quantum defects for ground state are found from the rydberg
		 * equation, and the ionization potential of the ion.	
		 * the lint needs to see assert and not ASSERT */
		assert(nelem>ipHYDROGEN && nelem<LIMELM );
		qd = 1.0 - nelem * sqrt(1/EionRYD[nelem]);
	}
	else
	{
		/* For levels with n > 5, the quantum defect	*/
		/* is approximately the same as if n equaled 5.	*/
		if ( n > 5L )
		{
			n = 5L;
		}
		/* For P singlets	*/
		if ( lqn==1L && s==0L )
		{
			qd = 1./(P1[n-2][0] + P1[n-2][1] * (nelem+1) * log((double)nelem+1.) );
		}
		/* Defects for orbitals with l>2 are approximately equal to zero.	*/
		else if ( lqn < 2L )
		{
			a = param[2*(lqn+1)-s-1][n-2][0];  
			b = param[2*(lqn+1)-s-1][n-2][1];  
			c = param[2*(lqn+1)-s-1][n-2][2];  
			qd = exp((a+c*(nelem+1))/(1.0+b*(nelem+1)));  
		}
		/* This fit is a simplification of table 11.9 from 
		 * >>refer Helike	defects	Drake, G.W.F., editor.  Atomic, Molecular & Optical Physics Handbook.
		 * >>refercon	Chapter 11, "High Precision Calculations for Helium", G.W.F. Drake. 
		 * >>refercon	AIP Press: Woodbury, New York, 1996
		 * This will give quasi-real energies for all transitions, allowing a reasonable
		 * determination of which decays are zeroed due to being below the plasma frequency.
		 * The 1/nelem dependence is arbitray.  	*/
		else
		{
			ASSERT( lqn >= 2L );
			qd = ( ( 0.0612/(double)nelem ) / pow((double)lqn, 4.44) );
		}
	} 

	return qd;
}

/* Function he_assign assigns, in order of increasing energy, 				*/
/* quantum numbers n,l, and s to each element of vector iso.quant_desig[ipHE_LIKE].			*/
static void he_assign( long int nelem)
{
	/*	in, il, and is are dummies for n, l, and s.		*/
	/*	i is the index in the vector iso.quant_desig[ipHE_LIKE].			*/
	long int in, il, is, ij, i = 0 , level;
	
	/* this loop is over quantum number n */
	for ( in = 1L; in <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] ; ++in )
	{
		for ( il = 0L; il < in; ++il )
		{
			for ( is = 1L; is >= 0 ; --is )
			{
				/* All levels except singlet P follow the ordering scheme:	*/
				/*	lower l's have lower energy	*/
				/* 	triplets have lower energy	*/
				if ( (il == 1L) && (is == 0L) )
					continue;
				/* n = 1 has no triplet, of course.	*/
				if ( (in == 1L) && (is == 1L) )
					continue;
				ij = 0;
				do 
				{
					iso.quant_desig[ipHE_LIKE][nelem][i].n = in;
					iso.quant_desig[ipHE_LIKE][nelem][i].s = is;
					iso.quant_desig[ipHE_LIKE][nelem][i].l = il;
					QuantumNumbers2Index[nelem][in][il][is] = i;
					++i;
					++ij;
				}	while ( (in == 2) && (il == 1) && (is == 1) && (ij < 3) );
			}
		}
		/*	Insert singlet P at the end of every sequence for a given n.	*/
		if ( in > 1L )
		{
			iso.quant_desig[ipHE_LIKE][nelem][i].n = in;
			iso.quant_desig[ipHE_LIKE][nelem][i].s = 0L;
			iso.quant_desig[ipHE_LIKE][nelem][i].l = 1L;
			QuantumNumbers2Index[nelem][in][1][0] = i;
			++i;
		}
	}
	/* now do the collapsed levels */
	in = iso.n_HighestResolved_max[ipHE_LIKE][nelem] + 1;
	for( level = i; level< iso.numLevels_max[ipHE_LIKE][nelem]; ++level)
	{
		iso.quant_desig[ipHE_LIKE][nelem][level].n = in;
		iso.quant_desig[ipHE_LIKE][nelem][level].s = -LONG_MAX;
		iso.quant_desig[ipHE_LIKE][nelem][level].l = -LONG_MAX;
		/* Point every l and s to same index for collapsed levels.	*/
		for( il = 0; il < in; ++il )
		{
			for( is = 0; is <= 1; ++is )
			{
				QuantumNumbers2Index[nelem][in][il][is] = level;
			}
		}
		++in;
	}
	--in;

	/* confirm that we did not overrun the array */
	ASSERT( i <= iso.numLevels_max[ipHE_LIKE][nelem] );

	ASSERT( (in > 0) && (in < (iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] + 1) ) );


	/* Verify iso.quant_desig[ipHE_LIKE] and QuantumNumbers2Index agree in all cases	*/
	for ( in = 2L; in <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem]  ; ++in )
	{
		for ( il = 0L; il < in; ++il )
		{
			for ( is = 1L; is >= 0 ; --is )
			{
				/* Ground state is not triplicate.	*/
				if ( (in == 1L) && (is == 1L) )
					continue;
				
				ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ QuantumNumbers2Index[nelem][in][il][is] ].n == in );
				if( in <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
				{
					/* Must only check these for resolved levels...
					 * collapsed levels in iso.quant_desig[ipHE_LIKE] have pointers for l and s that will blow if used.	*/
					ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ QuantumNumbers2Index[nelem][in][il][is] ].l == il );
					ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ QuantumNumbers2Index[nelem][in][il][is] ].s == is );
				}
			}
		}
	}
	
	return;
}


/* This is just a routine to print out all As in a particular format, nothing more. */
static void printCustomAs( void )
{
	long int i, ipLo, familyIndex, nHi, upperIndex, nelem;
	nelem = ipHELIUM;

	if ( (ofp = fopen("CustomAs.txt","w")) == NULL )
		return;
	fprintf(ofp,"ipLo\tnLo\tlLo\tsLo\tlHi\t" );
	for( i = 2; i<= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM]; i++ )
	{
		fprintf( ofp,"%li\t", i);
	}
	fprintf( ofp, "\n" );
		
	for( ipLo = 0; ipLo< iso.numLevels_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] - 1; ipLo++ )
	{
		if( N_(ipLo) < iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] )
		{
			for( familyIndex = N_(ipLo)*(N_(ipLo)+1)+1; familyIndex < (N_(ipLo)+1)*(N_(ipLo)+2)+1; familyIndex++ )
			{
				if( N_(familyIndex) <= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] )
				{
					if( ( abs(L_(ipLo)-L_(familyIndex)) == 1 ) && (S_(ipLo) == S_(familyIndex)) )
					{
						fprintf( ofp, "%li\t%li\t%li\t%li\t%li\t",
							ipLo, N_(ipLo), L_(ipLo), S_(ipLo), L_(familyIndex) );

						for( nHi = 2; nHi < N_(familyIndex); nHi++ )
						{
							fprintf( ofp, "\t" );
						}
			
						for( nHi = N_(familyIndex); nHi <= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM]; nHi++ )
						{
							upperIndex = QuantumNumbers2Index[ipHELIUM][nHi][L_(familyIndex)][S_(familyIndex)];

							fprintf( ofp, "%.4e\t", EmisLines[ipHE_LIKE][nelem][upperIndex][ipLo].Aul );
						}
						fprintf( ofp, "\n" );
					}
				}
			}
		}
	}
}

/* create he-like series, called by ContCreatePointers */
void HeCreate(void)
{
	double **energies, *n_effs, **SumAPerN, **RobbinsC, ***HeliumAs;

	long int i, i1;
	long int j, ipLo, ipHi, ipFirstCollapsed, nelem;

	static int nCalled = 0;

#	define chLine_LENGTH 1000
	char chLine[chLine_LENGTH] , 
		/* this must be longer than chDataPath, set in path.h */
		chFilename[FILENAME_PATH_LENGTH_2];

	FILE *ioDATA;
	int lgEOL, lgHIT;

	char chSpin[2]={'1','3'};
	char chL[6]={'S','P','D','F','G','H'};

#	ifdef DEBUG_FUN
	fputs( "<+>HeCreate()\n", debug_fp );
#	endif

	/* now MALLOC the space for the helium-like lines, 
	 * but this must be done exactly one time per coreload */
	if( nCalled > 0 )
	{
#		ifdef DEBUG_FUN
		fputs( " <->HeCreate()\n", debug_fp );
#		endif
		return;
	}
	++nCalled;

	if( helike.lgCompileRecomb )
	{
		helike.lgNoRecombInterp = FALSE;
	}

	/* This is the flag that says to simulate the work of Benjamin, Skillman, and Smits (1999). */
	if( helike.lgSetBenjamin )
	{
		max_num_levels = iso.numLevels_max[ipHE_LIKE][ipHELIUM];
		max_n = iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM];
		iso.lgInd2nu_On = FALSE;
		
		if( helike.lgCompileRecomb )
		{
			fprintf( ioQQQ, " Not compiling He recomb file...can not do with Benjamin command set.");
			helike.lgCompileRecomb = FALSE;
		}
		if( helike.lgFSM )
		{
			fprintf( ioQQQ, " Can not include f-s mixing with the Benjamin command!  I have turned it off.");
			helike.lgFSM = FALSE;
		}
	}
	else
	{
		/* first find the largest number of levels in any element in the he-like sequence */
		max_num_levels = 0;
		max_n = 0;
		for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
		{
			/* only check elements that are turned on */
			if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
			{
				max_num_levels = MAX2( max_num_levels, iso.numLevels_max[ipHE_LIKE][nelem] );
				max_n = MAX2( max_n, iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] );
			}
		}
	}

	if( helike.lgRandErrGen )
	{
		/* These are only needed if error generation is turned on */
		if( (helike.Error = (float****)MALLOC(sizeof(float***)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();

		if( (helike.ErrorFactor = (float****)MALLOC(sizeof(float***)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
	}

	/* collision strength for impact by protons and singly-ionized helium */
	/* for all species and between all levels */
	if( (helike.cs_proton = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();
	if( (helike.cs_heplus = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/* power law for collision strengths (for impact by electrons, protons,
	 * and singly-ionized helium) in neighborhood of current temperature. */
	/* for all species and between all levels */
	if( (helike.cs_elec_power = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();
	if( (helike.cs_prot_power = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();
	if( (helike.cs_heplus_power = (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/* In this array are stored the C values described in Robbins 68.	*/
	if( (HeliumAs = (double***)MALLOC(sizeof(double**)*(unsigned)(LIMELM) ))==NULL )
		BadMalloc();
	
	/* branching ratios from all levels to all lower levels for all species. */
	if( (helike.BranchRatio = (double***)MALLOC(sizeof(double**)*(unsigned)(LIMELM) ))==NULL )
		BadMalloc();

	/*char ***helike chLevel[lev][lev][str]*/
	if( (chLevel = (char **)MALLOC(sizeof(char *)*(unsigned)(max_num_levels) ) )==NULL )
		BadMalloc();

	/* In this array are stored the C values described in Robbins 68.	*/
	if( (RobbinsC = (double**)MALLOC(sizeof(double*)*(unsigned)(max_num_levels) ))==NULL )
		BadMalloc();
	
	/* In this array are stored the C values described in Robbins 68.	
	 * this is a permanent one for helium only.	*/
	if( (helike.CascadeProb = (double**)MALLOC(sizeof(double*)*(unsigned)(max_num_levels) ))==NULL )
		BadMalloc();

	/* Effective recombination into all levels for all species. */
	if( (helike.RadEffec = (double*)MALLOC(sizeof(double)*(unsigned)(max_num_levels) ))==NULL )
		BadMalloc();

	if( helike.lgRandErrGen && helike.lgHugeCaseB )
	{
		/* Uncertainty in cascade probability and effective recombination,
		 * only when error generation is turned on. */
		if( (helike.SigmaCascadeProb = (double**)MALLOC(sizeof(double*)*(unsigned)(max_num_levels) ))==NULL )
			BadMalloc();
		if( (helike.SigmaRadEffec = (double*)MALLOC(sizeof(double)*(unsigned)(max_num_levels) ))==NULL )
			BadMalloc();
	}

	/* Wavenumber level energies for each element are stored here.	*/
	if( (energies = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/* The lifetime of individual levels.  helike.Lifetime[LIMELM][iso.numLevels_max[ipHE_LIKE][nelem]]	*/
	if( (helike.Lifetime = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	if( helike.lgRandErrGen && helike.lgHugeCaseB )
	{
		if( (helike.SigmaAtot = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
	}
	
	/* The sum of all A's coming out of a given n...used to make sure trend is monotonic	*/
	if( (SumAPerN = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/* the reverse of iso.quant_desig[ipHE_LIKE]...QuantumNumbers2Index[nelem][n][l][s] 
	 * gives corresponding index in iso.quant_desig[ipHE_LIKE].	*/
	if( (QuantumNumbers2Index = (long****)MALLOC(sizeof(long***)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	for( ipLo=ipHe1s1S; ipLo < iso.numLevels_max[ipHE_LIKE][ipHELIUM] ;++ipLo )
	{
		if( (chLevel[ipLo] = (char*)MALLOC(sizeof(char)*(unsigned)10 ))==NULL )
			BadMalloc();
	}

	/* >>refer He triplets	Robbins, R.R. 1968, ApJ 151, 497R	*/
	/* >>refer He triplets	Robbins, R.R. 1968a, ApJ 151, 511R	*/
	for( ipLo=ipHe1s1S; ipLo < max_num_levels ;++ipLo )
	{
		if( (RobbinsC[ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)(ipLo+1) ))==NULL )
    		BadMalloc();

		if( (helike.CascadeProb[ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)(ipLo+1) ))==NULL )
			BadMalloc();
		
		if( helike.lgRandErrGen && helike.lgHugeCaseB )
		{
			if( (helike.SigmaCascadeProb[ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)(ipLo+1) ))==NULL )
				BadMalloc();
		}
	}
	
	for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
	{
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
		{
			if( (energies[nelem] = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) ))==NULL )
				BadMalloc();
			
			if( helike.lgRandErrGen )
			{
				/* Here we add one slot for rates involving the continuum and one additional for an electric fence. */
				if( (helike.Error[nelem] = (float***)MALLOC(sizeof(float**)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+2) ) )==NULL )
					BadMalloc();

				if( (helike.ErrorFactor[nelem] = (float***)MALLOC(sizeof(float**)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+2) ) )==NULL )
					BadMalloc();

				for( i = 1; i< iso.numLevels_max[ipHE_LIKE][nelem] + 1; i++ )
				{
					if( (helike.Error[nelem][i] = (float**)MALLOC(sizeof(float*)*(unsigned)(i+1) ) )==NULL )
						BadMalloc();

					if( (helike.ErrorFactor[nelem][i] = (float**)MALLOC(sizeof(float*)*(unsigned)(i+1) ) )==NULL )
						BadMalloc();

					for( j = 0; j<i+1; j++ )
					{
						if( (helike.Error[nelem][i][j] = (float*)MALLOC(sizeof(float)*(unsigned)(3) ) )==NULL )
							BadMalloc();

						if( (helike.ErrorFactor[nelem][i][j] = (float*)MALLOC(sizeof(float)*(unsigned)(3) ) )==NULL )
							BadMalloc();

						/* set each of these to negative one */
						for( i1=0; i1<3; i1++ )
						{
							helike.Error[nelem][i][j][i1] = -1.f;
							helike.ErrorFactor[nelem][i][j][i1] = -1.f;
						}
					}
				}

				if( helike.lgHugeCaseB )
				{
					if( (helike.SigmaAtot[nelem] = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) ))==NULL )
						BadMalloc();
				}
			}
			
			if( (helike.cs_proton[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) ) )==NULL )
				BadMalloc();

			if( (helike.cs_heplus[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) ) )==NULL )
				BadMalloc();


			if( (helike.cs_elec_power[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) ) )==NULL )
				BadMalloc();

			if( (helike.cs_prot_power[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) ) )==NULL )
				BadMalloc();

			if( (helike.cs_heplus_power[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) ) )==NULL )
				BadMalloc();

			if( (HeliumAs[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) ) )==NULL )
				BadMalloc();

			if( (helike.BranchRatio[nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]+1) ) )==NULL )
				BadMalloc();

			for( i = 1; i< iso.numLevels_max[ipHE_LIKE][nelem]; i++ )
			{
				if( (helike.cs_proton[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) ) )==NULL )
					BadMalloc();

				if( (helike.cs_heplus[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) ) )==NULL )
					BadMalloc();

				if( (helike.cs_elec_power[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) ) )==NULL )
					BadMalloc();

				if( (helike.cs_prot_power[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) ) )==NULL )
					BadMalloc();

				if( (helike.cs_heplus_power[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) ) )==NULL )
					BadMalloc();

				if( (HeliumAs[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) ) )==NULL )
					BadMalloc();

				if( (helike.BranchRatio[nelem][i] = (double*)MALLOC(sizeof(double)*(unsigned)(i+1) ) )==NULL )
					BadMalloc();

				for( j = 0; j< i; j++ )
                {
					helike.cs_proton[nelem][i][j] = -FLT_MAX;
					helike.cs_heplus[nelem][i][j] = -FLT_MAX;
					helike.cs_elec_power[nelem][i][j] = -FLT_MAX;
					helike.cs_prot_power[nelem][i][j] = -FLT_MAX;
					helike.cs_heplus_power[nelem][i][j] = -FLT_MAX;
					HeliumAs[nelem][i][j] = -DBL_MAX;
					helike.BranchRatio[nelem][i][j] = -DBL_MAX;
				}
			}

			if( (helike.Lifetime[nelem] = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) ))==NULL )
				BadMalloc();
			
			/* The sum of all A's coming out of a given n...used to make sure trend is monotonic	*/
			if( (SumAPerN[nelem] = (double*)MALLOC(sizeof(double)*
				(unsigned)(iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] + 1) ) )==NULL )
				BadMalloc();
			
			/* Allocate proper number of principal quantum number.	*/
			if( (QuantumNumbers2Index[nelem] = (long***)MALLOC(sizeof(long**)*
				(unsigned)(iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] + 1) ) )==NULL )
				BadMalloc();

			for( i = 0; i <= (iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem]); ++i)
			{
				/* Allocate proper number of angular momentum quantum number.	*/
				if( (QuantumNumbers2Index[nelem][i] = (long**)MALLOC(sizeof(long*)*(unsigned)MAX2(1,i) ) )==NULL )
					BadMalloc();

				for( i1 = 0; i1 < i; ++i1 )
				{
					/* Allocate 2 spaces for spin.	*/
					if( (QuantumNumbers2Index[nelem][i][i1] = (long*)MALLOC(sizeof(long)*(unsigned)2 ) )==NULL )
						BadMalloc();
				}
			}
		}
	}

	for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
	{
		/* only check elements that are turned on */
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
		{
			/* Function he_assign assigns, in order of increasing energy, 
			 * quantum numbers n,l,s, and j to each element of vector iso.quant_desig[ipHE_LIKE]. */
			he_assign(nelem);
	    }
	}

	/* >>allocate one more than needed to allow space for end of vector sentinal
	 * in case where no collapsed levels */
	if ( (n_effs= (double *)MALLOC((unsigned)(max_num_levels+1)*sizeof(double)) ) == NULL )
		BadMalloc();
	{
		/* option to print particulars of some line when called
		 * a prettier print statement is near where chSpin is defined below*/
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			for( nelem = ipHELIUM; nelem <LIMELM; ++nelem )
			{
				/* first print internal representation */
				for( i=0; i < iso.numLevels_max[ipHE_LIKE][nelem]; ++i )
				{
					fprintf(ioQQQ,"nelem %li i %li n %li s %li l %li\n", 
						nelem, i, N_(i), S_(i), L_(i) );
				}
				if( DEBUG_LOC )
					cdEXIT(EXIT_SUCCESS);
			}
		}
	}

	/* This routine reads in transition probabilities from a file. */ 
	HelikeTransProbSetup();

	/* main helium-like arrays, fill with sane values */
	for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
	{
		/* charge to 4th power, needed for scaling laws for As*/
		double z4 = POW2((double)nelem);
		z4 *= z4;

		/* must always do helium even if turned off */
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] ) 
		{
			for( i=0; i <= ( iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] ); ++i)
			{
				SumAPerN[nelem][i] = 0.;
			}

			/* this is the number of resolved levels, so first collapsed level is [ipFirstCollapsed] */
			ipFirstCollapsed = iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; 
			
			/* Initialize some ground state stuff, easier here than in loops.	*/
			RobbinsC[ipHe1s1S][ipHe1s1S] = 1.;
			/*HeliumAs[nelem][ipHe1s1S][ipHe1s1S] = 0.;*/
			/*helike.BranchRatio[nelem][ipHe1s1S][ipHe1s1S] = 0.;*/
			helike.Lifetime[nelem][ipHe1s1S] = 0.;
			helike.CascadeProb[ipHe1s1S][ipHe1s1S] = 1.;
			if( helike.lgRandErrGen && helike.lgHugeCaseB )
			{
				helike.SigmaAtot[nelem][ipHe1s1S] = 0.;
				helike.SigmaCascadeProb[ipHe1s1S][ipHe1s1S] = 0.;
			}
			iso.stat[ipHE_LIKE][nelem][ipHe1s1S] = 1.f;
			n_effs[ipHe1s1S] = N_(ipHe1s1S) - defect(nelem,ipHe1s1S);
			energies[nelem][ipHe1s1S] = EionWN[nelem];
			iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHe1s1S] = energies[nelem][ipHe1s1S]*WAVNRYD;
			ASSERT( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHe1s1S] > 0. );
			
			/*************************************************************/
			/**  Statistical weights, Line and level energies, and		**/
			/**  transition probabilities are calculated in this loop.	**/
			/*************************************************************/
			for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem] ; ipHi++ )
			{
				/* the effective quantum number */
				n_effs[ipHi] = N_(ipHi) - defect(nelem,ipHi);
				ASSERT( ( L_(ipHi)==1 && S_(ipHi)==0 ) || 
					( N_(ipHi) - n_effs[ipHi] >= 0. ) );

				HeliumAs[nelem][ipHi][ipHi] = 0.;
				
				if (ipHi >= ipFirstCollapsed )
				{
					iso.stat[ipHE_LIKE][nelem][ipHi] = 4.f*POW2((float)N_(ipHi));
				}
				else if ( ipHi>=ipHe2p3P0 && ipHi<=ipHe2p3P2 )
				{
					iso.stat[ipHE_LIKE][nelem][ipHi] = (2.f*J_(ipHi) + 1.f);
				}
				else
				{
					iso.stat[ipHE_LIKE][nelem][ipHi] = (2.f*L_(ipHi) + 1.f) * (2.f*S_(ipHi) + 1.f);
				}

				for( ipLo=ipHe1s1S; ipLo<ipHi ; ipLo++ )
				{
					double Enerwn, Aul, Aul1;
					
					if( helike.lgHugeCaseB == FALSE )
					{
						/* atomic number or charge and stage:			  */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].nelem = (int)(nelem+1);
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].IonStg = (int)(nelem);

						/* store statistical weights in EmisLines */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo = iso.stat[ipHE_LIKE][nelem][ipLo];
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi = iso.stat[ipHE_LIKE][nelem][ipHi];
						ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi> 0.);
						ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo> 0.);
					}
					
					/************************************************/
					/**********  Level and line energies ************/
					/************************************************/
					{
						
						/* energies (in wavenumbers) that correspond to quantum defect */
						/* for atomic helium use experimental energies */
						if( nelem==ipHELIUM && ipHi<NHE1LEVELS &&
							N_(ipHi)<=iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
						{
							energies[nelem][ipHi] = EionWN[ipHELIUM] - He1Energies[ipHi];
						} 
						else if( N_(ipHi) > iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
							energies[nelem][ipHi] = RYD_INF * POW2((double)nelem/(double)N_(ipHi));  
						else
							energies[nelem][ipHi] = he_energy(n_effs[ipHi],nelem,ipHi);
						
						/* convert wavenumbers to Ryd.	*/
						EthRyd = energies[nelem][ipHi] * WAVNRYD;
						
						/* Store threshold energy in larger cloudy array	*/
						iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHi] = EthRyd;
						ASSERT( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHi] > 0. );

						Enerwn = energies[nelem][ipLo] - energies[nelem][ipHi];

						/* This is the minimum line energy we will allow.  */
						if ( Enerwn == 0 )
							Enerwn = 0.0001;

						if( Enerwn < 0. )
							Enerwn = -1.0 * Enerwn;

						/* This is necessary because of the funny nature of the 2^3P sublevels in
						 * the helium atom.  Our bookkeeping has 2^3P0 as the lowest level, and 2^3P2
						 * as the lowest, but in actuality, 2^3P0 can decay to 2^3P1.  We negate so
						 * that the line energy is positive are properly treated later.  */
						/* if( ipHi==4 && ipLo==3 && nelem==ipHELIUM )
							Enerwn = -1.0 * Enerwn;
						*/

 						if( helike.lgHugeCaseB == FALSE )
						{
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN = (float)Enerwn;
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg = (float)(Enerwn*WAVNRYD*EN1RYD);
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK = (float)(Enerwn*WAVNRYD*TE1RYD );
						}
					}
					
					/************************************************/
					/**********  Transition probabilities ***********/
					/************************************************/
					{
						if ( ipHi >= ipFirstCollapsed )
						{
							if ( ipLo >= ipFirstCollapsed )
							{
								/* Neither upper nor lower is resolved...Aul's are easy.	*/
								Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4;
								putError(nelem,ipHi,ipLo,IPRAD,0.001f);

								ASSERT( Aul > 0.);
							}
							else 
							{
								/* Lower level resolved, upper not. First calculate Aul
								 * from upper level with ang mom one higher.	*/
								Aul = he_1trans( nelem, Enerwn, n_effs[ipHi], L_(ipLo)+1,
									S_(ipLo), -1, n_effs[ipLo], L_(ipLo), S_(ipLo), ipLo-3);
								Aul *= (2.*L_(ipLo)+3.) * (2.*S_(ipLo)+1.) / 
									(4.*(double)N_(ipHi)*(double)N_(ipHi));
								if ( L_(ipLo) != 0 )
								{
									/* For all l>0, add in transitions from upper level with ang mom one lower.	*/
									Aul1 = he_1trans( nelem, Enerwn, n_effs[ipHi], L_(ipLo)-1,
										S_(ipLo), -1, n_effs[ipLo], L_(ipLo), S_(ipLo), ipLo-3);
									/* now add in this part after multiplying by stat weight for lHi = lLo-1.	*/
									Aul += Aul1*(2.*L_(ipLo)-1.) * (2.*S_(ipLo)+1.) /
										(4.*(double)N_(ipHi)*(double)N_(ipHi));
								}
								putError(nelem,ipHi,ipLo,IPRAD,0.01f);
								ASSERT( Aul > 0.);
							}
							if( N_(ipHi) > N_(ipLo) )
								SumAPerN[nelem][N_(ipHi)] += Aul;
						}
						else
						{
							/* Both levels are resolved...do the normal bit.	*/
							if( Enerwn < 0. )
							{
								Aul = he_1trans( nelem, -1.*Enerwn,  n_effs[ipLo],
									L_(ipLo), S_(ipLo), ipLo-3, n_effs[ipHi], L_(ipHi), S_(ipHi), ipHi-3);
							}
							else
							{
								Aul = he_1trans( nelem, Enerwn,  n_effs[ipHi],
									L_(ipHi), S_(ipHi), ipHi-3, n_effs[ipLo], L_(ipLo), S_(ipLo), ipLo-3);
							}

							if( N_(ipHi) > N_(ipLo) )
							{
								SumAPerN[nelem][N_(ipHi)] += Aul*(2.*L_(ipHi) + 1.)*(2.*S_(ipHi) + 1.)/
									(4. * (double)N_(ipHi) * (double)N_(ipHi)) ;
							}
						}

						if( helike.lgHugeCaseB == FALSE )
						{
							/* set the transition probability */
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul = Aul;

							ASSERT(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul > 0.);
						}

						HeliumAs[nelem][ipHi][ipLo] = Aul;
						ASSERT(HeliumAs[nelem][ipHi][ipLo] > 0.);
					}
				}
			}

			/************************************************/
			/**********  Lifetimes of levels   **************/
			/************************************************/
			for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem] ; ipHi++ )
			{
				helike.Lifetime[nelem][ipHi] = 0.;
				for( ipLo=ipHe1s1S; ipLo<ipHi ; ipLo++ )
				{	
					/* This is not actually the lifetime...must invert when all A's added.	*/
					/* No stat. weight is needed because this is the lifetime of an individual level,
					 * Collapsed or resolved is irrelevant.	*/
					if( helike.lgHugeCaseB == FALSE )
					{
						helike.Lifetime[nelem][ipHi] += EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul;
					}
					else
					{
						helike.Lifetime[nelem][ipHi] += HeliumAs[nelem][ipHi][ipLo];
					}
				}
				
				/* convert sum of As into proper lifetime in seconds */
				if( opac.lgCaseB )
				{
					/*lint -e777 float test equality */
					if( helike.Lifetime[nelem][ipHi] == HeliumAs[nelem][ipHi][ipHe1s1S] )
					{
						/* Put something huge in */
						helike.Lifetime[nelem][ipHi] = 1./iso.SmallA;
					}
					/*lint +e777 float test equality */
					else
					{
						helike.Lifetime[nelem][ipHi] = 1./(helike.Lifetime[nelem][ipHi]-HeliumAs[nelem][ipHi][ipHe1s1S]);
					}
				}
				else
				{
					helike.Lifetime[nelem][ipHi] = 1./helike.Lifetime[nelem][ipHi];
				}
			}

			/************************************************/
			/**********  Fine Structure Mixing - FSM ********/
			/************************************************/
			if( helike.lgFSM == TRUE && helike.lgHugeCaseB == FALSE )
			{
				for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem] ; ipHi++ )
				{
					for( ipLo=ipHe1s1S; ipLo<ipHi ; ipLo++ )
					{
						DoFSMixing( nelem, ipLo, ipHi );
					}
				}
			}

			/************************************************/
			/**********  stuff gf, opacity, etc.   **********/
			/************************************************/
			for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem] ; ipHi++ )
			{
				RobbinsC[ipHi][ipHi] = 1.;
				helike.CascadeProb[ipHi][ipHi] = 1.;
				if( helike.lgRandErrGen && helike.lgHugeCaseB )
				{
					helike.SigmaAtot[nelem][ipHi] = 0.;
					helike.SigmaCascadeProb[ipHi][ipHi] = 0.;
				}
			
				for( ipLo=ipHe1s1S; ipLo<ipHi ; ipLo++ )
				{
					RobbinsC[ipHi][ipLo] = 0.;
					helike.CascadeProb[ipHi][ipLo] = 0.;

					if( helike.lgHugeCaseB == FALSE )
					{
						/* there are some negative energy transitions, where the order
						 * has changed, but these are not optically allowed, these are
						 * same n, different L, forbidden transitions */
						ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN > 0. ||
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul <= iso.SmallA );

						/* very forbidden transitions may have negative energies */
						if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul <= iso.SmallA || 
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN <= 0.)
						{
							/* this branch - bogus line */
							/* transition energy in various units:*/

							/* take abs value since some have order switched */
							/* make following an air wavelength */
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng = 1e6f;

							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf = 1e-20f;

							/* derive the abs coef, call to function is gf, wl (A), g_low */
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity = 0.;
						}
						else
						{
							/* take abs value since some have order switched */
							/* make following an air wavelength */
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng = 
								(float)fabs(1.0e8/
							  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN/
						  RefIndex( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN ));
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf = 
								(float)(GetGF(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul,
							  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
							  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi));
							ASSERT(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf > 0.);

							/* derive the abs coef, call to function is gf, wl (A), g_low */
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity = 
								(float)(abscf(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf,
							  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
							  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo));
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity = 
								MAX2(0.f, EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity );
						}
						{
							/* option to print particulars of some line when called
							 * a prettier print statement is near where chSpin is defined below
							 * search for "pretty print" */
							/*@-redef@*/
							enum {DEBUG_LOC=FALSE};
							/*@+redef@*/
							if( DEBUG_LOC  )
							{
								if( (nelem == ipHELIUM) && (EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul> 1e-19) &&
									N_(ipLo) <= 5 && N_(ipHi)<=10 )
								{
									fprintf(ioQQQ,"Z %li Lo %li n %li s %li l %li \t ", 
										nelem+1, ipLo, N_(ipLo), S_(ipLo), L_(ipLo) );
									fprintf(ioQQQ," Hi %li n %li s %li l %li \t", 
										ipHi, N_(ipHi), S_(ipHi), L_(ipHi) );
									fprintf(ioQQQ,"%.4e\t%.4e\tcs\t%.4e\n",
										/*EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,*/
										EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
										EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul ,
										EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs);
								}
							}
						}

						/* redistribution function = partial or complete? */
						if( ipLo == ipHe1s1S )
						{
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].iRedisFun = ipPRD;
						}
						else
						{
							/* >>chng 01 feb 27, had been -1, crd with core only,
							 * change to crd with wings as per discussion with Ivan Hubeny */
							/*EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].iRedisFun = -1;*/
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].iRedisFun = ipCRDW;
						}
						/* create array index that will blow up */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ipCont = INT_MIN;
					}

					iso.Pop2Ion[ipHE_LIKE][nelem][ipLo] = 0.;

					if( helike.lgRandErrGen && helike.lgHugeCaseB )
					{
						ASSERT( helike.Error[nelem][ipHi][ipLo][IPRAD] >= 0. );
						/* Uncertainties in A's are added in quadrature, square root is taken below. */
						helike.SigmaAtot[nelem][ipHi] += 
							pow( helike.Error[nelem][ipHi][ipLo][IPRAD] * HeliumAs[nelem][ipHi][ipLo], 2. );
					}
				}

				helike.BranchRatio[nelem][ipHi][ipHe1s1S] = 0.;
				for( ipLo=ipHe2s3S; ipLo<ipHi ; ipLo++ )
				{
					helike.BranchRatio[nelem][ipHi][ipLo] = HeliumAs[nelem][ipHi][ipLo] * helike.Lifetime[nelem][ipHi];
					/*fprintf(ioQQQ," DEBUG He %li %li %li %e %e %e\n", 
						nelem,ipHi,ipLo,helike.BranchRatio[nelem][ipHi][ipLo] ,
						HeliumAs[nelem][ipHi][ipLo] ,
						helike.Lifetime[nelem][ipHi]); */
					/*>>chng 05 jun 18, eps > 1 needed to pass assert */
					ASSERT( helike.BranchRatio[nelem][ipHi][ipLo] <= 1.0000001 );
				}

				if( helike.lgRandErrGen && helike.lgHugeCaseB )
				{
					helike.SigmaAtot[nelem][ipHi] = sqrt( helike.SigmaAtot[nelem][ipHi] );
				}

				if( nelem == ipHELIUM )
				{
					for( ipLo=ipHe1s1S; ipLo<ipHi ; ipLo++ )
					{
						double SigmaA, SigmaCul = 0.;
						/*fprintf(ioQQQ,"DEBUG %li %li \n", ipLo , ipHi );*/
						/*ASSERT( ipLo < ipHi );*/
						for ( i=ipLo /* + 1*/; i<ipHi; i++ )
						{
							RobbinsC[ipHi][ipLo] += helike.BranchRatio[nelem][ipHi][i] * RobbinsC[i][ipLo];
							if( helike.lgRandErrGen && helike.lgHugeCaseB )
							{
								/* Uncertainties in A's and cascade probabilities */
								SigmaA = helike.Error[nelem][ipHi][i][IPRAD]*HeliumAs[nelem][ipHi][i];
								SigmaCul += 
									pow(SigmaA*RobbinsC[i][ipLo]*helike.Lifetime[nelem][ipHi], 2.) +
									pow(helike.SigmaAtot[nelem][ipHi]*helike.BranchRatio[nelem][ipHi][i]*RobbinsC[i][ipLo]*helike.Lifetime[nelem][ipHi], 2.) +
									pow(helike.SigmaCascadeProb[i][ipLo]*helike.BranchRatio[nelem][ipHi][i], 2.);
							}
							
						}
						helike.CascadeProb[ipHi][ipLo] = RobbinsC[ipHi][ipLo];
						if( helike.lgRandErrGen && helike.lgHugeCaseB )
						{
							/*ASSERT( ipLo >= ipHi || SigmaCul >= 0. );*/
							SigmaCul = sqrt(SigmaCul);
							helike.SigmaCascadeProb[ipHi][ipLo] = SigmaCul;
						}
						/*fprintf(ioQQQ,"DEBUG %li %li \n", ipLo , ipHi );*/
					}
				}
			}

			/**********************************************************************/
			/** Allowed decay conversion probabilities. See Robbins68b, Table 1. **/
			/** These are not used by the code.  They're simply calculated for   **/
			/** output, (so Ryan can make unbelieveably pretty tables.)			 **/
			/**********************************************************************/
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
					
				if( DEBUG_LOC && (nelem == ipHELIUM ) )
				{
					/* To output Bm(n,l; ipLo), set ipLo, lo_l, and lo_s accordingly.	*/
					long int hi_l,hi_s;
					double Bm;
					
					/* tripS to 2^3P	
					hi_l = 0;
					hi_s = 1;
					ipLo = ipHe2p3P0;*/

					/* tripD to 2^3P	
					hi_l = 2;
					hi_s = 1;
					ipLo = ipHe2p3P0;*/

					/* tripP to 2^3S	*/
					hi_l = 1;
					hi_s = 1;
					ipLo = ipHe2s3S;	
					
					ASSERT( hi_l != iso.quant_desig[ipHE_LIKE][nelem][ipLo].l );
					
					fprintf(ioQQQ,"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo);
					fprintf(ioQQQ,"m\t2\t\t3\t\t4\t\t5\t\t6\n");

					for( ipHi=ipHe2p3P2; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem] ; ipHi++ )
					{
						/* Pick out excitations from metastable 2tripS to ntripP.	*/
						if ( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].l == 1) && (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s == 1) )
						{
							fprintf(ioQQQ,"\n%ld\t",iso.quant_desig[ipHE_LIKE][nelem][ipHi].n);
							j = 0;
							Bm = 0;
							for ( i = ipLo; i<=ipHi; i++)
							{
								if ( (iso.quant_desig[ipHE_LIKE][nelem][i].l == hi_l) && (iso.quant_desig[ipHE_LIKE][nelem][i].s == hi_s)  )
								{
									if ( (ipLo == ipHe2p3P0) && (i > ipHe2p3P2) )
									{
										Bm += RobbinsC[ipHi][i] * ( helike.BranchRatio[nelem][i][ipHe2p3P0] + 
											helike.BranchRatio[nelem][i][ipHe2p3P1] + helike.BranchRatio[nelem][i][ipHe2p3P2] );
									}
									else
										Bm += RobbinsC[ipHi][i] * helike.BranchRatio[nelem][i][ipLo];
										
									if ( (i == ipHe2p3P0) || (i == ipHe2p3P1) || (i == ipHe2p3P2) )
									{
										j++;
										if (j == 3)
										{
											fprintf(ioQQQ,"%2.4e\t",Bm);
											Bm = 0;
										}
									}
									else
									{
										fprintf(ioQQQ,"%2.4e\t",Bm);
										Bm = 0;
									}
								}
							}
						}
					}
					fprintf(ioQQQ,"\n\n");
				}
			}

			/**************************************/
			/*********  Extra Lyman lines    ******/
			/**************************************/
			ipLo = ipHe1s1S;
			/* This ipHi does NOT represent ipHe2s1S, but the n in n^1P. */
			for( ipHi=2; ipHi <iso.nLyman[ipHE_LIKE]; ipHi++ )
			{
				double Enerwn , Aul /*,RI2,npstar */ ;
				float gHi , gLo;

				/* pointers to transitions, ipLo, ipHi
					* actual quantum numbers, nLo nHi
					* atomic number or charge and stage: */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].nelem = (int)(nelem+1);
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].IonStg = (int)(nelem);

				/* stat weight of all n ^1P is 3 */
				gHi = 3.f;

				gLo = 1.f;

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gLo = gLo;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gHi = gHi;

				/* hydrogenic energy, first ryd, convert to wn */
				Enerwn = EionWN[nelem] - RYD_INF*POW2((double)nelem)/POW2((double)ipHi);
				
				if ( nelem == ipHELIUM )
				{
					/* A simple fit for the calculation of Helium lyman Aul's.	*/
					Aul = (1.508e10) / pow((double)ipHi,2.975);
				}
				else 
				{
					/* Fit to values given in 
					 * >>refer	He-like	As	Johnson, W.R., Savukov, I.M., Safronova, U.I., & 
					 * >>refercon	Dalgarno, A., 2002, ApJS 141, 543J	*/
					/* originally astro.ph. 0201454  */
					Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)ipHi,3.1);

				}

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Aul = Aul;

				/* the energy of the transition in wavenumbers, this can be negative
					* for highly forbidden upper levels where energy switches */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN = (float)Enerwn;

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyErg = 
					(float)(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN * WAVNRYD*
					EN1RYD);

				/* the damping constant times the velocity */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].dampXvel = (float)(Aul/
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN/PI4);

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyK = 
					(float)(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN * WAVNRYD*
					TE1RYD );

				/* transition energy in various units:*/

				/* take abs value since some have order switched */
				/* make following an air wavelength */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].WLAng = 
					(float)fabs(1.0e8/
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN/
					RefIndex( iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN));

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gf = 
					(float)(GetGF(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Aul,
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN,
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gHi));

				/* derive the abs coef, call to function is gf, wl (A), g_low */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].opacity = 
					(float)(abscf(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gf,
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN,
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gLo));

					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].iRedisFun = iso.ipResoRedist[ipHE_LIKE];

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Pesc = 1.0;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Pelec_esc = 0.0;

				/* destruction probability*/
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Pdest = 0.0;

				/* upper and lower level populations */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].PopHi = 0.0;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].PopLo = 0.0;

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].xIntensity = 0.0;

				/* following three will eventually be set by model atoms, but
					* must have reasonable value at start */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].pump = 0.;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].cool = 0.;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].heat = 0.;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].ColOvTot = 0.;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].PopOpc = 0.;
				{
					/* option to print particulars of some line when called
						* a prettier print statement is near where chSpin is defined below
						* search for "pretty print" */
					/*@-redef@*/
					enum {DEBUG_LOC=FALSE};
					/*@+redef@*/
					if( DEBUG_LOC )
					{
						fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",
							nelem+1,
							ipHi,
							iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Aul , 	
							iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].opacity
							);
					}
				}
			}

			/* This is above the test below so that, in the case that one of those asserts blows,
			 * we will still have this complete set of A's for diagnostics. */
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC && (nelem == ipHELIUM) )
				{
					printCustomAs();
				}
			}

			/******************************************************/
			/***  Lifetimes should increase monotonically with  ***/
			/***  increasing n...Make sure the A's decrease.	***/
			/******************************************************/
			{
				for( i=2; i < iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] - 1; ++i)
				{
					ASSERT( SumAPerN[nelem][i] > SumAPerN[nelem][i+1] );
				}
			}
			
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC /* && (nelem == ipHELIUM) */)
				{
					for( i = 2; i<= (iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem]); ++i)
					{
						fprintf(ioQQQ,"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[nelem][i]);
					}
				}
			}

			if( helike.lgHugeCaseB == FALSE )
			{
				/* this is the two photon transition, designated by wl of 0 */
				EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].WLAng = 0.;

				/* opacity in two-photon continuum is not correct since treated as line,
				 *TODO	2	 this is part of the induce 2nu problem, which must be fixed for both
				 * he and h sequences. */
				EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].opacity /= 1e4f;
				
				/* we now have lifetimes for everything - set up true damping constants */
				for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
				{
					for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
					{
						/* the damping constant times the velocity */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].dampXvel = (float)(
							1./helike.Lifetime[nelem][ipHi]/
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN/PI4);
						/*fprintf(ioQQQ," DEBUG helike %li %li %li %.3e\n",
							nelem , ipHi , ipLo, 
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].dampXvel );*/
					}
				}
			}
		}
	}

	/* This routine reads in precompiled recombination coefficients for an
	 * array of temperatures.  Interpolation is later necessary to arrive at 
	 * coefficients for a given temperature. */
	HelikeRecombSetup();
	
/*	printstuff(printflag); */

	/* free the memory	*/
	free(n_effs );

	/* create spectroscopic designation of labels */
	for( ipLo=ipHe1s1S; ipLo < iso.numLevels_max[ipHE_LIKE][ipHELIUM]-iso.nCollapsed_max[ipHE_LIKE][ipHELIUM]; ++ipLo )
	{
		nelem = ipHELIUM;
		sprintf( chLevel[ipLo] , "%li %c%c", N_(ipLo), chSpin[S_(ipLo)], chL[MIN2(5,L_(ipLo))] );
	}

	/* option to create trace printout */
	if( (trace.lgTrace && trace.lgHeBug && !helike.lgHugeCaseB) )
	{

		/* "pretty print" only helium, O, and Fe 
		 * comments elsewhere say to search for "pretty print" */
		/* next print designation of levels */
		fprintf(ioQQQ,"These indices may not represent the same levels for each element!\n");
		fprintf(ioQQQ,"index\tdesig\tE(wn,He)\tE(wn,O)\tE(wn,Fe)\n");
		for( ipHi=1; ipHi<( MAX3(iso.numLevels_max[ipHE_LIKE][ipHELIUM],
			iso.numLevels_max[ipHE_LIKE][ipIRON],iso.numLevels_max[ipHE_LIKE][ipOXYGEN]) - iso.nCollapsed_max[ipHE_LIKE][ipOXYGEN] ); ++ipHi )
		{
			float enHe=0. , enO=0., enFe=0.;
			nelem = -1;
			if( dense.lgElmtOn[ipHELIUM] &&(ipHi < iso.numLevels_max[ipHE_LIKE][ipHELIUM]) )
			{
				enHe = EmisLines[ipHE_LIKE][ipHELIUM][ipHi][0].EnergyWN;
				nelem = ipHELIUM;
			}
			else
			{
				enHe = 0.;
			}
			if( dense.lgElmtOn[ipOXYGEN] &&(ipHi < iso.numLevels_max[ipHE_LIKE][ipOXYGEN]) )
			{
				enO = EmisLines[ipHE_LIKE][ipOXYGEN][ipHi][0].EnergyWN;
				nelem = ipOXYGEN;
			}
			else
			{
				enO = 0.;
			}
			if( dense.lgElmtOn[ipIRON] &&(ipHi < iso.numLevels_max[ipHE_LIKE][ipIRON]) )
			{
				enFe = EmisLines[ipHE_LIKE][ipIRON][ipHi][0].EnergyWN;
				nelem = ipIRON;
			}
			else
			{
				enFe = 0.;
			}
			ASSERT( nelem >0 && nelem < LIMELM );
			if( ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] )
			{
				fprintf(ioQQQ,"%li\t%li %c%c\t%.5e\t%.5e\t%.5e\n", 
					ipHi, N_(ipHi), chSpin[S_(ipHi)], chL[MIN2(5,L_(ipHi))], 
					enHe, enO, enFe  );
			}
		
		}

		fprintf(ioQQQ,"lo\thi\t He states \n");
		nelem = ipHELIUM;
		for( ipHi=ipHe2s3S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ; ++ipHi )
		{
			for( ipLo=ipHe1s1S; ipLo<ipHi ; ++ipLo )
			{
				fprintf(ioQQQ," %li\t%li\t%li %c%c\t%li %c%c\n",
					ipLo, ipHi, 
					N_(ipLo), chSpin[S_(ipLo)], chL[MIN2(5,L_(ipLo))], 
					N_(ipHi), chSpin[S_(ipHi)], chL[MIN2(5,L_(ipHi))] );
			}
		}
	}

	if( !helike.lgHugeCaseB )
	{
		if( iso.numLevels_max[ipHE_LIKE][ipHELIUM] > ipHe2p1P )
		{
			/* some wavelengths of HeI are expected to be exact */
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe2p3P0][ipHe2s3S].WLAng = 10830.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe2p3P1][ipHe2s3S].WLAng = 10830.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe2p3P2][ipHe2s3S].WLAng = 10830.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe2p1P][ipHe2s1S].WLAng = 20600.;

			if(  iso.numLevels_max[ipHE_LIKE][ipHELIUM] > ipHe3d1D )
			{
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P0].WLAng = 7065.;
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P1].WLAng = 7065.;
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P2].WLAng = 7065.;
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe3p3P][ipHe2s3S].WLAng = 3889.;
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P0].WLAng = 5876.;
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P1].WLAng = 5876.;
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P2].WLAng = 5876.;
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d1D][ipHe2p1P].WLAng = 6678.;

				if(  iso.numLevels_max[ipHE_LIKE][ipHELIUM] > ipHe4d3D )
				{
					EmisLines[ipHE_LIKE][ipHELIUM][ipHe4d3D][ipHe2p3P0].WLAng = 4471.;
					EmisLines[ipHE_LIKE][ipHELIUM][ipHe4d3D][ipHe2p3P1].WLAng = 4471.;
					EmisLines[ipHE_LIKE][ipHELIUM][ipHe4d3D][ipHe2p3P2].WLAng = 4471.;
				}
			}
		}

		{
			/* prints table of wavelengths of lines to ground */
			/*@-redef@*/
			enum {PRINTIT=FALSE};
			/*@+redef@*/
			if( PRINTIT )
			{
				for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
				{
					if( dense.lgElmtOn[nelem] )
					{
						fprintf(ioQQQ,"%li\t%s\t",nelem+1,elementnames.chElementSym[nelem] );

						ipLo = ipHe1s1S;
						for( ipHi=ipLo+1; ipHi<=ipHe2p1P; ++ipHi )
						{
							if( ipHi==ipHe2s1S || ipHi==ipHe2p3P0 ) continue;
							fprintf(ioQQQ,"%.4f\t", EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng);
						}
						fprintf(ioQQQ,"\n");
					}
				}
				cdEXIT(EXIT_SUCCESS);
			}
		}
	

		/* get the collision strength data for the He 1 lines */

		/* check on path if file not here and path set */
		/* path was parsed in getset */
		if( lgDataPathSet == TRUE )
		{
			/*path set, so look only there */
			strcpy( chFilename , chDataPath );
			strcat( chFilename , "he1_cs.dat" );
		}
		else
		{
			/* path not set, check local space only */
			strcpy( chFilename , "he1_cs.dat" );
		}

		if( trace.lgTrace )
			fprintf( ioQQQ," HeCreate opening he1_cs.dat:");

		if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
		{
			fprintf( ioQQQ, " HeCreate could not open he1_cs.dat\n" );
			if( lgDataPathSet == TRUE )
				fprintf( ioQQQ, " even tried path\n" );

			if( lgDataPathSet == TRUE )
			{
				fprintf( ioQQQ, " HeCreate could not open he1_cs.dat\n");
				fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
				fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
			}

			puts( "[Stop in HeCreate]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* check that magic number is ok */
		if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
		{
			fprintf( ioQQQ, " HeCreate could not read first line of he1_cs.dat.\n");
			puts( "[Stop in HeCreate]" );
			cdEXIT(EXIT_FAILURE);
		}
		i = 1;
		i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
		/*i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
		i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);*/

		/* the following is to check the numbers that appear at the start of he1_cs.dat */
		if( i1 !=COLLISMAGIC )
		{
			fprintf( ioQQQ, 
				" HeCreate: the version of he1_cs.dat is not the current version.\n" );
			fprintf( ioQQQ, 
				" HeCreate: I expected to find the number %i and got %li instead.\n" ,
				COLLISMAGIC, i1 );
			fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
			puts( "[Stop in HeCreate]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* get the array of temperatures */
		lgHIT = FALSE;
		while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
		{
			/* only look at lines without '#' in first col */
			if( chLine[0] != '#')
			{
				lgHIT = TRUE;
				helike.nCS = 0;
				lgEOL = FALSE;
				i = 1;
				while( !lgEOL && helike.nCS < HE1CSARRAY)
				{
					helike.CSTemp[helike.nCS] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
					++helike.nCS;
				}
				break;
			}
		}
		--helike.nCS;
		if( !lgHIT )
		{
			fprintf( ioQQQ, " HeCreate could not find line in CS temperatures in he1_cs.dat.\n");
			puts( "[Stop in HeCreate]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* create space for array of CS values, if not already done */
		if( (helike.HeCS = (float ****)MALLOC(sizeof(float ***)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();

		for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
		{
			/* TODO	2	- this structure is currently only used for helium itself...
			 * stuff numbers in for other elements, or drop the [nelem] dimension off
			 * of helike.HeCS	*/
			if( nelem != ipHELIUM )
				continue;

			/* only grab core for elements that are turned on */
			if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
			{
				if( (helike.HeCS[nelem] = (float***)MALLOC(sizeof(float**)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]- iso.nCollapsed_max[ipHE_LIKE][nelem]) ))==NULL )
					BadMalloc();

				/* avoid allocating 0 bytes, some OS return NULL pointer, PvH */
				helike.HeCS[nelem][0] = NULL;
				for( ipHi=ipHe2s3S; ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem];++ipHi )
				{
					if( (helike.HeCS[nelem][ipHi] = (float**)MALLOC(sizeof(float*)*(unsigned)ipHi ))==NULL )
						BadMalloc();

					for( ipLo=ipHe1s1S; ipLo<ipHi; ++ipLo )
					{
						if( (helike.HeCS[nelem][ipHi][ipLo] = (float*)MALLOC(sizeof(float)*(unsigned)helike.nCS ))==NULL )
							BadMalloc();

						for( i=0; i<helike.nCS; ++i )
						{
							helike.HeCS[nelem][ipHi][ipLo][i] = 0.f;
						}
					}
				}
			}
		}

		/* now read in the CS data */
		lgHIT = FALSE;
		nelem = ipHELIUM;
		while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
		{
			char *chTemp;
			/* only look at lines without '#' in first col */
			if( (chLine[0] == ' ') || (chLine[0]=='\n') )
				break;
			if( chLine[0] != '#')
			{
				lgHIT = TRUE;

				/* get lower and upper level index */
				j = 1;
				ipLo = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
				ipHi = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
				ASSERT( ipLo < ipHi );
				if( ipHi >= iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] )
					continue;
				else
				{
					chTemp = chLine;
					/* skip over 4 tabs to start of cs data */
					for( i=0; i<3; ++i )
					{
						if( (chTemp = strchr( chTemp, '\t' )) == NULL )
						{
							fprintf( ioQQQ, " HeCreate no init cs\n" );
							puts( "[Stop in HeCreate]" );
							cdEXIT(EXIT_FAILURE);
						}
						++chTemp;
					}

					for( i=0; i<helike.nCS; ++i )
					{
						float a;
						if( (chTemp = strchr( chTemp, '\t' )) == NULL )
						{
							fprintf( ioQQQ, " HeCreate not scan cs, current indices: %li %li\n", ipHi, ipLo );
							puts( "[Stop in HeCreate]" );
							cdEXIT(EXIT_FAILURE);
						}
						++chTemp;
						sscanf( chTemp , "%e" , &a );
						helike.HeCS[nelem][ipHi][ipLo][i] = a;
					}
				}
			}
		}

		/* close the data file */
		fclose( ioDATA );

		/* option to print cs data for AGN */
		/* create spectroscopic designation of labels */
		{
			/* option to print particulars of some line when called
			 * a prettier print statement is near where chSpin is defined below*/
			/*@-redef@*/
			enum {AGN=FALSE};
			/*@+redef@*/
			if( AGN )
			{
#				define NTEMP 6
				double te[NTEMP]={6000.,8000.,10000.,15000.,20000.,25000. };
				double telog[NTEMP] ,
					cs ,
					ratecoef;
				nelem = ipHELIUM;
				fprintf( ioQQQ,"trans");
				for(i=0; i<NTEMP; ++i)
				{
					telog[i] = log10( te[i] );
					fprintf( ioQQQ,"\t%.3e",te[i]);
				}
				for(i=0; i<NTEMP; ++i)
				{
					fprintf( ioQQQ,"\t%.3e",te[i]);
				}
				fprintf(ioQQQ,"\n");

				for( ipHi=ipHe2s3S; ipHi< iso.numLevels_max[ipHE_LIKE][ipHELIUM]; ++ipHi )
				{
					for( ipLo=ipHe1s1S; ipLo < ipHi; ++ipLo )
					{
						/* print the designations of the lower and upper levels */
						fprintf( ioQQQ,"%s - %s",
							chLevel[ipLo] , chLevel[ipHi] );

						/* print the interpolated collision strengths */
						for(i=0; i<NTEMP; ++i)
						{
							phycon.alogte = telog[i];
							/* print cs */
							cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
							fprintf(ioQQQ,"\t%.2e", cs );
						}

						/* print the rate coefficients */
						for(i=0; i<NTEMP; ++i)
						{
							phycon.alogte = telog[i];
							phycon.te = pow(10.,telog[i] );
							tfidle(FALSE);
							cs = HeCSInterp( nelem , ipHi , ipLo, ipELECTRON );
							/* collisional deexcitation rate */
							ratecoef = cs/sqrt(phycon.te)*COLL_CONST/iso.stat[ipHE_LIKE][nelem][ipLo] *
								sexp( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK / phycon.te );
							fprintf(ioQQQ,"\t%.2e", ratecoef );
						}
						fprintf(ioQQQ,"\n");
					}
				}
				cdEXIT(EXIT_FAILURE);
			}
		}
	}
	for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
	{
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
		{
			for( ipHi=ipHe2s3S; ipHi < iso.numLevels_max[ipHE_LIKE][nelem] ;++ipHi )
			{
				free( HeliumAs[nelem][ipHi] );
			}
			free( HeliumAs[nelem] );
		}
	}
	free( HeliumAs );

	/* Free up some memory that will not be used again. */
	for( ipLo=ipHe1s1S; ipLo < max_num_levels ;++ipLo )
	{
		free( RobbinsC[ipLo] );
	}
	for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
	{
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
		{
			free( energies[nelem] );
			free( SumAPerN[nelem] );
		}
	}
	free( energies );
	free( SumAPerN );
	
	free( RobbinsC );

	if( helike.lgHugeCaseB && helike.lgRandErrGen )
	{
		/* These are only done if we're doing a simple
		 * one zone cascade calculation.  We do error
		 * propagation as well */
		HeLikeError(ipHELIUM);
		HeRecom(ipHELIUM);
	}

#	ifdef DEBUG_FUN
	fputs( " <->HeCreate()\n", debug_fp );
#	endif
	return;
}

