/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*CoolNitr evaluate total cooling due to nitrogen */
#include "cddefines.h"
#include "taulines.h"
#include "dense.h"
#include "phycon.h"
#include "ligbar.h"
#include "thermal.h"
#include "lines_service.h"
#include "embesq.h"
#include "atoms.h"
#include "cooling.h"
#include "nitro.h"
/* xNI_coll_stren: Author: Terry Yun 2006 */
/* This routine interpolates the collision strength of NI 
 * We are only interested in the first 5 levels, therefore there are 10 transitions.
 * Used polynomial fit, eqn: CS = a*t+b*t^1.5+c*t^0.5
 * The range of temp is 0 - 50,000 K */
/* #define PRINT */
static const int N1_SIZE = 10;

static const double aNI[N1_SIZE] = {2.755e-5,4.123e-5,7.536e-6,1.486e-5,4.516e-5,
				    -2.935e-6,4.000e-6,3.751e-6,-2.176e-6,1.024e-5};
static const double bNI[N1_SIZE] = {-8.150e-8,-1.220e-7,-2.226e-8,-4.390e-8,-1.130e-7,
				    8.000e-9,-1.1447e-8,-1.061e-8,5.610e-9,-3.227e-8};
static const double cNI[N1_SIZE] = {2.140e-4,3.272e-4,-4.944e-6,3.473e-6,-8.772e-4,
				    1.654e-3,1.675e-3,1.123e-3,3.867e-3,3.376e-4};

static double NI[6][6][3]; /* coefficients a b c into array */

/*xNI_coll_stren interpolate the collision strength of NI  */
STATIC double xNI_coll_stren(int ipLo, int ipHi)
{
	/* incorporating the N0 collision strengths give in
	 *>>refer	N1	cs	Tayal, S.S., 2006, ApJS, 163, 207 */

	double CS; /* collision strength */

	// Tayal's energy index IS NOT IN INCREASING ENERGY ORDER
	// call routine with proper energy order, swap it around
	// to Tayal order
	int ipLoTayal=-1, ipHiTayal=-1;

	/* Invalid entries returns '-1':the initial indices are smaller than the final indices */
	if(ipLo >= ipHi)
		return( -1. );

	/* Invalid returns '-1': the indices are greater than 5 or smaller than 0 */
	if(ipLo < 1 || ipHi > 5)
		return( -1. );

	int ipTayalOrder[6]={0,1,3,2,4,5};
	ipLoTayal = ipTayalOrder[ipLo];
	ipHiTayal = ipTayalOrder[ipHi];
	// special case of 2-3, must swap order again
	if( ipLo==2 && ipHi==3 )
	{
		ipLoTayal = 2;
		ipHiTayal = 3;
	}

	static bool lgMustInit = true;
	if( lgMustInit )
	{
		lgMustInit = false;
		/* assigning array initially to zero */ 
		for(long i=0; i<6; i++)
		{
			for(long j=0; j<6; j++)
			{
				set_NaN(NI[i][j][0]); /* coeff a */
				set_NaN(NI[i][j][1]); /* coeff b */
				set_NaN(NI[i][j][2]); /* coeff c */
			}
		}

		/* reading coeffs into 3D array */   
		/* The index is in physics scale */
		int index = 0;
		for(int i=1; i<6; i++)
		{
			for(int j=i+1; j<6; j++)
			{
				NI[i][j][0] = aNI[index];
				NI[i][j][1] = bNI[index];
				NI[i][j][2] = cNI[index];
				index++;
			}
		}
	}

#	ifdef PRINT
	/* physical index goes from 1 to 5 */
	for(long i=1; i<6; i++)
	{
		for(long j=i+1; j<6; j++)
		{
			printf("NI %i%i a:%e ", i, j, NI[i][j][0]);
			printf("%i%i b:%e\n", i, j, NI[i][j][1]);
			/* printf("%i%i c:%lf\n", i, j, NI[i][j][2]); */
		}
	}
#	endif

	double temp_max = 50e3;
	double temp_min = 0;
	double temp = phycon.te;
	temp = MAX2(temp, temp_min); /* return lager temp */
	temp = MIN2(temp, temp_max); /* return smaller temp */
	/* eqn: CS = a*t+b*t^1.5+c*t^0.5 */
	CS = NI[ipLoTayal][ipHiTayal][0]*phycon.te + NI[ipLoTayal][ipHiTayal][1]*phycon.te32 + 
		NI[ipLoTayal][ipHiTayal][2]*phycon.sqrte;

	return CS;
}


