# Trims NEC-4.1 and NEC-2D output file <filename>.NOU, sorting the output into the
# following ASCII text files:
#    <filename>.CM -- contains comment cards.
#    <filename>.AEG -- azimuth, elevation, and total gain table.
#    <filename>.CON -- drive conditions
#    <filename>.DAT -- gain statistics.
# 
# Call:
#    gawk -f noutrim.awk <target.xy> <filename.NOU>
#
# The target zones used to calculate statistics are defined by the file target.xy
# This file has two lines each containing 4 integers and a word:
#    Target 1: az-left to az-right, el-lo to el-hi, name (one word)
#    Target 2: az-left to az-right, el-lo to el-hi, name (one word)
# If a target zone is not needed, set the azimuths or elevations equal to each other.
# Both lines are REQUIRED, even if target zone(s) are not needed.
#
# NOTE: Points in Target zone 1 are EXCLUDED from the statistics for Target zone 2.
#
# NOTE: The variable LowestGainForAveraging is set to the lowest gain that will be used in averaging.
# This recognizes that real-world construction can not maintain tolerances and stability necessary for very low gain levels,
# e.g., more than 30 dB below the main beam of a large yagi.  Environmental scattering will also fill in deep nulls in patterns.
#
# NOTE: Only looks for Sommerfeld ground conditions.
#
# There are several states in processing a .NOU file:
#    Init -- First pass through NEC output file; used to initialize noutrim's output files.
#    FindVersion -- Looking for NEC version type.
#    PreComment -- Skipping lines to reach the comment cards.
#    Comment -- Copying comment cards to the .CM file.
#    PreAEGHeader -- Searching for the start of az-el gain data table.
#    PreAEGData -- Found the header; skipping lines to start the data.
#    AEGData -- Now copying the data.
# Finally, the END routine writes out the accumulated statistics.
#
# 02 Sep 23: Added logic to process NEC-2D as well as NEC-4.1 files.
# 02 Sep 23: Located NEC version number.  Annotate comment file with NEC version.
# 02 Sep 23: Find comments correctly for NEC-4.1 files.  Still not correct for NEC-2D.  No Comment case also not correct.
# 02 Dec 30: v2.0 - Corrected drive impedance and current calculations for a current source.
# 02 Dec 30: v2.0 - Annotated comments with NouTrim version number.
# 03 Jan 24: v2.1 - Added error warning if azimuth data wraps past 360 degrees.
# 03 May 2:  v2.2 - Prints version number on screen.
# 03 Jul 1:  v2.3 - Corrected deviation calculations to use deviation from mean power.
# 03 Jul 2:  v2.4 - Fixed all known errors in processing NEC-2D and -4.1 files:
#                   -- Handles blank comment cards in the middle of the CM deck properly (see bug 02 Sep 23).
#                   -- To save space in the plot annotation, multiple consecutive blank CM lines collapse into one blank line.
#                   -- Handles case of no comments correctly.
#                   -- No longer plots the word "NOUTrim" twice in succession in the left annotation column.
# 03 Jul 7:  v3.0 - Added gain statistics based on medians.
# 03 Jul 9:  v3.1 - Fixed bug in putting drivepoint impedances into Conditions file.
# 03 Jul 10: v3.2 - Added mean of median, and max/min stats for zones 1+2 combined.
# 03 Jul 11: v3.3 - Removed medians within each target zone from PostScript map (commented out).
# 03 Jul 13: v3.4 - Fixed bug so that combined target stats do not print when one target is missing.
# 03 Jul 16: v3.5 - Fixed bug in NEC-4.1 processing that skipped over the first comment line.
# 03 Jul 17: v3.6 - Clarified description of the combined target 1+2 statistics.
# 03 Jul 17: v3.7 - Corrected references from "mode" to "median".
# ===================================================================================================================
BEGIN {
   NoutrimVersion = "NOUTrim v3.7 2003 Jul 21 Mon"
   print NoutrimVersion

   # Logical constants
   TRUE = 1
   FALSE = 0

   # First round through NOU file will be for initialization.
   state = "Init"

   # Initialize variables for statistics outputs
   Gain = 0.0
   Azimuth = 0.0
   Elevation = 0.0
   LowestGainForAveraging = -15.0
   LowestPowerForAveraging = 10 ^ ( LowestGainForAveraging / 10)   # -15 dBi is the lowest gain that will be used for statistics.
   MaxGain = -100.00
   NrDataPoints = 0.0
   Power = 0.0

   # Initialize tag number to avoid having string forced into it.
   TagNr = 0.0

   # Get the target zones
      # Initialize variables to avoid having strings stuffed into the numeric ones.
      EXISTS = 1
      LEFT = 2
      RIGHT = 3
      LOW = 4
      HIGH = 5
      INCLUDESNORTH = 6
      NRDATAPOINTS = 7
      NAME = 8
      SUMPOWERS = 9
      MAXDBI = 10
      MAXDBIAZ = 11
      MAXDBIEL = 12
      MINDBI = 13
      MINDBIAZ = 14
      MINDBIEL = 15
      MEDIANPOWER = 16
      StepSize = 0.0

      # Initialize array for two target zones and a third, untargeted region.
      for (i=1; i<=3; i++)
      {
         Target[ i, EXISTS] = FALSE
         Target[ i, LEFT] = 0.0
         Target[ i, RIGHT] = 0.0
         Target[ i, LOW] = 0.0
         Target[ i, HIGH] = 0.0
         Target[ i, INCLUDESNORTH] = 0.0
         Target[ i, NRDATAPOINTS] = 0.0
         Target[ i, SUMPOWERS] = 0.0
         Target[ i, MAXDBI] = -100.00		# NEC never generates more than +99.99 dBi nor less than -99.99 dBi gain.
         Target[ i, MAXDBIAZ] = 0.0
         Target[ i, MAXDBIEL] = 0.0
         Target[ i, MINDBI] = 100.00		# NEC never generates more than +99.99 dBi nor less than -99.99 dBi gain.
         Target[ i, MINDBIAZ] = 0.0
         Target[ i, MINDBIEL] = 0.0
         Target[ i, MEDIANPOWER] = 0.0
      }

      # Now get the zone boundaries from the first input file (target.xy).
      for (i=1; i<=2; i++)
      {
         getline
         left = $1
         right = $2
         low = $3
         high = $4	
         if ( (left != right) && (low != high) )
         {
            Target[ i, EXISTS] = TRUE
            Target[ i, LEFT] = left
            Target[ i, RIGHT] = right
            Target[ i, LOW] = low
            Target[ i, HIGH] = high
            Target[ i, NAME] = $5
            if ( left > right ) Target[ i, INCLUDESNORTH] = TRUE

            # Write a description of the target borders for GMT to plot.
            # Vertical lines are great circle lines that follow longitude (azimuth).
            printf "%5.1f %4.1f\n%5.1f %4.1f\n>\n", left, low, left, high >"targetvertical.xy"    # great circle route up left longitude.
            printf "%5.1f %4.1f\n%5.1f %4.1f\n>\n", right, low, right, high >"targetvertical.xy"    # great circle route up right longitude.

            # Horizontal lines are non-great circles lines that follow latitude (elevation).  A separate GMT psxy statement is needed for these.
            printf "%5.1f %4.1f\n%5.1f %4.1f\n>\n", left, high, right, high >"targethorizontal.xy"    # straight across the top.
            printf "%5.1f %4.1f\n%5.1f %4.1f\n>\n", left, low, right, low >"targethorizontal.xy"    # straight across the bottom.
         }
      }
   # All finished with target zone definitions.
   # Any extra lines in the target zone file may cause trouble with the start of NEC .NOU file processing.
}

