/*WidthLine compute line width (cm/sec), using optical depth array information */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "physconst.h"
#include "doppvel.h"
#include "redis.h"
#include "phycon.h"
#include "pressure.h"
#include "touton.h"
#include "wind.h"
#include "widthline.h"

double WidthLine(EmLine * t)
{
	double WidthLine_v, 
	  a, 
	  aa, 
	  atau, 
	  b, 
	  r, 
	  tau, 
	  tauin, 
	  tauout, 
	  therm, 
	  vth;

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

	/* lgTauOutOn=.TRUE. if optical depth in outer direction is defined
	 * smaller of inner and outer optical depths is chosen for esc prob
	 * thermal broadening alone follows, lvg at 99
	 * following uses line width from Bonilha et al. Ap.J. (1979) 233 649
	 * return value is half width*(1-ESC PROB)
	 * TAU is used for width of source function
	 * T is used for escape probability
	 * this assumes incompl redis, a.tau^1/3 width */
	tauin = t->TauIn;
	tauout = t->TauTot;
	a = t->damp;
	vth = DoppVel.doppler[ t->nelem -1 ];

	if( touton.lgTauOutOn )
	{
		/* find optical depth to nearest boundary */
		/*tau = tauout/2.;*/
		/* >>chng 99 Apr 22, go back to smallest optical depth */
		tau = MIN2( tauin , tauout );
	}
	else
	{
		tau = tauin;
	}

	/* max optical depth is thermalization length */
	therm = 5.3e16/phycon.eden;

	if( tau > therm )
	{
		pressure.lgPradDen = TRUE;
		tau = therm;
	}

	/* !Wind means static */
	if( strcmp(redis.chRedis,"WIND") != 0 )
	{
		/* static geometry */
		if( tau< 0.0001 )
		{
			/* >>chng 99 Apr 22, this branch added for very optically thin lines */
			aa = 14.;
			b = 9.21;
			/* in very optically thin limt assume 1-escp is tau */
			WidthLine_v = vth*0.8862*aa/b*tau ;
		}
		else if ( tau <= 20. )
		{
			atau = log(MAX2(0.0001,tau));
			aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau;
			b = 6.5*tau - atau;
			WidthLine_v = vth*0.8862*aa/b*(1. - t->Pesc);
		}
		else
		{
			atau = log(MAX2(0.0001,tau));
			aa = 1. + 2.*atau/pow(1. + 0.3*a*tau,0.6667) + pow(6.5*
			  a*tau,0.333);
			b = 1.6 + 1.5/(1. + 0.20*a*tau);
			WidthLine_v = vth*0.8862*aa/b*(1. - t->Pesc);
		}


		/* we want full width, not half width */
		WidthLine_v *= 2.;

	}
	else
	{
		/* wind */
		r = a*tau/PI;
		if( r <= 1. )
		{
			WidthLine_v = vth*sqrt(log(MAX2(1.,tau))*.2821);
		}
		else
		{
			WidthLine_v = 2.*wind.windv0;
			if( r*vth <= WidthLine_v )
			{
				WidthLine_v = vth*r*log(WidthLine_v/(r*vth));
			}
		}
	}


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