#!/usr/bin/perl -w
# $Id: satfilter.pl, V1.3 djl Exp $
#
# Copyright (c) 2004 Douglas Larson, All Rights Reserved.
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 2 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
#
# =====================================================================
# dbg = 0, no messages.
# dbg = 1, print debug messages
# dbg = 2, print copious debug messages
use lib "$ENV{HOME}/perl";
use strict;
use warnings;
# support options+arguments
use Getopt::Long qw(:config permute pass_through);
use IO::File;
use File::Basename;
use Pod::Usage; # avoid redundant &Usage()
use Data::Dumper; # Used for debugging
use DateConversions;
my $progname = "satfilter";
my $prog_version = 1.3;
my $opt_version;
my @fileList=();
my $opt_help;
my $opt_man;
my $opt_verbose=0;
my $opt_interactive=0;
my $opt_filelist=0;
my $opt_debug=0;
my $opt_sat="";
my $opt_style="ggcm"; # [ggcm|dg]
# User provided start and stop times
my $opt_tstart;
my $opt_tend;
GetOptions(
'debug=i' => \$opt_debug,
'help|?' => \$opt_help,
'man!' => \$opt_man,
'version!' => \$opt_version,
'verbose!' => \$opt_verbose,
'sat=s' => \$opt_sat,
'i' => \$opt_interactive,
'f' => \$opt_filelist,
'tstart=s' => \$opt_tstart,
'tend=s' => \$opt_tend,
'style=s' => \$opt_style
);
pod2usage(-verbose => 1) && exit if ($opt_debug !~ /^[012]$/);
pod2usage(-verbose => 1) && exit if defined $opt_help;
pod2usage(-verbose => 2) && exit if defined $opt_man;
$main::DateConversions::DebugLevel = $opt_debug;
if( $opt_debug > 0 ) {
$opt_verbose = 1;
print "Unprocessed by Getopt::Long" if $ARGV[0];
foreach (@ARGV) {
print "\n\t$_";
}
}
if( $opt_interactive ) {
print STDOUT ("\t>>>Interactive Mode <<<\n");
push( @fileList, &getfilename("Filename") );
$opt_filelist=1;
}
# Time windowing feature
my $Window = 0;
my %tstart;
my %tend;
if($opt_tstart) {
%tstart=buildTimeWindow($opt_tstart);
} else {
my $str="01/01/1900 00:00:00";
%tstart=buildTimeWindow($str);
}
if($opt_tend) {
%tend=buildTimeWindow($opt_tend);
} else {
my $str="31/12/2100 23:59:59";
%tend=buildTimeWindow($str);
}
if( $opt_tstart || $opt_tend) {
$Window = 1;
}
if( length($opt_sat) == 0 ) {
print STDOUT ("No dataset selected! You *must* specify:".
"\n\tACE GEOTAIL, GOES12, OMNI, POLAR, or WIND\n\n");
pod2usage(-verbose => 1) && exit;
}
if( ($opt_filelist) || ($ARGV[0]) ) {
foreach (@ARGV) {
push(@fileList, $_);
}
}
my $fileCount = @fileList;
foreach (@fileList) {
if ( &checkFileExists() == -1 ) {
print STDERR "Aborting due to problems with input file(s).\n";
exit(1);
}
}
foreach my $fullFileName (@fileList) {
if (! open (INFILE,"<$fullFileName") ) {
print STDERR "Can't open input file $fullFileName\n";
next;
}
(my $file , my $dir, my $ext ) = fileparse($fullFileName, qr/\.[^.]*/);
$opt_debug and print STDOUT ("Filename= $fullFileName\n\tBasename= $file\n");
print STDOUT ("Processing= $fullFileName\n");
my $mode="columns"; # Options are [serial | columns]
my $maxstr;
my $str_start;
my @str_missing = ();
my @timecodeProvided = ();
my @timecodeDesired = ("yr", "mo", "dd", "hh", "mm", "ss", "ms", "doy", "fp_year", "fp_doy");
my @timecodeCutKey = ();
my @satset = ();
my @outputSet= ();
if ( $opt_sat =~ /^ACE$/i ) {
if ( $opt_style =~ "ggcm" ) { $file = "ac" }
$mode="columns";
$maxstr = 100;
# Keep in mind they keep changing the file formats!
$str_start="^year day hr min sec fp_year fp_doy ACEepoch";
$str_start="^BEGIN DATA";
@str_missing=('-9999\.90','-999\.900');
# Column Names that the satellite team uses.
# This commented region is for an earlier version of the ACE data.
# my $strColumnNames="year,day,hr,min,sec,fp_year,fp_doy,ACEepoch,
# Np,Tp,Alpha_ratio,Vp,V_rtn_r,V_rtn_t,V_rtn_n,
# V_gse_x,V_gse_y,V_gse_z,V_gsm_x,V_gsm_y,V_gsm_z,
# B_rtn_r,B_rtn_t,B_rtn_n,B_gse_x,B_gse_y,B_gse_z,
# pos_gsm_x,pos_gsm_y,pos_gsm_z,pos_hs_x,pos_hs_y,pos_hs_z,MAG_pts";
my $strColumnNames="year,day,hr,min,sec,fp_year,fp_doy,ACEepoch,
Np,Tp,Alpha_ratio,Vp,V_rtn_r,V_rtn_t,V_rtn_n,
V_gse_x,V_gse_y,V_gse_z,V_gsm_x,V_gsm_y,V_gsm_z,
B_rtn_r,B_rtn_t,B_rtn_n,B_gse_x,B_gse_y,B_gse_z,
B_gsm_x,B_gsm_y,B_gsm_z,Bmag,Lambda,Delta,dBrms,
pos_gse_x,pos_gse_y,pos_gse_z,pos_gsm_x,pos_gsm_y,pos_gsm_z,
pos_hs_x,pos_hs_y,pos_hs_z,MAG_pts";
my @ColumnNames = split(/,/, $strColumnNames);
foreach ( @ColumnNames ) {
$_ =~ s/\s*//;
push( @outputSet , [$_] );
}
# push(@outputSet, ["pp"]);
push(@outputSet, ["pp"]); #re-added WDC
push(@outputSet, ["rr"]); # copy from np later
# my @ColWidth=( 4, 3, 2, 2, 6,13,10,13,9,11,
# 10,10,10,10,10,10,10,10,10,10,
# 10,10,10,10,10,10,10,10,10,10,
# 10,10,10,10,10,10,10,10,10,10,
# 10,10,10,10);
my @ColWidth=( 4, 3, 2, 2, 6,13,10,13,10,11,10,
9,9,9,9,9,9,9,9,9,9,
8,8,8,8,8,8,8,8,8,8,
8,8,8,
11,11,11,11,11,11,11,11,11,
7);
my @expected=(
# 0, 5, 9, 12, 15, 22, 36, 48, 62, 72,
# 84, 96, 106, 116, 126, 136, 146, 156, 166, 176,
# 186, 198, 206, 216, 224, 234, 243, 251, 261, 270,
# 279, 286, 296, 306, 313, 325, 337, 349, 361, 373,
# 385, 396, 408, 426);
0, 5, 9, 12, 15, 22, 36, 47, 61, 72,
84, 95, 105, 115, 125, 135, 145, 155, 165, 175,
185, 195, 204, 213, 222, 231, 240, 249, 258, 267,
276, 285, 294, 303, 312, 324, 336, 348, 360, 372,
384, 396, 408, 420);
@timecodeProvided=("yr", "doy", "hh", "mm", "ss", "fp_year", "fp_doy");
# [ value, starting column, width]
for( my $i = 0; $i < scalar(@timecodeProvided) ; $i++ ) {
push (@timecodeCutKey, [$timecodeProvided[$i], $expected[$i], $ColWidth[$i]]);
}
# [ value, starting column, width]
#my @ColStart = &satsetGuess(\*INFILE,$str_start,@expected); #WDC: bad
my @ColStart = @expected;
close(INFILE);
open (INFILE,"<$fullFileName");
if( $opt_debug > 0 ) {
print Dumper(\@ColStart);
}
my $NumberOfColumns = scalar(@ColumnNames);
for( my $i=0; $i < $NumberOfColumns; $i++) {
my $ColumnName = $ColumnNames[$i];
my @record = ( $ColumnNames[$i], $ColStart[$i], $ColWidth[$i] );
$satset[$i]=[ @record ];
}
if( $opt_debug > 0 ) {
print STDOUT ("satset Array\n");
print STDOUT Dumper(\@satset);
print STDOUT ("End of ACE init section\n");
}
} # End of ACE section
if ( $opt_sat =~ /^ACE_ftpweb$/i ) {
if ( $opt_style =~ "ggcm" ) { $file = "ac" }
$mode="columns";
$maxstr = 130;
$str_start="^YYYY DAY HR MN";
@str_missing=('-9999\.90');
# Column Names that the satellite team uses.
my $strColumnNames="year,doy,hr,min,
pos_gse_x,pos_gse_y,pos_gse_z,
Btot,B_gse_x,B_gse_y,B_gse_z,
Np,Tp,Alpha_ratio,
Vtot,V_gse_x,V_gse_y,V_gse_z";
my @ColumnNames = split(/,/, $strColumnNames);
foreach ( @ColumnNames ) {
$_ =~ s/\s*//;
push( @outputSet , [$_] );
}
my @ColWidth=( 4, 4, 3, 3, 9, 9, 9,
9, 9, 9, 9,
7, 9, 7,
9, 9, 9, 9);
my @expected=(
0, 5, 9, 12, 15, 24, 33, 42,
51, 60, 69, 78,
85, 94, 101,
110, 119, 128, 137);
@timecodeProvided=("yr", "doy", "hh", "mm");
# [ value, starting column, width]
for( my $i = 0; $i < scalar(@timecodeProvided) ; $i++ ) {
push (@timecodeCutKey, [$timecodeProvided[$i], $expected[$i], $ColWidth[$i]]);
}
# [ value, starting column, width]
my @ColStart = &satsetGuess(\*INFILE,$str_start,@expected);
close(INFILE);
open (INFILE,"<$fullFileName");
if( $opt_debug > 0 ) {
print STDOUT ("Column Start Array\n");
print STDOUT Dumper(\@ColStart);
}
my $NumberOfColumns = scalar(@ColumnNames);
for( my $i=0; $i < $NumberOfColumns; $i++) {
my $ColumnName = $ColumnNames[$i];
my @record = ( $ColumnNames[$i], $ColStart[$i], $ColWidth[$i] );
$satset[$i]=[ @record ];
}
if( $opt_debug > 0 ) {
print STDOUT ("satset Array\n");
print Dumper(\@satset);
print STDOUT ("End of ACE_ftpweb init section\n");
}
} # End of ACE_ftpweb section
if ( $opt_sat =~ /GEOTAIL/i ) {
if ( $opt_style =~ "ggcm" ) { $file = "ge" }
$maxstr = 100;
$str_start="^dd-mm-yyyy hh:mm:ss.ms";
@str_missing=('0\.10000E\+35.*');
@timecodeProvided=("dd", "mo", "yr", "hh", "mm", "ss", "ms");
# [ value, starting column, width]
@timecodeCutKey = ( ["yr", 6, 4], ["mo", 3, 2], ["dd", 0, 2],
["hh", 11, 2], ["mm", 14, 2], ["ss", 17, 2], ["ms", 20, 3] );
# [ value, starting column, width]
@satset = (
["xgse", 95, 14], ["ygse", 110, 14], ["zgse", 125, 14],
["bxgse", 31, 6], ["bygse", 38, 6], ["bzgse", 45, 6],
["bmag", 24, 6] );
} # End of GEOTAIL section
if ( $opt_sat =~ /GOES12/i ) {
if ( $opt_style =~ "ggcm" ) { $file = "go" }
$mode = "serial";
$maxstr = 30;
$str_start="^#yyyy-mm-dd hh:mm";
@str_missing=('1\.0E33.*');
@timecodeProvided=("yr", "mo", "dd", "hh", "mm");
# [ value, starting column, width]
@timecodeCutKey = ( ["yr", 0, 4], ["mo", 5, 2], ["dd", 8, 2],
["hh", 11, 2], ["mm", 14, 2] );
# [ value, starting column, width]
@satset = ( ["he", 16, 7], ["hn", 16, 7], ["hp", 16, 7], ["ht", 16, 7] );
} # End of GOES12 section
if ( $opt_sat =~ /^OMNI$/i ) {
if ( $opt_style =~ "ggcm" ) { $file = "om" }
$mode="columns";
$maxstr = 100;
$str_start="^dd-mm-yyyy hh:mm:ss.ms";
#WDC: catches B, V and X; T, N error when V does
#WDC: problem: T error on 7/12/2013 11:29 when no N, V error
@str_missing=('999\.9','1\.00000E\+07');
@timecodeProvided=("dd", "mo", "yr", "hh", "mm", "ss", "ms");
# [ value, starting column, width]
@timecodeCutKey = ( ["yr", 6, 4], ["mo", 3, 2], ["dd", 0, 2],
["hh", 11, 2], ["mm", 14, 2], ["ss", 17, 2], ["ms", 20, 3] );
# [ value, starting column, width]
my @ColumnNames = ( "mf_id", "pl_id",
"btot", "bxgse", "bygse", "bzgse", "bygsm", "bzgsm",
"vtot", "vxgse", "vygse", "vzgse",
"np", "temp", "xgse", "ygse", "zgse" );
my $ColPos = 23; #initial position past time field
#my @ColWidth=(12, 15, 17, 14, 14, 14, 16, 17, 17, 17, 15, 14, 14, 14,14);
my @ColWidth=(12,15,17,14,14,14,14,14,16,17,17,17,15,14,14,14,14);
foreach ( @ColumnNames ) {
$_ =~ s/\s*//;
push( @outputSet , [$_] );
}
my $NumberOfColumns = scalar(@ColumnNames);
for ( my $i=0; $i < $NumberOfColumns; $i++) {
my $ColumnName = $ColumnNames[$i];
my @record = ( $ColumnNames[$i], $ColPos, $ColWidth[$i] );
$satset[$i]=[ @record ];
#if ( $ColumnNames[$i] =~ /np/i ) {
#my @record = ( "rr", $ColPos, $ColWidth[$i] );
#$satset[$i]=[ @record ];
#push(@outputSet, ["rr"]);
#}
$ColPos += $ColWidth[$i];
}
push(@outputSet, ["pp"]); # calculate later
push(@outputSet, ["rr"]); # copy from np later
if( $opt_debug > 0 ) {
print STDOUT ("satset Array\n");
print STDOUT Dumper(\@satset);
print STDOUT ("End of OMNI init section\n");
}
} # End of OMNI section
if ( $opt_sat =~ /POLAR/i ) {
if ( $opt_style =~ "ggcm" ) { $file = "po" }
$maxstr = 100;
$str_start="^DATA:";
@str_missing=('0\.10000E\+35.*');
@timecodeProvided=("yr", "mo", "dd", "hh", "mm");
# [ value, starting column, width]
@timecodeCutKey = ( ["yr", 0, 4], ["mo", 6, 2], ["dd", 9, 2],
["hh", 12, 2], ["mm", 14, 2], ["ss", 18, 2], ["ms", 21, 3] );
# [ value, starting column, width]
@satset = (
["xgsm", 25, 13], ["ygsm", 38, 13], ["zgsm", 51, 13],
["bxgsm", 65, 13], ["bygsm", 78, 13], ["bzgsm", 91, 13],
["bmag", 104, 13] );
} # End POLAR
if ( $opt_sat =~ /WIND/i ) {
if ( $opt_style =~ "ggcm" ) { $file = "wi" }
$maxstr = 100;
$str_start="^DATA:";
@str_missing=('0\.10000E\+35.*');
@timecodeProvided=("yr", "mo", "dd", "hh", "mm");
# [ value, starting column, width]
@timecodeCutKey = ( ["yr", 0, 4], ["mo", 5, 2], ["dd", 8, 2],
["hh", 11, 2], ["mm", 14, 2], ["ss", 17, 2], ["ms", 20, 3] );
# [ value, starting column, width]
@satset = (
["xgsm", 77, 13], ["ygsm", 89, 13], ["zgsm", 103, 13],
["bxgsm", 25, 13], ["bygsm", 37, 13], ["bzgsm", 51, 13],
["bmag", 64, 13] );
} # End WIND Section
$opt_debug and print STDOUT ("Dataset = $opt_sat, mode = $mode\n");
my $header=1;
my $footer=1;
my $BadDataCount=0;
my $GoodDataCount=0;
my $Element;
my $MainLoop = 0;
while(<INFILE>) {
chomp;
if($MainLoop == 1) { last; }
my $line = $_;
if( $header == 1 ) {
if( $mode =~ /serial/i ) {
$opt_debug and print STDOUT $mode." ".$line."\n";
if( $line =~ /^\#Element\:.*/ ) {
$Element = $line;
$Element =~ s/^\#Element\://g;
# Blow away leading and trailing whitespace
$Element =~ s/^\s+|\s*\n$//g;
$opt_debug and print STDOUT (">>>> ".$Element."\n");
# Skip rest of header until the timecode string.
until($line =~ $str_start ) {
$line = <INFILE>;
}
$header and print STDOUT ("Beginning of data record: ".$Element."\n");
$header=0;
$footer=1;
#foreach (@satset) {
foreach (@outputSet) {
if( $_->[0] =~ /$Element/ ) {
$opt_debug and print STDOUT ("$_->[1] $_->[2] $_->[0]\n");
my $fh = IO::File->new("> $file.$_->[0]")
or die "Could not open file: $!";
$_->[3]=*$fh;
}
}
}
} else {
# Skip rest of header until the timecode string.
until($line =~ $str_start ) {
$line = <INFILE>;
}
if( $line =~ /$str_start/ ) {
$header and print STDOUT ("Beginning of data record.\n");
$header=0;
$footer=1;
#foreach (@satset) {
if( $opt_debug > 0 ) {
print STDOUT ("<outputSet>\n");
print STDOUT Dumper(\@outputSet);
print STDOUT ("</outputSet>\n");
}
foreach (@outputSet) {
if( $opt_debug >= 2 ) {
print STDOUT ($_->[0]."\n");
}
my $fh = IO::File->new("> $file.$_->[0]")
or die "Could not open file: $!";
$_->[1]=*$fh;
}
}
} # Element:
} else {
if( $mode =~ /serial/i ) {
if( $_ =~ /\#\>/ ) {
$footer and print STDOUT ("End of record reached.\n");
$header=1;
$footer=1;
}
} else {
if( length($_) < $maxstr ) {
$footer and print STDOUT ("End of data record reached.\n");
$footer=0;
$MainLoop=1;
}
}
if( ( $header == 0 ) && ( $footer == 1 ) ) {
my %timecode = &buildTimecodeArray($line,\@timecodeProvided,\@timecodeDesired,\@timecodeCutKey);
my $str_timecode;
if( $opt_style =~ "dg" ) {
# Deltagraph time format
my $time=0.0;
if( $timecode{fp_year} > 0 ) {
$time = $timecode{fp_year};
} else {
my $doy = $timecode{doy};
my $dayfrac = &HMSToFraction($timecode{hh},$timecode{mm},$timecode{ss});
my $time = $doy + $dayfrac;
}
$str_timecode = sprintf( "%0f", $time);
} elsif ( $opt_style =~ "ggcm" ) {
$str_timecode = sprintf("%4d %2d %2d %2d %2d %6.3f",
$timecode{yr},$timecode{mo},$timecode{dd},
$timecode{hh},$timecode{mm},$timecode{ss});
}
my $GoodData=1;
foreach (@str_missing) {
if( $line =~ /$_/ ) { $GoodData=0; }
}
if( $GoodData == 0 ) {
$BadDataCount++;
} else {
$GoodDataCount++;
}
#if( $line =~ /$str_missing/ ) {
# $BadDataCount++;
# $GoodData=0;
#} else {
# $GoodDataCount++;
# $GoodData=1;
#}
if( $Window ) {
if( $timecode{fp_year} < $tstart{fp_year} ) {
if( $opt_debug > 0 ) {
print STDOUT ($Window." ".$timecode{fp_year}." < ".$tstart{fp_year}."\n");
}
$GoodData=0;
}
if( $timecode{fp_year} > $tend{fp_year} ) {
if( $opt_debug > 0 ) {
print STDOUT ($Window." ".$timecode{fp_year}." > ".$tend{fp_year}."\n");
}
$GoodData=0;
$footer=1;
$MainLoop=1;
}
}
if( $GoodData ) {
if( $opt_debug > 0 ) {
print STDOUT ($Window." ".$tstart{fp_year}." < ".$timecode{fp_year}." < ".$tend{fp_year}."\n");
}
# Now spit out the processed file(s).
if( $mode =~ /columns/i ) {
# First we read the whole line and create a line hash.
my %linearray=();
foreach (@satset) {
my $data = substr($line,$_->[1],$_->[2]);
# Delete leading blanks
$data =~ s/^\ *//g;
if( $data =~ /\ / ) {
# We only do this if the field width was too big
# Delete from end to blank
$data =~ s/\ *[\w|\-|\.]*$//g;
}
$linearray{$_->[0]}=$data;
if( $opt_debug > 0 ) {
print STDOUT ("offset=".$_->[1].":: width=".$_->[2]."::\n");
print STDOUT ("$str_timecode ::".$linearray{$_->[0]}."::\n");
if( $opt_debug > 2 ) {
print STDOUT ("\n$line\n");
exit(0);
}
}
}
if ( $opt_sat =~ /^ACE$/i ) {
$linearray{Vtot}=&VectorMagnitude($linearray{V_gse_x},
$linearray{V_gse_y}, $linearray{V_gse_z});
# Now do add Solar Wind Plasma Pressure.
# Psw = (Solar wind density) * (Solar wind velocity)**2
# $linearray{pp}= 1.6726E-6 * $linearray{Np} *
# $linearray{Vtot} * $linearray{Vtot};
# WDC: above is dyn pressure, we need plasma pressure
# Psw = nkT = (density) * k * (temperature)
$linearray{pp}= $linearray{Np} * 1.3806488E-23 *
$linearray{Tp} * 1.E18;
# copy np to rr file
$linearray{rr}= $linearray{Np};
# convert satellite position from km to RE
$linearray{pos_gse_x}= $linearray{pos_gse_x} / 6371.;
$linearray{pos_gse_y}= $linearray{pos_gse_y} / 6371.;
$linearray{pos_gse_z}= $linearray{pos_gse_z} / 6371.;
}
if ( $opt_sat =~ /^OMNI$/i ) {
# Psw = nkT = (density) * k * (temperature)
$linearray{pp}= $linearray{np} * 1.3806488E-23 *
$linearray{temp} * 1.E18;
# WRONG! - convert bow shock position from RE to km
#$linearray{xgse}= $linearray{xgse} * 6378.;
#$linearray{ygse}= $linearray{ygse} * 6378.;
#$linearray{zgse}= $linearray{zgse} * 6378.;
# copy np to rr file
$linearray{rr}= $linearray{np};
}
foreach (@outputSet) {
if( $opt_debug > 0 ) {
print STDOUT ($_->[0]."\n");
}
my $str_data = sprintf("%s %14.8f",$str_timecode,$linearray{$_->[0]});
print {$_->[1]} ("$str_data\n");
#print STDOUT ("$str_data\n");
}
} elsif ( $mode =~ /serial/i ) {
# need to assign $_ to the current dataset.
foreach (@satset) {
if( $_->[0] =~ /$Element/ ) {
my $data = substr($line,$_->[1],$_->[2]);
print {$_->[3]} ("$str_timecode $data\n");
}
}
}
} # GoodData
} # Header && footer
} # Header == 1
} # Loop through the input file
close (INFILE);
print STDOUT ("Records found: Good= $GoodDataCount, Bad=$BadDataCount\n");
} # Input filelist loop
exit(0);
# ************************************************************************
sub satsetGuess {
local(*INPUTFILE) = shift;
my $str_start = shift;
my @expected = @_;
print STDOUT "str=$str_start\n";
my @retArray=();
my $line = <INPUTFILE>;
until($line =~ $str_start ) {
$line = <INPUTFILE>;
}
# This next read should be a data line.
$line = <INPUTFILE>;
# Tokenize
my $columnNumber=0;
while ( $line =~ m/((\+|\-)*\w[\w\.(\+|\-)]*)/g ) {
my $strBefore=$`;
my $w = length($strBefore);
$retArray[$columnNumber]=$w;
if( $opt_debug > 0 ) {
print STDOUT "$columnNumber\texpected=$expected[$columnNumber]\tfound=$w\t$1\n";
}
$columnNumber++;
}
return(@retArray);
}
# ************************************************************************
sub getfilename {
#
# Pupose :: Interactively ask a user for a valid filename.
#
my $prompt = $_[0];
my $done = 0;
my $retval;
while ($done == 0) {
print STDOUT ($prompt," > ");
$_ = scalar(<STDIN>);
chomp;
$opt_debug and print STDOUT ("getfilename ::",$_,"::\n");
if ( &checkFileExists() == -1 ) {
$done = 0;
} else {
$retval=$_;
$done = 1;
}
}
return ($retval);
}
# ************************************************************************
sub checkFileExists {
#
# Purpose : If a file is empty, or does not exist, return error -1.
#
my $filename = $_;
my $retval=1;
$opt_debug and print STDOUT ("checkFileExists ::",$filename,"::\n");
if ( -d $filename ) {
print STDOUT ("File $filename is actually a directory\n");
$retval = -1;
}
if ( !( -e $filename ) ) {
print STDOUT ("File $filename does NOT exist!\n");
$retval = -1;
}
if ( -z $filename ) {
print STDOUT ("File $filename is an empty file\n");
$retval = -1;
}
return ($retval);
}
# ************************************************************************
sub buildTimecodeArray {
#
# Purpose : All modes must munch the timecode.
#
my $line = shift;
my $refTimecodeDesired = shift;
my $refTimecodeProvided = shift;
my $refTimecodeCutKey = shift;
my %timecode=();
# Initialize the timecode array.
# keys = yr, mo, dd, hh, mm, ss, ms, doy, fp_year, fp_doy
if( $opt_debug == 2 ) {
print STDOUT ("timecodeDesired:\n");
print Dumper($refTimecodeDesired);
print STDOUT ("timecodeProvided:\n");
print Dumper($refTimecodeProvided);
}
$timecode{yr} = -1;
$timecode{mo} = -1;
$timecode{dd} = -1;
$timecode{hh} = -1;
$timecode{mm} = -1;
$timecode{ss} = -1;
$timecode{ms} = -1;
$timecode{fp_year} = -1;
$timecode{fp_doy} = -1;
$timecode{doy} = -1;
# Load what we are provided into the timecode array.
foreach (@$refTimecodeCutKey) {
my $str_temp = substr($line,$_->[1],$_->[2]);
# Blow away leading and trailing whitespace
$str_temp =~ s/^\s+|\s*\n$//g;
$timecode{$_->[0]} = scalar($str_temp);
if( $opt_debug >= 2 ) {
print STDOUT ("timecode{".$_->[0]."}=".$timecode{$_->[0]}."\n");
}
}
# Now we must fill in any missing bits.
if( $timecode{yr} =~ -1 ) {
if( $timecode{fp_year} > 0 ) {
$timecode{yr} = int($timecode{fp_year});
}
}
if( $timecode{fp_year} =~ -1 ) {
if( $timecode{yr} > 0 ) {
$timecode{fp_year} = 1.0 * $timecode{yr};
} else {
print STDOUT ("Ambiguous: yr=". $timecode{yr}." fp_year=".$timecode{fp_year}."\n");
exit(2);
}
}
my $isleap = &f_isleap($timecode{yr});
if( $timecode{mo} =~ -1 ) {
if( $timecode{fp_year} > 0 ) {
my $year = int($timecode{fp_year});
my $fp_doy=0;
if( $isleap) {
$fp_doy = 366.0 * ($timecode{fp_year} - $year);
} else {
$fp_doy = 365.0 * ($timecode{fp_year} - $year);
}
#$timecode{doy} = int($fp_doy);
#print STDOUT ("timecode{fp_doy}=".$timecode{fp_doy}." ".$doy."\n");
$timecode{mo} = &doytomonth($year,$timecode{doy});
}
}
#if( $timecode{dd} =~ -1 ) {
# $timecode{dd} = &doytodom($timecode{yr},$timecode{doy});
#}
if( $opt_sat =~ /ACE/i ) {
# ACE actually uses DOY in the day column.
if( ( $timecode{doy} > 0 ) && ( $timecode{fp_doy} =~ -1 ) ) {
$timecode{dd} = &doytodom($timecode{yr},$timecode{doy});
$timecode{fp_doy}= 1.0 * $timecode{doy};
} else {
$timecode{dd} = &doytodom($timecode{yr},$timecode{fp_doy});
}
print STDOUT ("doytodom: ".$timecode{dd}."\n");
}
if( $timecode{ss} =~ -1 ) {
$timecode{ss} = 0;
}
if( $timecode{ms} =~ -1 ) {
$timecode{ms} = 0;
}
if( $opt_debug >= 2 ) {
print STDOUT ("timecode:\n");
print Dumper(\%timecode);
if( $opt_debug >= 3 ) { exit(0); }
}
return(%timecode);
}
# ************************************************************************
sub buildTimeWindow {
my $str_time = shift;
my %Time;
if( $str_time ) {
my $retcode = &DateConversions::validateDate($str_time);
if( $retcode == 0 ) {
if( $opt_debug > 0 ) {
print STDOUT ("\nbuildTimeWindow:$retcode (0-bad, 1-OK)\n");
}
$Time{dd} = 31;
$Time{mo} = 12;
$Time{yr} = 2100;
$Time{hh} = 23;
$Time{mm} = 59;
$Time{ss} = 59;
$Time{ms} = 0;
} else {
# DD/MO/YEAR HH:MM:SS
$Time{dd} = int(substr($str_time, 0, 2));
$Time{mo} = int(substr($str_time, 3, 2));
$Time{yr} = int(substr($str_time, 6, 4));
$Time{hh} = int(substr($str_time, 11, 2));
$Time{mm} = int(substr($str_time, 14, 2));
$Time{ss} = int(substr($str_time, 17, 2));
$Time{ms} = 0;
}
$Time{doy} = &datetodoy($Time{yr}, $Time{mo},$Time{dd},
$Time{hh}, $Time{mm},$Time{ss},$Time{ms});
my $tmp = &HMSToFraction($Time{hh},$Time{mm},$Time{ss},$Time{ms});
$Time{fp_doy} = $Time{doy} + $tmp;
my $isleap = &f_isleap($Time{yr});
if( $isleap) {
$Time{fp_year} = $Time{yr} + ($Time{fp_doy} / 366.0);
} else {
$Time{fp_year} = $Time{yr} + ($Time{fp_doy} / 365.0);
}
}
if( $opt_debug > -1 ) {
print Dumper(\%Time);
}
return(%Time);
}
# ************************************************************************
sub VectorMagnitude {
my (@vector) = @_;
my $vecmag= sqrt( $vector[0] * $vector[0] +
$vector[1] * $vector[1] +
$vector[2] * $vector[2] );
return($vecmag);
}
# ************************************************************************
# POD
# ************************************************************************
BEGIN {
print "\n";
} END {
if(defined $opt_version){
print
"\nModules, Perl, OS, Program info:\n",
" Pod::Usage $Pod::Usage::VERSION\n",
" Getopt::Long $Getopt::Long::VERSION\n",
" strict $strict::VERSION\n",
" Perl $]\n",
" OS $^O\n",
" $0 $prog_version\n",
"\n\n";
}
}
=pod
=head1 NAME
satfilter.pl - A filter for magnetometer data.
=head1 SYNOPSIS
satfilter.pl [--help] [--man] --sat=[ACE|ACE_ftpweb|GEOTAIL|GOES12|OMNI|POLAR|WIND] file(s)
=head1 EXAMPLE WORKFLOW
=over 4
=item B<Run the script:>
=item I<satfilter.pl --verbose ~/Some_Magnetometer_Data>
=item I<satfilter.pl --sat=ACE --style=ggcm --debug=2 ACE_MAGSWE_Data.txt>
=back
=head1 DESCRIPTION
ACE http://www.srl.caltech.edu/ACE/ASC/level2/lvl2DATA_SWEPAM.html
GEOTAIL, GOES12, POLAR, WIND magnetometer data from online sources is
converted into a {timcode value} format, with one value per output file.
POLAR, WIND, GEOTAIL data - ftp://pwgdata.gsfc.nasa.gov/pub/
Various Spacecraft (out of date) - ftp://nssdcftp.gsfc.nasa.gov/spacecraft_data/
Space Environment Center (SPIDR service) - http://www.sec.noaa.gov/
Goddard Space Physics Data Facility http://sscweb.gsfc.nasa.gov/
=head1 ARGUMENTS
=over 4
=item --help Print Options and Arguments.
=item --man Print complete man page.
=item --mode=[serial|columns]
=back
=head1 OPTIONS
=over 4
=item --versions Print Modules, Perl, OS, Program info
=item --debug=0 Do NOT print debugging information [default]
=item --debug=1 Print debugging information
=item --fmt=[ggcm|dg] The default is 6 columns with year, month, day, hour, minute, seconds. dg=Timecode in Decimal Day Format (good for plotting)
=item --ts="DD/MO/YEAR HH:MM:SS" Start time [default is beginning of file]
=item --te="DD/MO/YEAR HH:MM:SS" End time [default is the end of file]
=back
=head1 AUTHOR
Dr. Douglas Larson
University of New Hampshire
Space Science Center
Durham, NH
03824
=head1 TESTED
Getopt::Long
Perl 5.8.1
=head1 SUBROUTINES
getfilename(string) - Interactively ask for a filename.
checkFileExists() - Dies various existance checks on a filename.
datetodoy - Convert a date to the Day Of Year.
f_isleap -Return 0 if a year is NOT leap year and a 1 otherwise.
satsetGuess - Try and guess which columns the data reside.
buildTimecodeArray - Given whatever time is provided, make a uniform time array.
buildTimeWindow - This is used to limit the times saved. Good for event studies.
doytomonth - Day Of Year to month number.
HMSToFraction - Convert wallclock to fractional time.
validateDate - Make sure the user's command line time selection is OK.
=head1 BUGS
None that I know of. ;-)
=head1 TODO
CDF and/or HDF.
=head1 UPDATES
=over 4
=item Wednesday, 12 May, 2004
Satellite magnetometer data filter.
=item Friday, 20 May, 2005
Extensive modifications. Added a time windowing capability. Added GGCM
output option. Added ACE full year data file. Many new subroutines associated
with dates and time.
=item Wednesday, 22 May, 2007
Added the ACE_ftpweb option and a minor update to the start string on ACE option. The ACE team changed their software such that the data has 'BEGIN METADATA; and 'BEGIN DATA', which makes finding the header much easier. They also added a half dozen new columns of commonly needed info like GSE position.
=item Wednesday, 4 December, 2013
Added OMNI option. Data must come from CDAWeb with the following items included: source IMF s/c, source Plasma s/c, B, Bx(GSE), By(GSE), Bz(GSE), V, Vx(GSE), Vy(GSE), Vz(GSE), proton density, temperature, bow shock location X(GSE). Added calculation of plasma pressure ("pp") for ACE and OMNI. Changed to output to file with GGCM prefix ("ac", etc.) when GGCM option selected. Still need to have output files properly named for ggcm.
=back
=cut