/* 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 */
/*PressureChange called by ConvPresTempEdenIoniz
 * evaluate the current pressure, change needed to get it to converge,
 * the global static variable pressure_change_factor
 * aplies this correction factor to all gas constituents,
 * sets conv.lgConvPres true if good pressure, false if pressure change capped */
/*lgConvPres finds amount by which density should change to move towards pressure equilibrium
 * returns true if pressure is converged */
#include "cddefines.h"
#include "physconst.h"
#include "abund.h"
#include "hmi.h"
#include "trace.h"
#include "wind.h"
#include "thermal.h"
#include "dense.h"
#include "geometry.h"
#include "radius.h"
#include "mole.h"
#include "dynamics.h"
#include "pressure.h"
#include "colden.h"
#include "rt.h"
#include "tabden.h"
#include "fabden.h"
#include "conv.h"

/* this is the pressure change pressure_change_factor, needed to move current pressure to correct pressure */
static double pressure_change_factor;

/*PressureChange evaluate the current pressure, and change needed to
 * get it to PresTotlInit,
 * return value is true is density was changed, false if no changes were necessary */
int PressureChange( 
	/* this is change factor, 1 at first, becomes smaller as oscillations occur */
	double dP_chng_factor )
{
	long int ion, 
	  nelem,
	  lgChange,
	  mol;

	double abun, 
	  edensave, 
	  hold;

	/* biggest multiplicative change in the pressure that we will allow */
	/* allowed error in the pressure is conv.PressureErrorAllowed*/
	double pdelta;

	static double FacAbun, 
	  FacAbunSav, 
	  OldAbun;

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

	edensave = dense.eden;

	/* first evaluate total pressure for this location, and current conditions
	 * CurrentPressure is just sum of gas plus local line radiation pressure */
	/* this sets values of pressure.PresTotlCurr, also calls tfidle */
	PresTotCurrent();

	/* this is special case where we are working on first zone, at
	 * illuminated face of the cloud.  Here we want to remember the
	 * initial pressure in case constant pressure is needed */
	/* >>chng 05 jan 10, chng from nzone==1 to nzone<=1 so pressure not changed
	 * during search phase */
	if( nzone <= 1 )
	{
		/* this is first zone, lock onto pressure */
		pressure.PresTotlInit = pressure.PresTotlCurr;
		pressure.PresRamInit = pressure.PresRamCurr;
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
			  " PressureChange called, this is first zone, so reseting pressure to %e\n", 
			  pressure.PresTotlInit );
		}
	}

	/* this is equal to zero upon first call in this iteration,
	 * init some vars if this is first call */
	if( !conv.nTotalIoniz )
	{
		/* one time initialization of variables */
		conv.hist_pres_npres = -1;
		conv.hist_pres_nzone = -1;
		conv.hist_pres_limit = 0;
	}

	/* this will save the  history of density - pressure relationship
	 * for the current zone */
	if( nzone!=conv.hist_pres_nzone )
	{
		/* first time in this zone - reset counters */
		conv.hist_pres_nzone = nzone;
		/* this counter will be updated after vars are saved so will be
		 * total number of saved values */
		conv.hist_pres_npres = 0;
	}
	/* do we need to allocate, or reallocate, memory for the vectors */
	if( conv.hist_pres_limit==0 )
	{
		conv.hist_pres_limit = 200;
		conv.hist_pres_density = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) );
		conv.hist_pres_current = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) );
		conv.hist_pres_correct = (double *)MALLOC( (unsigned)(conv.hist_pres_limit)*sizeof(double) );
	}
	/* ran out of space - allocate more */
	if( (conv.hist_pres_npres+1) >= conv.hist_pres_limit )
	{
		conv.hist_pres_limit *= 3;
		if( (conv.hist_pres_density = (double *)realloc( conv.hist_pres_density, (unsigned)(conv.hist_pres_limit)*sizeof(double)))==NULL )
			BadMalloc();
		if( (conv.hist_pres_current = (double *)realloc( conv.hist_pres_current, (unsigned)(conv.hist_pres_limit)*sizeof(double)))==NULL )
			BadMalloc();
		if( (conv.hist_pres_correct = (double *)realloc( conv.hist_pres_correct, (unsigned)(conv.hist_pres_limit)*sizeof(double)))==NULL )
			BadMalloc();
	}

	/* >>chng 04 feb 11, add option to remember current density and pressure */
	conv.hist_pres_density[conv.hist_pres_npres] = dense.gas_phase[ipHYDROGEN];
	conv.hist_pres_current[conv.hist_pres_npres] = pressure.PresTotlCurr;
	conv.hist_pres_correct[conv.hist_pres_npres] = pressure.PresTotlCorrect;
	++conv.hist_pres_npres;

	/* remember old hydrogen density */
	hold = dense.gas_phase[ipHYDROGEN];

	/* this will be set true if density or abundances change in this zone */
	lgChange = FALSE;

	/* this evaluates current pressure, sets pressure_change_factor and 
	 * pressure.PresTotlCorrect, updates velocity,
	 * and returns TRUE if pressure has converged */
	conv.lgConvPres = lgConvPres();

	/* if convergence is OK at present state, so no change reqd, simply return 
	 * >>chng 05 feb 04, cannot do this test here since variable abundances, test has
	 * not been done, so variable abundances did not work,
	 * caught by Marcelo Castellanos and David Valls-Gabaud 
	if( conv.lgConvPres )
		return FALSE;*/

	/* >> chng 02 dec 13 rjrw: short-circuit if nothing changes */
	if( pressure_change_factor != 1. )
	{
		lgChange = TRUE;
	}

	/* allow 3 percent changes,
	 * dP_chng_factor is initially 1, becomes smaller if sign of pressure change, changes */
	pdelta = 0.03 * dP_chng_factor;

	/* make sure that change is not too extreme */
	pressure_change_factor = MIN2(pressure_change_factor,1.+pdelta);
	pressure_change_factor = MAX2(pressure_change_factor,1.-pdelta);

	{
		/*@-redef@*/
		enum{DEBUG_LOC=FALSE};
		static long int nsave=-1;
		/*@+redef@*/
		if( DEBUG_LOC /*&& nzone > 150 && iteration > 1*/ )
		{
			if( nsave-nzone ) fprintf(ioQQQ,"\n");
			nsave = nzone;
			fprintf(ioQQQ,"nnzzone\t%li\t%.2f%%\t%.3f\t%.2e\t%.2e\t%.2e\t%.2e\n", 
				nzone,
				pressure_change_factor,
				/* when this is negative we need to raise the density */
				(pressure.PresTotlCorrect-pressure.PresTotlCurr)*100./pressure.PresTotlCorrect, 
				pressure.PresTotlCorrect,
			    pressure.PresTotlCurr, 
				pressure.PresGasCurr,
				pressure.PresRadCurr ) ;
		}
	}

	/* >>chng 97 jun 03, added variable abundances for element table command
	 * option to get abundances off a table with element read command */
	if( abund.lgAbTaON )
	{
		lgChange = TRUE;
		for( nelem=1; nelem < LIMELM; nelem++ )
		{
			if( abund.lgAbunTabl[nelem] )
			{
				abun = AbundancesTable(radius.Radius,radius.depth,nelem+1)*
				  dense.gas_phase[ipHYDROGEN];

				hold = abun/dense.gas_phase[nelem];
				dense.gas_phase[nelem] = (float)abun;

				for( ion=0; ion < (nelem + 2); ion++ )
				{
					dense.xIonDense[nelem][ion] *= (float)hold;
				}
			}
		}
	}

	/* this is set false if fluctuations abundances command entered,
	 * when density variations are in place this is true,
	 * and dense.chDenseLaw is "SINE" */
	if( !dense.lgDenFlucOn )
	{
		/* abundance variations are in place */
		lgChange = TRUE;
		if( nzone <= 1 )
		{
			OldAbun = 1.;
			FacAbun = 1.;
			if( dense.lgDenFlucRadius )
			{
				/* cycle over radius */
				FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
				dense.flcPhase) + dense.csecnd;
			}
			else
			{
				/* cycle over column density */
				FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
				dense.flcPhase) + dense.csecnd;
			}
		}
		else
		{
			OldAbun = FacAbunSav;
			/* rapid abundances fluctuation */
			if( dense.lgDenFlucRadius )
			{
				/* cycle over radius */
				FacAbunSav = dense.cfirst*cos(radius.depth*dense.flong+
				dense.flcPhase) + dense.csecnd;
			}
			else
			{
				/* cycle over column density */
				FacAbunSav = dense.cfirst*cos( colden.colden[ipCOL_HTOT]*dense.flong+
				dense.flcPhase) + dense.csecnd;
			}
			FacAbun = FacAbunSav/OldAbun;
		}
	}
	else
	{
		/* abundance variations are NOT in place */
		FacAbun = 1.;
	}

	/* chng 02 dec 11 rjrw -- remove test, saves marginal time could generate nasty intermittent bug when 
	 * ( pressure_change_factor*FacAbun == 1 ) && (pressure_change_factor != 1) */
	if( lgChange )
	{
		/* H, He not affected by abundance fluctuations, so only change these
		 * by the pressure change factor */
		for( nelem=ipHYDROGEN; nelem <= ipHELIUM; ++nelem )
		{
			dense.xMolecules[nelem] *= (float)pressure_change_factor;
			dense.gas_phase[nelem] *= (float)pressure_change_factor;
			for( ion=0; ion < (nelem + 2); ion++ )
			{
				/* >>chng 96 jul 12 had only multiplied total abun, not ions */
				dense.xIonDense[nelem][ion] *= (float)pressure_change_factor;
			}
		}

		/* Li on up is affect by both pressure and abundance variations,
		 * so multiply by both factors */
		for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
		{
			dense.xMolecules[nelem] *= (float)(pressure_change_factor*FacAbun);
			dense.gas_phase[nelem] *= (float)(pressure_change_factor*FacAbun);
			for( ion=0; ion < (nelem + 2); ion++ )
			{
				/* >>chng 96 jul 12 had only multiplied total abun, not ions */
				dense.xIonDense[nelem][ion] *= (float)(pressure_change_factor*FacAbun);
			}
		}

		dense.eden *= pressure_change_factor;

		/* >>chng 96 oct 25 update some te-ne variables */
		tfidle(FALSE);

		/* molecules done in hmole, only change pressure, not abundances */
		for (mol = 0; mol < N_H_MOLEC; mol++) 
		{
			hmi.Hmolec[mol] *= (float)pressure_change_factor;
		}
		hmi.H2_total *= (float)pressure_change_factor;

		/* molecules done in comole */
		/* >>chng 03 sep 15, upper limit had not included the C, O atoms/ions */
		/*for( ion=0; ion < NUM_HEAVY_MOLEC; ion++ )*/
		for( ion=0; ion < NUM_COMOLE_CALC; ion++ )
		{
			co.hevmol[ion] *= (float)(pressure_change_factor*FacAbun);
		}
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
		  " PressureChange called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n", 
		  hold, dense.gas_phase[ipHYDROGEN], geometry.FillFac );

		if( trace.lgNeBug )
		{
			fprintf( ioQQQ, " EDEN change PressureChange from to %10.3e %10.3e %10.3e\n", 
			  edensave, dense.eden, edensave/dense.eden );
		}
	}

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (nzone>215)/**/ )
		{
			fprintf( ioQQQ, 
				"%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%c\n", 
			  radius.depth, 
			  pressure.PresTotlCurr, 
			  pressure.PresTotlInit + pressure.PresInteg, 
			  pressure.PresTotlInit, 
			  pressure.PresGasCurr, 
			  pressure.PresRamCurr, 
			  pressure.PresRadCurr, 
			  /* subtract continuum rad pres which has already been added on */
			  pressure.PresInteg - pressure.pinzon, 
			  wind.windv/1e5,
			  sqrt(5.*pressure.PresGasCurr/3./dense.xMassDensity)/1e5,
			  TorF(conv.lgConvPres) );
		}
	}

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

	return lgChange;
}

