/* 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 */
/*ConvIterCheck check whether model has converged or whether more iterations
 * are needed - implements the iter to converg comnd */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "phycon.h"
#include "mole.h"
#include "elementnames.h"
#include "dynamics.h"
#include "stopcalc.h"
#include "dense.h"
#include "iterations.h"
#include "colden.h"
#include "punch.h"
#include "rt.h"
#include "conv.h"

/*ConvIterCheck check whether model has converged or whether more iterations
 * are needed - implements the iter to converg comnd */
void ConvIterCheck( void )
{
	int lgConverged;
	long int nelem, 
		i,
		ipISO,
		ipHi, ipLo;

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

	/* =======================================================================*/
	/* this is an option to keep iterating until it converges
	 * iterate to convergence option
	 * autocv is percentage difference in optical depths allowed,
	 * =0.20 in block data
	 * checking on Ly and Balmer lines */
	/*>>chng 04 oct 19, promote loop to do all iso-electronic series */
	lgConverged = TRUE;
	strcpy( conv.chNotConverged, "Converged!" );
	if( iteration > 1 && conv.lgAutoIt )
	{
		for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
		{
			for( nelem=ipISO; nelem < LIMELM; nelem++ )
			{
				if( dense.lgElmtOn[nelem] )
				{
					/* check both H-alpha and Ly-alpha for all h-like species - 
					 * only if balmer lines thick 
					 *  so check if Ha optical deth significant */
					if( TRUE && EmisLines[ipISO][nelem][3][2].TauIn > 0.5 )
					{
						/* test if Lya converged, nLyaLevel is upper level of Lya for iso seq */
						if( fabs(EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].TauTot/
							(EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].TauIn*rt.DoubleTau)-1.) > 
							conv.autocv )
						{
							/* not converged to within AUTOCV, normally 15 percent */
							lgConverged = FALSE;

							/* for iterate to convergence, print reason why it was not converged 
							* on 3rd and higher iterations */
							sprintf( conv.chNotConverged, "%s-like Lya",elementnames.chElementSym[ipISO] );

							if( punch.lgPunConv )
							{
								fprintf( punch.ipPunConv, " %s-like Lya, nelem= %li iteration %li old %.3e new %.3e \n",
									elementnames.chElementSym[ipISO] ,
									nelem, 
									iteration,
									EmisLines[ipISO][nelem][ipH2p][ipH1s].TauTot ,
									EmisLines[ipISO][nelem][ipH2p][ipH1s].TauIn	);
							}
						}

						/* now check if major subordinate line is converged - for H-like this will
						 * be Ha, and for He-like, the 23P - 23S transition - this will not work for
						 * NISO > 2 so must check against this */
						if(ipISO==ipH_LIKE )
						{
							ipHi = 3;
							ipLo = 2;
						}
						else if( ipISO==ipHE_LIKE )
						{
							ipHi = ipHe2p3P2;
							ipLo = ipHe2s3S;
						}
						else
							/* fails when niso increased, add more transitions */
							TotalInsanity();

						if( fabs(EmisLines[ipISO][nelem][ipHi][ipLo].TauTot/
							(EmisLines[ipISO][nelem][ipHi][ipLo].TauIn*rt.DoubleTau)-1.) > conv.autocv )
						{
							/* not converged to within AUTOCV, normally 15 percent */
							lgConverged = FALSE;

							/* for iterate to convergence, print reason why it was not converged 
							* on 3rd and higher iterations */
							sprintf( conv.chNotConverged, "%s-like subord",elementnames.chElementSym[ipISO] );

							if( punch.lgPunConv )
							{
								fprintf( punch.ipPunConv, " %s-like subord, nelem= %li iteration %li old %.3e new %.3e \n" ,
									elementnames.chElementSym[ipISO],
									nelem, iteration,
									EmisLines[ipISO][nelem][ipHi][ipLo].TauTot ,
									EmisLines[ipISO][nelem][ipHi][ipLo].TauIn	);
							}
						}
					}
				}
			}
		}

		/* >>chng 03 sep 07, add this test */
		/* check on changes in major column densities */
		for( i=0; i<NCOLD; ++i )
		{
			double differ = fabs(colden.colden_old[i]-colden.colden[i]) /
				SDIV(colden.colden[i]) ;
			/* was the species column density significant relative to
			 * the total H column density, and was its abundance changing? */
			if( (colden.colden[i]/colden.colden[ipCOL_HTOT] > 1e-5) &&
				(differ > conv.autocv) )
			{
				/* not converged to within conv.autocv, normally 20 percent */
				lgConverged = FALSE;

				/* for iterate to convergence, print reason why it was not converged 
				 * on 3rd and higher iterations */
				strcpy( conv.chNotConverged, "H mole col" );

				if( punch.lgPunConv )
				{
					fprintf( punch.ipPunConv, " H mole col species %li iteration %li old %.2e new %.2e H col den %.2e\n",
						i,iteration,
						colden.colden_old[i],
						colden.colden[i],
						colden.colden[ipCOL_HTOT] );
				}
			}
		}

		/* >>chng 03 sep 07, add this test */
		/* check on changes in major column densities */
		for( i=0; i<NUM_HEAVY_MOLEC; ++i )
		{
			double differ = fabs(co.hevcol_old[i]-co.hevcol[i]) /
				SDIV(co.hevcol[i]) ;
			/* was the species abundance and changing? */
			if( (co.hevcol[i]/colden.colden[ipCOL_HTOT] > 1e-5) &&
				(differ > conv.autocv) )
			{
				/* not converged to within conv.autocv, normally 20 percent */
				lgConverged = FALSE;

				/* for iterate to convergence, print reason why it was not converged 
				 * on 3rd and higher iterations */
				strcpy( conv.chNotConverged, "CO mol col" );
				/*fprintf(ioQQQ,"debugggreset\t CO mole %li %li %.2e %.2e\n",
					i,iteration,co.hevcol_old[i],co.hevcol[i]);*/

				if( punch.lgPunConv )
				{
					fprintf( punch.ipPunConv, "CO mol col, old:%.3e new:%.3e\n" ,
						co.hevcol_old[i] ,
						co.hevcol[i] );
				}
			}
		}

		/* check on dynamical convergence in wind model with negative velocity */
		if( dynamics.lgAdvection )
		{
			double error1 = dynamics.convergence_error / SDIV(dynamics.error_scale2);
			double error2 = dynamics.discretization_error / SDIV(dynamics.error_scale2);
			/* >>chng 02 nov 29, as per Will Henney email */
			if( MAX2(error1/dynamics.convergence_tolerance, error2 ) > conv.autocv )
			/*if( MAX2(error1, error2 ) > conv.autocv )*/
			{
				lgConverged = FALSE;
				/* for iterate to convergence, print reason why it was not converged 
				 * on 3rd and higher iterations */
				strcpy( conv.chNotConverged, "Dynamics  " );
				if( punch.lgPunConv )
				{
					fprintf( punch.ipPunConv, " Dynamics\n" );
				}
			}
		}

		if( punch.lgPunConv && lgConverged )
		{
			fprintf( punch.ipPunConv, " converged\n" );
		}

		/* lower limit to number of iterations if converged */
		if( lgConverged )
			iterations.itermx = MIN2(iterations.itermx,iteration);

		/* >>chng 96 dec 20, moved following to within if on lgAutoIt
		 * this is test for stopping on first zone */
		if( phycon.te < StopCalc.tend && nzone == 1 )
		{
			lgConverged = TRUE;
			strcpy( conv.chNotConverged, "          " );
			iterations.itermx = MIN2(iterations.itermx,iteration);
		}
	}

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

}
