/*ligbar obtain collision strength for any Li-sequence line */
#include "cddefines.h"
#include "physconst.h"
#include "phycon.h"
#include "taulines.h"
#include "ionfracs.h"
#include "showme.h"
#include "ligbar.h"

void ligbar(long int ized, 
  EmLine * t2s2p, 
  EmLine * t2s3p, 
  double *cs2s2p, 
  double *cs2s3p)
{
	double a, 
	  b, 
	  c, 
	  excit, 
	  gbar;

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

	/* compute collision strength for Li-seq g-bar approx
	 * summarized in Cochrane and McWhirter Physica Scripta 28, 25.
	 * ized is nuclear charge, can be anything larger than 2
	 * Kirk Korista
	 *
	 * ized is nuclear charge, 6 for carbon, etc
	 * t2s2p is tau array for stronger memeber of 2s 2p multiplet,
	 * which is treated as two separate lines
	 * t2s3p is next transition up, treated as an averaged multiplet
	 *
	 * cs2s2p is the cs for the single 2s2p line that comes in
	 * cs2s3p is the multiplet cs for that transition, which is
	 * treated as  a multiplet average.  If t2s3p is ever separated
	 * (as t2s2p was) then cs2s3p will be the single line not the multiplet
	 *
	 * T2S2P, T2S3P are line information array, defined in block data
	 */

	/* no need to evaluate coll strength if population is zero */
	if( xIonFracs[ t2s2p->IonStg ][ t2s2p->nelem -1 ] == 0. )
	{
		*cs2s2p = 1.;
		*cs2s3p = 1.;
		
#		ifdef DEBUG_FUN
		fputs( " <->ligbar()\n", debug_fp );
#		endif
		return;
	}

	if( ized < 3 )
	{
		/* this is a sanity check */
		fprintf( ioQQQ, " LIGBAR called with insane charge, ized=%4ld\n", 
		  ized );
		ShowMe();
		puts( "[Stop in LIGBAR]" );
		exit(1);
	}

	else if( ized == 6 )
	{
		/* CIV 1549 */
		a = 0.292;
		b = 0.289;
		c = 2.67;
	}

	else if( ized == 7 )
	{
		/* NV 1240 */
		a = 0.387;
		b = 0.247;
		c = 3.93;
	}

	else if( ized == 8 )
	{
		/* OVI 1035 -- values interpolated */
		a = 0.40;
		b = 0.256;
		c = 4.12;
	}

	else if( ized == 10 )
	{
		/* NeVIII 774 */
		a = 0.426;
		b = 0.273;
		c = 4.50;
	}

	else if( ized == 12 )
	{
		/* Mg 10 615 -- these values are general */
		a = 0.45;
		b = 0.27;
		c = 5.0;
	}

	else if( ized == 18 )
	{
		/* Ar 16 365 */
		a = 0.311;
		b = 0.294;
		c = 6.65;
	}

	else if( ized == 26 )
	{
		/* Fe 24 213 */
		a = 0.435;
		b = 0.314;
		c = 6.92;
	}

	else
	{
		/* use general formula for all other cases */
		a = 0.6 - 1.5/((float)(ized) - 2.);
		b = 0.27;
		c = 5.;
	}

	/* evaluate expression in terms of coefficients
	 * tarray(ipLnBolt) = line energy in degrees kelvin */
	excit = t2s2p->EnergyK/phycon.te;

	/* excit = e1/(te * 1.380622e-16) */
	gbar = a + b*log(1./excit+c);

	/* tarray(ipLnGF) = gf; tarray(ipLnBolt) excit temp kelvin */
	*cs2s2p = gbar*197.47*EVDEGK*t2s2p->gf/t2s2p->EnergyK;

	/* Small fudge factors to CMcW83 2s-2p fits:
	 * fits 3.57% too small compared to R-matrix calc. for Mg X.
	 * scaled all, initially, by this constant. Pradhan & Peng (1994)
	 * compilation cites a pc with Burgess, which further scales
	 * cs(C IV) by 1.0429 and cs(N V) by 0.9691, approximately. 
	 * The scaled cs(OVI) matched well with Burgess, so no further
	 * scaling was done for more highly ionized species. */

	if( ized == 6 )
	{
		*cs2s2p *= 1.08013;
	}

	else if( ized == 7 )
	{
		*cs2s2p *= 1.00370;
	}

	else
	{
		*cs2s2p *= 1.0357;
	}


	/* use general formula for 2s3p */
	a = -0.244;
	b = 0.25;
	c = 4.;

	/* excit = e2/(te * 1.380622e-16) */
	excit = t2s3p->EnergyK/phycon.te;
	gbar = a + b*log(1./excit+c);
	/* tarray(ipLnGF) = gf */
	*cs2s3p = gbar*197.47*EVDEGK*t2s3p->gf/t2s3p->EnergyK;
	/* cs2s3p = gbar * 197.47*eVdegK *  GF2/(e2/1.60184e-12)
	 * */

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