#!/afs/sipb/project/perldev/p -w
# convert "return" output into plottable data

use strict;
use IO::File;

use vars qw( $dir );

$dir = shift(@ARGV)
  or die "usage: $0 directory\n";
$dir .= '/' unless $dir =~ m@/$@;

use vars qw( @xlim @ylim $xstep $ystep @ilim @jlim @angles );
use vars qw( $i $j $x $y );
use vars qw( %files $OUT $max_rp );

@xlim = (-.72,1.38);
@ylim = (-.42,.6);
$xstep = $ystep = .06;

@angles = map sprintf('%03d', 10*$_), 0..17;

@ilim = map round($_/$xstep), @xlim;
@jlim = map round($_/$ystep), @ylim;

$max_rp = 0;

for ($i = $ilim[0];  $i <= $ilim[1];  $i++) {
  print "col $i\n";
  $x = $i * $xstep;
  for ($j = $jlim[0];  $j <= $jlim[1];  $j++) {
    $y = $j * $ystep;

    my @fns = map "${dir}return_$x,$y"."_$_", @angles;
    if (! scalar grep( (! -e $_) || (-s _ != 2853), @fns )) {
      # All files OK; .
      my @fds = map IO::File->new($_, 'r'), @fns;

      my $first = 1;
    LINE:
      while (1) {
	my @lines = map scalar <$_>, @fds;
	last LINE unless defined $lines[0];
	($first = 0, next LINE) if $first;
	# error checking? ha!

	my $T = (split(/\s+/, $lines[0], 2))[0];
	my $name = "ret_$T";  $name =~ s/\./_/g;  $name .= '.m';

	if (exists $files{$name}) {
	  $OUT = IO::File->new($name, 'a');

	} else {
	  $OUT = IO::File->new($name, 'w');
	  $OUT->print(<<"EndOfStartup");
%%% do not edit by hand

theta = [@angles] * pi/180;
sinth = sin(theta);
costh = cos(theta);

global scale

title('T = $T');
hold on;

EndOfStartup
	  $files{$name} = 1;
	}

	my @rp = map( (split(/\s+/, $_, 3))[1], @lines );
	$max_rp = max($max_rp, @rp);

	print $OUT "retpoint($x, $y, scale*[@rp], sinth, costh);\n";
	$OUT->close;
      }
    }

  }
}

$OUT = IO::File->new('all_files.m', 'w');
$OUT->print(<<"EndOfFtext");
%%% list of files

global scale dist f_text
if isempty(scale),  get_scale;  end
if isempty(dist),   dist=.03;   end
f_text = sprintf('f=%g', scale.^2/dist);

EndOfFtext

foreach (sort keys %files) {
  IO::File->new($_, 'a')->print("\nhold off\naxis equal\n");
  s/\.m$//;
  s/_/ /;  # first one only!
  $OUT->print("$_\n");
}
$OUT->close;

IO::File->new('get_scale.m', 'w')->print(<<"EndOfScaleCalc");
%%% default for scale
global scale dist

maximum = $max_rp;
dist = 0.06/2;

scale = sqrt(dist) / maximum;
EndOfScaleCalc

sub round ($) { return int($_[0] + ($_[0]>0 ? +.5 : -.5)); }
sub max (@) {
  my $m = shift;
  foreach (@_) { $m = $_ if $m < $_; }
  $m;
}