#!/usr/bin/perl -w

print('# sample of vector field data: flow around Rankine vortex' . "\n");

($a,$b) = (1.00, 0.12);  # see below

print('# (u,v) =  (-b*y,    b*x   )  [r2<=1]' . "\n",
      '#       or (-b*y/r2, b*x/r2)  [r2>1] ' . "\n",
      '#       where r2 = (x*x+y*y)/(a*a), ' . "a=$a, b=$b" . "\n\n");

($xmin,$xmax) = (-2.5, 2.5);
($ymin,$ymax) = (-2.0, 2.0);
$dx = $dy = ($xmax-$xmin)/30;

$eps = 1.e-6;  # epsilon, very small

print('# format: x y u v' . "\n");

for( $xx = $xmin ;  $xx <= $xmax+$eps ;  $xx += $dx ) {
    for( $yy = $ymin ;  $yy <= $ymax+$eps ;  $yy += $dy ) {
	$x = $xx + (0.2*$dx)*(2.0*rand()-1.0);
	$y = $yy + (0.2*$dy)*(2.0*rand()-1.0);
	$r2 = ($x*$x+$y*$y)/($a*$a);
	if( $r2 > 1.0 ) {
	    ($u,$v) = (-$b*$y/$r2, $b*$x/$r2);
	} else {
	    ($u,$v) = (-$b*$y,     $b*$x    );
	}
	printf("%f %f %f %f\n", $x,$y, $u, $v );
    }
}
