#!/export/opt1/perl/bin/perl -w
#
# SPLITTER, version 2.0
#
#
# For command line input:
#
# a = archive temporary files
# e = eigenvalue method only
# i = STOP output images
# o = STOP ORACLE output
# s = Silver & Chan error methods only (STOP BOOTSTRAPPING);
# t = transverse method only

### CURRENT NOTES ###


use Oraperl;
use strict;

do '/export/home/matt/wfars/shared.pl';
do '/export/home/matt/wfars/stalta_picker.pl';



my (@filename);
my ($file, $newfile, $e_file, $n_file, $z_file, $one_file, $two_file, $all_files, $rf_out);
my ($reqid, $sta);
my (@temp, @sac_out, $line);
my ($hypid, $net);
my ($evdate, $evmsec, $evla, $evlo, $evdp);
my ($stla, $stlo);
my ($temp, $sac, $command);
my ($sks_start, $skks_start, $sks_pick);
my (@method, $method, $hz, $req_hz, $dec);
my ($pick, $sks_end);
my ($seaz, $delta, $sacdate);
my ($eig_or_ratio, $eig_cor_ratio, $tra_or_ratio, $tra_cor_ratio, $eig_win, $tra_win);

#my ($i, $j, $target, $size, $gmt);
#my (@radial, @transv, @epmon, @epmoe, @epmcn, @epmce, @tpmon, @tpmoe, @tpmcn, @tpmce);
#my (@eigbootang, @eigbootdt, @trabootang, @trabootdt);
#my (@index, @eigang, @eigdt, @traang, @tradt);

my ($eang, $eangmin, $eangmax, $edt, $edtmin, $edtmax, $tang, $tangmin, $tangmax, $tdt, $tdtmin, $tdtmax);

my ($num_boot, $ar_comp, $ar_total);
my ($a_flag, $i_flag, $s_flag, $o_flag);

my ($phid, $rdid, $pol_dir, $err_method);
my ($reporter, $request, $date, $name, $reqidmin, $reqidmax);

my ($eigscang, $eigscdt, $trascang, $trascdt);

my ($cut1, $cut2);

my ($net_name);

### Environment variables ###

$net_name = "CN";

#$ENV{QWFARS}      = "/export/home/matt/dummy_wfars";
$ENV{ARCHIVE}     = "/space/matt/archive";

$ENV{SPLIT}       = "MATT.".$net_name."_REQUEST";
$ENV{RESULT}      = "MATT.".$net_name."_RESULT";


$ENV{HYPOCENTER}  = "ISCDBADM.HYPOCENTER";
$ENV{SITE}        = "MATT.STATION";
$ENV{PHASE}       = "MATT.PHASE";
$ENV{ASSOCIATION} = "MATT.ASSOCIATION";
$ENV{REPORT}      = "MATT.REPORT";
$ENV{AUTHOR}      = "ULE_SKS";

$ENV{FOLDER}      = "$ENV{QINPUT}/$ENV{AUTHOR}";

#($ENV{LDA})       = oracle_login();
($ENV{TEMP_DIR})  = make_temp_dir();

### DETERMINE DEFAULT VALUES ###

$req_hz   = 10;

### Variables for STALTA picker ###



$num_boot = 200;

@method = qw(e t);
$a_flag = "0";
$i_flag = "1";
$s_flag = "0";
$o_flag = "1";
$ar_total = "";

### Perform command line tests ###

foreach $ar_comp (@ARGV) {

    if (substr($ar_comp,0,1) eq "-") {
        $ar_total .= $ar_comp;
    }

}

### Loop over each option ###

if ($ar_total =~ /e/) { @method = "e";print "Eigenvalue method only\n" }
if ($ar_total =~ /t/) { @method = "t";print "Transverse method only\n" }
if ($ar_total =~ /a/) { $a_flag = "1";print "Archive temporary files\n" }
if ($ar_total =~ /i/) { $i_flag = "0";print "Stop images\n" }
if ($ar_total =~ /s/) { $s_flag = "1";print "S&C only\n" }
if ($ar_total =~ /o/) { $o_flag = "0";print "No ORACLE output\n" }


### Find list of 'Z' files ###

@filename = `ls $ENV{QWFARS}/data/SPLIT/*Z`;

$reqidmin = "";
$reqidmax = "";

($ENV{LDA}) = oracle_login();

### Assign new reportor ID ###

if ($o_flag == 1) {

    $request = "SELECT iscdbadm.repid.NEXTVAL FROM DUAL";
    ($reporter) = oracle_request($request);
    #oracle_logout();
    $name = $reporter.".log";

    open (SUMMARY, ">$ENV{FOLDER}/$name");

} else {
    open (SUMMARY, ">./log.split");
} 

