tidy & vectorise examples 8,11,22, replace bespoke f2mnmax with PDL::…

…minmax
PDLPorters · Mar 31, 2024 · d03805f · d03805f
1 parent 6349b6b
commit d03805f
Show file tree

Hide file tree

Showing 5 changed files with 52 additions and 145 deletions.
diff --git a/t/x08.pl b/t/x08.pl
@@ -61,66 +61,46 @@ sub cmap1_init {
     $s = pdl [0.0,    # minimum saturation
               0.0];   # minimum saturation
   } else {
-    $h = pdl [240,    # blue -> green -> yellow -> */
-              0];     # -> red
-
-    $l = pdl [0.6, 0.6];
-    $s = pdl [0.8, 0.8];
+    ($h, $l, $s) = (pdl(240, 0), pdl(0.6, 0.6), pdl(0.8, 0.8));
   }
 
   plscmap1n (256);
   plscmap1l (0, $i, $h, $l, $s, pdl []);
 }
 
-my $LEVELS = 10;
-
 # Parse and process command line arguments
-
 my $rosen;
 plParseOpts (\@ARGV, PL_PARSE_SKIP | PL_PARSE_NOPROGRAM);
 GetOptions ("rosen" => \$rosen);
 
-my $nlevel = LEVELS;
-my $indexxmin = 0;
-my $indexxmax = XPTS;
+my ($indexxmin, $indexxmax) = (0, XPTS);
 # parameters of ellipse (in x, y index coordinates) that limits the data.
 # x0, y0 correspond to the exact floating point centre of the index range.
-my $x0 = 0.5 * ( XPTS - 1 );
-my $a  = 0.9 * $x0;
-my $y0 = 0.5 * ( YPTS - 1 );
-my $b  = 0.7 * $y0;
-my $sombrero = !$rosen;
+my ($x0, $y0) = (0.5 * ( XPTS - 1 ), 0.5 * ( YPTS - 1 ));
+my ($a, $b) = (0.9 * $x0, 0.7 * $y0);
 
 # Initialize plplot
-
 plinit ();
 
-my $x = (sequence (XPTS) - int(XPTS / 2)) / int(XPTS / 2);
-$x *=  1.5
-  if $rosen;
-my $xx = $x->dummy(1,YPTS);
-
-my $y = (sequence (YPTS) - int(YPTS / 2)) / int(YPTS / 2);
-$y += 0.5
-  if $rosen;
-my $yy = $y->dummy(0,XPTS);
+my ($x, $y) = map +(sequence($_) - int($_ / 2)) / int($_ / 2), XPTS, YPTS;
+$x *= 1.5 if $rosen;
+$y += 0.5 if $rosen;
+my ($xx, $yy) = ($x->dummy(1,YPTS), $y->dummy(0,XPTS));
 
 my $z;
 if ($rosen) {
   $z = (1 - $xx) ** 2 + 100 * ($yy - ($xx ** 2)) ** 2;
   # The log argument may be zero for just the right grid.
   $z = log ($z);
-}
-else {
+} else {
   my $r = sqrt ($xx * $xx + $yy * $yy);
   $z = exp (-$r * $r) * cos (2.0 * pi * $r);
 }
 $z->inplace->setnonfinitetobad;
 $z->inplace->setbadtoval(-5); # -MAXFLOAT would mess-up up the scale
 
 my (@indexymin, @indexymax);
-my $i = sequence( XPTS );
-my $square_root = sqrt( 1. - hclip( ( ( $i - $x0 ) / $a ) ** 2, 1 ) );
+my $square_root = sqrt(1. - hclip(( (sequence(XPTS) - $x0) / $a ) ** 2, 1));
 # Add 0.5 to find nearest integer and therefore preserve symmetry
 # with regard to lower and upper bound of y range.
 my $indexymin = lclip( 0.5 + $y0 - $b * $square_root, 0 )->indx;
@@ -133,8 +113,8 @@ sub cmap1_init {
   $zlimited->index($i)->slice($j) .= $z->index($i)->slice($j);
 }
 
