#!/usr/local/bin/perl
#  dlt.pl
#  produce tables and plots of average density, lapse rate, and temperature
#  v0.1   Keith Lofstrom   2014 Oct 27

use strict;
use List::Util qw(sum min max);

my  $s         = 4 ; # half altitude span

my  @densfile  = ( "dens070", "dens120", "dens200" ) ; # density      g/cm^3
my  @tempfile  = ( "temp070", "temp120", "temp200" ) ; # temperature  K
my  @f107      = ( 70, 120, 200                    ) ; # e-22 W/Hz-m^2

my  @alt ;   # single dimensional altitude array
my  @den ;   # density array of arrays   
my  @lap ;   # lapse rate array of arrays
my  @tem ;   # temperature array of arrays

foreach my $i ( 0..$#f107 ) {
   my @alt1 = () ;  # altitude km
   my @den1 = () ;  # density kg/m^3
   my @tem1 = () ;  # temperature K
   
   # create temperature by averaging over a day
   open( IN, $tempfile[$i] ) ;
   while (<IN>) {
      next if /^#/ ;
      my @innum = split ;
      shift @innum                         ; # first value is altitude
      push @tem1 , sum( @innum ) / @innum  ; # average of temperatures
   }

   # create density and altitude arrays by averaging over a day
   open( IN, $densfile[$i] ) ;
   while (<IN>) {
      next if /^#/ ;
      my @innum = split ;
      push @alt1 , shift @innum            ; # first value is altitude
      push @den1 , sum( @innum ) / @innum  ; # average of densities
   }
   close( IN );

   if ( $#alt1 !=  $#den1 ) {
      die( "density ", $f107[$i]," ", $#alt1," ",$#den1  )
   }

   if ( $#alt1 !=  $#tem1 ) {
      die( "temperature ", $f107[$i]," ", $#alt1," ",$#tem1  )
   }

   # compute lapse rate array and copy altitude (redundant)
   foreach my $j ( 0..$#alt1 ) {
      my $n0 = max( 0,      $j-$s );
      my $n1 = min( $#alt1, $j+$s );
      $lap[$i][$j] = ($alt1[$n1]-$alt1[$n0])/log($den1[$n0]/$den1[$n1]);
      $alt[ $j]    =  $alt1[ $j]   ; # probably can do this with a  push
      $tem[$i][$j] =  $tem1[ $j]   ;
      $den[$i][$j] =  $den1[ $j] * 1000.0 ; # convert from g/cm3 to kg/m3
   }
}

print "#f107";
foreach my $i ( 0..$#f107 ) { printf "|%15.0f         ", $f107[ $i ] ; }
print "\n# alt";
foreach my $i ( 0..$#f107 ) { printf "|  density  lapse   temp "; }
print "\n" ;


foreach my $j ( 0..$#alt  ) {
   printf "%4.0f", $alt[$j] ;
   foreach my $i( 0..$#f107 ) {
      printf "%11.3e%7.2f%7.1f", $den[$i][$j], $lap[$i][$j], $tem[$i][$j] ;
   }
   print "\n" ;
}

exit ;