#!/usr/bin/perl

#======================================================================

#                    A N T S I N T E G R A T E . P L 

#                    doc: Fri Feb 28 22:54:04 1997

#                    dlm: Fri Oct 15 23:14:11 2010

#                    (c) 1997 Andreas Thurnherr

#                    uE-Info: 27 62 NIL 0 0 72 2 2 4 NIL ofnI

#======================================================================

# Integrator Library; used by [integrate] and [transport]

#	([transport] is deprecated...)

# HISTORY:

#	Nov 14, 2000: - divorced from [integrate]

#	Nov 16, 2000: - added $opt_i

#				  - BUG: $opt_c forgot to output last record

#				  - BUG: copying of buffer had not worked

#	Nov 17, 2000: - changed $opt_i to $opt_l

#	Nov 24, 2000: - aimless changes

#	Mar 10, 2002: - added $dx output if $dfnr defined

#	Mar 30, 2002: - cosmetics

#	May 24, 2002: - added -b)ox to -f

#	May 25, 2002: - BUG: -c did not output records with missing y values

#				  - changed to dying on missing x values

#	Dec 20, 2005: - BUG: -cf produced short 1st record (dfnr not set)

#   Jul  1, 2006: - Version 3.3 [HISTORY]

#	Oct 15, 2010: - removed warning about integrand set to nan

# NOTES:

#

#	&integrate() can run in several moments:

#		0:	integral

#		1:	integral weighted by (signed) distance from zero

#		2:	integral weighted by square of distance from zero

#

#	&integrate() returns the single-value sum; without -f @sum is

#	set as a side-effect

#  

#	unless xmin,xmax are given as args, the integral is taken over the

#	entire data

#

#	records with missing y values have the integrand(s) set to NaN

#

#	because this routine was ones an integral part of [integrate] it

#	uses a number of global variables (ug-lee!):

#		$xfnr					y = f(x)

#		$opt_c					output sum after each step (cummulative)

#		$opt_f					integrate only given field

#		  $yfnr   (if $opt_f)	... this one

#		  $opt_b				use box rule

#		  $iScale (if $opt_f)	scale integrated values

#		  $opt_l  (if $opt_f)	output individual summands

sub integrate ($@)										# ret integral, set $m