void CoolNitr(void)
{
	realnum p2;
	double a21, 
	  a31, 
	  a32, 
	  cs, 
	  cs2s2p, 
	  cs2s3p, 
	  cs21, 
	  cs31, 
	  cs32, 
	  p3;
	long int i;

	double a12, 
	  a13, 
	  a14, 
	  a15, 
	  a23, 
	  a24, 
	  a25, 
	  a34, 
	  a35, 
	  a45, 
	  cs12, 
	  cs13, 
	  cs14, 
	  cs15, 
	  cs23, 
	  cs24, 
	  cs25, 
	  cs34, 
	  cs35, 
	  cs45;

	double pop[5];

	DEBUG_ENTRY( "CoolNitr()" );

	static const double gN1[5]={4.,6.,4.,2.,4.};
	static const double exN1[4]={19224.464,8.713,9605.74,0.386};

	/* trans prob from 
	 * >>refer	n1	as	Froese Fischer, C. & Tachiev, G. 2004, ADNDT, 87, 1
	 */
	// 5200, g_up = 6, cs=0.337 @ 1e4K in Tayal06, 
	// NB - his J levels in ^2D are not in energy order
	cs12 = xNI_coll_stren(1, 2);
	a12 = 7.57e-6;
	double a12Tot = a12 + atoms.d5200r;

	// 5198, g_up = 4, cs=0.224 @ 1e4K in Tayal06
	// NB - his J levels in ^2D are not in energy order
	cs13 = xNI_coll_stren(1, 3);
	a13 = 2.03e-5;
	double a13Tot = a13 + atoms.d5200r;

	// 3466.543 0.055 @ 1e4K, g_up = 2
	cs14 = xNI_coll_stren(1, 4);
	a14 = 2.61e-3;

	// 3466.497 0.109 @ 1e4K, g_up = 4
	cs15 = xNI_coll_stren(1,5);
	a15 = 6.50e-3;

	/* this is fraction of 2D excitations that emit a photon in the 5200 multiplet
	 * this is used for collisional quenching in prt lines routine
	 * true line emission is in numerator while all loss terms are in
	 * denominator */
	nitro.quench_5200 = (a12+a13) / ((a12Tot+a13Tot) + (cs12 + cs13) * dense.cdsqte);

	// 0.257 @ 1e4K, g_up = 4
	cs23 = xNI_coll_stren(2, 3);
	a23 = 1.07e-8;

	// 10398.2 0.139 @ 1e4K, g_up = 2
	cs24 = xNI_coll_stren(2, 4);
	a24 = 3.45e-2;

	// 10397.7 0.366 @ 1e4K, g_up = 4
	cs25 = xNI_coll_stren(2, 5);
	a25 = 6.14e-2;

	// 10407.6 0.141 @ 1e4K, g_up = 2
	cs34 = xNI_coll_stren(3, 4);
	a34 = 5.27e-2;

	// 10407.2 0.195 @ 1e4K, g_up = 4
	cs35 = xNI_coll_stren(3, 5);
	a35 = 2.75e-2;

	// 0.123 @ 1e4K, g_up = 4
	cs45 = xNI_coll_stren(4, 5);
	a45 = 9.50e-17;

	// FUV allowed and intercombination lines that pump [NI]
	// the branching ratios used below are based on
	// >>refer	n1	as	Froese Fischer, C. & Tachiev, G. 2004, ADNDT, 87, 1
	//
	// the colission strengths are taken from
	// >>refer	n1	cs	Tayal, S. S., 2006, ApJS, 163, 207
	// the fits are based on data between 500K and 50 kK
	//
	// update rates for the driving transitions
	cs = exp(-11.3423 + 0.8379*min(phycon.alnte,10.82));
	PutCS( cs, TauLines[ipNI_pumpIndirect] );
	atom_level2( TauLines[ipNI_pumpIndirect] );

	static const double csN1_a[NI_NDP] = {
		-12.3982, -9.4523, -12.5580, -10.8813, -13.6532, -9.9035, -11.4470, -5.4776, -6.3304
	};
	static const double csN1_b[NI_NDP] = {
		0.7458, 0.3865, 0.7330, 0.6853, 0.7712, 0.3919, 0.6734, 0.1789, 0.1966
	};

  	for( i=0; i < NI_NDP; ++i )
	{
		cs = exp(csN1_a[i] + csN1_b[i]*min(phycon.alnte,10.82));
		PutCS( cs, TauLines[ipNI_pumpDirect[i]] );
		atom_level2( TauLines[ipNI_pumpDirect[i]] );
	}
	// units of pump rate s-1
	double RPI = TauLines[ipNI_pumpIndirect].Emis().pump();
	double RPD0 = TauLines[ipNI_pumpDirect[0]].Emis().pump();
	double RPD1 = TauLines[ipNI_pumpDirect[1]].Emis().pump();
	double RPD2 = TauLines[ipNI_pumpDirect[2]].Emis().pump();
	double RPD3 = TauLines[ipNI_pumpDirect[3]].Emis().pump();
	double RPD4 = TauLines[ipNI_pumpDirect[4]].Emis().pump();
	double RPD5 = TauLines[ipNI_pumpDirect[5]].Emis().pump();
	double RPD6 = TauLines[ipNI_pumpDirect[6]].Emis().pump();
	double RPD7 = TauLines[ipNI_pumpDirect[7]].Emis().pump();
	double RPD8 = TauLines[ipNI_pumpDirect[8]].Emis().pump();
	// correct for removal of photons that go straight back to ground
	// beta is the branching ratio directly back to ground
	// for pumping lines not listed below, beta is less than 0.01
	double beta = 0.7955;
	double Premove = TauLines[ipNI_pumpIndirect].Emis().Pesc() +
		TauLines[ipNI_pumpIndirect].Emis().Pelec_esc() + TauLines[ipNI_pumpIndirect].Emis().Pdest();
	RPI *= (1.-beta)/(1.-beta*(1.-Premove));
	beta = 0.1384;
	Premove = TauLines[ipNI_pumpDirect[0]].Emis().Pesc() +
		TauLines[ipNI_pumpDirect[0]].Emis().Pelec_esc() + TauLines[ipNI_pumpDirect[0]].Emis().Pdest();
	RPD0 *= (1.-beta)/(1.-beta*(1.-Premove));
	// the branching ratios to the lowest 4 excited levels are calculated in the low-density
	// limit (i.e. no collisions are treated) and exclude transitions straight back to ground
	// from the upper level of the driving line
	double pump14 = 0.0113*RPI + 0.0239*RPD0 + 0.0090*RPD1 + 0.1617*RPD2 + 0.0167*RPD3 + 0.4404*RPD4;
	double pump15 = 0.0112*RPI + 0.0265*RPD0 + 0.6253*RPD1 + 0.5108*RPD2 + 0.0824*RPD3 + 0.2588*RPD4;
	double pump13 = 0.3441*RPI + 0.8621*RPD0 + 0.0233*RPD1 + 0.0895*RPD2 + 0.1068*RPD3 + 0.1644*RPD4;
	double pump12 = 0.0417*RPI + 0.0468*RPD0 + 0.3408*RPD1 + 0.2328*RPD2 + 0.7937*RPD3 + 0.1338*RPD4;

	pump14 += 0.4881*RPD5 + 0.0876*RPD6 + 0.0450*RPD7 + 0.1777*RPD8;
	pump15 += 0.1569*RPD5 + 0.0484*RPD6 + 0.2240*RPD7 + 0.0854*RPD8;
	pump13 += 0.2908*RPD5 + 0.8397*RPD6 + 0.0694*RPD7 + 0.7369*RPD8;
	pump12 += 0.0623*RPD5 + 0.0238*RPD6 + 0.6615*RPD7 + 0.0000*RPD8;

	// estimate of pumping contribution to intensity of [NI] 5199
	nitro.pump5199 = (pump12+pump13+0.9710*pump14+0.9318*pump15) * nitro.quench_5200 * 
		dense.xIonDense[ipNITROGEN][0] * 3.82e-12;

	/**** Terry's addition, recombination from N+ **************/
	/* rate coefficient (cm3 s-1) from Table 3 of
	 *>>refer	NI	rec	Pequignot, D., Petijean, P. & Boisson, C. 1991, A&A, 251, 680 */
	double eff_recrate_2D = 1.108e-13 * pow(phycon.te*1e-4, -0.6085) /
		(1. - 0.0041 * pow(phycon.te*1e-4, -0.3975) );
	double eff_recrate_2P = 0.659e-13 * pow(phycon.te*1e-4, -0.6158);
	double fac_n1;
	if( dense.xIonDense[ipNITROGEN][0] > 0. )
		fac_n1 = dense.eden * dense.xIonDense[ipNITROGEN][1] / dense.xIonDense[ipNITROGEN][0];
	else
		fac_n1 = 0.;

	// assume levels are populated according to statistical weight
	double rec14 = eff_recrate_2P * fac_n1 * 2./6.;
	double rec15 = eff_recrate_2P * fac_n1 * 4./6.;
	double rec13 = eff_recrate_2D * fac_n1 * 4./10.;
	double rec12 = eff_recrate_2D * fac_n1 * 6./10.;

	// estimate of recombination contribution to intensity of [NI] 5199
	nitro.rec5199 = (rec12+rec13+0.9710*rec14+0.9318*rec15) * nitro.quench_5200 * 
		dense.xIonDense[ipNITROGEN][0] * 3.82e-12;

	pump14 += rec14;
	pump15 += rec15;
	pump13 += rec13;
	pump12 += rec12;

	double Cooling , CoolingDeriv;
	atom_pop5(gN1,exN1,cs12,cs13,cs14,cs15,cs23,cs24,cs25,cs34,cs35,
		cs45,a12Tot,a13Tot,a14,a15,a23,a24,a25,a34,a35,a45,pop,
		dense.xIonDense[ipNITROGEN][0],	&Cooling , &CoolingDeriv ,
		pump12 , pump13 , pump14 , pump15 );

	nitro.xN5200 = pop[1]*a12*3.818e-12;
	nitro.xN5198 = pop[2]*a13*3.821e-12;
	nitro.xN3467 = pop[3]*a14*5.729e-12;
	nitro.xN3466 = pop[4]*a15*5.729e-12;
	nitro.xN10398 = pop[3]*a24*1.910e-12;
	nitro.xN10397 = pop[4]*a25*1.910e-12;
	nitro.xN10408 = pop[3]*a34*1.908e-12;
	nitro.xN10407 = pop[4]*a35*1.908e-12;

	// population of excited NI - used in ionization balance routines
	atoms.p2nit = (realnum)(pop[1]+pop[2])/SDIV(dense.xIonDense[ipNITROGEN][0]);

	// total cooling from 5-level system
	thermal.dCooldT += CoolingDeriv;
	CoolAdd("N  1",5199, Cooling );

	/* N I 1200, cs from trans probability */
	PutCS(4.1,TauLines[ipT1200]);
	atom_level2(TauLines[ipT1200]);

	/* N II 1084, cs from trans prob */
	PutCS(5.5,TauLines[ipT1085]);
	atom_level2(TauLines[ipT1085]);

	/* [N II] coll data from 
	 * >>referold	n2	cs	Stafford, R.P., Bell, K.L, Hibbert, A. & Wijesundera, W.P.,
	 * >>rereroldcon 1994, MNRAS 268, 816,
	 * at 10,000K (v weak T dep)
	 * >>chng 00 dec 11, to Lennon & Burke
	 * >>referold	n2	cs	Lennon, D.J., & Burke, V.M., 1994, A&AS, 103, 273-277
	 * transit prob from 
	 * >>referold	n2	as	Nussbaumer, H., & Rusca, C. 1979, A&A, 72, 129
	 * >>chng 00 dec 11, to Lennon & Burke cs
	 *>>refer	n2	as	Froese Fischer, C., & Tachiev, G. 2004, At. Data Nucl. Data Tables, 87, 1
	 */
	a21 = 3.889e-3;
	a31 = 0.03140;
	a32 = 1.136;
	/* >>chng 02 may 02, put in option to switch between Lennon Burke and Stafford et al. ,
	 * default state of this var is true */
	/** \todo	1	update cs these to following reference:
	>>refer	n2	cs	Hudson, C.E. & Bell, K.L. 2004, MNRAS, 348, 1275 and A&A, 430, 725 
	 * they agree with Lennon & Burke
	 * >>chng 10 feb 24 ML: cs values updated with Hudson & Bell. Not sure why they were wrong.*/
	cs21 = 2.722;
	cs31 = 0.3162;
	cs32 = 0.806;

	/* POP3(   G1,G2,G3,   O12,  O13,  O23,  A21, A31, A32,*/
	p3 = atom_pop3(  9.,5.,1.,  cs21, cs31, cs32, a21, a31, a32 ,
	  /*   E12,E23,P2,ABUND,GAM2) */
	  21955.,24982.,&p2,dense.xIonDense[ipNITROGEN][1],0.,0.,0.);

	nitro.c5755 = p3*a32*3.46e-12;

	nitro.c6584 = p2*a21*3.03e-12;
	thermal.dCooldT += nitro.c6584*(2.2e4*thermal.tsq1 - thermal.halfte);
	CoolAdd("N  2",5755,nitro.c5755);
	CoolAdd("N  2",6584,nitro.c6584);

	/* nitro.xN2_A3_tot is fraction of excitations that produce a photon 
	 * and represents the correction for collisiona deexcitation */
	nitro.xN2_A3_tot = (a31+a32) /(a31+a32 + (cs31+cs32)/1.*dense.cdsqte );
	ASSERT( nitro.xN2_A3_tot <= 1. );

	/* N II fine structure lines, */
	/* >>chng 00 dec 11, to Lennon & Burke CS
	* >>refer	n2	cs	Hudson, C.E. & Bell, K.L. 2004, MNRAS, 348, 1275 and A&A, 430, 725
	 */
	cs21 = 0.431;
	cs32 = 1.15;
	cs31 = 0.273;

	PutCS(cs21,TauLines[ipT205]);
	PutCS(cs32,TauLines[ipT122]);
	PutCS(cs31,*TauDummy);
	atom_level3(TauLines[ipT205],TauLines[ipT122],*TauDummy);

	/* N II 2140, data 
	 * >>referold	n2	cs	Stafford, R.P., Bell, K.L, Hibbert, A. & Wijesundera, W.P. 1994, MNRAS, 268, 816
	 * >>referold	n2	cs	Lennon, D.J., & Burke, V.M., 1994, A&AS, 103, 273-277
	 * A from 
	 * >>referold	n2	as	Brage, T., Hibbert, A., Leckrone, D.S. 1997, ApJ, 478, 423
	 * >>refer	n2	as	Froese Fischer, C., & Tachiev, G. 2004, At. Data Nucl. Data Tables, 87, 1
	 * >>refer	n2	cs	Hudson, C.E. & Bell, K.L. 2004, MNRAS, 348, 1275 and A&A, 430, 725
	 */
	cs21 = 1.125;

	PutCS(cs21,TauLines[ipT2140]);
	atom_level2(TauLines[ipT2140]);

	/* N III 989.8, cs from 
	 * >>refer	N3	cs	Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425 */
	PutCS(7.12,TauLines[ipT990]);
	atom_level2(TauLines[ipT990]);

	/* 57 micron N III, A=
	 * >>refer	n3	as	Froese Fischer, C. 1983, J.Phys. B, 16, 157
	 * collision strength from 
	 * >>refer	n3	cs	Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425 */
	cs = MIN2(1.90,0.2291*phycon.te10*phycon.te10);
	PutCS(cs,TauLines[ipT57]);

	static vector< pair<TransitionList::iterator,double> > N3Pump;
	N3Pump.reserve(32);

	/* one time initialization if first call */
	if( N3Pump.empty() )
	{
		// set up level 1 pumping lines
		pair<TransitionList::iterator,double> ppa( TauLines.begin()+ipT990, 1./6. ); 
		N3Pump.push_back( ppa );
		pair<TransitionList::iterator,double> ppb( TauLines.begin()+ipT374g, 1./6. ); 
		N3Pump.push_back( ppb );
		pair<TransitionList::iterator,double> ppc( TauLines.begin()+ipT315, 1./6. ); 
		N3Pump.push_back( ppc );
		pair<TransitionList::iterator,double> ppd( TauLines.begin()+ipT324, 1./2. ); 
		N3Pump.push_back( ppd );
		pair<TransitionList::iterator,double> ppe( TauLines.begin()+ipT333, 2./3. ); 
		N3Pump.push_back( ppe );
		// set up level 2 pumping lines
		for( i=0; i < nWindLine; ++i )
		{
			/* don't test on nelem==ipIRON since lines on physics, not C, scale */
			if( (*TauLine2[i].Hi()).nelem() == 7 && (*TauLine2[i].Hi()).IonStg() == 3 )
			{
#				if	0
				DumpLine( TauLine2.begin()+i );
#				endif
				double branch_ratio;
				// the branching ratios used here ignore cascades via intermediate levels
				// usually the latter are much slower, so this should be reasonable
				if( fp_equal( (*TauLine2[i].Hi()).g(), realnum(2.) ) )
					branch_ratio = 2./3.; // 2S upper level
				else if( fp_equal( (*TauLine2[i].Hi()).g(), realnum(6.) ) )
					branch_ratio = 1./2.; // 2P upper level
				else if( fp_equal( (*TauLine2[i].Hi()).g(), realnum(10.) ) )
					branch_ratio = 1./6.; // 2D upper level
				else
					TotalInsanity();
				pair<TransitionList::iterator,double> pp2( TauLine2.begin()+i, branch_ratio ); 
				N3Pump.push_back( pp2 );
			}
		}
	}

	/* now sum pump rates */
	double pump_rate = 0.;
	vector< pair<TransitionList::iterator,double> >::const_iterator n3p;
	for( n3p=N3Pump.begin(); n3p != N3Pump.end(); ++n3p )
	{
		const TransitionList::iterator t = n3p->first;
		double branch_ratio = n3p->second;
		pump_rate += (*t).Emis().pump()*branch_ratio;
#		if	0
		dprintf( ioQQQ, "N III %.3e %.3e\n",
			 (*t).WLAng , (*t).Emis().pump()*branch_ratio );
#		endif
	}

	/* N III] N 3, N  3, 1765 multiplet */
	/*atom_level2(TauLines[ipT57]);*/
	/*AtomSeqBoron compute cooling from 5-level boron sequence model atom */
	AtomSeqBoron(TauLines[ipT57], 
	  TauLines[ipN3_1749], 
	  TauLines[ipN3_1747], 
	  TauLines[ipN3_1754], 
	  TauLines[ipN3_1752], 
	  TauLines[ipN3_1751], 
	  0.201 , 1.088 , 0.668 , 2.044 , pump_rate,"N  3");
	/*fprintf(ioQQQ," n4 %.3e\n",  (
		TauLines[ipN3_1749].xIntensity() + 
		TauLines[ipN3_1747].xIntensity() + 
		TauLines[ipN3_1754].xIntensity()+ 
		TauLines[ipN3_1752].xIntensity()+ 
		TauLines[ipN3_1751].xIntensity()) / dense.xIonDense[ipNITROGEN][2] );*/

	/* N IV 1486, N 4, N  4,collisions within 3P just guess
	 * cs to ground from 
	 * >>refer	n4	cs	Ramsbottom, C.A., Berrington, K.A., Hibbert, A., Bell, K.L. 1994,
	 * >>refercon Physica Scripta, 50, 246 */
	if( phycon.te > 1.584e4 )
	{
		cs = 21.346/(phycon.te10*phycon.te10*phycon.te10*phycon.te02);
	}
	else
	{
		cs = 75.221/(phycon.sqrte/phycon.te03/phycon.te02);
	}
	/* >>chng 01 sep 09, AtomSeqBeryllium will reset this to 1/3 so critical density correct */
	cs = MAX2(0.01,cs);
	PutCS(cs,TauLines[ipT1486]);
	AtomSeqBeryllium(.9,.9,3.0,TauLines[ipT1486],.0115);
	embesq.em1486 = (realnum)(atoms.PopLevels[3]*0.0115*1.34e-11);
	/*fprintf(ioQQQ," n4 %.3e\n", (TauLines[ipT1486].xIntensity() + embesq.em1486 ) / dense.xIonDense[ipNITROGEN][3] );*/

	/* N IV 765, cs from 
	 * >>refer	n4	cs	Ramsbottom, C.A., Berrington, K.A., Hibbert, A., Bell, K.L. 1994,
	 * >>refercon Physica Scripta, 50, 246 */
	/* >>refer	n4	as	Flemming, J., Brage, T., Bell, K.L., Vaeck, N., Hibbert, A., 
	 * >>refercon	Godefroid, M., & Froese Fischer, C., 1995, ApJ, 455, 758*/
	cs = MIN2(4.0,1.864*phycon.te03*phycon.te03);
	PutCS(cs,TauLines[ipT765]);
	atom_level2(TauLines[ipT765]);

	/* N V 1240
	 * >>refer	n5	cs	Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
	ligbar(7,TauLines[ipT1239],TauLines[ipT209],&cs2s2p,&cs2s3p);
	PutCS(cs2s2p,TauLines[ipT1239]);
	PutCS(cs2s2p*0.5,TauLines[ipT1243]);
	PutCS(1.0,*TauDummy);
	atom_level3(TauLines[ipT1243],*TauDummy,TauLines[ipT1239]);
	/*fprintf(ioQQQ," n5 %.3e\n", (TauLines[ipT1243].xIntensity() + TauLines[ipT1239].xIntensity() ) / dense.xIonDense[ipNITROGEN][4] );*/

	/* N V 209 */
	PutCS(cs2s3p,TauLines[ipT209]);
	atom_level2(TauLines[ipT209]);
	return;
}