#/usr/local/bin/perl

require "/data1/WWW/Tropical/Bin/RungeKutta.pl";
require "/data1/WWW/Tropical/Bin/StrikeProb.pl";

$i = 9;
$r_over_e = 25;
$s_over_e = 30;

$piover2 = atan2(1,0);
$pi = $piover2*2;

$a1 = &AVGINT($i,$r_over_e,$r_over_e**2,$s_over_e,$s_over_e**2);
$a2 = &AVGINTRK4($i,$r_over_e,$r_over_e**2,$s_over_e,$s_over_e**2);
$a3 = &AVGINT2($i,$r_over_e,$r_over_e**2,$s_over_e,$s_over_e**2);


print "$a1 $a2 $a3\n";

$kh = 0;
$k = 0;
$k0 = 0;
$k00 = 0;
for ($i=0; $i<=1000; ++$i) {
   $gp = &StrikeProb'GAMMP($i+1,$s_over_e*$s_over_e);
   $kh += $gp*&StrikeProb'AVGINT2($i,$r_over_e,$r_over_e*$r_over_e,$s_over_e,$s_over_e*$s_over_e);
   $k += exp(-($r_over_e*$r_over_e)+log($r_over_e)*(2*$i)-&StrikeProb'LNFACTORIAL($i))*$gp;
   $k0 += $piover2*$gp*$gp;
   $k00 += $gp/($s_over_e*$s_over_e);
}
$kh = $kh / ($piover2*$s_over_e*$s_over_e);
$k0 = $k0 / ($piover2*$s_over_e*$s_over_e);
print "$r_over_e $s_over_e k: $k, kh: $kh, k0: $k0 k00: $k00 \n";

exit;

sub AVGINTRK4 {
   my($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) = @_;
   my($rho,$rho2,$h,$integral,$step,$lowbound,$upbound,$term,$nterms,$facti);

   $lowbound = $r_over_e - $s_over_e;
   if ($lowbound < 0) {$lowbound = 0;}
   $upbound = ($r_over_e+$s_over_e);
   @RK4PARAMS = ($i,&LNFACTORIAL($i),$r_over_e,$s_over_e,$r_over_e2,$s_over_e2);   # must have global scope
   my $t = $lowbound;
   my $nsteps = 35;
   my $dt = ($upbound-$lowbound)/($nsteps);
   my (@y,@errors,$epsilon);
   $errors[0] = 0.01;
   $y[0] = 0;
   my $call = 0;

   $NCALLS = 0;
   while ($t+$dt < $upbound) {
      ($t, $dt, @y) = rk4_auto(\@y,\&dydt,$t,$dt,\@errors);
      #print "call ",$call++," ",$t," ",$dt,"\n";
   }
   $dt = $upbound - $t;
   if ($dt) {($t, @y) = &rk4( \@y, \&dydt, $t, $dt );}
   print "auto used ",$NCALLS," function calls\n";

   my @ystep = (0);
   my $tstep = $lowbound;
   my $dtstep = ($upbound-$lowbound)/$nsteps;

   $NCALLS = 0;
   for ($step=0; $step<=$nsteps-1; ++$step) {
      ($tstep, @ystep) = &rk4( \@ystep, \&dydt, $tstep, $dtstep );   
   }
   print "stepwise got ",$tstep," ",$ystep[0],"\n";
   print "stepwise used ",$NCALLS," function calls\n";

   print "Final: ",$t," ",$dt," ",$lowbound," ",$upbound,"\n";
   $y[0];
}


sub dydt {
   my ($t,@y) = @_;
   my ($rho,$rho2,$term);
   $rho = $t;
   $rho2 = $rho*$rho;
   my ($i,$LNFACTI,$r_over_e,$s_over_e,$r_over_e2,$s_over_e2) = @RK4PARAMS;

   if ($rho <= 0) {
      $term = 0;
   } else {
      $term = exp(log($rho)*(2*$i+1) - $LNFACTI - $rho2);
   }
   if (($rho+$r_over_e) <= $s_over_e || $r_over_e == 0) {
      $term *= $pi;
   } else {
      $term *= &ACOS(($r_over_e2-$s_over_e2+$rho2)/(2*$r_over_e*$rho));
   }
   ++$NCALLS;
   return ($term);
}
                                                                                                                
# ---------------------------

@gammacof = (76.18009173,
             -86.50532033,
             24.01409822,
             -1.231739516,
             0.120858003e-2,
             -0.536382e-5);

sub GAMMALN {

   # ln of the gamma function.

   my($z) = @_;
   my($answer,$x,$y,$tmp,$j,$ser);
   $y = $z;
   $x = $z;
   $tmp = $x+5.5;
   $tmp -= ($x+0.5)*log($tmp);
   $ser = 1.000000000190015;
   for ($j=0; $j<=5; $j++) {$ser += $gammacof[$j]/(++$y);}
   -$tmp+log(2.5066282746310005*$ser/$x);
}

sub GAMMA {

   # The gamma function

   my($z) = @_;
   exp(&GAMMALN($z));
}


sub LNFACTORIAL {

   # log factorial function
   # Positive integers only

   my($x) = @_;

   if ($x > 12) {&GAMMALN($x+1);}
   else {
      log(&FACTORIAL($x));
   }
}

sub FACTORIAL {

   # factorial function
   # Positive integers only

   my($x) = @_;

   if ($x > 12) {&GAMMA($x+1);}
   else {
      if ($x <= 1) { 1;} else {$x*&FACTORIAL($x-1);}
   }
}

sub ACOS {

# Compute the arc cosine
# Input: Cosine (forced to -1 <= cosine <= 1)
# Output: Angle in radians (0<angle<pi)

   my($cosine) = @_;
   my($argument);
   $argument = 1.0-$cosine**2;
   if ($argument<=0) {
      if ($cosine > 0) {0;} else {$pi;}
   } else {
      $piover2 - atan2($cosine/(sqrt($argument)),1);
   }
}

sub tan { sin($_[0]) / cos($_[0]) }

sub AVGINT2 {

    # Trapezoidal rule integration

    my($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) = @_;
    my($rho,$rho2,$h,$integral,$step,$lowbound,$upbound,$term,$nterms,$lnfacti,$facti);
    my($lastterm,$lastrho);
    
    $integral = 0;
    $nterms = 30;  # if you change this, change $step below too
    $lowbound = $r_over_e - $s_over_e;
    $upbound = $r_over_e + $s_over_e;
    if ($lowbound < 0) {$lowbound = 0;}
    # use a variable step size which is small for
    # small rho but gets bigger for larger rho
    # 210 is 1+2+3+4+...+19+20
    my $nsum = 0;
    my $ns;
    for ($ns=0; $ns<=$nterms; ++$ns) {$nsum += $ns;}
    $step = ($upbound-$lowbound)/$nsum;  # 210 only works for nterms = 20
    print "nterms: $nterms, nsum: $nsum, step: $step\n";
    $lnfacti = &LNFACTORIAL($i);
    $lastrho = $lowbound;
    $lastterm = 0;
    for ($h=0; $h<=$nterms; ++$h) {
       $rho = $step*$h + $lastrho;
       $rho2 = $rho*$rho;
       if ($rho <= 0) {
          $term = 0;
       } else {
          $term = exp(log($rho)*(2*$i+1) - $lnfacti - $rho2);
       }
       if (($rho+$r_over_e) <= $s_over_e || $r_over_e == 0) {
	  $term *= $pi;
       } else {
          $term *= &ACOS(($r_over_e2-$s_over_e2+$rho2)/(2*$r_over_e*$rho));
       }
       if ($h) {$integral += (0.5*($term+$lastterm)*($rho-$lastrho));}
       $lastrho = $rho;
       $lastterm = $term;
    }
    print "check these are equal $rho $upbound\n";
    $integral;

}



sub AVGINT {

    # Trapezoidal rule integration

    my($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) = @_;
    my($rho,$rho2,$h,$integral,$step,$lowbound,$term,$nterms,$lnfacti,$facti);
    
    $integral = 0;
    $nterms = 20;
    $lowbound = $r_over_e - $s_over_e;
    if ($lowbound < 0) {$lowbound = 0;}
    $step = ($r_over_e+$s_over_e-$lowbound)/$nterms;
    $lnfacti = &LNFACTORIAL($i);
    for ($h=0; $h<=$nterms; ++$h) {
       $rho = $step*$h + $lowbound;
       $rho2 = $rho*$rho;
       if ($rho <= 0) {
          $term = 0;
       } else {
          $term = exp(log($rho)*(2*$i+1) - $lnfacti - $rho2);
       }
       if (($rho+$r_over_e) <= $s_over_e) {
	  $term *= $pi;
       } else {
          $term *= &ACOS(($r_over_e2-$s_over_e2+$rho2)/(2*$r_over_e*$rho));
       }
       if ($h == 0 || $h == $nterms) {
          $integral += $term/2;
       } else {
          $integral += $term;
       }
    }
    $integral * $step;

}

sub GAMMP {
                                                                                                                               
   # Incomplete Gamma Function, Gamma(a)*P(a,x)
   # $ x must be positive or zero, $a must be positve
                                                                                                                               
   my($a,$x) = @_;
   my($gammp);
                                                                                                                               
   if ($x < ($a + 1)) {
      $gammp = &GSER($a,$x);
   } else {
      $gammp = 1 - &GCF($a,$x);
   }
                                                                                                                               
   $gammp;
                                                                                                                               
}

$FPMIN = 1.0e-30;
                                                                                                                               
sub GCF {
                                                                                                                               
   # Incomplete Gamma Function, Q(a,x)
   # $ x must be positive or zero, $a must be positve
                                                                                                                               
   my($a,$x) = @_;
   my($gln,$b,$c,$d,$h,$i,$an,$absd,$del,$absdel1);
                                                                                                                               
   $gln = &GAMMALN($a);
   $b = $x + 1.0 - $a;
   $c = 1.0/1.0e-30;
   $d = 1.0/$b;
   $h = $d;
   for ($i=1; $i<=$ITMAX; $i++) {
      $an = -$i*($i-$a);
      $b += 2.0;
      $d = $an*$d+$b;
      $absd = $d;
      if ($absd < 0.0) {$absd = -$absd;}
      if ($absd < $FPMIN) {$d = $FPMIN;}
      $c = $b + $an/$c;
      $absc = $c;
      if ($absc < 0.0) {$absc = -$absc;}
      if ($absc < $FPMIN) {$c = $FPMIN;}
      $d = 1.0/$d;
      $del = $d*$c;
      $h *= $del;
      $absdel1 = $del - 1.0;
      if ($absdel1 < 0) {$absdel1 = -$absdel1;}
      if ($absdel1 < $EPS) {last;}
   }
   exp(-$x+$a*log($x)-$gln)*$h;
     
}

$ITMAX = 100;
$EPS = 3.0e-7;
                                                                                                                               
sub GSER {
                                                                                                                               
   # Incomplete Gamma Function, P(a,x)
   # $ x must be positive or zero, $a must be positve
                                                                                                                               
   my($a,$x) = @_;
   my($gln,$gser,$ap,$del,$sum,$n,$absdel,$abssum);
                                                                                                                               
   $gln = &GAMMALN($a);
   if ($x <= 0.0) {
      $gser = 0.0;    # Input Error
   } else {
      $ap = $a;
      $del = $sum = 1.0/$a;
      for ($n=1; $n<=$ITMAX; $n++) {
         ++$ap;
         $del *= $x/$ap;
         $sum += $del;
         $absdel = $del;
         if ($absdel < 0.0) {$absdel = -$absdel;}
         $abssum = $sum;
         if ($abssum < 0.0) {$abssum = -$abssum;}
         if ($absdel < $abssum*$EPS) {
            last;
         }
      }
      $gser = $sum * exp(-$x+$a*log($x)-$gln)
   }
   $gser;
}

1