require "/data1/WWW/Tropical/Bin/RungeKutta.pl";
                                                                                                                                        
sub AVGINTRK4 {

   # Do the integral using Runge-Kutta with auto step scaling.

   my($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) = @_;

   my $lowbound = $r_over_e - $s_over_e;
   if ($lowbound < 0) {$lowbound = 0;}
   my $upbound = ($r_over_e+$s_over_e);
   @RK4PARAMS = ($i,&LNFACTORIAL($i));   # must have global scope
   my $t = $lowbound;
   my $nsteps = 20;
   my $dt = ($upbound-$lowbound)/$nsteps;
   my @y,@errors,$step;
   $errors[0] = 0.00001;
   $y[0] = 0.0;

   while ($t+$dt < $upbound) {
      ($t, $dt, @y) = rk4_auto(\@y,\&dydt,$t,$dt,\@errors);
   }
   $dt = $upbound - $t;
   if ($dt) {($t, @y) = &rk4( \@y, \&dydt, $t, $dt );}

   $y[0];
}


sub dydt {
   my ($t,@y) = @_;
   my $rho,$rho2,$term;
   $rho = $t;
   $rho2 = $rho*$rho;
   my ($i,$LNFACTI) = @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));
   }
   return ($term);
}