/*widlin determine half width of any line with known optical depths */
#include <string.h>
#include "cddefines.h"
#include "physconst.h"
#include "redis.h"
#include "phycon.h"
#include "pressure.h"
#include "touton.h"
#include "rtescprob.h"
#include "wind.h"
#include "widlin.h"

double widlin(double tauin, 
	  double tauout, 
	  double a, 
	  double vth)
{
	double aa, 
	  atau, 
	  b, 
	  r, 
	  t, 
	  tau, 
	  therm, 
	  widlin_v;

#	ifdef DEBUG_FUN
	fputs( "<+>widlin()\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
	 * */
	if( touton.lgTauOutOn )
	{
		tau = tauout/2.;
		t = MIN2(tauin,tauout-tauin);
		t = MAX2(tauin/100.,t);
	}
	else
	{
		tau = tauin;
		t = tauin;
	}

	/* max optical depth is thermalization length */
	therm = 5.3e16/phycon.eden;
	if( tau > therm )
	{
		pressure.lgPradDen = TRUE;
		tau = therm;
	}

	if( (strcmp(redis.chRedis,"COMP") == 0 ) ||
		 (strcmp(redis.chRedis,"INCO") == 0 ) )
	{
		/* converted from go-to 15 june 95
		 * */
		atau = log(MAX2(0.0001,tau));
		if( tau <= 20. )
		{
			aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau;
			b = 6.5*tau - atau;
		}
		else
		{
			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);
		}

		/* this is half width of line */
		widlin_v = vth*0.8862*aa/b*(1. - escinc(t,a));

	}
	else if( strcmp(redis.chRedis,"WIND") == 0 )
	{
		/* WIND */
		r = a*tau/PI;
		if( r <= 1. )
		{
			widlin_v = vth*sqrt(log(MAX2(1.,tau))*.2821);
			
#	ifdef DEBUG_FUN
			fputs( " <->widlin()\n", debug_fp );
#	endif
			return( widlin_v );
		}

		widlin_v = 2.*wind.windv0;
		if( r*vth > widlin_v )
			{ 
#	ifdef DEBUG_FUN
			fputs( " <->widlin()\n", debug_fp );
#	endif
			return( widlin_v );}

		widlin_v = vth*r*log(widlin_v/(r*vth));
	}

	else
	{
		fprintf( ioQQQ,"widlin called with insane redistribution function ==%s==\n",
			redis.chRedis );
		puts( "[Stop in widlin]" );
		exit(1);
	}


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