/*lgConvPres finds amount by which density should change to move towards pressure equilibrium
 * returns true if pressure is converged */
int lgConvPres(void)
{

	static double pold = 0.;

	double windnw,
	  dnew, 
	  term;
	int lgRet;

	/* make sure this is set by one of the following branches - set to zero here, 
	 * then assert that it is greater than zero at end */
	pressure.PresTotlCorrect = 0.;

	/* evaluate a series of possible pressure options, and set the file static variable
	 * pressure_change_factor */
	/* inside out globule */
	if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
	{
		/* GLBDST is distance from globule, or glbrad-DEPTH */
		if( radius.glbdst < 0. )
		{
			fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured,  sorry.\n" );
			fprintf( ioQQQ, " This is routine PressureChange, GLBDST is%10.2e\n", 
			  radius.glbdst );
			puts( "[Stop in PressureChange]" );
			cdEXIT(EXIT_FAILURE);
		}
		pressure_change_factor = (radius.glbden*pow(radius.glbrad/(radius.glbdst),radius.glbpow))/
		  dense.gas_phase[ipHYDROGEN];
		pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
	}

	/* this is positive wind velocity the outflowing wind beyond sonic point */
	else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv > 0. )
	{

		static double rsave = 0.;
		/* this is logic for "old" wind solution,
		 * which assumes positive velocity, well above sonic point.
		 * following makes sure wind v only updated once per zone */
		/*lint -e777 float test equality */
		if( radius.depth != rsave )
		/*lint +e777 float test equality */
		{
			rsave = radius.depth;
			/* evaluate radiative acceleration */
			RT_radiative_acceleration();
			/* Wind model
			 * G, COMASS = mass of star in solar masses */
			wind.agrav = (float)((6.67e-8*wind.comass)*(SOLAR_MASS/radius.Radius)/
			  radius.Radius);

			/* is this the first zone? */
			if( nzone > 2 )
			{
				/* acceleration due to grad P(rad), xMassDensity is gm per cc */
				wind.AccelPres = (float)(-(pressure.PresRadCurr - pold)/radius.drad/
				  dense.xMassDensity);
			}
			else
			{
				wind.AccelPres = 0.;
			}

			/* this is numerically unstable */
			wind.AccelPres = 0.;

			/* AccelXXX is computed in radinc, is continuum and line acceleration
			 * AccelCont is continuum accel 
			 * AccelLine is line accel 
			 * AccelPres is gradient of local gas pressure */
			wind.AccelTot = wind.AccelCont + wind.AccelLine + wind.AccelPres;

			/* remember largest radiative acceleration */
			wind.AccelMax = (float)MAX2(wind.AccelMax,wind.AccelTot);

			/* keep track of average acceleration */
			wind.AccelAver += wind.AccelTot*(float)radius.drad_x_fillfac;
			wind.acldr += (float)radius.drad_x_fillfac;

			/* following is form of energy equation, only used to check if neg vel */
			term = POW2(wind.windv) + 2.*(wind.AccelTot - wind.agrav)*
			  radius.drad;

			/* remember the current wind velocity
			vold = wind.windv; */
			if( term <= 1e3 )
			{
				/* wind vel is way below sonic point, give up, do not change vel */
				wind.lgVelPos = FALSE;
			}
			else
			{
				/* square of new wind velocity for outer edge of this zone */
				windnw = (double)(wind.windv*wind.windv) + (double)(2.*
				  (wind.AccelTot-wind.agrav))*radius.drad;

				/* WINDV is velocity at OUTER edge of this zone */
				wind.windv = (float)sqrt(windnw);
				wind.lgVelPos = TRUE;
			}

			/* following needed for expansion cooling */
			dynamics.dDensityDT = (float)(wind.AccelTot/wind.windv + 2.*wind.windv/
			  radius.Radius);

			/* following used to control new DRAD
			wind.dVeldRad = (float)(fabs(wind.windv-vold)/wind.windv/radius.drad); */

		}
		else
		{
			pressure_change_factor = 1.;
		}

		pressure_change_factor = wind.emdot/(wind.windv*dense.gas_phase[ipHYDROGEN])/radius.r1r0sq;
		pold = pressure.PresRadCurr;
		pressure.PresTotlCorrect = pressure.PresTotlCurr * pressure_change_factor;
	}

	/* this is negative wind velocity the new dynamics */
	else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. )
	{
		/* sets pressure.PresTotlCorrect  */
		pressure_change_factor = DynaPresChngFactor();
	}

	/* this is impossible - wind with identically zero velocity */
	else if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv == 0. )
	{
		/* declare insanity  */
		fprintf( ioQQQ, " PROBLEM WIND called with zero velocity - this is impossible.\n Sorry.\n" );
		/* TotalInsanity announces fatal problem, ShowMe, then cdEXIT with failure */
		TotalInsanity();
	}

	else if( strcmp(dense.chDenseLaw,"SINE") == 0 )
	{
		/* rapid density fluctuation */
		if( dense.lgDenFlucRadius )
		{
			pressure_change_factor = (dense.cfirst*cos(radius.depth*dense.flong+dense.flcPhase) + 
			dense.csecnd)/dense.gas_phase[ipHYDROGEN];
		}
		else
		{
			pressure_change_factor = (dense.cfirst*cos(colden.colden[ipCOL_HTOT]*dense.flong+dense.flcPhase) + 
			dense.csecnd)/dense.gas_phase[ipHYDROGEN];
		}
		pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
	}

	else if( strcmp(dense.chDenseLaw,"POWR") == 0 )
	{
		/* power law function of radius */
		dnew = dense.den0*pow(radius.Radius/radius.rinner,(double)dense.DensityPower);
		pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
		pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
	}

	else if( strcmp(dense.chDenseLaw,"POWD") == 0 )
	{
		/* power law function of depth */
		dnew = dense.den0*pow(1. + radius.depth/dense.rscale,(double)dense.DensityPower);
		pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
		pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
	}

	else if( strcmp(dense.chDenseLaw,"POWC") == 0 )
	{
		/* power law function of column density */
		dnew = dense.den0*pow(1.f + colden.colden[ipCOL_HTOT]/
		  dense.rscale,dense.DensityPower);
		pressure_change_factor = dnew/dense.gas_phase[ipHYDROGEN];
		pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
	}

	else if( strcmp(dense.chDenseLaw,"CPRE") == 0 )
	{
		/* constant pressure */
		if( pressure.lgContRadPresOn )
		{
			/* >>chng 01 oct 31, replace pneed with CorretPressure */
			/* this has pressure due to incident continuum */
			pressure.PresTotlCorrect = pressure.PresTotlInit + pressure.PresInteg;
		}
		else
		{
			/* this does not have pressure due to incident continuum*/
			pressure.PresTotlCorrect = pressure.PresTotlInit*
				/* following term normally unity, power law set with option par on cmmnd*/
				pow(radius.Radius/radius.rinner,(double)pressure.PresPowerlaw);
		}

		/* ratio of correct to current pressures */
		pressure_change_factor = pressure.PresTotlCorrect / pressure.PresTotlCurr;
	}

	else if( strncmp( dense.chDenseLaw ,"DLW" , 3) == 0 )
	{
		if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
		{
			/* call ACF sub */
			pressure_change_factor = fabden(radius.Radius,radius.depth)/dense.gas_phase[ipHYDROGEN];
		}
		else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
		{
			/* call table interpolation subroutine
			 * >>chng 96 nov 29, added tabden */
			pressure_change_factor = tabden(radius.Radius,radius.depth)/dense.gas_phase[ipHYDROGEN];
		}
		else
		{
			fprintf( ioQQQ, " Insanity, PressureChange gets chCPres=%4.4s\n", 
			  dense.chDenseLaw );
			puts( "[Stop in PressureChange]" );
			cdEXIT(EXIT_FAILURE);
		}
		pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
		if( pressure.PresTotlCorrect > 3. || pressure.PresTotlCorrect< 1./3 )
		{
			static int lgWARN2BIG=FALSE;
			if( !lgWARN2BIG )
			{
				lgWARN2BIG = TRUE;
				fprintf(ioQQQ,"\n\n >========== Warning!  The tabulated or functional change in density as a function of depth was VERY large. This is zone %li.\n",nzone);
				fprintf(ioQQQ," >========== Warning!  This will cause convergence problems. \n");
				fprintf(ioQQQ," >========== Warning!  The current radius is %.3e. \n",radius.Radius);
				fprintf(ioQQQ," >========== Warning!  The current depth is %.3e. \n",radius.depth);
				fprintf(ioQQQ," >========== Warning!  Consider using a more modest change in density vs radius. \n\n\n");
			}
		}
	}

	else if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
	{
		/* this is the default, constant density */
		pressure_change_factor = 1.;
		pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
	}

	else
	{
		fprintf( ioQQQ, " Unknown pressure law=%s= This is a major internal error.\n", 
		  dense.chDenseLaw );
		puts( "[Stop in PressureChange]" );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	/* one of the branches above must have reset this variable,
	 * and it was init to 0 at start.  Confirm that non-zero */
	ASSERT( pressure.PresTotlCorrect > FLT_MIN );

	/* now see whether current pressure is within error bounds */
	if( pressure_change_factor > 1. + conv.PressureErrorAllowed || pressure_change_factor < 1. - conv.PressureErrorAllowed )
	{
		lgRet = FALSE;
		conv.lgConvPres = FALSE;
	}
	else
	{
		lgRet = TRUE;
		conv.lgConvPres = TRUE;
	}

	return lgRet;
}
