/* 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 */
/*iso_continuum_lower - limit max prin. quan. no. due to continuum lowering processes	*/
#include "cddefines.h"
#include "phycon.h"
#include "dense.h"
#include "iso.h"
/* #include "helike.h" */
#include "hydrogenic.h"

void iso_continuum_lower( long ipISO, long nelem )
{
	double a;
	long int np, nd, ns, nc;
	long eff_charge;
	
	/* size of rate matrices will be defined according to the n calculated here	*/
	
	ASSERT( nelem < LIMELM );
	/* this may change at a future date. */
	ASSERT( ipISO <= 1 );

	if( ipISO == ipH_LIKE )
		eff_charge = nelem + 1;
	else
		eff_charge = nelem;
	
	/* Particle packing - the density here should be density of all nuclei in the plasma */
	/* This one is just nuclear charge, which is independent of iso, and always nelem+1. */
	a = sqrt( 1.8887E8 * (nelem+1.) / pow((double)dense.xNucleiTotal, 0.333) );
	if( a > (double)iso.n_HighestResolved_max[ipISO][nelem]+(double)iso.nCollapsed_max[ipISO][nelem] )
	{
		np = iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem] + 1;
	}
	else
		np = (long)a;

	/* Debye shielding - the density here is electron density	*/
	/* This one depends on effective charge. */
	a = 2.6E7 * eff_charge * eff_charge * pow( phycon.te/dense.eden, 0.25);
	if( a > (double)iso.n_HighestResolved_max[ipISO][nelem]+(double)iso.nCollapsed_max[ipISO][nelem] )
	{
		nd = iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem] + 1;
	}
	else
		nd = (long)a;
	
	/* Stark broadening - this should be the density of singly charged ions, 
	 * both positive and negative.  The sum of protons, electrons, and HeII should be
	 * good enough.	*/
	/* This one depends on effective charge. */
	a = 3171. * pow( (double)eff_charge, 0.8 ) * pow( dense.eden + (double)dense.xIonDense[ipHYDROGEN][1]
		+ (double)dense.xIonDense[ipHELIUM][1], -0.1333);
	if( a > (double)iso.n_HighestResolved_max[ipISO][nelem]+(double)iso.nCollapsed_max[ipISO][nelem] )
	{
		ns = iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem] + 1;
	}
	else
		ns = (long)a;
	
	ASSERT( np > 3 );
	ASSERT( nd > 3 );
	ASSERT( ns > 3 );

	nc = MIN3(np, nd, ns);
	
	/* I assert greater than three because the code depends upon having at least up to n=3, and
	 * because it would take EXTREMELY dense conditions to lower the continuum that much, and
	 * the code would probably get a wrong answer then anyway.  */
	ASSERT( nc > 3 );

	/* Here is the case where the critical n is less than what has been malloc'd for resolved levels */
	/* need to truncate and remove the collapsed levels. */
	if( ipISO == ipHE_LIKE )
	{
		if( nc < iso.n_HighestResolved_max[ipHE_LIKE][nelem])
		{
			iso.n_HighestResolved_local[ipHE_LIKE][nelem] = nc;
			iso.nCollapsed_local[ipHE_LIKE][nelem] = 0;
			iso.numLevels_local[ipHE_LIKE][nelem] = nc*nc + nc + 1;
			iso.lgFullSize[ipHE_LIKE][nelem] = TRUE;
		}
		/* Here is the case where the critical n is lies among the one or more collapsed levels */
		/* we just get rid of any that are too high. */
		else if( nc <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] + iso.nCollapsed_max[ipHE_LIKE][nelem] )
		{
			iso.n_HighestResolved_local[ipHE_LIKE][nelem] = iso.n_HighestResolved_max[ipHE_LIKE][nelem];
			iso.nCollapsed_local[ipHE_LIKE][nelem] = nc - iso.n_HighestResolved_local[ipHE_LIKE][nelem];
			iso.numLevels_local[ipHE_LIKE][nelem] = iso.n_HighestResolved_local[ipHE_LIKE][nelem]*iso.n_HighestResolved_local[ipHE_LIKE][nelem] +
				iso.n_HighestResolved_local[ipHE_LIKE][nelem] + 1 + iso.nCollapsed_local[ipHE_LIKE][nelem];
			iso.lgFullSize[ipHE_LIKE][nelem] = TRUE;
		}
		/* This is usually where control will flow, because in most conditions the continuum will not be lowered.
		* Nothing changes in this case. */
		else
		{
			iso.numLevels_local[ipHE_LIKE][nelem] = iso.numLevels_max[ipHE_LIKE][nelem];
			iso.nCollapsed_local[ipHE_LIKE][nelem] = iso.nCollapsed_max[ipHE_LIKE][nelem];
			iso.n_HighestResolved_local[ipHE_LIKE][nelem] = iso.n_HighestResolved_max[ipHE_LIKE][nelem];
			iso.lgFullSize[ipHE_LIKE][nelem] = FALSE;
		}

	}
	else
	{
		if( nc <= iso.numLevels_max[ipISO][nelem])
		{
			/* This flag tells us later if we have to reevaluate recombination 
			 * only a problem if eden has changed since last iteration, but temperature
			 * has not, and the eden change caused a change in the nc calculated here. */
			if( nc != iso.numLevels_local[ipISO][nelem] )
				hydro.lgReevalRecom = TRUE;
			else
				hydro.lgReevalRecom = FALSE;
		
			iso.numLevels_local[ipISO][nelem] = nc;
			iso.nCollapsed_local[ipISO][nelem] = 0;
			iso.lgFullSize[ipISO][nelem] = TRUE;
		}
		/* This is usually where control will flow, because in most conditions the continuum will not be lowered.
		* Nothing changes in this case. */
		else
		{
			hydro.lgReevalRecom = FALSE;
			iso.n_HighestResolved_local[ipISO][nelem] = iso.n_HighestResolved_max[ipISO][nelem];
			iso.nCollapsed_local[ipISO][nelem] = 0;
			iso.numLevels_local[ipISO][nelem] = iso.numLevels_max[ipISO][nelem];
			iso.lgFullSize[ipISO][nelem] = FALSE;
		}
	}

	/* None of these can be greater than that which was originally malloc'd. */
	ASSERT( iso.numLevels_local[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] );
	ASSERT( iso.nCollapsed_local[ipISO][nelem] <= iso.nCollapsed_max[ipISO][nelem] );
	ASSERT( iso.n_HighestResolved_local[ipISO][nelem] <= iso.n_HighestResolved_max[ipISO][nelem] );

	/* Lyman lines can not be greater than original malloc or critical pqn. */
	iso.nLyman[ipISO] = MIN2( nc, iso.nLyman[ipISO]);

	return;
}
