#!/usr/bin/perl

# program to spit out hstcp values

my ( $Wl, $Wh, $Pl, $Ph, $B ) = @ARGV;

#$w   current window size
#$Wl  Low_Window
#$Wh  High_Window
#$Pl  Low_Loss
#$Ph  High_loss
#$B   High Decrease (ie decrease when at $Wl


print "Low_Window=$Wl\nHigh_Window=$Wh\nLow_Loss=$Pl\nHigh_Loss=$Ph\nHigh Decrease=$B\n\n";

$S = ( log($Wh) - log($Wl) ) / ( log($Ph) - log($Pl) );

# make $S less accurate - just like Sally's
$S = sprintf( "%.3f", $S );

print "S=$S\n";

# set $w to be the minimum
$w = $Wl;


while ( $w <= $Wh ) {
       
    $b = ($B - 0.5) * ( ( log($w) - log($Wl) ) / ( log($Wh) - log($Wl) ) ) + 0.5;

    $p_ =  $Pl / $Wl**(1/$S);
    # make p less accurate
    $p_ = sprintf( "%.3f", $p_ );

#    print "  p_=$p_\n";

#    $p = ( $Pl * $w**(1/$S) ) / ($Wl**(1/$S));
#    $p = $p_ * $w**(1/$S);

    $p_prime = sprintf( "%.1f", 1 / $p_ );
    $S_prime = sprintf( "%.1f", -1*1/$S );

#    print "p'=$p_prime, s'=$S_prime\n";

#    $a = ( ($w**2) * 2.0 * $b * $p ) / ( 2.0 - $b );
    $a = ( ($w**2) * 2.0 * $b ) / ( (2.0 - $b) * $w**$S_prime * $p_prime);

    if ( $w eq $Wl || $a > $olda + 1) {
	printf ("%6d\t%.1f\t%3.2f\n", $w, $a, $b);

#    print ("$w\t$a\t$b\n");
	$olda = $a;
    }

    $w++;
}