# HG changeset patch
# User Marco Caliari <marco.caliari@univr.it>
# Date 1244662021 -7200
# Node ID 22a4433b41e3153f0fc573ba4d1cef84c550ae1a
# Parent  29563379fa9bdc3af9bb6045c3e4791e4afcb9cb
Better waypoint transform pour quadgk

diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,11 +1,15 @@
+2009-06-10  Marco Caliari <marco.caliari@univr.it>
+
+	* general/quadgk.m: Better waypoint transform.
+
 2009-06-09  David Bateman  <dbateman@free.fr>
 
-	* general/guadgk.m: Add test case and fixed doubly infinite 
+	* general/quadgk.m: Add test case and fixed doubly infinite 
 	waypoint transform for x = 0 case.
 
 2009-06-09  Marco Caliari <marco.caliari@univr.it>
 
-	* general/guadgk.m: Fix doubly infinite transformation to the finite
+	* general/quadgk.m: Fix doubly infinite transformation to the finite
 	interval.
 
 2009-06-08  Ben Abbott <bpabbott@mac.com>
diff --git a/scripts/general/quadgk.m b/scripts/general/quadgk.m
--- a/scripts/general/quadgk.m
+++ b/scripts/general/quadgk.m
@@ -193,18 +193,16 @@
       ##   g(t)  = t / (1 - t^2)
       ##   g'(t) =  (1 + t^2) / (1 - t^2) ^ 2
       ## waypoint transform is then
-      ##   t =  (-1 + sqrt(1 + 4 * g(t) .^ 2)) ./ (2 * g(t))
+      ##   t =  (2 * g(t)) ./ (1 + sqrt(1 + 4 * g(t) .^ 2))
       if (!isempty (waypoints))
-	trans = @(x) (-1 + sqrt(1 + 4 * x .^ 2)) ./ (2 * x);
+	trans = @(x) (2 * x) ./ ((1 + sqrt(1 + 4 * x .^ 2));
 	subs = [-1; trans(waypoints); 1];
-	## The waypoint transform gives NaN for x=0;
-	subs(isnan(subs)) = 0;
       else
 	subs = linspace (-1, 1, 11)'; 
       endif
       h = 2;
       trans = @(t) t ./ (1 - t.^2);
-      f = @(t) f (trans(t)) .* (1 + t .^ 2) ./ ((1 - t .^ 2) .^ 2);
+      f = @(t) f (t ./ (1 - t .^ 2)) .* (1 + t .^ 2) ./ ((1 - t .^ 2) .^ 2);
     elseif (isinf(a))
       ## Formula defined in Shampine paper as two separate steps. One to
       ## weaken singularity at finite end, then a second to transform to