# ===================================================================================================================
# Line-by-line NEC output file processing
# This part of a GAWK command program is looped through until the entire contents of all remaining files have been processed.
{
   # Initializing noutrim's output files.
   if (state == "Init")
   {
      # make filenames.
      split( FILENAME, file, ".")
      CommentFile = file[1] ".CM"
      ConditionsFile = file[1] ".CON"
      AEGFile = file[1] ".AEG"
      StatisticsFile = file[1] ".DAT"

      # Initialize the GMT paragraph format strings for pstext in output files.
      printf ">  +0.00 3.50 8 0 0 TL 9p 4.00 l\n" >CommentFile
      printf ">  +0.00 3.50 8 0 0 TL 9p 2.00 l\n" >ConditionsFile
      printf ">  +0.00 3.50 8 0 0 TL 9p 4.00 l\n" >StatisticsFile

      state="FindVersion"
   }

   # ===================================================================================================================
   # Find the NEC version used to produce these results.
   if (state == "FindVersion")
   {
      if( $0 ~ /NUMERICAL ELECTROMAGNETICS CODE/)
      {
         NECversion = "unknown"		# Zero is flag value indicating no NEC version was recognized.

         if ($0 ~ /NEC-2D/ )
         {
            NECversion = "2D"
            state="PreComment"   # Must now hunt down a line that says "COMMENTS"
         }

         if ($0 ~ /NEC-4.1/ )
         {
            NECversion = "4.1"
            for (i=1; i<=8; i++) getline	            # skip 9 lines to reach first comment line.
            state = "Comment"
         }

         if (NECversion == "unknown")
         {
            printf "NEC version not recognized!\n"
            exit   # bails out, executing the END statements.
         }
      }
   }

   # ===================================================================================================================
   # Search for the start of comment cards.
   if (state == "PreComment")
   {
      if (NECversion == "2D" &&  $0 ~ /COMMENTS/)
      {
         # skip two blank lines to reach the first comment.
         getline
         getline
         state = "Comment"
	 getline		# Now grab the first line with a comment.
      }
   }   

   # ===================================================================================================================
   # Copy the NEC comment cards to the comment file.
   if (state == "Comment")
   {
      NrConsecutiveBlankLines = 0.0		# Used to avoid printing multiple blank lines or lines with asterisks that NEC4 inserts.
      NoCommentCards = TRUE                     # Used to detect the no comment card case.
      getline
      for ( ; !($0 ~ /STRUCTURE SPECIFICATION/); getline)
      {
	 # Eliminate multiple blank lines or lines with all asterisks.
         if ((NF == 0) || ($0 ~ /\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*/ ) )
         {
            NrConsecutiveBlankLines++
            if ( NoCommentCards == FALSE && NrConsecutiveBlankLines == 1) printf "\n\n" >CommentFile		# Print just one consecutive blank line if it's not the first comment card.
         }
         else
         {
            NrConsecutiveBlankLines = 0.0
            NoCommentCards = FALSE		# Found at least one non-blank comment card.

            # NEC comment lines come with trailing spaces, which must be eliminated to plot correctly in GMT.
            LastChar = length()		# Find the last character number.
            LastNonSpaceChar = LastChar	# Starting at the right end of the line...
            while ( substr( $0, LastNonSpaceChar, 1) == " ") LastNonSpaceChar--	# Work towards the left until a non-space is found.
            $0 = substr( $0, 1, LastNonSpaceChar)		# Now set $0 to just the substring from start to last non-space char.

            sub( /\040\040/, "")		# First two chars are always a space and can be eliminated.
            printf "%s\n", $0		# Spit it onto the screen for the user to check.

            gsub( /\040/, "\\040")		# Substitute '\040' for all other spaces to keep GMT from deleting the starting spaces.
            gsub( /\#/, "\\#")		# Replace # with \# to keep AWK from deleting lines of text.

            printf "%s\n\n", $0 >CommentFile		# Spit out the cleaned-up comment line.
         }
      }
      # Seems like no more comments.
      if (NoCommentCards == TRUE) printf "(No comments.)\n\n" >CommentFile
      printf "NEC version %s\n", NECversion
      printf "NEC version %s\n\n", NECversion >CommentFile
      printf "\nK3NA %s\n\n", NoutrimVersion >CommentFile
      state = "PreEXCard"
   }

   # ===================================================================================================================
   # Searching for the quantity of EX cards, which will tell us how many feedpoints exist
   if (state == "PreEXCard")
   {
      NrEXCards = 0
      if( NECversion == "4.1" && $0 ~ /INPUT LINE/)        # true when we hit the start of the input card listing.
      {
         while ( $0 ~ /INPUT LINE/)
         {
            if ( $5 ~ /EX/) NrEXCards++
            getline    # grab the next line.  All the EX cards should be together.
         }
         state = "PreFrequency"            # No more INPUT LINE cards, so we are done counting EX cards.
      }
      if( NECversion == "2D" && $0 ~ /DATA CARD/)        # true when we hit the start of the input card listing.
      {
         while ( $0 ~ /DATA CARD/)
         {
            if ( $6 ~ /EX/) NrEXCards++
            getline    # grab the next line.  All the EX cards should be together.
         }
         state = "PreFrequency"            # No more INPUT LINE cards, so we are done counting EX cards.
      }
   }

   # ===================================================================================================================
   # Searching for the start of frequency data.
   if (state == "PreFrequency") {
      if ($0 ~ /FREQUENCY=/) {
         printf "%7.3f MHz\n\n\n", $2 >ConditionsFile
         state = "PreGround"
      }
   }

   # ===================================================================================================================
   # Searching for the start of ground information.
   if (state == "PreGround") {
      if ($0 ~ /FINITE GROUND/) {
         printf "%s\n\n", tolower( $0) >ConditionsFile
         getline
         printf "\\040\\040\\040%4.1f dielectric constant\n\n", $4  >ConditionsFile
         getline
         printf "\\040\\040\\040%6.1f mS/m\n\n\n", 1000 * $2  >ConditionsFile
         state = "PreImpedance"
      }
   }

   # ===================================================================================================================
   # Searching for drive point impedance results.
   if ( (state == "PreImpedance") && ($0 ~ /ANTENNA INPUT PARAMETERS/) )
   {
      # true when we hit the top of the antenna input parameters listing.
      # The actual data is fixed field size.  Minus signs caused fields to run together, so normal GAWK processing can't be used.
      # Use substr to extract the correct data.
      printf "Drive impedance; voltage or current:\n\n" >ConditionsFile
      getline    # skip next three lines
      getline
      getline

      for (i=1; i<=NrEXCards; i++) {
         getline    # grab a line of data and extract values

         # For a voltage source model, assume the tag number is <900.
         # For a current source model, assume the tag number is >=900.
         # Multiply substr by 1.0 to force it to be a number.
         TagNr = 1.0 * substr( $0, 1, 6)

         if ( TagNr < 900) {
             Vr = substr( $0, 13, 12)
             Vi = substr( $0, 25, 12)
             Zr = substr( $0, 61, 12)
             Zi = substr( $0, 73, 12)
             Vmag = sqrt( Vr^2 + Vi^2)
             Vphase = (180/3.14159) * atan2( Vi, Vr)

             # special printing in order to put "j" and plus/minus signs in the right places
             printf "\\#%i:  ", i >ConditionsFile
             if (Zi > 0) printf "%7.1f + j%7.1f @~\\127@~", Zr,  Zi >ConditionsFile
             else        printf "%7.1f - j%7.1f @~\\127@~", Zr, -Zi >ConditionsFile
             printf "\\040\\040\\040%7.3fV @~\\320@~ %5.1f\\217\n\n", Vmag, Vphase >ConditionsFile
	 }
         else {
             # Current source model.  The following is per W4RNL email dated 02 Dec 24/0435 EST.
             # The listed "admittance" for the dummy wire is the actual drive impedance.
             # The listed "voltage" for the dummy wire is the actual drive current + 90 degrees.
             Ir = substr( $0, 13, 12)
             Ii = substr( $0, 25, 12)
             Zr = substr( $0, 85, 12)
             Zi = substr( $0, 97, 12)
             Imag = sqrt( Ir^2 + Ii^2)
             Iphase = -90 + ( (180/3.14159) * atan2( Ii, Ir))

             # special printing in order to put "j" and plus/minus signs in the right places
             printf "\\#%i:  ", i >ConditionsFile
             if (Zi > 0) printf "%7.1f + j%7.1f @~\\127@~", Zr,  Zi >ConditionsFile
             else        printf "%7.1f - j%7.1f @~\\127@~", Zr, -Zi >ConditionsFile
             printf "\\040\\040\\040%7.3fA @~\\320@~ %5.1f\\217\n\n", Imag, Iphase >ConditionsFile
         }
      }
      state = "PreAEGHeader"
    }

   # ===================================================================================================================
   # Searching for the start of az-el gain data.
   if ( (state == "PreAEGHeader") && ( $0 ~ /RADIATION PATTERNS/) )   # Found AEG header. Skip next 4 lines to reach data.
   {
      getline	# NEC 2D & 4.1:  This line is blank.
      getline	# NEC 2D & 4.1:  - - ANGLES - -
      getline	# NEC 2D & 4.1:  THETA     PHI
      getline	# NEC 2D & 4.1:  DEGREES  DEGREES
      state = "AEGData"
      getline	# First line of actual data.
   }

   # ===================================================================================================================
   # In AEG data until a blank record is reached.
   if (state == "AEGData") {
      if ( NF > 0) {
         # Convert elevation and azimuth information
         Azimuth = 360-$2
         Elevation = 90-$1
         if ( $2 > 360) printf "Error: Az > 360 at phi=%5.2f, theta=%6.2f.\n", $2, $1
         else
         {
            # Grab gain and write to Az-El-Gain file.
            Gain = $5
            printf "%6.2f %5.2f %6.2f\n", Azimuth, Elevation, Gain >AEGFile

            # Is this the largest gain encountered to date?
            if ( Gain > MaxGain)
            {
               MaxGain = Gain
               MaxGainAzimuth = Azimuth
               MaxGainElevation = Elevation
            }

            # Accumulate sum of powers for later averaging in efficiency calculation.
            NrDataPoints++
            Power = 10^( Gain/10)
            SumPower = SumPower + Power

            # Is this point within target zone 1, 2 or the non-target area (zone 3)?
            InTarget = 3   # Begin by assuming point is not in a target.
            # Check first for Zone 2, then zone 1, since a point in zone 1 is not considered part of zone 2.
            for( i=2; i>=1; i--)
            {
               # Check if within range of elevations
               if ( (Target[ i, EXISTS] == TRUE) &&
                    ( Elevation >= Target[ i, LOW ] ) &&
                    ( Elevation <= Target[ i, HIGH] ) )
               {
                  # Check if within range of azimuths
                  if ( Target[ i, INCLUDESNORTH] == TRUE )
                  {
                     if ( ( Azimuth >= Target[ i, LEFT ] ) || ( Azimuth <= Target[ i, RIGHT] ) ) InTarget = i
                  }
                  else
                  {
                     if ( ( Azimuth >= Target[ i, LEFT ] ) && ( Azimuth <= Target[ i, RIGHT] ) ) InTarget = i
                  }
               }
            }

            # To save a pass through the gain bin counts later, accumulate sum of powers in each target for later averaging (mean).
            Target[ InTarget, NRDATAPOINTS]++
            if ( Power < LowestPowerForAveraging ) PowerOrFloor = LowestPowerForAveraging
            else PowerOrFloor = Power
            Target[ InTarget, SUMPOWERS] = Target[ InTarget, SUMPOWERS] + PowerOrFloor

            # See if the new datapoint gain is a maximum or minimum within the target zone.
            if ( Gain > Target[ InTarget, MAXDBI] )
            {
               Target[ InTarget, MAXDBI] = Gain
               Target[ InTarget, MAXDBIAZ] = Azimuth
               Target[ InTarget, MAXDBIEL] = Elevation
            }
            if ( Gain < Target[ InTarget, MINDBI] )
            {
               Target[ InTarget, MINDBI] = Gain
               Target[ InTarget, MINDBIAZ] = Azimuth
               Target[ InTarget, MINDBIEL] = Elevation
            }

            # To calculate median statistics, power or gain values in each zone must be in an ordered list.
            # Take advantage of NEC's granularity in calculating gain to a precision of 0.01 dBi.
            # Bin the gain values from the floor to the maximum gain in the target zone, counting the number of times each gain value appears.
            # This approach takes less space than an ordered list of 20,629 points.
            # From a floor of -15.00 dBi to a maximum of +21.00 dBi, just 3600 bins are required per zone if a value appears at each possible bin.
            # Maximum number of bins for 3 zones is 10,800 bins; usually there will be far less.
            # Thus, it is faster to use bin counts than to create three ordered lists.

            if (Gain < LowestGainForAveraging) GainOrFloor = LowestGainForAveraging
            else GainOrFloor = Gain

            # GAWK arrays and array indices are strings.  Thus, the initial value of an array position is "", not numeric zero.
            # To avoid confusion, force the index to be an integer that will always get converted into the same kind of string.
            GainIndex = int( 100 * GainOrFloor)

            # Bump up the count in the appropriate gain bin.
            # Use one dimension arrays for each target zone so that the convention "if (value in array)" can be used later.
            if (InTarget == 1) TargetGain1[ GainIndex]++
            if (InTarget == 2) TargetGain2[ GainIndex]++
            if (InTarget == 3) TargetGain3[ GainIndex]++
         }
      }
      else state = "PreAEGHeader"
   }
}

# All done parsing the .NOU file.
END {
   print NrDataPoints " data points found."

   printf "%+4.2f dBi maximum gain at %3i az  %2i el.\n", MaxGain, MaxGainAzimuth, MaxGainElevation
   printf "Maximum gain:  %+4.2f dBi at %3i\217 az  %2i\217 el.\n\n", MaxGain, MaxGainAzimuth, MaxGainElevation >StatisticsFile

   printf "%5.1f %% power efficiency.\n", 100 * (SumPower/NrDataPoints) / 2
   printf "%5.1f %% power efficiency.\n\n", 100 * (SumPower/NrDataPoints) / 2 >StatisticsFile

   # Summarize the target zone statistics.
   for (i=1; i<=4; i++)			# i=3 for calculating remainder of sky outside zones 1 and 2; i=4 for calculating zones 1+2 together.
   {
      if ( (Target[ i, EXISTS] == TRUE) \
           || (i == 3) \
           || ( (i == 4) && Target[ 1, EXISTS] == TRUE && Target[ 2, EXISTS] == TRUE) )
      {
         # Calculate mean power and gain, except for combination of zones 1+2 (where only median is of interest).
         if ( i!=4)
         {
            MeanPower = Target[ i, SUMPOWERS] / Target[ i, NRDATAPOINTS]
            MeanGain = 10 * log( MeanPower) / 2.302585			# Note: GAWK log function is base e.  Divide by 2.302585 to get base 10.
         }

         # Median is the middle figure in the ordered list.
         if (i==4) MidPoint = int( (Target[ 1, NRDATAPOINTS] + Target[ 2, NRDATAPOINTS]) / 2)
         else MidPoint = int( Target[ i, NRDATAPOINTS] / 2)

         # Now go through the existing bins once, starting at the floor gain value, accumulating the following data:
         #   -- sum of power deviation above and below mean (has to be done on both sides of mean, as some values may equal mean.)
         #   -- nr points below mean.
         #   -- nr points above mean.
         #   -- sum of power deviation above and below median (must be done on both sides of median).
         #   -- nr of points below median
         #   -- nr of points above median (can be different as some points may equal the median value).

         PointsFound = 0.0
         PointsBelowMean = 0.0
         PointsBelowMedian = 0.0
         PointsAboveMean = 0.0
         PointsAboveMedian = 0.0
         SumPowerDeviationBelowMean = 0.0
         SumPowerDeviationBelowMedian = 0.0
         SumPowerDeviationAboveMean = 0.0
         SumPowerDeviationAboveMedian = 0.0
         MedianFound = FALSE
         if ( i!=4) TotalPoints = Target[ i, NRDATAPOINTS]
         else TotalPoints = Target[ 1, NRDATAPOINTS] + Target[ 2, NRDATAPOINTS]

         for ( GainIndex = int( 100*LowestGainForAveraging); \
              PointsFound < TotalPoints; \
              GainIndex++)
         {
            # How many points have this gain value?
            NrPoints = 0
            if ( i==1 && GainIndex in TargetGain1) NrPoints = TargetGain1[ GainIndex]
            if ( i==2 && GainIndex in TargetGain2) NrPoints = TargetGain2[ GainIndex]
            if ( i==3 && GainIndex in TargetGain3) NrPoints = TargetGain3[ GainIndex]
            if ( i==4)
            {
               NrPoints=0
               if ( GainIndex in TargetGain1) NrPoints = NrPoints + TargetGain1[ GainIndex]
               if ( GainIndex in TargetGain2) NrPoints = NrPoints + TargetGain2[ GainIndex]
            }

            if (NrPoints != 0)
            {
               PointsFound = PointsFound + NrPoints	# Used to see if we have found all the points yet and can stop hunting.
               Gain = GainIndex/100
               Power = 10^( Gain/10)

               if ( Gain < MeanGain && i != 4 )		# Don't calculate means for zone 1+2.
               {
                  PointsBelowMean = PointsBelowMean + NrPoints
                  SumPowerDeviationBelowMean = SumPowerDeviationBelowMean + ( NrPoints * (Power - MeanPower))		# Accumulates a negative number.
               }

               if ( Gain > MeanGain && i != 4)		# Don't calculate means for zone 1+2.
               {
                  PointsAboveMean = PointsAboveMean + NrPoints
                  SumPowerDeviationAboveMean = SumPowerDeviationAboveMean + ( NrPoints * (Power - MeanPower))
               }

               if ( PointsFound < MidPoint)		# Have not reached median yet, so still accumulating below-median statistics.
               {
                  PointsBelowMedian = PointsBelowMedian + NrPoints
                  SumPowerDeviationBelowMedian = SumPowerDeviationBelowMedian + ( NrPoints * Power)			# Note: We will subtract the MedianPower later when we find out its value.
               } 

               if ( PointsFound >= MidPoint)
               {
                  if ( !MedianFound)		# Crossed mid-point of ordered list; grab the median values.
                  {
                     MedianFound = TRUE
                     MedianGain = Gain
                     MedianPower = Power
                     Target[ i, MEDIANPOWER] = MedianPower
                     SumPowerDeviationBelowMedian = SumPowerDeviationBelowMedian - (PointsBelowMedian * MedianPower)	# Back out the MedianPower, now that we know what it is.
                  }

                  if ( Gain > MedianGain)
                  {
                     PointsAboveMedian = PointsAboveMedian + NrPoints
                     SumPowerDeviationAboveMedian = SumPowerDeviationAboveMedian + ( NrPoints * (Power - MedianPower))
                  } 
               }
            }		# Done with this gain value.
         }		# Done looking at all gain values that exist in this target.

         # Calculate average deviation from mean gain
         if (i!=4)
         {
            AveragePowerDeviationAboveMean = SumPowerDeviationAboveMean / PointsAboveMean
            AverageGainDeviationAboveMean = 10 * log( (MeanPower + AveragePowerDeviationAboveMean) / MeanPower) / 2.302585   # Note: GAWK log function base e.  Divide by 2.302585 to get base 10.
            AveragePowerDeviationBelowMean = SumPowerDeviationBelowMean / PointsBelowMean
            AverageGainDeviationBelowMean = 10 * log( (MeanPower + AveragePowerDeviationBelowMean) / MeanPower) / 2.302585   # Note: GAWK log function base e.  Divide by 2.302585 to get base 10.
         }

         # Calculate average deviation from median gain
         # Note: the median may be at the gain floor, especially outside the two target areas.
         #       That will cause PointsBelowMedian to be zero.

         AveragePowerDeviationAboveMedian = SumPowerDeviationAboveMedian / PointsAboveMedian
         AverageGainDeviationAboveMedian = 10 * log( (MedianPower + AveragePowerDeviationAboveMedian) / MedianPower) / 2.302585		# Note: GAWK log function base e.  Divide by 2.302585 to get base 10.
         if ( PointsBelowMedian > 0)
         {
            AveragePowerDeviationBelowMedian = SumPowerDeviationBelowMedian / PointsBelowMedian
            AverageGainDeviationBelowMedian = 10 * log( (MedianPower + AveragePowerDeviationBelowMedian) / MedianPower) / 2.302585
         }
         else
         {
            AveragePowerDeviationBelowMedian = 0
            AverageGainDeviationBelowMedian = 0
         }

         # Screen and DAT file output
         # Print appropriate headline:
         if ( i <  3)
         {
            printf "\nTarget %1i - %s: %3i-%3i az  %2i-%2i el.\n",           i, Target[ i,NAME], Target[ i,LEFT], Target[ i,RIGHT], Target[ i,LOW], Target[ i,HIGH]
            printf "\nTarget %1i - %s: %3i-%3i\217 az  %2i-%2i\217 el.\n\n", i, Target[ i,NAME], Target[ i,LEFT], Target[ i,RIGHT], Target[ i,LOW], Target[ i,HIGH] >StatisticsFile
         }
         if ( i == 3)
         {
            printf "\nNon-targeted area statistics:\n"
            printf "\nNon-targeted area statistics:\n\n" >StatisticsFile
         }
         if ( i == 4)
         {
            printf "\nEqual weight mean of target 1 and 2 statistics:\n"
            printf "\nEqual weight mean of target 1 and 2 statistics:\n\n" >StatisticsFile
         }

         # Screen diagnostic: datapoint count
         if ( i != 4) printf "   %i data points.\n", Target[ i, NRDATAPOINTS]

         # Median statistics
         printf             "   Median gain*: %+5.2f dBi.  Avg dev: %+5.2f above, %+5.2f below.\n",   MedianGain, AverageGainDeviationAboveMedian, AverageGainDeviationBelowMedian
         printf "\\040\\040\\040Median gain*: %+5.2f dBi.  Avg dev: %+5.2f above, %+5.2f below.\n\n", MedianGain, AverageGainDeviationAboveMedian, AverageGainDeviationBelowMedian >StatisticsFile

         # Max and min dBi and locations
         if ( i != 4)
         {
            printf "   Max gain:  %+5.2f dBi at %3i az  %2i el.\n", Target[ i, MAXDBI], Target[ i, MAXDBIAZ], Target[ i, MAXDBIEL]
            printf "   Min gain:  %+5.2f dBi at %3i az  %2i el.\n", Target[ i, MINDBI], Target[ i, MINDBIAZ], Target[ i, MINDBIEL]
            printf "\\040\\040\\040Min: %+5.2f dBi at %3i\217 az  %2i\217 el.", Target[ i, MINDBI], Target[ i, MINDBIAZ], Target[ i, MINDBIEL] >StatisticsFile
            printf "\\040\\040\\040Max: %+5.2f dBi at %3i\217 az  %2i\217 el.\n\n", Target[ i, MAXDBI], Target[ i, MAXDBIAZ], Target[ i, MAXDBIEL] >StatisticsFile
         }
         else
         {
            # Did Target 1 or 2 have the max or min gains?
            if ( Target[ 1, MAXDBI] > Target[ 2, MAXDBI]) MaxTarget=1
            else MaxTarget = 2
            if ( Target[ 1, MINDBI] < Target[ 2, MINDBI]) MinTarget=1
            else MinTarget = 2

            # Print out correct one
            printf "   Max gain:  %+5.2f dBi at %3i az  %2i el.\n", Target[ MaxTarget, MAXDBI], Target[ MaxTarget, MAXDBIAZ], Target[ MaxTarget, MAXDBIEL]
            printf "\\040\\040\\040Max: %+5.2f dBi at %3i\217 az  %2i\217 el.", Target[ MaxTarget, MAXDBI], Target[ MaxTarget, MAXDBIAZ], Target[ MaxTarget, MAXDBIEL] >StatisticsFile
            printf "   Min gain:  %+5.2f dBi at %3i az  %2i el.\n", Target[ MinTarget, MINDBI], Target[ MinTarget, MINDBIAZ], Target[ MinTarget, MINDBIEL]
            printf "\\040\\040\\040Min: %+5.2f dBi at %3i\217 az  %2i\217 el.\n\n", Target[ MinTarget, MINDBI], Target[ MinTarget, MINDBIAZ], Target[ MinTarget, MINDBIEL] >StatisticsFile

            # Print out equal-weighted mean of the two targets' median, max and min
            MeanMedian = 10 * log( ( Target[ 1, MEDIANPOWER] + Target[ 2, MEDIANPOWER]) / 2 ) / 2.302585
            MeanMax = 10 * log( ( 10^( Target[ 1, MAXDBI]/10) + 10^( Target[ 2, MAXDBI]/10) ) / 2) / 2.302585
            MeanMin = 10 * log( ( 10^( Target[ 1, MINDBI]/10) + 10^( Target[ 2, MINDBI]/10) ) / 2) / 2.302585
            printf             "   Mean median: %+5.2f dBi.  Mean max:  %+5.2f dBi.  Mean min: %+5.2f dBi.\n", MeanMedian, MeanMax, MeanMin
            printf "\\040\\040\\040Mean median: %+5.2f dBi.  Mean max:  %+5.2f dBi.  Mean min: %+5.2f dBi.\n\n", MeanMedian, MeanMax, MeanMin >StatisticsFile
         }

         # Mean statistics, except for combined target 1+2. 
         if ( i != 4)
         {
            printf             "   Mean gain*: %+5.2f dBi.  Avg dev: %+5.2f above, %+5.2f below.\n",   MeanGain, AverageGainDeviationAboveMean, AverageGainDeviationBelowMean
#           # To save space, don't annotate map with mean data.
#           printf "\\040\\040\\040Mean gain*: %+5.2f dBi.  Avg dev: %+5.2f above, %+5.2f below.\n\n", MeanGain, AverageGainDeviationAboveMean, AverageGainDeviationBelowMean >StatisticsFile
         }
      }
      else
      {
         if ( i!=4)
         {
            printf "\nTarget %1i not used.\n", i
            printf "\nTarget %1i not used.\n\n", i >StatisticsFile
         }
      }
   }

   # Footnote the minimum gain used in averaging within zones.
   printf "\n* Gain floor for calculating median, mean and deviation is %+3.1f dBi.\n", LowestGainForAveraging
   printf "\n* Gain floor for calculating median, mean and deviation is %+3.1f dBi.\n\n", LowestGainForAveraging >StatisticsFile
   printf "NEC x-axis is North azimuth.\n" >StatisticsFile

   close( CommentFile)
   close( ConditionsFile)
   close( AEGFile)
   close( StatisticsFile)
}