/* PresChng called by ConvPresTempEdenIoniz
 * evaluate the current pressure, change needed to get it to converge,
 * aplies this correction factor to all gas constituents,
 * sets conv.lgConvPres true if good pressure, false if pressure change capped
 */
#include <string.h>
#include <float.h>
#include "cddefines.h"
#include "physconst.h"
#include "abuntabllg.h"
#include "zonecnt.h"
#include "hmi.h"
#include "trace.h"
#include "wind.h"
#include "timefc.h"
#include "fluct.h"
#include "filfac.h"
#include "denpow.h"
#include "radius.h"
#include "phycon.h"
#include "h.h"
#include "he.h"
#include "ionfracs.h"
#include "densty.h"
#include "globul.h"
#include "pressure.h"
#include "colden.h"
#include "dndt.h"
#include "struc.h"
#include "hehp.h"
#include "hevmolec.h"
#include "charexc.h"
#include "elecden.h"
#include "tfidle.h"
#include "radacl.h"
#include "tabden.h"
#include "fabden.h"
#include "tababun.h"
#include "converge.h"
/* biggest change in the pressure that we will allow */
#define	PDELTA	0.07
/* allowed error in the pressure */
#define	PERROR	0.02

void PresChng( void )
{
	long int ion, 
	  nelem;

	double abun, 
	  dnew, 
	  edensave, 
	  hold, 
	  pneed, 
	  term;

	static double FacAbun, 
	  FacAbunSav, 
	  OldAbun;

	static double vold = 0.;
	double windnw;

	static double factor = 1.;
	static double pold = 0.;
	static double rsave = 0.;

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

	edensave = phycon.eden;

	/* first evaluate total pressure for this location, and current conditions
	 * pnow is just sum of gas plus local line radiation pressure */
	pressure.pnow = (float)TotalPressure();

	/* make sure this is set by one of the following branches */
	pressure.pcorec = 0.;

	/* 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 */
	if( ZoneCnt.nzone == 1 )
	{
		/* this is first zone, lock onto pressure */
		pressure.PressureInit = pressure.pnow;
		pressure.RamPresInit = (float)RamPressure();

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
			  " PresChng called, this is first zone, so reseting pressure to%10.3e\n", 
			  pressure.PressureInit );
		}
	}

	/* remember old hydrogen density */
	hold = phycon.hden;

	/* inside out globule */
	if( strcmp(pressure.chCPres,"GLOB") == 0 )
	{
		/* GLBDST is distance from globule, or glbrad-DEPTH */
		if( globul.glbdst < 0. )
		{
			fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured,  sorry.\n" );
			fprintf( ioQQQ, " This is routine PresChng, GLBDST is%10.2e\n", 
			  globul.glbdst );
			puts( "[Stop in PresChng]" );
			exit(1);
		}
		factor = (globul.glbden*pow(globul.glbrad/(globul.glbdst),globul.glbpow))/
		  phycon.hden;
		pressure.pcorec = (float)(pressure.pnow*factor);
	}

	else if( (strcmp(pressure.chCPres,"WIND") == 0) && wind.windv > 0. )
	{
		/* this is logic for wind solution
		 * following makes sure wind v only updated once per zone */
		if( radius.depth != rsave )
		{
			rsave = radius.depth;
			/* evaluate radiative acceleration */
			radacl();
			/* 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( ZoneCnt.nzone > 2 )
			{
				/* acceleration due to grad P(rad), xMassDensity is gm per cc */
				wind.AccelPres = (float)(-(pressure.RadPres - pold)/radius.drad/
				  densty.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.dReff;
			wind.acldr += (float)radius.dReff;

			/* 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 */
			dndt.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
		{
			/* this does not make sense - why zero?*/
			/*fixit(); logic appears to need pcorec be defined below - why?*/
			factor = 1.;
		}

		factor = wind.emdot/(wind.windv*phycon.hden)/radius.r1r0sq;
		pold = pressure.RadPres;
		pressure.pcorec = pressure.pnow * (float)factor;
	}

	else if( (strcmp(pressure.chCPres,"WIND") == 0) && wind.windv < 0. )
	{
		double error , slope , DenNew , RamNeeded;
		/* this is logic for d-critical flow solution, developed with Wil Henney*/

		/* following makes sure old wind v only updated once per zone,
		 * difference between old and new velocity is used to control zoning in nextdr */
		if( radius.depth != rsave )
		{
			rsave = radius.depth;
			/* evaluate radiative acceleration */
			radacl();

			/* remember the current wind velocity */
			vold = wind.windv;
		}

		/* this is the desired current ram pressure */
		RamNeeded = pressure.PressureInit + pressure.RamPresInit -
			pressure.GasPres - pressure.PresInteg;

		/* this is the current ram pressure, at this location */
		pressure.RamPres = densty.xMassDensity*
			POW2(wind.windv0*struc.DenMass[0]/densty.xMassDensity);

		/* now form the difference between current and desired pressure */
		/* this is error, the amount we need to make close to zero */
		error = RamNeeded - pressure.RamPres;

		/* this is the derivative of the pressure, 
		 * negative if ram pressure dominated, positive if gas pressure dominates */
		slope = (pressure.GasPres - pressure.RamPres)/phycon.hden;

		/* this is updated density */
		DenNew = phycon.hden + error / slope;

		/* this is the correction factor, to multipy densities of species, to make error go to zero, 
		 * it can be quite large and we want only small changes to keep things linear*/
		factor = DenNew / phycon.hden;
		factor = MAX2( 0.99 , factor );
		factor = MIN2( 1.01 , factor );

		/* remember the current total pressure */
		pressure.pnow = (float)pressure.RamPres + pressure.GasPres + pressure.PresInteg;

		/* update the current desired pressure, 
		 * the differences between pnow and pcorec will be
		 * used in ConvPresTempEdenIoniz to determine whether the 
		 * current zone has converged  */
		pressure.pcorec = pressure.PressureInit + pressure.RamPresInit;

		/* update the wind velocity from mass conservation */
		wind.windv = wind.windv0*struc.DenMass[0]/densty.xMassDensity;

		fprintf(ioQQQ,
			"nz %li now %.2e corr %.2e error %.2e slope %.2e pgas %.2e pram %.2e hden %.2e \n",
			ZoneCnt.nzone,
			pressure.pnow , pressure.pcorec , error,slope,
			pressure.GasPres , pressure.RamPres, phycon.hden );

		/* this happens if no possible solun, */
		/*if( RamNeeded <=0. || 
			fabs(1.-pressure.GasPres/pressure.RamPres) < 0.05 )*/
		if( RamNeeded <=0. )
		{
			/* this will say to stop */
			wind.lgVelPos = FALSE;
			/* do not change the density */
			factor = 1.;
			/* trick the routine above us to think we are converged */
			pressure.pnow = pressure.pcorec;
		}

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

		/*fprintf(ioQQQ,
			"vold %.2e vnew %.2e dVdr %.2e \n", vold , wind.windv ,wind.dVeldRad);*/
	}

	else if( strcmp(pressure.chCPres,"TIME") == 0 )
	{
		/* called by TIMER for N = N(T) */
		fprintf( ioQQQ, " TIME does not now work.  Sorry.\n" );
		factor = timefcCom.timefc;
		puts( "[Stop in PresChng]" );
		exit(1);
	}

	else if( strcmp(pressure.chCPres,"SINE") == 0 )
	{
		/* rapid density fluctuation */
		factor = (fluct.cfirst*cos(radius.depth*fluct.flong+fluct.flcPhase) + 
		  fluct.csecnd)/phycon.hden;
		pressure.pcorec = (float)(pressure.pnow*factor);
	}

	else if( strcmp(pressure.chCPres,"POWR") == 0 )
	{
		/* power law function of radius */
		dnew = denpow.den0*pow(radius.Radius/radius.rinner,denpow.DensityPower);
		factor = dnew/phycon.hden;
		pressure.pcorec = (float)(pressure.pnow*factor);
	}

	else if( strcmp(pressure.chCPres,"POWD") == 0 )
	{
		/* power law function of depth */
		dnew = denpow.den0*pow(1. + radius.depth/denpow.rscale,denpow.DensityPower);
		factor = dnew/phycon.hden;
		pressure.pcorec = (float)(pressure.pnow*factor);
	}

	else if( strcmp(pressure.chCPres,"POWC") == 0 )
	{
		/* power law function of column density */
		dnew = denpow.den0*pow(1. + coldenCom.colden[IPCOLUMNDENSITY-1]/
		  denpow.rscale,denpow.DensityPower);
		factor = dnew/phycon.hden;
		pressure.pcorec = (float)(pressure.pnow*factor);
	}

	else if( strcmp(pressure.chCPres,"CPRE") == 0 )
	{
		/* constant pressure */
		if( pressure.lgConPres )
		{
			/* this has pressure due to incident continuum */
			pneed = pressure.PressureInit + pressure.PresInteg;
			factor = pneed/pressure.pnow;
		}
		else
		{
			/* this does not */
			pneed = pressure.PressureInit*
				pow(radius.Radius/radius.rinner,
				pressure.PresPowerlaw);
			factor = pneed/pressure.pnow;
		}

		pressure.pcorec = (float)(pressure.pnow*factor);
		pold = pressure.pnow;
	}

	else if( strncmp( pressure.chCPres ,"DLW" , 3) == 0 )
	{
		if( strcmp(pressure.chCPres,"DLW1") == 0 )
		{
			/* call ACF sub */
			factor = fabden(radius.Radius,radius.depth)/phycon.hden;
		}
		else if( strcmp(pressure.chCPres,"DLW2") == 0 )
		{
			/* call table interpolation subroutine
			 * >>chng 96 nov 29, added tabden */
			factor = tabden(radius.Radius,radius.depth)/phycon.hden;
		}
		else
		{
			fprintf( ioQQQ, " Insanity, PresChng gets chCPres=%4.4s\n", 
			  pressure.chCPres );
			puts( "[Stop in PresChng]" );
			exit(1);
		}
		pressure.pcorec = (float)(pressure.pnow*factor);
	}

	else if( strcmp(pressure.chCPres,"CDEN") == 0 )
	{
		/* this is the default, constant density */
		factor = 1.;
		pressure.pcorec = (float)(pressure.pnow*factor);
	}

	else
	{
		fprintf( ioQQQ, " Unknown pressure law=%4.4s STOP.\n", 
		  pressure.chCPres );
		puts( "[Stop in PresChng]" );
		exit(1);
	}

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

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

	/* also check that pressure itself is ok */
	if( fabs(1.-pressure.pnow/pressure.pcorec) > PERROR )
	{
		conv.lgConvPres = FALSE;
	};

	/* save correct pressure */
	pressure.pcorec = (float)(pressure.pnow*factor);

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

	/* >>chng 97 jun 3, added variable abundances for element table command
	 * option to get abundances off a table with element read command */
	if( AbunTablLG.lgAbTaON )
	{
		for( nelem=1; nelem < LIMELM; nelem++ )
		{
			if( AbunTablLG.lgAbunTabl[nelem] )
			{
				abun = tababun(radius.Radius,radius.depth,nelem+1)*
				  phycon.hden;

				hold = abun/xIonFracs[0][nelem];
				xIonFracs[0][nelem] = (float)abun;

				for( ion=1; ion <= (nelem + 2); ion++ )
				{
					xIonFracs[ion][nelem] *= (float)hold;
				}

			}
		}
	}

	/* this happens if fluctuations abundances command entered */
	if( !fluct.lgDenFluc )
	{
		if( ZoneCnt.nzone <= 1 )
		{
			OldAbun = 1.;
			FacAbun = 1.;
			FacAbunSav = fluct.cfirst*cos(radius.depth*fluct.flong+
			  fluct.flcPhase) + fluct.csecnd;
		}

		else
		{
			OldAbun = FacAbunSav;
			/* rapid abundances fluctuation */
			FacAbunSav = fluct.cfirst*cos(radius.depth*fluct.flong+
			  fluct.flcPhase) + fluct.csecnd;
			FacAbun = FacAbunSav/OldAbun;
		}
	}

	else
	{
		FacAbun = 1.;
	}

	if( factor*FacAbun != 1. )
	{
		/* nelem,0 is abundance of nelem
		 * H, He not affected by abundance fluctuations */
		for( nelem=0; nelem < 2; nelem++ )
		{
			for( ion=0; ion <= (nelem + 2); ion++ )
			{
				/* >>chng 96 july 12 had only multiplied total abun, not ions */
				xIonFracs[ion][nelem] *= (float)factor;
			}
		}

		/* Lit on up is, so with FacAbun */
		for( nelem=2; nelem < LIMELM; nelem++ )
		{
			for( ion=0; ion <= (nelem + 2); ion++ )
			{
				/* >>chng 96 july 12 had only multiplied total abun, not ions */
				xIonFracs[ion][nelem] *= (float)(factor*FacAbun);
			}
		}

		phycon.eden *= (float)factor;
		ElecDen.EdenMet *= (float)factor;
		ElecDen.EdenNotH *= (float)factor;

		/* >>chng 96 oct 25 update some te-ne variables */
		tfidle();
		phycon.hden *= (float)factor;
		h.hi *= (float)factor;
		h.hii *= (float)factor;
		he.ahe *= (float)factor;
		he.hei *= (float)factor;
		he.heii *= (float)factor;
		he.heiii *= (float)factor;
		he.hei1 *= (float)factor;
		he.hei3 *= (float)factor;

		/* charge transfer rates */
		CharExc.CTHion[0][0] *= (float)factor;
		CharExc.CTHion[0][1] *= (float)factor;
		CharExc.CTHrec[0][0] *= (float)factor;
		CharExc.CTHrec[0][1] *= (float)factor;
		CharExc.chrener *= (float)factor;

		/* molecules done in hmole */
		hmi.htwo *= (float)factor;
		hmi.h2plus *= (float)factor;
		hmi.hminus *= (float)factor;
		hehpCom.hehp *= (float)factor;

		/* molecules done in comole */
		for( ion=0; ion < NHEVML; ion++ )
		{
			hevmolec.hevmol[ion] *= (float)FacAbun;
		}
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
		  " PresChng called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n", 
		  hold, phycon.hden, filfac.FillFac );

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


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