#!/afs/athena/contrib/perl/perl
# plots output from radio-telescope interface

require "ctime.pl";
require "getopts.pl";

# multips -p -o c -j b,l -r # -c #

&parse_args;

($hgt,$wid) = (7 *72, 9.5 *72);  # inches * 72

@scale = ($wid / 30, $hgt / ($lim[1]-$lim[0]));

&ps_head($hgt,$wid,-$lim[0]*$scale[1],@lim);
foreach $i (@ARGV) {
    &plot_file($i,@scale);
}
&ps_foot;

print STDERR "done (total pages: $pg)\n";

sub plot_file {
    local($fn,@sc) = @_;
    local(@ARGV) = ($fn);
    local($cmd) = "moveto";
    local(@max) = undef;

    &ps_pre(++$pg,$fn);
    while (<>) {
	($x,@val) = split(/\s+/);
	if (&abs($val[1]-$val[0]-$val[2])>2e-6) {
	    $val[0]=$val[1]-$val[2];  # no warning here
	}
	print $x*$sc[0]," ",$val[$yval]*$sc[1]," $cmd\n";
	@max=($x,$val[$yval]) if (!defined(@max) || $max[1]<$val[$yval]);
	$cmd = "lineto";
    }
    &ps_delim(@max,@sc);
}

sub limits {
    local(@ARGV) = @_;
    local($x,@val,$tmp);
    local($min,$max) = (undef,undef);

    while (<>) {
	($x,@val) = split(/\s+/);
	$tmp = &abs($val[1]-$val[0]-$val[2]);
	if ($tmp>2e-6) {
	    $val[0]=$val[1]-$val[2];
	    warn("$0: miscalculated difference (fixed?)\n");
	}
	$min=$val[$yval] if (!defined($min) || $min>$val[$yval]);
	$max=$val[$yval] if (!defined($max) || $max<$val[$yval]);
    }
    ($min,$max);
}

sub abs { @_[0] >= 0 ? @_[0] : -@_[0]; }

sub ps_head {
    local($hgt,$wid,$yzero,$min,$max) = @_;

    print "%!PS-Adobe-1.0\n";
    print "%%Creator: aplot.pl\n";
    if (@ARGV && $#ARGV >= 1) {
	print "%%Title: ",@ARGV[0],"..",@ARGV[$#ARGV],"\n";
    } else {
	print "%%Title: ",@ARGV[0],"\n";
    }

    print "%%CreationDate: ",&ctime($^T),"\n";
    printf "/hgt $hgt def /wid $wid def\n";
    print "/yzero $yzero def\n";
    printf "/min (%.3f) def /max (%.3f) def\n\n",$min,$max;
    print <<"--END--";
/inch {72 mul} def
/cshow {dup stringwidth pop 2 div neg 0 rmoveto show} def
/rshow {dup stringwidth pop neg 0 rmoveto show} def
/circ {1 index 5 add 1 index moveto 5 0 360 arc} def

/setup { % (file)
  8.5 inch 0 translate 90 rotate
  /Helvetica findfont 30 scalefont setfont
  5.5 inch 25 moveto cshow
  /Helvetica findfont 15 scalefont setfont
  0 setlinewidth
  1 inch dup translate
  0 0 moveto 0 hgt lineto stroke
  -9 -10 moveto min rshow
  -9 hgt moveto max rshow
  0 yzero translate
  0 0 moveto wid 0 lineto stroke
  -9 -5 moveto (0.0) rshow
  /Helvetica findfont 20 scalefont setfont
  3 setlinewidth 1 setlinecap 1 setlinejoin
} def

%%EndProlog
--END--
}

sub ps_pre {
    print "\n%%Page: @_[0] @_[0]\n\n($label@_[1]) setup\n\n";
}

sub ps_delim {
    local($mx,$my,@sc) = @_;
    print "stroke\n";
    printf "%s %s moveto (%d, %.3f) cshow",$mx*$sc[0],$my*$sc[1]+5,$mx,$my;
    print "\nshowpage\n";
}

sub ps_foot {
    print "\n%%Trailer\n";
}

sub usage {
    
    die join("\n",(@_,<<"DIE_NOW"));
Usage:    aplot [-l num,num] [-c expr] [-d|-r] filename [filename...]

            -l   limits for Y-axis
            -c   circle if expr is true (doesn't work)
            -d   data spectra only
            -r   reference spectra only
DIE_NOW
}

sub parse_args {
    do Getopts('l:c:dr') || &usage;

    @ARGV || &usage("You must specify a filename.");
    $opt_d && $opt_r && &usage("-d and -r are mutually exclusive.");

    if (defined($opt_l)) {
	eval "\@lim = ($opt_l);" || &usage("Invalid Y axis limits.");
	($#lim = 1) || &usage("Need *two* Y axis limits.");
	($lim[0] < $lim[1]) || &usage("Bogus Y axis limits.");
    } else {
	@lim = &limits(@ARGV);
	$lim[0]=0 if $lim[0]>0;
    }

    if (defined($opt_c)) {
	$circ_expr = $opt_c;
    }

    ($yval,$label) = (0,"");
    ($yval,$label) = (1,"DATA: ") if $opt_d;
    ($yval,$label) = (2,"REF: ") if $opt_r;
}