#!/usr/bin/perl
#
# Program rotated_AN.pl
#
# Purpose: 
# --------
#  To plot rotated D-terms from all experiments on sample plot
#
# usage: rotated_AN.pl *dterms*
#
# dterm files must have format bl111x.dterms.xxx, where xxx is the evpa.rotation
#
# assume 4 IFs are in the experiment
$N_IF = 4;  print "Assuming 4 IFs in experiment.\n";

# limits for plotting DTERMS
$dt_size = 0.05; 

# plotting panels per page = $nx * $ny
$nx = 2;     
$ny = 2;    

# set up some variables
$n = 0;
$IF = 1;
$real = 1;
$newan  = MEDIAN_DTERMS;    #new antenna table

foreach $ANTABLE (@ARGV) {

    open(AN,$ANTABLE) || die ("cannot open $ANTABLE: $!\n");
    $TEMPLATE = $ANTABLE;
    print "Processing $ANTABLE...  plotted with symbol $n in PGPLOT\n";
    $evparot[$n] = 2*substr($ANTABLE,14,6)*3.14159/180.;


# Skip AN file header 
    until (<AN> =~ /BEGIN/) {
	next;
    }
    
# Read in DTERMS
    while (<AN>) {
	$line = $_;
	
        if ($line !~ /END/) {
	    $ant = substr($line, 6, 2);

	    $aname = substr($line,10,2);  #Assign antenna name to number

	    if ($aname =~ /\w/) {$antname[$ant] = $aname;}
            if($real == 1) {      
	     $rcp_re[$ant][$IF][$n] = substr($line ,130,15);
	     $lcp_re[$ant][$IF][$n] = substr($line ,168,15);
             $real = 0; 
	    } else {
	     $rcp_im[$ant][$IF][$n] = substr($line,130,15);
	     $lcp_im[$ant][$IF][$n] = substr($line,168,15);
	     $real = 1; 
             $IF++;  if($IF > $N_IF) { $IF = 1; }
	    }
	}
    }
 


    close AN;

    $line = "";
    $n++;
    $n_ant = $ant;
}


###############################################
#
#  Now do some plotting + find median D-terms
#
###############################################

use PGPLOT;  # Load PGPLOT module
print "PGPLOT module version $PGPLOT::VERSION\n\n";

pgqinf("VERSION",$val,$len);
print "PGPLOT $val library\n\n";

# "?" will prompt for device

$dev = "?" unless defined $dev; 

pgbegin(0,$dev,$nx,$ny);  # Open plot device 

pgscf(2);             # Set character font
pgslw(4);             # Set line width
pgsch(1.2);           # Set character height
pgsci(3);

