#!/usr/bin/perl -w

# Humlicek 1979 Quant. Spectrosc. Radiat. Transfer _21_ 309 in perl
# Calculate the "real" part of the Voigt function (ignore imaginery part)
# This can be substituted into Eqn 32 from lecture 1 to get the tau as
# a function of lambda.

# Need to define x and y for each absorption line as per assignment 
# definition of Voigt profile

$b1=0.850;
@t=(0.314240376,0.947788391,1.59768264,2.27950708,3.02063703,3.8897249);
@c=(1.01172805,-0.75197147,1.2557727e-2,1.00220082e-2,-2.42068135e-4,
    5.00848061e-7);
@s=(1.393237,0.231152406,-0.155351466,6.21836624e-3,9.19082986e-5, 
    -6.27525958e-7);

$u=0;

$y1 = $y + 1.500;
$y2 = $y1*$y1;
$b2 = 18.100*$y + 1.6500;
    
if ( $y > $b1 || abs($x) < $b2){
    for ($i=0;$i<6;$i++){
	$r  = $x - $t[$i];
	$d  = 1.00/($r*$r+$y2);
	$d1 = $y1*$d;
	$d2 = $r*$d;
	$r  = $x + $t[$i];
	$d  = 1.00/($r*$r+$y2);        
	$d3 = $y1*$d;
	$d4 = $r*$d;
	$u  = $u + $c[$i]*($d1+$d3) - $s[$i]*($d2-$d4);
    }
}

else{
    if (abs($x) < 12.00){$u = exp(-$x*$x)};
    $y3 = $y + 3;
    for($i=0;$i<6;$i++){
	$r  = $x - $t[$i];
	$r2 = $r*$r;
	$d  = 1.00/($r2+$y2);
	$d1 = $y1*$d;
	$d2 = $r*$d;        
	$u  = $u + $y*($c[$i]*($r*$d2-1.5*$d1)+$s[$i]*$y3*$d2)/($r2+2.25);
	$r  = $x + $t[$i];
	$r2 = $r*$r;
	$d  = 1.0/($r2+$y2);        
	$d3 = $y1*$d;
	$d4 = $r*$d;
	$u  = $u + $y*($c[$i]*($r*$d4-1.5*$d3)-$s[$i]*$y3*$d4)/($r2+2.25);
    }
}