-my $zmin = min ($z);
-my $zmax = max ($z);
+my ($zmin, $zmax) = (min($z), max($z));
+my $nlevel = LEVELS;
 my $step = ($zmax - $zmin) / ($nlevel + 1);
 my $clevel = $zmin + $step + $step * sequence ($nlevel);
 
@@ -160,21 +140,17 @@ sub cmap1_init {
             "bnstu", "x axis", "bnstu", "y axis", "bcdmnstuv", "z axis");
     plcol0 (2);
 
+    cmap1_init(($ifshade == 0) || 0);
     if ($ifshade == 0) {        # diffuse light surface plot
-      cmap1_init (1);
-      plsurf3d ($x, $y, $z, 0, pdl []);
+      plsurf3d($x, $y, $z, 0, pdl []);
     } elsif ($ifshade == 1) { # magnitude colored plot
-      cmap1_init (0);
-      plsurf3d ($x, $y, $z, MAG_COLOR, pdl []);
+      plsurf3d($x, $y, $z, MAG_COLOR, pdl []);
     }
     elsif ($ifshade == 2) {     # magnitude colored plot with faceted squares
-      cmap1_init (0);
-      plsurf3d ($x, $y, $z, MAG_COLOR | FACETED, pdl []);
+      plsurf3d($x, $y, $z, MAG_COLOR | FACETED, pdl []);
     } elsif ($ifshade == 3) {   # magnitude colored plot with contours
-      cmap1_init (0);
-      plsurf3d ($x, $y, $z, MAG_COLOR | SURF_CONT | BASE_CONT, $clevel);
-    } else {                  # magnitude colored plot with contours and index limits.
-      cmap1_init (0);
+      plsurf3d($x, $y, $z, MAG_COLOR | SURF_CONT | BASE_CONT, $clevel);
+    } else {          # magnitude colored plot with contours and index limits.
       plsurf3dl(
         $x, $y, $zlimited, MAG_COLOR | SURF_CONT | BASE_CONT, $clevel,
         $indexxmin, $indexxmax, $indexymin, $indexymax,

diff --git a/t/x09.pl b/t/x09.pl
@@ -94,17 +94,6 @@ sub polar {
 #  plFree2dGrid ($cgrid2);
 }
 
-# f2mnmx
-#
-# Returns min & max of input 2d array
-
-sub f2mnmx {
-  my $f = shift;
-  my $fmin = min ($f);
-  my $fmax = max ($f);
-  return ($fmin, $fmax)
-}
-
 # Shielded potential contour plot example
 
 sub potential {
@@ -122,8 +111,8 @@ sub potential {
 
   my $rmax = 0.5 + (PRPTS - 1);
 
-  my ($xmin, $xmax) = f2mnmx ($xg);
-  my ($ymin, $ymax) = f2mnmx ($yg);
+  my ($xmin, $xmax) = minmax($xg);
+  my ($ymin, $ymax) = minmax($yg);
 
   my $x0 = ($xmin + $xmax) / 2;
   my $y0 = ($ymin + $ymax) / 2;
@@ -161,7 +150,7 @@ sub potential {
   my $div2i = sqrt (($xg - $d2i) ** 2 + ($yg + $d2i) ** 2 + $eps ** 2);
   my $z = $q1 / $div1 + $q1i / $div1i + $q2 / $div2 + $q2i / $div2i;
 
-  my ($zmin, $zmax) = f2mnmx ($z);
+  my ($zmin, $zmax) = minmax($z);
 
   # Positive and negative contour levels
 

diff --git a/t/x11.pl b/t/x11.pl
@@ -72,25 +72,13 @@ sub cmap1_init {
 
 plinit ();
 
-my $z = zeroes (XPTS, YPTS);
-
-my $x = 3 * (sequence (XPTS) - int(XPTS / 2)) / int(XPTS / 2);
-my $y = 3 * (sequence (YPTS) - int(YPTS / 2)) / int(YPTS / 2);
-
-# The code below may be vectorized to improve speed
-for (my $i = 0; $i < XPTS; $i++) {
-  my $xx = $x->index ($i);
-  for (my $j = 0; $j < YPTS; $j++) {
-    my $yy = $y->index ($j);
-    $z->slice ("$i,$j") .=
-      3. * (1.-$xx)*(1.-$xx) * exp(-($xx*$xx) - ($yy+1.)*($yy+1.)) -
-      10. * ($xx/5. - pow($xx,3.) - pow($yy,5.)) * exp(-$xx*$xx-$yy*$yy) -
-      1./3. * exp(-($xx+1)*($xx+1) - ($yy*$yy));
-  }
-}
-
-my $zmax = max ($z);
-my $zmin = min ($z);
+my ($x, $y) = map 3*(sequence($_) - int($_ / 2)) / int($_ / 2), XPTS, YPTS;
+my ($xx, $yy) = ($x->dummy(1,YPTS), $y->dummy(0,XPTS));
+my $z =
+  3. * (1.-$xx)*(1.-$xx) * exp(-($xx*$xx) - ($yy+1.)*($yy+1.)) -
+  10. * ($xx/5. - pow($xx,3.) - pow($yy,5.)) * exp(-$xx*$xx-$yy*$yy) -
+  1./3. * exp(-($xx+1)*($xx+1) - ($yy*$yy));
+my ($zmin, $zmax) = (min($z), max($z));
 my $step = ($zmax - $zmin) / ($nlevel + 1);
 my $clevel = $zmin + $step + $step * sequence ($nlevel);
 

diff --git a/t/x15.pl b/t/x15.pl
@@ -37,7 +37,6 @@
 
 sub plot1 ();
 sub plot2 ();
-sub f2mnmx ($);
 sub cmap1_init1	();
 sub cmap1_init2	();
 
@@ -78,7 +77,7 @@ sub main {
 	}
     }
 
-    ($zmin, $zmax) = f2mnmx ($z);
+    ($zmin, $zmax) = minmax($z);
 
     plot1 ();
     plot2 ();
@@ -283,13 +282,4 @@ ()
 
 }
 
-# f2mnmx
-#
-# Returns min & max of input 2d array
-
-sub f2mnmx ($) {
-    my $f = shift;
-    return (min ($f), max ($f));
-}
-
 main ();
diff --git a/t/x22.pl b/t/x22.pl
@@ -121,18 +121,10 @@ sub transform {
 # Vector plot of flow through a constricted pipe
 # with a coordinate transform
 sub constriction2 {
-    my $nx = nx;
-    my $ny = ny;
-    my $nc   = 11;
-    my $nseg = 20;
-
-    my $dx = 1.0;
-    my $dy = 1.0;
-
-    my $xmin = -$nx / 2 * $dx;
-    my $xmax = $nx / 2 * $dx;
-    my $ymin = -$ny / 2 * $dy;
-    my $ymax = $ny / 2 * $dy;
+    my ($nx, $ny, $nc, $nseg) = (nx, ny, 11, 20);
+    my ($dx, $dy) = (1.0, 1.0);
+    my ($xmin, $xmax) = (-$nx / 2 * $dx, $nx / 2 * $dx);
+    my ($ymin, $ymax) = (-$ny / 2 * $dy, $ny / 2 * $dy);
 
     plstransform( \&transform, $xmax );
 
@@ -163,13 +155,6 @@ sub constriction2 {
     plstransform( undef, undef );
 }
 
-sub f2mnmx {
-  my $f = shift;
-  my $fmin = min ($f);
-  my $fmax = max ($f);
-  return ($fmin, $fmax);
-}
-
 # Vector plot of the gradient of a shielded potential (see example 9)
 
 sub potential {
@@ -178,67 +163,46 @@ sub potential {
   # Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
   # Also put in smoothing term at small distances.
 
-  my $rmax = nr;
-
-  my $eps = 2;
-
-  my $q1 = 1;
-  my $d1 = $rmax / 4;
-
-  my $q1i = - $q1 * $rmax / $d1;
-  my $d1i = $rmax * $rmax / $d1;
-
-  my $q2 = -1;
-  my $d2 = $rmax / 4;
-
-  my $q2i = - $q2 * $rmax / $d2;
-  my $d2i = $rmax * $rmax / $d2;
-
+  my ($rmax, $eps, $q1, $q2) = (nr, 2, 1, -1);
+  my ($d1, $d2) = ($rmax / 4, $rmax / 4);
+  my ($q1i, $d1i) = (- $q1 * $rmax / $d1, $rmax * $rmax / $d1);
+  my ($q2i, $d2i) = (- $q2 * $rmax / $d2, $rmax * $rmax / $d2);
   my $r = (0.5 + sequence (nr))->dummy (1, ntheta);
   my $theta = (2 * pi / (ntheta - 1) *
                (0.5 + sequence (ntheta)))->dummy (0, nr);
-  my $x = $r * cos ($theta);
-  my $y = $r * sin ($theta);
-  my $cgrid2 = plAlloc2dGrid ($x, $y);
-  my $div1 = sqrt (($x - $d1) ** 2 + ($y - $d1) ** 2 + $eps * $eps);
-  my $div1i = sqrt (($x - $d1i) ** 2 + ($y - $d1i) ** 2 + $eps * $eps);
-  my $div2 = sqrt (($x - $d2) ** 2 + ($y + $d2) ** 2 + $eps * $eps);
-  my $div2i = sqrt (($x - $d2i) ** 2 + ($y + $d2i) ** 2 + $eps * $eps);
+  my ($x, $y) = ($r * cos ($theta), $r * sin ($theta));
+  my ($div1, $div1i) =
+    map sqrt(($x - $_) ** 2 + ($y - $_) ** 2 + $eps * $eps), $d1, $d1i;
+  my ($div2, $div2i) =
+    map sqrt(($x - $_) ** 2 + ($y + $_) ** 2 + $eps * $eps), $d2, $d2i;
   my $z = $q1 / $div1 + $q1i / $div1i + $q2 / $div2 + $q2i / $div2i;
   my $u = -$q1 * ($x - $d1) / ($div1**3) - $q1i * ($x - $d1i) / ($div1i ** 3)
           -$q2 * ($x - $d2) / ($div2**3) - $q2i * ($x - $d2i) / ($div2i ** 3);
   my $v = -$q1 * ($y - $d1) / ($div1**3) - $q1i * ($y - $d1i) / ($div1i ** 3)
           -$q2 * ($y + $d2) / ($div2**3) - $q2i * ($y + $d2i) / ($div2i ** 3);
-
-  my ($xmin, $xmax) = f2mnmx ($x);
-  my ($ymin, $ymax) = f2mnmx ($y);
-  my ($zmin, $zmax) = f2mnmx ($z);
-
-  plenv ($xmin, $xmax, $ymin, $ymax, 0, 0);
-  pllab ("(x)", "(y)",
-         "#frPLplot Example 22 - potential gradient vector plot");
-
+  my ($xmin, $xmax, $ymin, $ymax, $zmin, $zmax) = map minmax($_), $x, $y, $z;
   # Plot contours of the potential
   my $dz = ($zmax - $zmin) / nlevel;
   my $clevel = $zmin + (sequence (nlevel) + 0.5) * $dz;
+  # the perimeter of the cylinder
+  $theta = (2 * pi / (nper - 1)) * sequence (nper);
+  my $px = $rmax * cos ($theta);
+  my $py = $rmax * sin ($theta);
 
+  plenv ($xmin, $xmax, $ymin, $ymax, 0, 0);
+  pllab ("(x)", "(y)",
+         "#frPLplot Example 22 - potential gradient vector plot");
   plcol0 (3);
   pllsty (2);
+  my $cgrid2 = plAlloc2dGrid ($x, $y);
   plcont ($z, 1, nr, 1, ntheta, $clevel, \&pltr2, $cgrid2);
   pllsty (1);
   plcol0 (1);
-
   # Plot the vectors of the gradient of the potential
   plcol0 (2);
   plvect ($u, $v, 25.0, \&pltr2, $cgrid2);
   plcol0 (1);
-
-  # Plot the perimeter of the cylinder
-  $theta = (2 * pi / (nper - 1)) * sequence (nper);
-  my $px = $rmax * cos ($theta);
-  my $py = $rmax * sin ($theta);
-  plline ($px , $py);
-
+  plline ($px , $py); # Plot the cylinder
 }
 
 #--------------------------------------------------------------------------