for($ant = 1; $ant < $n_ant+1;  $ant++) {
    for($IF = 1; $IF < $N_IF+1; $IF++) {
	pgenv(-$dt_size,$dt_size,-$dt_size,$dt_size,1,1);    
        pglabel("Real","Imaginary","Ant: $antname[$ant] -- IF: $IF"); # Labels
        $j = 0; $k = 0;
	for($i = 0; $i < $n; $i++) {
            pgsci(2);
	    pgsch(1.5);
	    $xx = $rcp_re[$ant][$IF][$i];
	    $yy = $rcp_im[$ant][$IF][$i];

	    $x = $xx*cos($evparot[$i])-$yy*sin($evparot[$i]);
	    $y = $xx*sin($evparot[$i])+$yy*cos($evparot[$i]);

            if(sqrt($x*$x) > $dt_size || sqrt($y*$y) > $dt_size) {
		printf "R point out of range for epoch %d .\n",$i+1;
                # note that this point is still 'plotted' 
                # (although it falls outside the viewport) 
                # and is used for median analysis
            }
            if($x == 0 && $y == 0) {
		print "R point identically zero, not plotted.\n";
                # such a point is not used for median analysis
            } else {
		pgpoint(1,$x,$y,$i+97); 
#		pgpoint(1,$x,$y,$i); 

                ########################################
                # Store plotted point for finding median
                ########################################
                $Mrcp_re[$j] = $x;
                $Mrcp_im[$j] = $y;
                $j++
	    }
            pgsci(4);
	    $xx= $lcp_re[$ant][$IF][$i];
	    $yy= $lcp_im[$ant][$IF][$i];

	    $x = $xx*cos(-$evparot[$i])-$yy*sin(-$evparot[$i]);
	    $y = $xx*sin(-$evparot[$i])+$yy*cos(-$evparot[$i]);
            if(sqrt($x*$x) > 0.1 || sqrt($y*$y) > 0.1) {
		print "L point out of range.\n";
                # note that this point is still 'plotted' 
                # (although it falls outside the viewport) 
                # and is used for median analysis
            }
            if($x == 0 && $y == 0) {
		print "L point identically zero, not plotted.\n";
                # such a point is not used for median analysis
            } else {
		pgpoint(1,$x,$y,$i+97);
#		pgpoint(1,$x,$y,$i);
                ########################################
                # Store plotted point for finding median
                ########################################
                $Mlcp_re[$k] = $x;
                $Mlcp_im[$k] = $y;
                $k++
	    }
        }
        #################################
        # 
        # Find median values for this IF 
        #
        #################################
        @R_re = sort {$a <=> $b} @Mrcp_re; 
        @R_im = sort {$a <=> $b} @Mrcp_im; 
        @L_re = sort {$a <=> $b} @Mlcp_re; 
        @L_im = sort {$a <=> $b} @Mlcp_im; 
        undef @Mrcp_re;    # these need to be undefined for the next antenna
        undef @Mrcp_im;
        undef @Mlcp_re;
        undef @Mlcp_im;

        if($j == 0) {
            $rx = 0; 
	    $ry = 0;
        } elsif(int($j/2)) { 
            $rx = $R_re[int($j/2)];
            $ry = $R_im[int($j/2)];
        } else {
	    $rx = ($R_re[int($j/2)] + $R_re[int($j/2)-1])/2;
	    $ry = ($R_im[int($j/2)] + $R_im[int($j/2)-1])/2;
        }
	pgsci(5); pgsch(4.8); 
#pgpoint(1,$rx,$ry,2); 
        if($k == 0) {
            $lx = 0; 
	    $ly = 0;
        } elsif(int($k/2)) { 
            $lx = $L_re[int($k/2)];
            $ly = $L_im[int($k/2)];
        } else {
	    $lx = ($L_re[int($k/2)] + $L_re[int($k/2)-1])/2;
	    $ly = ($L_im[int($k/2)] + $L_im[int($k/2)-1])/2;
        }
#        pgpoint(1,$lx,$ly,2);
 pgsch(1.2); pgsci(3);
	########################
        #0.9
        # Store median values
        #
	########################
        $Med_R_re[$ant][$IF] = $rx; 
        $Med_R_im[$ant][$IF] = $ry; 
        $Med_L_re[$ant][$IF] = $lx; 
        $Med_L_im[$ant][$IF] = $ly; 
    } 
}

pgend;

######################################
#
# Write out median D-terms
#  -- template based on last AN table
#     processed above
#
#######################################

open(OLDAN,$TEMPLATE) || die ("cannot open $ANTABLE: $!\n");  
open(OUTAN,">$newan") || die "cannot open $avgright: $!\n";

# Copy header to new AN table
until ($line =~ /BEGIN/) {
    $line = <OLDAN>; 
    print OUTAN $line;
    next;
}

# Output the rest of the new AN table 
$IF = 1;
$real = 1;

while (<OLDAN>)  {
    if (/END/) {
	print OUTAN $_;
	next;
    } else {
	$ant = substr($_, 6, 2);
        $dat[0] = substr($_,0,129);
	$dat[1] = substr($_, 147,20);
        if($real == 1) {      
	    $right = sprintf("%8.6E", $Med_R_re[$ant][$IF]);
	    $left = sprintf("%8.6E", $Med_L_re[$ant][$IF]);
            $real = 0; 
	} else {
	    $right = sprintf("%8.6E", $Med_R_im[$ant][$IF]);
	    $left = sprintf("%8.6E", $Med_L_im[$ant][$IF]);
	    $real = 1; 
	    $IF++;  if($IF > $N_IF) { $IF = 1; }
	}
	write(OUTAN);
    }
}
format OUTAN  =
@<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @>>>>>>>>>>>>>  @>>>>>>>>>>>>>>>>>>>>> @>>>>>>>>>>>>>
$dat[0], $right, $dat[1], $left;
.