{

	my($moment,$zero,$xmin,$xmax) = @_;					# optional args

	my($lastx,$cur,$curwt,$sum,$cury);

	my($dx,$r,$f,@out,@nan,$warned);

	for ($f=0; $f<$antsBufNFields; $f++) {				# prep nan output

		$nan[$f] = nan;

	}

	for ($r=0; $r<=$#ants_; $r++) {

		croak("$0: can't handle non-numeric x values\n")

			unless (numberp($ants_[$r][$xfnr]));

#	1st, find a valid x value, and a y value on -f (one field only).

#	On -c, generate output: 0 on open limits, nan (only possible

#	if -f is set) otherwise

		unless (defined($lastx)) {

			if (defined($opt_f)) {

				if ($opt_c) {

					@out = @{$ants_[$r]};

					$out[$yfnr] = numberp($ants_[$r][$yfnr]) ? 0 : nan;

					$out[$dfnr] = nan if (defined($dfnr));

					&antsOut(@out);

				}

				next unless (numberp($ants_[$r][$yfnr]));

				$lastx = $ants_[$r][$xfnr];

				croak("$0: lower limit ($xmin) < first valid x-value ($lastx)\n")

					if (defined($xmin) && $lastx > $xmin);  

				$lasty = $ants_[$r][$yfnr];

				next;

			}

			$lastx = $ants_[$r][$xfnr];

			if ($opt_c) {

				for ($f=0; $f<$antsBufNFields; $f++) {

					$out[$f] = ($f == $xfnr) ? $ants_[$r][$f] :

									(numberp($ants_[$r][$f]) ? 0 : nan);

				}

				&antsOut(@out);

			}

			next;

		}

#	next, update x&y while below lower limit for later interpolation

#	NB: only possible on -f!

		if (defined($xmin) && $ants_[$r][$xfnr] < $xmin) {

			if (numberp($ants_[$r][$yfnr])) {

				$lastx = $ants_[$r][$xfnr];

				$lasty = $ants_[$r][$yfnr];

            }

			if ($opt_c) {

				$ants_[$r][$yfnr] = nan;

				&antsOut(@{$ants_[$r]});

            }

			next;

		}

#	we have an x-value > min; is it valid?

#	NB: xmin is undefined once lower limit is handled

		croak("$0: Error: no data within integration limits\n")

			if (defined($xmin) && defined($xmax) && $ants_[$r][$xfnr] > $xmax);

#	finally! we have a valid x-value; if it's the first, interpolate

#	y at xmin if that's defined

#	undefined xmin at end so this code is not executed again

#	NB: xmin can only be defined on -f!

		if (defined($xmin)) {

			unless (numberp($ants_[$r][$yfnr])) {

				if ($opt_c) {

					$ants_[$r][$yfnr] = nan;

					&antsOut(@{$ants_[$r]});

				}

				next;

			}

			$lasty += ($xmin-$lastx) / ($ants_[$r][$xfnr]-$lastx)

						* ($ants_[$r][$yfnr]-$lasty);

			$lastx = $xmin; undef($xmin);

		}

#	it is also possible (though not on the first time round), that we

#	have just passed the upper limit (xmax); simulate normal code but

#	using xmax (and interpolated y value) instead of real data

#	NB: xmax can only be defined on -f!

		if (defined($xmax) && $ants_[$r][$xfnr] >= $xmax) {

			unless (numberp($ants_[$r][$yfnr])) {

				if ($opt_c) {

					$ants_[$r][$yfnr] = nan;

					&antsOut(@{$ants_[$r]});

				}

				next;

			}

			$dx = $xmax - $lastx;

			$cury = $lasty + $dx / ($ants_[$r][$xfnr]-$lastx)

								* ($ants_[$r][$yfnr]-$lasty);

			croak("$0: x-field must be monotonically increasing [xmax=$xmax, lastx=$lastx]\n")

				if ($dx < 0);

			$cur  = ($cury+$lasty)/2 * $dx;

			$cur *= (($xmax+$lastx)/2-$zero)**$moment if ($moment);

			$sum += $cur*$iScale;

			if ($opt_c) {								# cummulative

				@out = @{$ants_[$r]};					# copy everything

				$out[$yfnr] = $sum;

				$out[$dfnr] = $dx if (defined($dfnr));

				&antsOut(@out);

			} elsif ($opt_l) {							# individual summands

				@out = @{$ants_[$r]};					# copy everything

				$out[$yfnr] = $cur*$iScale;

				$out[$dfnr] = $dx if (defined($dfnr));

				&antsOut(@out);

            }

            $lastx = $ants_[$r][$xfnr];

			last;

		}

#	phoar! we are finally handling the normal case  

		$dx = $ants_[$r][$xfnr] - $lastx;				# calc dx

		croak("$0: x-field must be monotonically increasing [curx=$ants_[$r][$xfnr], lastx=$lastx]\n")

			if ($dx < 0);

		if (defined($opt_f)) {							# integrate single field

			unless (numberp($ants_[$r][$yfnr])) {

				if ($opt_c) {

					$ants_[$r][$yfnr] = nan;

					&antsOut(@{$ants_[$r]});

				}

				next;

			}

			$cur = $opt_b ?

				$ants_[$r][$yfnr] * $dx :				# box rule

				($ants_[$r][$yfnr]+$lasty)/2 * $dx; 	# interpolate

			$cur *= (($ants_[$r][$xfnr]+$lastx)/2-$zero)**$moment

				if ($moment);

			$sum += $cur*$iScale;

			$lasty = $ants_[$r][$yfnr];

			if ($opt_c) {								# cummulative

				@out = @{$ants_[$r]};					# copy everything

				$out[$yfnr] = $sum;

			} elsif ($opt_l) {							# individual summands

				@out = @{$ants_[$r]};					# copy everything

				$out[$yfnr] = $cur*$iScale;

            }			

		} else {										# integrate all

			for ($f=0; $f<$antsBufNFields; $f++) {

				next if ($f == $xfnr);					# except x-field

				if (numberp($ants_[$r][$f]))	{		# val found

					$sum[$f] += $ants_[$r][$f] * $dx	# box-rule (no interp)

						unless isnan($sum[$f]);			# had missing vals

				} else {								# val missing 

					$sum[$f] = nan; 					# mark for later

#					unless ($warned) {					# warn user

#						&antsInfo("Warning: integrand(s) set to nan due to missing vals");

#						$warned = 1;

#					}

				}

				$out[$f] = $sum[$f] if ($opt_c);

			}

		}

		if ($opt_c || $opt_l) {

			$out[$xfnr] = $ants_[$r][$xfnr];

			$out[$dfnr] = $dx if (defined($dfnr));

			&antsOut(@out);

		}

		$lastx = $ants_[$r][$xfnr];						# last good x

	}

	croak("$0: empty input!\n")							# never initialized

	    unless (defined($lastx));

	croak("$0: upper limit > last valid x-value!\n")

	    if (defined($xmax) && $xmax>$lastx);

	return $sum;

}

1;