foreach $file (@filename) { 
    
    ### Undefine potential problem variables ###

    #undef $gmt;
    #undef @index;
    #undef @eigang;
    #undef @eigdt;
    #undef @traang;
    #undef @tradt;
    #undef @radial;
    #undef @transv;
    #undef @epmon;
    #undef @epmoe;
    #undef @epmcn;
    #undef @epmce;
    #undef @tpmon;
    #undef @tpmoe;
    #undef @tpmcn;
    #undef @tpmce;
    #undef @eigbootang;
    #undef @eigbootdt;
    #undef @trabootang;
    #undef @trabootdt;
    
    #($ENV{LDA}) = oracle_login();

    
    chomp($file);
    
    ### Find REQID, STA ###
    
    @temp = split(/\//,$file);
    @temp = split(/\./,$temp[-1]);
    
    $reqid = $temp[0];
    $sta   = $temp[1];
    
    print "Attempting $reqid at $sta\n";
    
    ### Sort outmin/max REQID values for log file ###
    
    if ($reqidmin eq "") {
    
        $reqidmin = $reqid;
        $reqidmax = $reqid;
        
    } else {
    
        if ($reqid < $reqidmin) { $reqidmin = $reqid; }
        if ($reqid > $reqidmax) { $reqidmax = $reqid; }
    
    }
    
    ### Find 'E' and 'N' files ###
    
    $z_file = $file;
    chop($file);
    $e_file    = $file."E";
    $n_file    = $file."N";
    $one_file  = $file."1";
    $two_file  = $file."2";
    $all_files = $file."*";
    
    $rf_out = 0;
    
    ### Check E and N component exist ###
    
    @temp = stat($n_file);
    if ($temp[0] ne "") { 
        ### N file found
        @temp = stat($e_file);
        if ($temp[0] ne "") { 
            ### N file found, E file found
            ### ALL OK, but check anyway ###
            ($rf_out) = rotate_files($n_file, $e_file, $n_file, $e_file);
        } else {
            ### N file found, E file NOT found
            ### Something weid, dump!
            print "$e_file not found, skipping\n";
            print SUMMARY "$e_file not found, skipping\n";
            remove_files($all_files);
            next;
        }
    } else {
        ### N file NOT found
        @temp = stat($one_file);
        if ($temp[0] ne "") { 
            ### 1 file found
            @temp = stat($two_file);
            if ($temp[0] ne "") { 
                ### 1 file found, 2 file found
                ### Rotate into N and E files
                print "Not North & East Components, correcting....\n";
                print SUMMARY "Not North & East Components, correcting....\n";

                ($rf_out) = rotate_files($one_file, $two_file, $n_file, $e_file);
            } else {
                ### 1 file found, 2 file NOT found
                ### Something weird, dump!
                print "$two_file not found, skipping\n";
                print SUMMARY "$two_file not found, skipping\n";
                remove_files($all_files);
                next;
            }
        } else {
        ### N file NOT found, 1 file NOT found
        print "Other components missing, skipping\n";
        print SUMMARY "Other components missing, skipping\n";
        remove_files($all_files);
        next;
        }
    }
    
    if ($rf_out eq "-1") {
        remove_files($all_files);
        next;
    }
    
    #@temp = stat($e_file);
    #if ($temp[0] eq "") { 
    #    remove_files($z_file, $n_file);
    #    print "$e_file not found, skipping\n";
    #    print SUMMARY "$e_file not found, skipping\n";
    #    next;
    #}
    
    #@temp = stat($n_file);
    #if ($temp[0] eq "") { 
    #    remove_files($z_file, $e_file);
    #    print "$n_file not found, skipping\n";
    #    print SUMMARY "$n_file not found, skipping\n";
    #    next;
    #}    
    
    ### Identify event/station details ###
    
    ($hypid, $net, $seaz) = identify_request_details($reqid, $sta);
    ($evdate, $evmsec, $evla, $evlo, $evdp) = identify_event_details($hypid);
    ($stla, $stlo) = identify_station_details($net, $sta);
    

    ($sacdate) = get_sac_time($evdate);
    
   
    ####################################    
    ###                              ###
    ### Change header values         ###
    ###   sta lat, lon,depth, height ###
    ###   CMPAZ, CMPINC              ###
    ###   EVLA, EVLO, EVDP           ###
    ###                              ###
    ####################################    
    
    change_sac_header($z_file, $e_file, $n_file, $stla, $stlo, 
        $evla, $evlo, $evdp, $sacdate, $evmsec);
    
    ### Archive files (make copies) ###
   
    archive_original_files($z_file, $e_file, $n_file, $net, $sta);
     
    ### Do STA/LTA pick process ###
    
    
    ### Find onset (using IASPEI temporarily) ###
    
    ($delta, $sks_start) = get_start_time($reqid, $sta, $evdp);
    
    $delta = sprintf("%6.2f",$delta);
    
    
    print "$sks_start\n";
    
    ### LOGOFF ORACLE (prevents timeout errors) ###
    
    #oracle_logout();
    
    $cut1 = $sks_start-60;
    $cut2 = $sks_start+60;

    $command =  "sac <<EOF\n";
    $command .= "cut $cut1 $cut2\n";
    $command .= "r $n_file $e_file\n";
    $command .= "rmean;rtrend\n";
    $command .= "w $n_file $e_file\n";
    $command .= "EOF\n";
    
    #print "$command\n";
    
    @sac_out = `$command`;
    
    #print "@sac_out\n";
    

    ($pick, $sks_pick) = stalta_picker($n_file, $e_file, $sks_start);
    
    
    
    ### Decide if measurement to be made ###
    
    if ($pick == 1) {
    
        if ($o_flag == 1) {
    
            #($ENV{LDA}) = oracle_login();
            ($phid, $rdid) = assign_phid($reqid, $sta, $net);
            #oracle_logout();
            
            print "Accepted -> phid: $phid\n";
        } else {
        
            $phid = 0;
            
        
        }
        
        
            
        $sks_start = $sks_pick - 5;
        $sks_end = $sks_start + 52;
            
            
        ### Filter files ###
    
        filter_files($z_file, $e_file, $n_file);
    
        ### Move to temporary folder ###
    
        foreach $file ($z_file, $e_file, $n_file) {
    
            @temp = split(/\//,$file);
            $newfile = "$ENV{TEMP_DIR}/$temp[-1]";
        
            `mv $file $newfile`;
        
            $file = $newfile;
    
        }
            
                
        ### Prepare N and E components ###
        
        $sac =  "r $n_file\n";
        $sac .= "lh delta";
        
        $command = "sac <<EOF\n$sac\nEOF\n";
        @sac_out = `$command`;
        
        foreach $line (@sac_out) {
            chomp ($line);
            if ($line =~ /delta/) {
                @temp = split(/ +/,$line);
                $hz   = 1/$temp[3];
                $dec  = $hz / $req_hz;
            }
        }
        
    
        $command =  "sac <<EOF\n";
        $command .= "cut $sks_start $sks_end\n";
        $command .= "r $n_file $e_file\n";
        $command .= "decimate $dec\n";
        $command .= "w $ENV{TEMP_DIR}/temp1 $ENV{TEMP_DIR}/temp2\n";
        $command .= "EOF\n";
        
        
        print "$command\n";
        
        @sac_out = `$command`;
        
        ### Perform measurements ###

        $command =  "/export/home/matt/autosplit/autosplit <<EOF\n";
        $command .= "$ENV{TEMP_DIR}/temp1\n";
        $command .= "$ENV{TEMP_DIR}/temp2\n";
        $command .= "$num_boot\n";
        $command .= "3\n";
        $command .= "1\n";
        $command .= "EOF\n";
        
        @sac_out = `$command`;
        
        open (TEMP, ">$ENV{TEMP_DIR}/dump.out");
        print TEMP "@sac_out\n\n";
        close TEMP;
        

        ($eang, $eangmin, $eangmax, $edt, $edtmin, $edtmax, $eig_win, $eig_or_ratio, $eig_cor_ratio, $pol_dir,
         $tang, $tangmin, $tangmax, $tdt, $tdtmin, $tdtmax, $tra_win, $tra_or_ratio, $tra_cor_ratio,
         $eigscang, $eigscdt, $trascang, $trascdt) = get_results(@sac_out);
        
        create_gifs($phid, $net, $sta, $evdate, $evla, $evlo, $evdp, $delta, $seaz, @sac_out);
        
        #
        ##
        ### Enter values into database here!!!
        ##
        #
        
        #$err_method = "BOOTS";
       
        if ($o_flag == 1) {
    
            #$ENV{LDA}) = oracle_login();
            
            values_to_tables($hypid, $phid, $rdid, $net, $sta, $delta, $seaz, $sks_start, $eang, $eangmin, $eangmax, $edt, $edtmin, $edtmax, $eig_win, $eig_or_ratio, $eig_cor_ratio, $pol_dir,
                $tang, $tangmin, $tangmax, $tdt, $tdtmin, $tdtmax, $tra_win, $tra_or_ratio, $tra_cor_ratio, $reporter, $eigscang, $eigscdt, $trascang, $trascdt);
  
            
            ora_commit($ENV{LDA});
            #oracle_logout();
            
        }
      

        
        printf "EIGEN %6s %5s %3.1f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n",$reqid, $sta, $eig_win, $eang, $eangmin, $eangmax, $edt, $edtmin, $edtmax, $eig_or_ratio, $eig_cor_ratio, $pol_dir;
        printf SUMMARY "EIGEN %6s %5s %3.1f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n",$reqid, $sta, $eig_win, $eang, $eangmin, $eangmax, $edt, $edtmin, $edtmax, $eig_or_ratio, $eig_cor_ratio, $pol_dir;
        
        printf "TRANS %6s %5s %3.1f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n",$reqid, $sta, $tra_win, $tang, $tangmin, $tangmax, $tdt, $tdtmin, $tdtmax, $tra_or_ratio, $tra_cor_ratio;
        printf SUMMARY "TRANS %6s %5s %3.1f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n",$reqid, $sta, $tra_win, $tang, $tangmin, $tangmax, $tdt, $tdtmin, $tdtmax, $tra_or_ratio, $tra_cor_ratio;
       
        print SUMMARY "SILVER & CHAN RESULTS: EIGEN $eang +/- $eigscang, $edt +/- $eigscdt\n";
        print SUMMARY "SILVER & CHAN RESULTS: TRANS $tang +/- $trascang, $tdt +/- $trascdt\n";

        print "SILVER & CHAN RESULTS: EIGEN $eang +/- $eigscang, $edt +/- $eigscdt\n";
        print "SILVER & CHAN RESULTS: TRANS $tang +/- $trascang, $tdt +/- $trascdt\n";
    
    } else { # related to 'if ($pick == 1)'

        
        print SUMMARY "$n_file -> no pick\n";
        print "$n_file -> no pick\n";
        
        remove_files($all_files);
        
    }   
    
    

    ### Clear temporary directory ###
            
    if ($a_flag != 1) { unlink <$ENV{TEMP_DIR}/*>; }

    if ($o_flag == 1) {
    
        #($ENV{LDA}) = oracle_login();
        
        $command = "UPDATE $ENV{SPLIT}
                    SET PROCESSED = 'P'
                    WHERE net = '$net'
                    AND sta = '$sta'
                    AND reqid = '$reqid'";
                
        oracle_command($command);
    
        #oracle_logout();
        
    }


} # End of @filename foreach loop


oracle_logout();

close SUMMARY;

#
##
### Do REPORTER stuff here ###
##
#

if ($o_flag == 1) {

    ($ENV{LDA}) = oracle_login();

    $temp = "REQIDs $reqidmin-$reqidmax";

    $command = "INSERT INTO $ENV{REPORT} (REPID, REPORTER, PRODUCT, REPORTER_ID, DIRNAME, FILENAME, COLLECTOR, LDDATE)
              VALUES ('$reporter', '$ENV{AUTHOR}', 'Shear splitting', '$temp', '$ENV{AUTHOR}', '$name', '$ENV{USER}', SYSDATE)";

    print "$command\n";

    oracle_command($command);

    oracle_logout();

}

send_mail("matt\@isc.ac.uk","SPLIT COMPLETE!");

exit (0);

#
##
#

#DOC#################################################
#DOC
#DOC subroutine archive_original_files
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Assign new phase ID
#DOC
#DOC input : $z_file - Vertical component filename
#DOC         $e_file - East component filename
#DOC         $n_file - North component filename
#DOC         $net - Network code
#DOC         $sta - Station code
#DOC
#DOC output : none
#DOC
#DOC internal : $file - Used to loop over each file
#DOC            $newfile - New filename
#DOC            @temp - Used in split command
#DOC
#DOC environment variables : $ENV{ARCHIVE} - Path to archive
#DOC
#DOC calls subroutines : none
#DOC
#DOC#################################################

sub archive_original_files {

    my ($z_file, $e_file, $n_file, $net, $sta) = @_;
    my ($file, $newfile, @temp);
    
    foreach $file ($z_file, $e_file, $n_file) {
        
        @temp = split(/\//,$file);
        $newfile = "$ENV{ARCHIVE}/$net/$sta/$temp[-1]";
        `cp $file $newfile`;
        chmod(0666,"$newfile");
        
    }
    
}

#DOC#################################################
#DOC
#DOC subroutine assign_phid
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Assign new phase ID
#DOC
#DOC input : $reqid - Request ID number
#DOC         $sta - Station code
#DOC
#DOC output : $phid - Phase ID
#DOC          $rdid - Reading ID
#DOC
#DOC internal : $request - SQLplus request
#DOC            $command - SQLplus command
#DOC
#DOC environment variables : $ENV{SPLIT} - Oracle table
#DOC
#DOC calls subroutines : oracle_request
#DOC                     oracle_command
#DOC
#DOC#################################################

sub assign_phid {

    my ($reqid, $sta) = @_;

    my ($phid, $rdid, $request, $command);

    $request = "SELECT iscdbadm.phid.NEXTVAL 
                FROM DUAL";
                
    ($phid) = oracle_request($request);
    
    $request = "SELECT iscdbadm.rdid.NEXTVAL 
                FROM DUAL";
                
    ($rdid) = oracle_request($request);
    
    $command = "UPDATE $ENV{SPLIT}
                SET sks_phid = '$phid'
                WHERE sta = '$sta'
                AND reqid = '$reqid'";
                
    oracle_command($command);
    
    return ($phid, $rdid);

}

#DOC#################################################
#DOC
#DOC subroutine change_sac_header
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Sets relevant SAC header variables (for archived files)
#DOC
#DOC input : $z_file - Vertical component file
#DOC         $e_file - East component file
#DOC         $n_file - North component file
#DOC         $stla - Station latitude (degrees)
#DOC         $stlo - Station longitude (degrees)
#DOC         $evla - Event latitude (degrees)
#DOC         $evlo - Event longitude (degrees)
#DOC         $evdp - Event depth (km)
#DOC         $evdate - Event date
#DOC         $evmsec - Milliseconds after event date
#DOC
#DOC output : none
#DOC
#DOC internal : $file - Loop over each file name
#DOC            $sac - SAC commands in SAC format
#DOC            $cmpaz - Component azimuth
#DOC            $cmpinc - Component incident angle
#DOC            $kcmpnm - Component name
#DOC            $command - Executable command
#DOC            @sac_out - Output from SAC
#DOC            %o_marker - Hash contains 'o' offset values (therefore can differ)
#DOC            @temp - Temporary array
#DOC            $temp - Key for '%o_marker'
#DOC            $line - Used to loop over '@sacout'
#DOC            $key - Used to loop over each key in '%o_marker'
#DOC
#DOC environment variables : none
#DOC
#DOC calls subroutines : sac (independant: /usr/local/sac10.6f/bin/sac)
#DOC
#DOC#################################################


sub change_sac_header {

    my ($z_file, $e_file, $n_file, $stla, $stlo, 
        $evla, $evlo, $evdp, $evdate, $evmsec) = @_;

    my ($file, $sac, $cmpaz, $cmpinc, $kcmpnm);
    my ($command, @sac_out, %o_marker, @temp, $temp, $line, $key);
    
    foreach $file ($z_file, $e_file, $n_file) {
        
        $temp = substr($file,-1,1);
        
        if ($temp eq "E") {
            $cmpaz  = 90;
            $cmpinc = 90;
            $kcmpnm = substr($file,-3,2)."E     ";
        } elsif ($temp eq "N") {
            $cmpaz  = 0;
            $cmpinc = 90;
            $kcmpnm = substr($file,-3,2)."N     ";
        } elsif ($temp eq "Z") {
            $cmpaz  = 0;
            $cmpinc = 0;
            $kcmpnm = substr($file,-3,2)."Z     ";
        }
        
        $sac .= "r $file\n";
        $sac .= "rmean; rtrend\n";
        $sac .= "w $file\n";
        $sac .= "r $file\n";
        #$sac .= "chnhdr CMPAZ $cmpaz CMPINC $cmpinc\n";
        $sac .= "chnhdr STLA $stla STLO $stlo\n";
        $sac .= "chnhdr EVLA $evla EVLO $evlo EVDP $evdp\n";
        $sac .= "chnhdr o GMT $evdate $evmsec\n";
        $sac .= "chnhdr IZTYPE 9 KCMPNM $kcmpnm\n";
        $sac .= "chnhdr LCALDA TRUE\n";

        $sac .= "wh\n";
        $sac .= "lh o\n";
        
    }
   
    $command = "sac <<EOF\n$sac\nEOF\n";
    @sac_out = `$command`;

    ### Use output to find o value and modify correctly ###
    
    undef %o_marker;
        
    foreach $line (@sac_out) {
    
        chomp($line);
        
        if ($line =~ /FILE: /) {
            @temp = split(/ +/,$line);
            $temp = $temp[2];
        } elsif ($line =~ /o = /) {
            @temp = split(/ +/,$line);
            $o_marker{$temp} = $temp[3];
        }
    
    }
    
    undef $sac;

    foreach $key (keys (%o_marker)) {
        $sac .= "r $key\n";
        if ($o_marker{$key} < 0) {
            $o_marker{$key} = $o_marker{$key} * -1;
            $sac .= "chnhdr allt $o_marker{$key}\n";
        } else {
            $sac .= "chnhdr allt -$o_marker{$key}\n";
        }
        
        $sac .= "w over\n";

    }
    
    $command = "sac <<EOF\n$sac\nEOF\n";
    @sac_out = `$command`;
    
}

#DOC#################################################
#DOC
#DOC subroutine create_gifs
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Make gif files, (using convert and gmt)
#DOC
#DOC input : $reqid - Request ID number
#DOC         $net - Network code
#DOC         $sta - Station code
#DOC         $evdate - Event date
#DOC         $evla - Event latitude
#DOC         $evlo - Event longitude
#DOC         $evdp - Event depth (km)
#DOC         $delta - Angular distance (degrees)
#DOC         $seaz - Station to event azimuth (degrees)
#DOC         @sac_out - STDOUT from splitting program
#DOC
#DOC output : none
#DOC
#DOC internal : @index - Array containing window lengths
#DOC            @eigang - Array containing eigenvalue angle results
#DOC            @eigdt - Array containing eigenvalue lag times
#DOC            @traang - Array containing transverse energy angle results
#DOC            @tradt - Array contianing transverse energy lag times
#DOC            $eig_or_ratio - Ratio of eigenvalues in original traces, (eigenvalue method)
#DOC            $eig_cor_ratio - Ratio of eigenvalues in corrected traces, (eigenvalue method)
#DOC            $tra_or_ratio - Ratio of eigenvalues in original traces, (transverse method)
#DOC            $tra_cor_ratio - Ratio of eigenvalues in corrected traces, (transverse method)
#DOC            $eig_win - Eigenvalue window length
#DOC            $tra_win - Transverse energy window length
#DOC            $i - General loop variable
#DOC            @radial - Radial waveform component
#DOC            @transv - Transverse waveform component
#DOC            @epmon - Eigenvalue method original particle motion (north component)
#DOC            @epmoe - Eigenvalue method original particle motion (east component)
#DOC            @epmcn - Eigenvalue method original corrected motion (north component)
#DOC            @epmce - Eigenvalue method original corrected motion (east component)
#DOC            @tpmon - Transverse energy method original particle motion (north component)
#DOC            @tpmoe - Transverse energy method original particle motion (east component)
#DOC            @tpmcn - Transverse energy method original corrected motion (north component)
#DOC            @tpmce - Transverse energy method original corrected motion (east component)
#DOC            $eang - Eigenvalue angle result
#DOC            $eange - Eigenvalue angle error
#DOC            $edt - Eigenvalue lag result
#DOC            $edte - Eigenvalue lag error
#DOC            $tang - Transverse energy angle result
#DOC            $tange - Transverse energy angle error
#DOC            $tdt - Transverse energy lag result
#DOC            $tdte - Transverse energy lag error
#DOC            $line - Used to loop over '@sac_out'
#DOC            $gmt - Input passed into 'gmt'
#DOC            $size - Size of waveforms in number of points, (used to create size of window in gmt)
#DOC            $size2 - Size of waveforms in seconds, (used to create size of window in gmt)
#DOC            $file - Output filename, (gmt)
#DOC            $target - GIF filename for output from 'convert'
#DOC            @eigbootang - Eigenvalue bootstrap angle results
#DOC            @eigbootdt - Eigenvalue bootstrap lag time results
#DOC            @trabootang - Transverse energy bootstrap angle results
#DOC            @trabootdt - Transverse energy bootstrap lag time results
#DOC            $temp - Temporary value
#DOC
#DOC environment variables : $ENV{TEMP_DIR} - Temporary directory
#DOC
#DOC calls subroutines : gmt (independent: /distrib/gmt3.4/bin/<file e.g. 'psxy'>)
#DOC                     convert (independent: /usr/local/bin/convert)
#DOC
#DOC#################################################

sub create_gifs {
    
    my ($phid, $net, $sta, $evdate, $evla, $evlo, $evdp, $delta, $seaz, @sac_out) = @_;
    my (@index, @eigang, @eigdt, @traang, @tradt);
    my ($eig_or_ratio, $eig_cor_ratio, $tra_or_ratio, $tra_cor_ratio, $eig_win, $tra_win);
    my ($i);
    my (@radial, @transv, @epmon, @epmoe, @epmcn, @epmce, @tpmon, @tpmoe, @tpmcn, @tpmce);
    my ($eang, $eangmin, $eangmax, $edt, $edtmin, $edtmax, $tang, $tangmin, $tangmax, $tdt, $tdtmin, $tdtmax);
    my ($line, $gmt, $size, $size2, $file, $target);
    my (@eigbootang, @eigbootdt, @trabootang, @trabootdt);
           
    foreach $line (@sac_out) {
    
        chomp ($line);
        
        if ($line =~ /Seismograms/) {
            @temp = split(/ +/,$line);
            push (@radial, $temp[2]);
            push (@transv, $temp[3]);
        } elsif ($line =~ /Window length/) {
            @temp = split(/ +/,$line);
            push (@index, $temp[3]);
            push (@eigang, $temp[4]);
            push (@eigdt, $temp[5]);
            push (@traang, $temp[6]);
            push (@tradt, $temp[7]);
        } elsif ($line =~ /eigenvalue original ratio/) {
            @temp = split(/ +/,$line);
            $eig_or_ratio = $temp[4];
        } elsif ($line =~ /eigenvalue corrected ratio/) {
            @temp = split(/ +/,$line);
            $eig_cor_ratio = $temp[4];
        } elsif ($line =~ /Eigenvalue window/) {
            @temp = split(/ +/,$line);
            $eig_win = $temp[3];
        } elsif ($line =~ /transverse original ratio/) {
            @temp = split(/ +/,$line);
            $tra_or_ratio = $temp[4];
        } elsif ($line =~ /transverse corrected ratio/) {
            @temp = split(/ +/,$line);
            $tra_cor_ratio = $temp[4];
        } elsif ($line =~ /Transverse window/) {
            @temp = split(/ +/,$line);
            $tra_win = $temp[3];
        } elsif ($line =~ /EIGENVALUE:/) {
            @temp  = split(/ +/,$line);
            $eang  = $temp[4];
            #$eange = $temp[6];
            $edt   = $temp[9];
            #$edte  = $temp[11];
        } elsif ($line =~ /TRANSVERSE:/) {
            @temp = split(/ +/,$line);
            $tang  = $temp[4];
            #$tange = $temp[6];
            $tdt   = $temp[9];
            #$tdte  = $temp[11];
        } elsif ($line =~ /PM PLOTS  1/) {
            @temp = split(/ +/,$line);
            push (@epmon, $temp[4]);
            push (@epmoe, $temp[5]);
            push (@epmcn, $temp[6]);
            push (@epmce, $temp[7]);
        } elsif ($line =~ /PM PLOTS  2/) {
            @temp = split(/ +/,$line);
            push (@tpmon, $temp[4]);
            push (@tpmoe, $temp[5]);
            push (@tpmcn, $temp[6]);
            push (@tpmce, $temp[7]);
        } elsif ($line =~ /Bootstrap values  1/) {
            @temp = split(/ +/,$line);
            push (@eigbootang, $temp[4]);
            push (@eigbootdt, $temp[5]);
        } elsif ($line =~ /Bootstrap values  2/) {
            @temp = split(/ +/,$line);
            push (@trabootang, $temp[4]);
            push (@trabootdt, $temp[5]);
        } elsif ($line =~ /EIGENVALUE ANGLE PERCENTILES/) {
            @temp = split(/ +/,$line);
            $eangmin = $temp[4];
            $eangmax = $temp[5];
        } elsif ($line =~ /EIGENVALUE LAG PERCENTILES/) {
            @temp = split(/ +/,$line);
            $edtmin = $temp[4];
            $edtmax = $temp[5];
        } elsif ($line =~ /TRANSVERSE ANGLE PERCENTILES/) {
            @temp = split(/ +/,$line);
            $tangmin = $temp[4];
            $tangmax = $temp[5];
        } elsif ($line =~ /TRANSVERSE LAG PERCENTILES/) {
            @temp = split(/ +/,$line);
            $tdtmin = $temp[4];
            $tdtmax = $temp[5];
        }
        
    }
       
    $size  = @radial;
    $size2 = $size/10;
    
    
    ### Create GMT output
    $gmt =  "/export/opt1/gmt3.4/bin/gmtset ";
    $gmt .= "ANOT_MIN_ANGLE		  20 ";
    $gmt .= "ANOT_MIN_SPACING	  0 ";
    $gmt .= "ANOT_FONT		  Helvetica ";
    $gmt .= "ANOT_FONT_SIZE		  11p ";
    $gmt .= "ANOT_OFFSET		  0.2c ";
    $gmt .= "BASEMAP_AXES		  WESN ";
    $gmt .= "BASEMAP_FRAME_RGB	  0/0/0 ";
    $gmt .= "BASEMAP_TYPE		  fancy ";
    $gmt .= "COLOR_BACKGROUND	  0/0/0 ";
    $gmt .= "COLOR_FOREGROUND	  255/255/255 ";
    $gmt .= "COLOR_NAN		  128/128/128 ";
    $gmt .= "COLOR_IMAGE		  adobe ";
    $gmt .= "COLOR_MODEL		  rgb ";
    $gmt .= "D_FORMAT		  %lg ";
    $gmt .= "DEGREE_FORMAT		  0 ";
    $gmt .= "DOTS_PR_INCH		  300 ";
    $gmt .= "ELLIPSOID		  WGS-84 ";
    $gmt .= "FRAME_PEN		  1.25p ";
    $gmt .= "FRAME_WIDTH		  0.2c ";
    $gmt .= "GLOBAL_X_SCALE		  1 ";
    $gmt .= "GLOBAL_Y_SCALE		  1 ";
    $gmt .= "GRID_CROSS_SIZE	  0c ";
    $gmt .= "GRID_PEN		  0.25p ";
    $gmt .= "GRIDFILE_SHORTHAND	  FALSE ";
    $gmt .= "HEADER_FONT		  Helvetica ";
    $gmt .= "HEADER_FONT_SIZE	  36p ";
    $gmt .= "HSV_MIN_SATURATION	  1 ";
    $gmt .= "HSV_MAX_SATURATION	  0.1 ";
    $gmt .= "HSV_MIN_VALUE		  0.3 ";
    $gmt .= "HSV_MAX_VALUE		  1 ";
    $gmt .= "INTERPOLANT		  akima ";
    $gmt .= "IO_HEADER		  FALSE ";
    $gmt .= "N_HEADER_RECS		  1 ";
    $gmt .= "LABEL_FONT		  Helvetica ";
    $gmt .= "LABEL_FONT_SIZE	  12p ";
    $gmt .= "LINE_STEP		  0.025c ";
    $gmt .= "MAP_SCALE_FACTOR	  0.9996 ";
    $gmt .= "MAP_SCALE_HEIGHT	  0.2c ";
    $gmt .= "MEASURE_UNIT		  cm ";
    $gmt .= "N_COPIES		  1 ";
    $gmt .= "OBLIQUE_ANOTATION	  1 ";
    $gmt .= "PAGE_COLOR		  255/255/255 ";
    $gmt .= "PAGE_ORIENTATION	  landscape ";
    $gmt .= "PAPER_MEDIA		  a4 ";
    $gmt .= "PSIMAGE_FORMAT		  hex ";
    $gmt .= "TICK_LENGTH		  0.2c ";
    $gmt .= "TICK_PEN		  0.5p ";
    $gmt .= "UNIX_TIME		  FALSE ";
    $gmt .= "UNIX_TIME_POS		  -2c/-2c ";
    $gmt .= "VECTOR_SHAPE		  0 ";
    $gmt .= "VERBOSE		  FALSE ";
    $gmt .= "WANT_EURO_FONT		  TRUE ";
    $gmt .= "X_AXIS_LENGTH		  25c ";
    $gmt .= "Y_AXIS_LENGTH		  15c ";
    $gmt .= "X_ORIGIN		  2.5c ";
    $gmt .= "Y_ORIGIN		  2.5c ";
    $gmt .= "XY_TOGGLE		  FALSE ";
    $gmt .= "Y_AXIS_TYPE		  hor_text \n";    
        
    
        
    ### Eigenvalue page ###

    $file = "$ENV{TEMP_DIR}/eigen.ps";
    
    ### Eigenvalue corrected PM
        
    $gmt .= "psbasemap -R-1/1/1/3 -JX3/3 -Bf10a100::/f1a4:\"Corrected PM\":Wsen -K -X2 -Y2 > $file\n";
    $gmt .= "psxy -R-1/1/-1/1 -JX3/3 -O -K -A << END >> $file\n";
    $i=0;
    foreach (@epmce) { $gmt .= "$epmce[$i] $epmcn[$i]\n"; $i++ }
    $gmt .= "END\n";
        
    ### Eigenvalue original PM
                
    $gmt .= "psbasemap -R-1/1/1/3 -JX3/3 -Bf10a100::/f1a4:\"Original PM\":Wsen -O -K -Y3 >> $file\n";
    $gmt .= "psxy -R-1/1/-1/1 -JX3/3 -O -K  -A << END >> $file\n";
    $i=0;
    foreach (@epmoe) { $gmt .= "$epmoe[$i] $epmon[$i]\n"; $i++ }
    $gmt .= "END\n";
        
    ### Eigenvalue results text
 
    $gmt .= "pstext -R120/410/0/1 -JX10/3 -O -K -N -P -X4 -Y-2 << END >> $file\n";
        
    @temp = ("Window length = $eig_win", "Angle : $eang", "Range: $eangmin -> $eangmax", "Lag : $edt","Range: $edtmin -> $edtmax",
             "Original ratio : $eig_or_ratio", "Corrected ratio : $eig_cor_ratio");
    $i=1.4;
    foreach (@temp) { $gmt .= "110 $i 12 0 4 LT $_\n"; $i -=0.2 }
    $gmt .= "END\n";
        
    ### Eigenvalue angle bootstrap histogram
        
    $gmt .= "psbasemap -R-90/90/10/60 -JX5/6 -Bf15a45:\"Angle (\217)\":/f10a100:\"Bootstrapping Results\":WSen  -O -K -X7 -Y-1 >> $file\n";
    $gmt .= "pshistogram -R-90/90/0/50 -JX5/6 -O -K -C -W1 -G128/128/128 << END >> $file\n";
    foreach (@eigbootang) { 
        if ($_ > 90) {$_ -= 180;} 
        elsif ($_ < -90) {$_ += 180;} 
        $gmt .= "$_\n"; }
    $gmt .= "END\n";
    
    ### Eigenvalue time bootstrap histogram
    
    $gmt .= "psbasemap -R-12/12/0/50 -JX5/6 -Bf1a4:\"Lag (s)\":/f10a100::wSen -O -K -X6 >> $file\n";
    $gmt .= "pshistogram -R-12/12/0/50 -JX5/6 -O -K -C -W0.1 -G128/128/128 << END >> $file\n";
    foreach (@eigbootdt) { $gmt .= "$_\n"; }
    $gmt .= "END\n";
    
    ### DIVIDER ####

    $gmt .= "psxy -R0/$size/0/1 -JX22/1 -O -K -W6/0/0/0  -A -X-17 -Y6 << END >> $file\n";
    $gmt .= "0 0.5\n";
    $gmt .= "520 0.5\n";
    $gmt .= "END\n";

    ### Eigenvalue dt results
        
    $gmt .= "psbasemap -R0/$size2/0/12 -JX22/2 -Bf1a10:\"Window length (s)\":/f1a4:\"Lag (s)\":WSen -O -K -Y2.5 >> $file\n";
    $gmt .= "psxy -R0/$size/0/12 -JX22/2 -O -K -Sc0.15c -G255/0/0  -A << END >> $file\n";
    $i=0;
    foreach (@index) { $gmt .= "$index[$i] $eigdt[$i]\n"; $i++; }
    $gmt .= "END\n";
    
    ### Eigenvalue angle results
    
    $gmt .= "psbasemap -R0/$size/-90/90  -JX22/2 -Bf10a100::/f15a90:\"Angle (\217)\":Wsen -O -K -Y2.5 >> $file\n";
    $gmt .= "psxy -R0/$size/-90/90 -JX22/2 -O -K -Sc0.15c -G255/0/0  -A << END >> $file\n";
    $i=0;
    foreach (@index) { $gmt .= "$index[$i] $eigang[$i]\n"; $i++; }
    $gmt .= "END\n";
    
    ### Eigenvalue window selection
    
    $temp = $eig_win * 10;
    
    $gmt .= "psxy -R0/$size/-90/90 -JX22/2 -O -K -Sc0.05c -N -G0/0/0  -A << END >> $file\n";
    for ($i=-315;$i<90;$i+=1) { $gmt .= "$temp $i\n"; }
    $gmt .= "END\n";
    
    ### Eigenvalue transverse component
    
    $gmt .= "psbasemap -R0/$size/1/3 -JX22/2 -Bf10a100:\"Seismogram\":/f4a4:\"Trans\":Wsen -O -K -Y2 >> $file\n";
    $gmt .= "psxy -R0/$size/-1/1 -JX22/2 -O -K  -A << END >> $file\n";
    $i=0;
    foreach (@transv) { $gmt .= "$i $transv[$i]\n"; $i++; }
    $gmt .= "END\n";
    
    ### Eigenvalue radial component
        
    $gmt .= "psbasemap -R0/$size/1/3 -JX22/2 -Bf10a100:\"Seismogram\":/f4a4:\"Radial\":Wsen -O -K -Y2 >> $file\n";
    $gmt .= "psxy -R0/$size/-1/1 -JX22/2 -O -K -A << END >> $file\n";
    $i=0;
    foreach (@radial) { $gmt .= "$i $radial[$i]\n"; $i++; }
    $gmt .= "END\n";
       
    ### Eigenvalue heading text
    
    $gmt .= "pstext -R120/410/0/1 -JX22/2 -O -N -P -Y1.5 << END >> $file\n";
        
    @temp = ("Eigenvalue method -    Station : $sta    Network : $net    Phase ID: $phid",
             "Event : $evdate    Latitude : $evla    Longitude :  $evlo    Depth : $evdp km",
             "Distance : $delta    Back Azimuth : $seaz");
    $i=1;
    foreach (@temp) { $gmt .= "130 $i 12 0 4 LT $_\n"; $i -=0.2; }
    $gmt .= "END\n";
    
    ### Transverse page
    
    $file = "$ENV{TEMP_DIR}/trans.ps";
    
    ### Transverse corrected PM
        
    $gmt .= "psbasemap -R-1/1/1/3 -JX3/3 -Bf10a100::/f1a4:\"Corrected PM\":Wsen -K -X2 -Y2 > $file\n";
    $gmt .= "psxy -R-1/1/-1/1 -JX3/3 -O -K -A << END >> $file\n";
    $i=0;
    foreach (@tpmce) { $gmt .= "$tpmce[$i] $tpmcn[$i]\n"; $i++; }
    $gmt .= "END\n";
    
    ### Transverse original PM
                
    $gmt .= "psbasemap -R-1/1/1/3 -JX3/3 -Bf10a100::/f1a4:\"Original PM\":Wsen -O -K -Y3 >> $file\n";
    $gmt .= "psxy -R-1/1/-1/1 -JX3/3 -O -K  -A << END >> $file\n";
    $i=0;
    foreach (@tpmoe) { $gmt .= "$tpmoe[$i] $tpmon[$i]\n"; $i++; }
    $gmt .= "END\n";
    
    ### Transverse results text
    
    $gmt .= "pstext -R120/410/0/1 -JX10/3 -O -K -N -P -X4 -Y-2 << END >> $file\n";
    
    @temp = ("Window length = $tra_win", "Angle : $tang", "Range: $tangmin -> $tangmax", "Lag : $tdt","Range: $tdtmin -> $tdtmax",
             "Original ratio : $tra_or_ratio", "Corrected ratio : $tra_cor_ratio");
    $i=1.4;
    foreach (@temp) { $gmt .= "110 $i 12 0 4 LT $_\n"; $i -=0.2; }
    $gmt .= "END\n";
    
    ### Transverse angle bootstrap histogram
    
    $gmt .= "psbasemap -R-90/90/10/60 -JX5/6 -Bf15a45:\"Angle (\217)\":/f10a100:\"Bootstrapping Results\":WSen -O -K -X7 -Y-1 >> $file\n";
    $gmt .= "pshistogram -R-90/90/0/50 -JX5/6 -O -K -C -W1 -G128/128/128 << END >> $file\n";
    foreach (@trabootang) { 
        if ($_ > 90) {$_ -= 180;} 
        elsif ($_ < -90) {$_ += 180;} 
        $gmt .= "$_\n"; }
    $gmt .= "END\n";
    
    ### Transverse dt bootstrap histogram
    
    $gmt .= "psbasemap -R-12/12/0/50 -JX5/6 -Bf1a4:\"Lag (s)\":/f10a100::wSen -O -K -X6 >> $file\n";
    $gmt .= "pshistogram -R-12/12/0/50 -JX5/6 -O -K -C -W0.1 -G128/128/128 << END >> $file\n";
    foreach (@trabootdt) { $gmt .= "$_\n"; }
    $gmt .= "END\n";
    
    ### DIVIDER ####

    $gmt .= "psxy -R0/$size/0/1 -JX22/1 -O -K -W6/0/0/0  -A -X-17 -Y6 << END >> $file\n";
    $gmt .= "0 0.5\n";
    $gmt .= "520 0.5\n";
    $gmt .= "END\n";
    
    ### Transverse dt results
    
    $gmt .= "psbasemap -R0/$size2/0/12 -JX22/2 -Bf1a10:\"Window length (s)\":/f1a4:\"Lag (s)\":WSen -O -K -Y2.5  >> $file\n";
    $gmt .= "psxy -R0/$size/0/12 -JX22/2 -O -K -Sc0.15c -G255/0/0  -A << END >> $file\n";
    $i=0;
    foreach (@index) { $gmt .= "$index[$i] $tradt[$i]\n"; $i++; }
    $gmt .= "END\n";
    
    ### Transverse angle results
    
    $gmt .= "psbasemap -R0/$size/-90/90  -JX22/2 -Bf10a100::/f15a90:\"Angle (\217)\":Wsen -O -K -Y2.5 >> $file\n";
    $gmt .= "psxy -R0/$size/-90/90 -JX22/2 -O -K -Sc0.15c -G255/0/0  -A << END >> $file\n";
    $i=0;
    foreach (@index) { $gmt .= "$index[$i] $traang[$i]\n"; $i++; }
    $gmt .= "END\n";
        
    ### Transverse window selection
    
    $temp = $tra_win * 10;
        
    $gmt .= "psxy -R0/$size/-90/90 -JX22/2 -O -K -Sc0.05c -N -G0/0/0  -A << END >> $file\n";
    for ($i=-315;$i<90;$i+=1) { $gmt .= "$temp $i\n"; } $gmt .= "END\n";
    
    ### Transverse transverse component ###
    
    $gmt .= "psbasemap -R0/$size/1/3 -JX22/2 -Bf10a100::/f4a4:\"Trans\":Wsen -O -K -Y2 >> $file\n";
    $gmt .= "psxy -R0/$size/-1/1 -JX22/2 -O -K  -A << END >> $file\n";
    $i=0;
    foreach (@transv) { $gmt .= "$i $transv[$i]\n"; $i++; }
    $gmt .= "END\n";
        
    ### Transverse radial component ###
    
    $gmt .= "psbasemap -R0/$size/1/3 -JX22/2 -Bf10a100::/f4a4:\"Radial\":Wsen -O -K -Y2 >> $file\n";
    $gmt .= "psxy -R0/$size/-1/1 -JX22/2 -O -K -A << END >> $file\n";
    $i=0;
    foreach (@radial) { $gmt .= "$i $radial[$i]\n"; $i++; }
    $gmt .= "END\n";
    
    ### Transverse heading text
    
    $gmt .= "pstext -R120/410/0/1 -JX21/2 -O -N -P -Y1.5 << END >> $file\n";
    @temp = ("Transverse energy method -    Station : $sta    Network : $net    Phase ID: $phid",
             "Event : $evdate    Latitude : $evla    Longitude :  $evlo    Depth : $evdp km",
             "Distance : $delta    Back Azimuth : $seaz");
    $i=1;
    foreach (@temp) { $gmt .= "130 $i 12 0 4 LT $_\n"; $i -=0.2; }
    $gmt .= "END\n";
    
    ### Execute GMT script
    
    ### Archive file
    
    open (TEMPGMT, ">$ENV{TEMP_DIR}/gmt.pm");
    print TEMPGMT "$gmt";
    close TEMPGMT;
    
    chmod(0774,"$ENV{TEMP_DIR}/gmt.pm");
    
    system ("cd $ENV{TEMP_DIR};gmt.pm");

    
    
    ### Convert to gif and move ###
    
    $file = "$ENV{TEMP_DIR}/eigen.ps";
    $target = "/export/home/matt/public_html/results/new/$net/$sta/$phid.eigen.gif";
    
    $command = "/export/opt1/im/bin/convert $file -rotate 90 $target";
    
    `$command`;
    
    $file = "$ENV{TEMP_DIR}/trans.ps";
    $target = "/export/home/matt/public_html/results/new/$net/$sta/$phid.trans.gif";
    
    $command = "/export/opt1/im/bin/convert $file -rotate 90 $target";
    
    `$command`;

}

#DOC#################################################
#DOC
#DOC subroutine filter_files
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Filter SAC files
#DOC
#DOC input : @files - List of files to be filtered
#DOC
#DOC output : none
#DOC
#DOC internal : $sac - SAC commands in SAC format
#DOC            $command - Executable command
#DOC            @sac_out - Output from SAC
#DOC            $file - Loop variable
#DOC
#DOC environment variables : none
#DOC
#DOC calls subroutines : sac (independant: /usr/local/sac10.6f/bin/sac)
#DOC
#DOC#################################################

sub filter_files {

    my (@files) = @_;
    my ($sac, $command, @sac_out, $file);
    
    foreach $file (@files) {
    
        $sac .= "r $file\n";
        $sac .= "rmean; rtrend\n";
        $sac .= "bp bu co 0.01 0.3 np 2 p 2\n";
        $sac .= "w $file\n";
    
    }

    $command = "sac <<EOF\n$sac\nEOF\n";
    @sac_out = `$command`;
    
}

#DOC#################################################
#DOC
#DOC subroutine find_phase_delay
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Identify phase delay from 'ttimes' output (WILL NOT WORK FOR 'P')
#DOC
#DOC input : $phase - Phase to be found (e.g. "SKSac")
#DOC         @data - Output from 'ttimes'
#DOC
#DOC output : $delay - Time between event origin and predicted arrival in seconds
#DOC
#DOC internal : @temp - Temporary array used in split
#DOC            $i - Loop variable
#DOC            $size - Size of '@data'
#DOC
#DOC environment variables : none
#DOC
#DOC calls subroutines : none
#DOC
#DOC#################################################

sub find_phase_delay {

    my ($phase, @data) = @_;
    my (@temp, $i, $size, $delay);
    
    $size = @data;
    
    
    for ($i=28;$i<$size;$i++) {
    
        if ($data[$i] =~ $phase) {
        
            @temp  = split(/ +/,$data[$i]);
            $delay = $temp[3];
            
            last;
        
        }
    
    
    }
    
    return $delay;    

}

#DOC#################################################
#DOC
#DOC subroutine get_results
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Identify splitting results from output
#DOC
#DOC input : @input - STDOUT from splitting program
#DOC
#DOC output : $eang - Eigenvalue fast polarization direction (degrees from N)
#DOC          $eange - Associated error
#DOC          $edt - Eigenvalue lag time (seconds)
#DOC          $edte - Associated error
#DOC          $eig_win - Length of window (in data points NOT seconds)
#DOC          $eig_or_ratio - Ratio of eigenvalues in original traces
#DOC          $eig_cor_ratio - Ratio of eigenvalues in rotated and time lagged traces
#DOC          $pol_dir - Polorization direction
#DOC          $tang - Transverse energy fast polarization direction (degrees from N)
#DOC          $tange - Associated error
#DOC          $tdt - Transverse error lag time (seconds)
#DOC          $tdte - Associated error
#DOC          $tra_win - Length of windo (in data points NOT seconds)
#DOC          $tra_or_ratio - Ratio of eigenvalues in original traces
#DOC          $tra_cor_ratio - Ratio of eigenvalues in rotated and time lagged traces
#DOC
#DOC internal : @temp - Temporary array used in split
#DOC            $line - Used to loop over each line in '@input'
#DOC
#DOC environment variables : none
#DOC
#DOC calls subroutines : none
#DOC
#DOC#################################################


sub get_results {

    my (@input) = @_;
    my ($eang, $eangmin, $eangmax, $edt, $edtmin, $edtmax, $eig_win, $eig_or_ratio, $eig_cor_ratio, $pol_dir);
    my ($tang, $tangmin, $tangmax, $tdt, $tdtmin, $tdtmax, $tra_win, $tra_or_ratio, $tra_cor_ratio);
    my ($line, @temp);
    my ($eigscang, $eigscdt, $trascang, $trascdt);
    
    foreach $line (@input) {
    
        #print "$line";
    
        chomp ($line);
        
        if ($line =~ /eigenvalue original ratio/) {
            @temp = split(/ +/,$line);
            $eig_or_ratio = $temp[4];
        } elsif ($line =~ /eigenvalue corrected ratio/) {
            @temp = split(/ +/,$line);
            $eig_cor_ratio = $temp[4];
        } elsif ($line =~ /Eigenvalue window/) {
            @temp = split(/ +/,$line);
            $eig_win = $temp[3];
        } elsif ($line =~ /transverse original ratio/) {
            @temp = split(/ +/,$line);
            $tra_or_ratio = $temp[4];
        } elsif ($line =~ /transverse corrected ratio/) {
            @temp = split(/ +/,$line);
            $tra_cor_ratio = $temp[4];
        } elsif ($line =~ /Transverse window/) {
            @temp = split(/ +/,$line);
            $tra_win = $temp[3];
        } elsif ($line =~ /EIGENVALUE:/) {
            @temp  = split(/ +/,$line);
            $eang  = $temp[4];
            #$eange = $temp[6];
            $edt   = $temp[9];
            #$edte  = $temp[11];
        } elsif ($line =~ /TRANSVERSE:/) {
            @temp = split(/ +/,$line);
            $tang  = $temp[4];
            #$tange = $temp[6];
            $tdt   = $temp[9];
            #$tdte  = $temp[11];
        } elsif ($line =~ /EIGENVALUE ANGLE PERCENTILES/) {
            @temp = split(/ +/,$line);
            $eangmin = $temp[4];
            $eangmax = $temp[5];
        } elsif ($line =~ /EIGENVALUE LAG PERCENTILES/) {
            @temp = split(/ +/,$line);
            $edtmin = $temp[4];
            $edtmax = $temp[5];
        } elsif ($line =~ /TRANSVERSE ANGLE PERCENTILES/) {
            @temp = split(/ +/,$line);
            $tangmin = $temp[4];
            $tangmax = $temp[5];
        } elsif ($line =~ /TRANSVERSE LAG PERCENTILES/) {
            @temp = split(/ +/,$line);
            $tdtmin = $temp[4];
            $tdtmax = $temp[5];
        } elsif ($line =~ /POL. DIRECTION/) {
            @temp = split(/ +/,$line);
            $pol_dir = $temp[3];
        } elsif ($line =~ /Eigenvalue S&C results/) {
            @temp = split(/ +/,$line);
            $eigscang = $temp[7];
            $eigscdt  = $temp[10];
        } elsif ($line =~ /Transverse S&C results/) {
            @temp = split(/ +/,$line);
            $trascang = $temp[7];
            $trascdt  = $temp[10];
        }
            
    }
    
    return ($eang, $eangmin, $eangmax, $edt, $edtmin, $edtmax, $eig_win, $eig_or_ratio, $eig_cor_ratio, $pol_dir,
            $tang, $tangmin, $tangmax, $tdt, $tdtmin, $tdtmax, $tra_win, $tra_or_ratio, $tra_cor_ratio,
            $eigscang, $eigscdt, $trascang, $trascdt);


}

#DOC#################################################
#DOC
#DOC subroutine get_sac_time
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Converts date from "DD-MM-YYYY HH24:MI:SS" format to "YYYY DDD HH24 MI SS"
#DOC
#DOC input : $evdate - Event date (in "DD-MM-YYYY HH24:MI:SS" format)
#DOC
#DOC output : $evdate - Event date (in "YYYY DDD HH24 MI SS" format)
#DOC
#DOC internal : $request - SQLplus request
#DOC
#DOC environment variables : $ENV{SPLIT} - Oracle table
#DOC
#DOC calls subroutines : oracle_request (in 'shared.pl')
#DOC
#DOC#################################################

sub get_sac_time {

    my ($evdate) = @_;
    my ($request);
    
    $request = "SELECT to_char(to_date('$evdate','DD-MM-YYYY HH24:MI:SS') , 'YYYY DDD HH24 MI SS')
                FROM dual";
    
    ($evdate) = oracle_request($request);

    return ($evdate);

}

#DOC#################################################
#DOC
#DOC subroutine get_start_time
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Obtain predicted phase onset time
#DOC
#DOC input : $reqid - Request ID number
#DOC         $sta - Station code
#DOC         $evdp - Event depth (km)
#DOC
#DOC output : $delta - Angular distance (degrees)
#DOC          $sks_delay - Arrival time in seconds after origin time
#DOC
#DOC internal : $command - Executable command
#DOC            @data - Output from executable command
#DOC            $request - SQLplus request
#DOC
#DOC environment variables : $ENV{SPLIT} - Oracle table
#DOC
#DOC calls subroutines : find_phase_delay (local)
#DOC                     oracle_request (in 'shared.pl')
#DOC                     ttimes (independant: /export/dev/ops/bin/ttimes)
#DOC
#DOC#################################################

sub get_start_time {

    my ($reqid, $sta, $evdp) = @_;
    my ($command, $delta, @data, $sks_delay, $request);
    
    $request = "SELECT delta
                FROM $ENV{SPLIT} 
                WHERE sta = '$sta' 
                AND reqid = '$reqid'";
    
    
    ($delta) = oracle_request($request);
   
    
    $command = "/export/dev/ops/bin/ttimes <<EOF\nALL\n\n$evdp\n$delta\n-1\n-1\nEOF\n";

    @data = `$command`;
    
    ($sks_delay)  = find_phase_delay("SKSac",@data);
    #($skks_delay) = find_phase_delay("SKKSac",@data);
    


    return ($delta, $sks_delay);

}

#DOC#################################################
#DOC
#DOC subroutine identify_event_details
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Identify event infomation
#DOC
#DOC input : $hypid - Hypocenter table key
#DOC
#DOC output : $evdate - Event origin time
#DOC          $evmsec - Milliseconds component of origin time
#DOC          $evla - Event latitude (degrees)
#DOC          $evlo - Event longitude (degrees)
#DOC          $evdp - Event depth (km)
#DOC
#DOC internal : $request - SQLplus request
#DOC            $temp - Output from SQLplus
#DOC
#DOC environment variables : $ENV{HYPOCENTER} - Oracle table
#DOC
#DOC calls subroutines : oracle_request (in 'shared.pl')
#DOC
#DOC#################################################

sub identify_event_details {

    my ($hypid) = @_;
    my ($evdate, $evmsec, $evla, $evlo, $evdp);
    my ($request, $temp);
    
    $request = "SELECT day, msec, lat, lon, depth
                FROM $ENV{HYPOCENTER}
                WHERE hypid = '$hypid'";
                
    ($temp) = oracle_request($request);
    
    ($evdate, $evmsec, $evla, $evlo, $evdp) = split(/::/,$temp);
    
    if ($evdp eq "") {
        
        $request = "select DEPDP from $ENV{HYPOCENTER}
                    where HYPID = '$hypid'";
                        
        ($evdp) = oracle_request($request);
            
        print "Depth substitution performed\n";
        print SUMMARY "Depth substitution performed\n";
        
    }
    
    return ($evdate, $evmsec, $evla, $evlo, $evdp);

}

#DOC#################################################
#DOC
#DOC subroutine identify_request_details
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Identify precalculated request details, (see output)
#DOC
#DOC input : $reqid - Request ID number
#DOC         $sta - Station code
#DOC
#DOC output : $hypid - Hypocenter table key
#DOC          $net - Station network
#DOC          $seaz - Station-to-event azimuth (back azimuth)
#DOC
#DOC internal : $request - SQLplus request
#DOC            $temp - Output from SQLplus
#DOC
#DOC environment variables : $ENV{SPLIT} - Oracle table
#DOC
#DOC calls subroutines : oracle_request (in 'shared.pl')
#DOC
#DOC#################################################

sub identify_request_details {


    my ($reqid, $sta) = @_;
    my ($hypid, $net, $seaz);
    my ($request, $temp);
    
    $request = "SELECT hypid, net, seaz
                FROM $ENV{SPLIT}
                WHERE sta = '$sta'
                AND reqid = '$reqid'";
                
    ($temp) = oracle_request($request);
    
    ($hypid, $net, $seaz) = split(/::/,$temp);
    
    return ($hypid, $net, $seaz);

}

#DOC#################################################
#DOC
#DOC subroutine identify_station_details
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Identify station details from ORACLE
#DOC
#DOC input : $net - FDSN network code
#DOC         $sta - Station code
#DOC
#DOC output : $stla - Station latitude (degrees)
#DOC          $stlo - Station longitude (degrees)
#DOC
#DOC internal : $request - SQLplus request
#DOC            $temp - Output from SQLplus
#DOC
#DOC environment variables : $ENV{SITE} - Oracle site table
#DOC
#DOC calls subroutines : oracle_request (in 'shared.pl')
#DOC
#DOC#################################################

sub identify_station_details {

    my ($net, $sta) = @_;
    my ($stla, $stlo);
    my ($request, $temp);
    
    $request = "SELECT lat, lon
                FROM $ENV{SITE}
                WHERE sta = '$sta'
                AND   net = '$net'";
                
    ($temp) = oracle_request($request);
    
    ($stla, $stlo) = split(/::/,$temp);
    
    return ($stla, $stlo);

}

#DOC#################################################
#DOC
#DOC subroutine remove_files
#DOC
#DOC author - M Evans
#DOC date - August 2003
#DOC
#DOC purpose : Removes unwanted files
#DOC
#DOC input : @files - List of filenames to be removed
#DOC
#DOC output : none
#DOC
#DOC internal : $file - Variable looped over in 'foreach' statement
#DOC
#DOC environment variables : none
#DOC
#DOC calls subroutines : none
#DOC
#DOC#################################################

sub remove_files {


    my (@files) = @_;
    my ($file);
    
    foreach $file (@files) {
        system ("rm $file");
    }
    

}

#DOC#################################################
#DOC
#DOC subroutine rotate_files
#DOC
#DOC author - M Evans
#DOC date - January 2004
#DOC
#DOC purpose : Rotates BH1, BH2 into BHN, BHE channels
#DOC
#DOC input : $one_file - Input channel
#DOC         $two_file - Input channel
#DOC
#DOC output : $n_file - North component
#DOC          $e_file - East component
#DOC
#DOC internal : $sac - SAC command
#DOC
#DOC environment variables : none
#DOC
#DOC calls subroutines : SAC
#DOC
#DOC#################################################

sub rotate_files {


    my ($one_file, $two_file, $n_file, $e_file) = @_;
    my ($command, $sac, $line, @sac_out, @bval, @eval, @temp, $bmax, $emin, $val, @az, @n1, @n2);
    
    undef @az;
    undef @bval;
    undef @eval;
        
    @n1 = split(/\//, $one_file); # 'Name' is now $n1[-1]
    @n2 = split(/\//, $two_file); # 'Name' is now $n2[-1]
    
    ### Check for orthogonality ###
    
    $sac =  "r $one_file $two_file\n";
    $sac .= "lh cmpaz\n";
        
    $command = "sac <<EOF\n$sac\nEOF\n";
    @sac_out = `$command`;
    
    foreach $line (@sac_out) {
        chomp ($line);
        if ($line =~ /cmpaz =/) {
            @temp = split(/ +/,$line);
            push (@az, $temp[3]);
        }
    }
    
    #print "@sac_out\n\n\n";
    
    print "Azimuths detected $n1[-1]: $az[0]   $n2[-1]: $az[1]\n";
    print SUMMARY "Azimuths detected $n1[-1]: $az[0]   $n2[-1]: $az[1]\n";
    
    if (abs($az[0]-$az[1]) != "90.") {
        if (abs($az[0]-$az[1]) != "270.") {
            print "Not othogonal, skipping....\n";
            print SUMMARY "Not orthogonal, skipping....\n";
        
            return ("-1");
        }
    }
    
    ### Synchronise and identify b and e values ###
    
    $sac =  "r $one_file $two_file\n";
    $sac .= "synchronise round on\n";
    $sac .= "w $one_file $two_file\n";
    $sac .= "lh b e\n";
        
    $command = "sac <<EOF\n$sac\nEOF\n";
    @sac_out = `$command`;
    
    #print "$sac\n\n@sac_out\n";
    
    foreach $line (@sac_out) {
        chomp ($line);
        if ($line =~ /b =/) {
            @temp = split(/ +/,$line);
            push (@bval, $temp[3]);
        } elsif ($line =~ /e =/) {
            @temp = split(/ +/,$line);
            push (@eval, $temp[3]);
        }
    }
    
    ### Identify cut values ###
    
    $bmax = $bval[0];
    
    foreach $val (@bval) {
        if ($val > $bmax) {
            $bmax = $val;
        }
    }
    

    $emin = $eval[0];
    
    foreach $val (@eval) {
        if ($val < $emin) {
            $emin = $val;
        }
    }
    
    $emin --;
    
    ### Cut and rotate

    $sac =  "cut $bmax $emin\n";
    $sac .= "r $one_file $two_file\n";
    $sac .= "rotate to 0 normal\n";
    $sac .= "w $one_file $two_file\n";

    $command = "sac <<EOF\n$sac\nEOF\n";
    @sac_out = `$command`;
    
    #print "$sac\n\n@sac_out\n";

    
    if ($one_file ne $n_file) {
        system ("mv $one_file $n_file");
        system ("mv $two_file $e_file");
    }
    
    return ("1");
   
}

######################################################

sub values_to_tables {


    my ($hypid, $phid, $rdid, $net, $sta, $delta, $seaz, $sks_start, $eang, $eangmin, $eangmax, $edt, $edtmin, $edtmax, $eig_win, $eig_or_ratio, $eig_cor_ratio, $pol_dir,
         $tang, $tangmin, $tangmax, $tdt, $tdtmin, $tdtmax, $tra_win, $tra_or_ratio, $tra_cor_ratio, $reporter, $eigscang, $eigscdt, $trascang, $trascdt) = @_;
    my ($command,  $eigscangmin, $eigscdtmin, $trascangmin, $trascdtmin, $eigscangmax, $eigscdtmax, $trascangmax, $trascdtmax);
    
    ### DO S&C CALCULATIONS
    
    $eigscangmin = $eang - $eigscang;
    $eigscangmax = $eang + $eigscang;
    $eigscdtmin  = $edt - $eigscdt;
    $eigscdtmax  = $edt + $eigscdt;
    
    $trascangmin = $tang - $trascang;
    $trascangmax = $tang + $trascang;
    $trascdtmin  = $tdt - $trascdt;
    $trascdtmax  = $tdt + $trascdt;
     
    ### LIMITS ###
    
    my ($or_max, $cor_min, $ang_lim, $dt_limit);
    my ($qc, $d_flag);
    
    #if ($net eq 'IR') {$net = "";}
   
    $ang_lim  = 45;
    $dt_limit = 5;
    $or_max   = 1/20;
    $cor_min  = 1/25;
    
    $qc = "";
    
    ### Angle ###
    
    if (($eangmax - $eangmin) > 45) {
        $qc .= "A"
    } else {
        $qc .= " "
    }
    
    ### Time ### (set if value > limit(5), or dt outside range)
    
    if ($edt > $dt_limit)  {
        $qc .= "T"
    } elsif ($edt < $edtmin) {
        $qc .= "T"
    } elsif ($edt > $edtmax) {
        $qc .= "T"
    } else {
        $qc .= " "
    }
    
    ### Time Error ###
    
    if ($edtmin < 0)  {
        $qc .= "E"
    } else {
        $qc .= " "
    }
    
    ### Original Eigenvalues ###
    
    if ($eig_or_ratio < $or_max)  {
        $qc .= "O"
    } else {
        $qc .= " "
    }
    
    ### Corrected Eigenvalues ###
    
    if ($eig_cor_ratio > $cor_min)  {
        $qc .= "C"
    } else {
        $qc .= " "
    }
    
    if ($qc eq "     ") {
        $qc = "";
        $d_flag = "";
    } else {
        $d_flag = "D";
    }
   
    $command = "INSERT INTO $ENV{RESULT} VALUES ('$hypid', '$phid', '$net', '$sta','EIGEN', 'BOOTS', 
      '$sks_start', '$eig_win', '$eang', '$eangmin', '$eangmax', '$edt', '$edtmin', '$edtmax', '$pol_dir','$d_flag', '$qc',
      '$eig_or_ratio', '$eig_cor_ratio', '$ENV{AUTHOR}', '$reporter',SYSDATE,'')";
    
    print "$command\n";
    
    oracle_command($command);

    $command = "INSERT INTO $ENV{RESULT} VALUES ('$hypid', '$phid', '$net', '$sta','EIGEN', 'SILVER', 
      '$sks_start', '$eig_win', '$eang', '$eigscangmin', '$eigscangmax', '$edt', '$eigscdtmin', '$eigscdtmax', '$pol_dir','$d_flag', '$qc',
      '$eig_or_ratio', '$eig_cor_ratio', '$ENV{AUTHOR}', '$reporter',SYSDATE,'')";
    
    print "$command\n";
    
    oracle_command($command);
    
    
    ####################################################################
    
    
    
    $qc = "";
    
    ### Angle ### (set if range > 45 degrees)
    
    if (($tangmax - $tangmin) > 45) {
        $qc .= "A"
    } else {
        $qc .= " "
    }
    
    ### Time ### (set if value > limit(5), or dt outside range)
    
    if ($tdt > $dt_limit)  {
        $qc .= "T"
    } elsif ($tdt < $tdtmin) {
        $qc .= "T"
    } elsif ($tdt > $tdtmax) {
        $qc .= "T"
    } else {
        $qc .= " "
    }
    
    ### Time Error ### (set if limits include zero)
    
    if ($tdtmin < 0)  {
        $qc .= "E"
    } else {
        $qc .= " "
    }
    
    ### Original Eigenvalues ### (set if PM too linear)
    
    if ($tra_or_ratio < $or_max)  {
        $qc .= "O"
    } else {
        $qc .= " "
    }
    
    ### Corrected Eigenvalues ### (set if corrected not linear enough)
    
    if ($tra_cor_ratio > $cor_min)  {
        $qc .= "C"
    } else {
        $qc .= " "
    }
    
    if ($qc eq "     ") {
        $qc = "";
        $d_flag = "";
    } else {
        $d_flag = "D";
    }
       
    $command = "INSERT INTO $ENV{RESULT} VALUES ('$hypid', '$phid', '$net', '$sta','TRANS', 'BOOTS', 
      '$sks_start', '$tra_win', '$tang', '$tangmin', '$tangmax', '$tdt', '$tdtmin', '$tdtmax', '','$d_flag', '$qc',
      '$tra_or_ratio', '$tra_cor_ratio', '$ENV{AUTHOR}', '$reporter',SYSDATE,'')";
    
    print "$command\n";
    
    oracle_command($command);
    
    $command = "INSERT INTO $ENV{RESULT} VALUES ('$hypid', '$phid', '$net', '$sta','TRANS', 'SILVER', 
      '$sks_start', '$tra_win', '$tang', '$trascangmin', '$trascangmax', '$tdt', '$trascdtmin', '$trascdtmax', '','$d_flag', '$qc',
      '$tra_or_ratio', '$tra_cor_ratio', '$ENV{AUTHOR}', '$reporter',SYSDATE,'')";
    
    print "$command\n";
    
    oracle_command($command);

    #
    ##
    ### General stuff ###
    ##
    #
    
    $command = "INSERT INTO $ENV{PHASE} (PHID, RDID,  PHASE, NET, STA, AUTHOR, REPORTER, LDDATE, CHAN)
                  VALUES ('$phid', '$rdid', 'SKSac', '$net', '$sta','$ENV{AUTHOR}', '$reporter', SYSDATE, 'BHR')";
    
    oracle_command($command);
    
 
    $command = "INSERT INTO $ENV{ASSOCIATION} (HYPID, PHID, PHASE, NET, STA, DELTA, SEAZ, AUTHOR, REPORTER, LDDATE)
                  VALUES ('$hypid', '$phid', 'SKSac', '$net', '$sta', '$delta', '$seaz','$ENV{AUTHOR}', '$reporter', SYSDATE)";

    oracle_command($command);

}

#
## end - values_to_table
#