Skip to content

Commit

Permalink
tidy example 22, PDL version
Browse files Browse the repository at this point in the history
  • Loading branch information
mohawk2 committed Mar 30, 2024
1 parent ff2ce74 commit cf1c5c0
Showing 1 changed file with 18 additions and 39 deletions.
57 changes: 18 additions & 39 deletions t/x22.pl
Original file line number Diff line number Diff line change
Expand Up @@ -163,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 $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) = minmax($x);
my ($ymin, $ymax) = minmax($y);
my ($zmin, $zmax) = minmax($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
}

#--------------------------------------------------------------------------
Expand Down

0 comments on commit cf1c5c0

Please sign in to comment.