#!/usr/bin/perl 
#
# For extending the unit-conversion part, have a look into
# GWR UIPC_Unit_Conversions.txt at http://www.oso.chalmers.se/~hgs/SG-docs/
#
sub sidiode {
    my ($x) = @_;
    my $y,$r,$xm,$t,$i;
    my @vl = ( 1.32412, 1.11732, 0.923142, 0.079767 );
    my @vu = ( 1.69812, 1.42013, 1.139350, 0.999614 );
    my @nd = ( 10, 11, 12, 11 );
    my @a = ( [  7.556358,   -5.917261,  0.237238, -0.334636, -0.058642,
                -0.019929,   -0.020715, -0.014814, -0.008789, -0.008554,  
                 0.000000,    0.000000 ],
              [ 17.304227,   -7.894688,  0.453442,  0.002243,  0.158036,
                -0.193093,    0.155717, -0.085185,  0.078550, -0.018312,  
                 0.039255,    0.000000 ],
              [ 71.818025,  -53.799888,  1.669931,  2.314228,  1.566635,
                 0.723026,   -0.149503,  0.046876, -0.388555,  0.056889, 
                 0.116823,    0.058580 ],
             [ 287.756797, -194.144823, -3.837903, -1.318325, -0.109120, 
                -0.393265,    0.146911, -0.111192,  0.028877, -0.029286,  
	       0.015619,    0.000000 ] );
    use Math::Trig;

    if    ( $x > $vu[0]                 ){ $r = 0 }
    elsif ( $x > $vu[1] && $x <= $vu[0] ){ $r = 0 }
    elsif ( $x > $vu[2] && $x <= $vu[1] ){ $r = 1 }
    elsif ( $x > $vu[3] && $x <= $vu[2] ){ $r = 2 }
    else                                 { $r = 3 }

    $xm  = (($x-$vl[$r])-($vu[$r]-$x))/($vu[$r]-$vl[$r]);
#    print STDERR "$x in range $r $vu[$r], mapped to $xm\n";
    if ( $xm < -1 || $xm > 1 )
    { if ( $qsderr )
      {
	  print STDERR "Sidiode in ch. $ch[$c]: Range $r $vu[3] $vl[3], signal $x out of range. $xm\n".
                       "Sidiode in ch. $ch[$c]: Probably sensor disconnected\n" ; 
      }
      $qsderr=0;
      return $uv
    }
    $t = 0.0;
    for( $i=0; $i < $nd[$r]; $i++ )
    {
      $t += $a[$r][$i]*cos($i*acos($xm))
    }
#    $rs = sprintf("%10.3f K ",$t);
#    return $rs

    return $t
}

sub uconv {
    my ($x,$n) = @_;
    my $a,$b,$c,$r,$y,$yt;
    if    ( $n =~ m/GEPLHe19/ ){ $y = 100 - 12.3817*$x }
    elsif ( $n =~ m/PSI/ )     { $y = ($x-6)/4 }
    elsif ( $n =~ m/PT100V/   ){ $y = $x/0.385 - 1.0 }
    elsif ( $n =~ m/factor,(-*\d+\.\d+)/) { $y = $1 * $x }
    elsif ( $n =~ m/TDVOhm/ )
    {
# HGS did it: TDVOhm (based on 1000*signal[V]) and TDOhm (1*)
#
      $y = ( -7.11936 + sqrt( 50.6853 - 0.036635*(1854 - 1000*$x)) )/0.0183175
#       $r=1854; $c=$r-1000*$x; $b=3.84 * 0.001 * $r; $a=4.94 * 0.000001 * $r;
#       $y = ((-$b) + sqrt(($b*$b) - (4.0*$a*$c))) / (2.0*$a)
    }
    elsif ( $n =~ m/TDOhm/ )
    {
      $y = ( -7.11936 + sqrt( 50.6853 - 0.036635*(1854 - $x)) )/0.0183175
#       $r=1854; $c=$r-$x; $b=3.84 * 0.001 * $r; $a=4.94 * 0.000001 * $r;
#       $y = ((-$b) + sqrt(($b*$b) - (4.0*$a*$c))) / (2.0*$a)
    }
    elsif ( $n =~ m/SiDiode-3/ ){ $yt = $x/4.98 }
    elsif ( $n =~ m/SiDiode-4/ || $n =~ m/SiDiode-2/ )
    { $yt = $x/4.98; 
      if ( $yt > 1.0918 ){ $yt=$yt - (0.2363 * ( $yt - 1.25 ))}
    }
    else  
    {
      if ( $qv && ( $opt_v || $opt_V )) 
      { if ( $opt_v)
	{ print "\n"; print "Conversion $n unknown\n" }
	else 
	{ $qv = 0 }
      }
      $y=$x 
    };
    if    ( $n =~ m/SiDiode/ ){ $y = sidiode($yt) }
    return $y;
} 

use Getopt::Std;
$q=0; $d=0; $kc=0; $xc=-1; $rc=1000000; $mr=1; $uv=-9999999.999; $co=0;
$qsderr=1;
getopts('C:cH:hI:l:m:N:n:qS:s:uvVX:x:');

if ( $opt_h ){
print <<"EOF";
USAGE\:
           xchan-units [opt] file

PURPOSE\:
           See http\://holt.oso.chalmers.se/hgs/hgs.man/xchan-units.html

           Extracts time and measurement columns from GWR files
           of type A1* A2* G1* (and *.TSF ?)
           and prints to stdout.
           Columns may be addressed by number or name.
           Name has the advantage that it is immune to changes in
           the column order. (If/when a name changes spelling,
           we'll have a new headache). 
          
           For .GGP files there is a sister script xggp.
           Options that apply identically to xggp are indicated
           by a   +   below.
                  |
                  V
OPTIONS\:
 -c               + Prints a table of the channels in the file.
 -q               - Like -c, all in one row, no decorations. 

 -C cset[,cset..] - units conversion presets:
                    cset = signal:factor  e.g. G1-Sig:-776.0

 -H #c[,#c...]    + Prints the header of channel number c.
                    Each c may be a range c1..c2
                  - Works also with channel names.

 -x #c[,#c...]    + The requested channel numbers to retrieve
                    Each c may be a range c1..c2

 -N c[,c...]      - The requested channel names to retrieve,
                    name must be spelled exactly. An alternative to -x

 -I #c[,#c...]    - Signal information for specified channel number.
 -S c[,c...]      - Signal information for specified channel name. 
                    Verbose mode, no data out.
                    prints   CONVERSION #chn S-name Conversion Units
                    like     CONVERSION 31 T2-Ext TDVOhm Celsius
                    upon xchan-units -S T2-Ext RAW_o054/A1121013.054
                    or   xchan-units -I 31 RAW_o054/A1121013.054

 -m m             + Print only each m\'th line
 -n r             + The number of records for output              [100000]
                    A value < 0 will suppress printing;
                    0 doesn\'t work, "Look at that stupid pe-herl!"
                    (Jagger&Richards)

 -u               - Do unit conversion (still incomplete)
                    Available\: factors, He-lvl, SiDiodes 2 3 4,
                    Pt100V, TDOhm, TDVOhm

 -s s             - Skip s records

 -l l             - Output l records (same function as -n but different
                    logic - maybe useful after a differenciated change. 

 -v               - Verbose mode, helpful for debug.
 -V               - reduced verbose. Shows selected channels
                    and unit conversions.
EOF
exit
}

%conv = (  "T1-Ext" => "TDVOhm",
           "T2-Ext" => "TDVOhm", 
           "T3-Ext" => "TDVOhm", 
          "Temp-G1" => "TDVOhm", 
          "Temp-G2" => "TDVOhm", 
          "Temp-TX" => "TDVOhm", 
          "Temp-TY" => "TDVOhm", 
          "Temp-TE" => "TDVOhm", 
          "Temp-AX" => "TDVOhm", 
          "Temp-CH" => "TDVOhm",
          "Dwr-Htr" => "factor,31.2",
          "HtrCrnt" => "factor,4.56"
         );

%units = ( "T1-Ext" => "Celsius",
           "T2-Ext" => "Celsius", 
           "T3-Ext" => "Celsius", 
          "Temp-G1" => "Celsius", 
          "Temp-G2" => "Celsius", 
          "Temp-TX" => "Celsius", 
          "Temp-TY" => "Celsius", 
          "Temp-TE" => "Celsius", 
          "Temp-AX" => "Celsius", 
          "Temp-CH" => "Celsius",
          "Dwr-Htr" => "mA",
          "HtrCrnt" => "mA"
         );
# above hard-wired conversion definitions,
# since GWR did not put in conversion directives
# for the environmental sensors into A1 and A2 file headers. 
#
if ( $opt_s ){ $s=$opt_s+1 } else { $s=0 }
if ( $opt_l ){ $o=$opt_l } else { $o=100000 }
if ( $opt_C )
{
    @convs = split(',',$opt_C);
    foreach $c ( @convs )
    { ( $cn, $cv ) = split(':',$c);
      $cv = "factor," . $cv;
      %conv = ( %conv, $cn => $cv );
      if ( $opt_v ){ print "Preset conversion $cn -> $cv\n" }
    }
}

$kr=-1; $qv=1; $qrc=1;
if ( $opt_m ){ $mr=$opt_m }
if ( $opt_n ){ $rc=$opt_n; if ( $rc < 0 ){ $rc=1; $qrc=0} }
if ( $opt_S || $opt_I ){ $rc=-1; $opt_x="" }

LINE: while (<>)
{ if ( $d ){ if ( $s > 0 ){$s+=-1}}
  if ( ( $o > 0 && $s <= 0 ) || ! $d ){
      if ( $d && $s <= 0 ){ $o+=-1 }
	chomp; 
	if ( $d && $xc>=0 )
        {  $rc+=-1; $kr++; $mm = $kr % $mr; 
	   if ( $rc ge 0 )
	   {
#          print "$_\n";
	     @rec = split(/ +/,$_);
	     if ( $mm eq 0 )
	     { if ( $qrc)
               { printf "%4d %2.2d %2.2d %2.2d %2.2d %2.2d ",
			   @rec[0,1,2,3,4,5]
               }
	       for ( $c=0; $c<=$xc; $c++ )
	       { 
		   if ( $opt_u ) # unit conversion
		   {
		     $xv=$rec[$ch[$c]+5];

		     if ( $opt_v )
{ print "-> $c :: $ch[$c] :: $name[$ch[$c]-1] :: $conv{$name[$ch[$c]-1]}\n";
  printf "%14.8f -> ",$xv 
}
                     if ( $rec[$ch[$c]+5] != $uv )
                     { 
		       $rec[$ch[$c]+5]=uconv($rec[$ch[$c]+5],$conv{$name[$ch[$c]-1]})
		     }
                     if ( $rec[$ch[$c]+5] == $uv )
                     {
			 $rec[$ch[$c]+5]=9999999.9999999;
# $rec[$ch[$c]+5]="miss"
                     }
		   }
                   if ($qrc){ printf "%15.7f ",$rec[$ch[$c]+5] }
	       } 
               if ( $qrc){ print "\n" }
	     }
	   } else { last LINE }
	} else {
# we've put in an 'else {' here, hoping for a speed-up
#
        if ( $q )
	{ if ( $_ =~ /^\[/ )
          { 
	    $q=0; $d=0;
            #print "End at $_\n"; 
          } else
	  { #print "look at $_\n";
	    if ( $_ =~ /.*:.*:(.*)/ ){ $a=$1 }
            $kc++;
            @rec = ( @rec , substr($a,0,-1));
            if ( $opt_v ){ print "Q $kc ->" . substr($a,0,-1) . "<-\n" }
          }
        }
	if ( $co )
        { @recc = split(/ +/,$_); 
          { $recc[$#recc]=substr($recc[$#recc],0,-1)}
	  if ( $opt_v ){ print "COMMENTS $#recc: @recc \n" }
	  if ( $recc[0]=~ m/\d+/ )
	  { if ( $recc[7] =~ m/(\d+\.\d+)/ ){ $recc[7]="factor,$1" }
         elsif ( $recc[7] =~ m/none/ ){ $recc[7]="factor,1.0" }
#
# Correcting some typos in the header:
            if ( $recc[1] =~ m/Temp_(\w\w)/ ){ $recc[1]="Temp-$1" }
         elsif ( $recc[1] =~ m/Dwr_Htr/ ){ $recc[1]="Dwr-Htr" }
         elsif ( $recc[1] =~ m/FB_Mod/ ){ $recc[1]="FB-Mod" }
         elsif ( $recc[1] =~ m/(\w+)_DC/ ){ $recc[1]="$1-DC" }
            if ( $recc[6] =~ m/Volts/ ){ $recc[6]="Volt" }
#
# Add this only if  $recc[1] => $recc[7] is new
#
# The preset can be specified from the command line
#
            if ( ! $conv{$recc[1]} )
            {  %conv = ( %conv,  $recc[1] => $recc[7] );
              %units = ( %units, $recc[1] => $recc[6] );
              if ( $opt_v )
	      { print "Add conversion: $recc[1] -> $recc[7] $recc[6]\n" } 
	    }
	  }
	}
#
# [DATA] part:
#
        if ( $_=~ /\[DATA\]/ )
        { $d=1; $co=0; 
#          print "DATA section $opt_c\n";
	  for ( $k=0; $k<$kc; $k++)
	  {
            $rec[$k] =~ s/_/-/g; 
            $rec[$k] =~ s/ //g; $name[$k]=$rec[$k]; 
            if ( $opt_v ){ print "NAME $k ->$name[$k]<-\n" }
	  }
	  if ($opt_H)
          { @loop = ( split(",",$opt_H) ); 
            foreach $l ( @loop )
            { if ( $l =~ /(\d+)\.\.(\d+)/ )
              { @iloop = ( $1..$2 );
                foreach $il ( @iloop )
                { printf "%2.2d:%s ",$il,$rec[$il-1] }
	      }
              else
              { if ( $l =~ m/^\d+/)
                {printf "%2.2d:%s ",$l,$rec[$l-1]}
                else
                { for ($ll=1; $ll<$#rec;$ll++)
                  {if ( $rec[$ll-1] eq $l ){printf "%2.2d:%s ",$ll,$rec[$ll-1]}}
		}
	      }
	    }
            print "\n";
            last LINE
	  }
          if ($opt_q)
          { for ( $j=0; $j<$kc; $j++)
	    { printf "%s ", $rec[$j] };
	    print "\n"; 
	    last LINE
	  }
          if ($opt_c)
          { print "Channels\n";
            print "        1        2        3        4        5        6        7        8        9        10\n";
	    for ( $j=0; $j<$kc; $j+=10)
	    { 
	      printf "%3d ",$j;
	      for ( $k=0; $k<10; $k++)
	      { 
	        printf "%9s", "$rec[$j+$k] " 
	      } 
	      print "\n\n" 
	    }
#            print "@rec \n"; 
	    last LINE;
          }
	  if ( $opt_x | $opt_I )
          { 
	    if ( $opt_x.$opt_I eq "H"){ 
		print "$kc\n";
		exit
            }
	    @cht=split(/,/,$opt_x.$opt_I); @chp= (); @ch= ();
            foreach $c ( @cht )
            { if ( $c =~ /(\d+)\.\.(\d+)/ )
              { @chp = ( @chp, ( $1..$2 ) ) } 
              else 
              { if ( $c =~ m/^\d+/ )
                { @chp = ( @chp, $c ) }
                else
                { for ($ll=1; $ll<$#rec;$ll++)
                  { $cc=$rec[$ll-1];
                    if ( $cc =~ m/$c/ ){ @chp = ( @chp, $ll ) }
                  }
		}
              }
            }
            if ($opt_v || $opt_V )
            { print "SELECTION REQUESTED: @chp\n" }
            $xc=$#chp;
            for ( $i=0; $i<=$xc; $i++ ) 
            { if ( $chp[$i] > $kc )
              { if ($opt_v || $opt_V )
                { print "Channel $chp[$i] not in file\n" }
                @ch = ( @ch, $ch[$#ch] )
	      }
              else 
              { @ch = ( @ch, $chp[$i] ) }
            }
            $xc=$#ch;
            if ($opt_v || $opt_V )
            { print "SELECTION GRANTED:   @ch\n" }
            
	  }
	  if ($opt_N || $opt_S )
          { 
	    if ( $opt_N.$opt_S eq "H"){ 
		print "$kc\n";
		exit
            }
	    @cht=split(/,/,$opt_N.$opt_S); @chp= (); @ch= ();
            foreach $c ( @cht )
	    { $in=0;
	      foreach $n ( @name )
		{ $in ++;
                  if ( $c eq $n )
                  { @chp = ( @chp, ( $in ) ) }
		}
            }
            if ($opt_v || $opt_V )
            { print "SELECTION REQUESTED: @chp\n" }
            $xc=$#chp;
            for ( $i=0; $i<=$xc; $i++ ) 
            { if ( $chp[$i] > $kc )
              { if ($opt_v || $opt_V )
                { print "Channel $chp[$i] not in file\n" }
                @ch = ( @ch, $ch[$#ch] )
	      }
              else 
              { @ch = ( @ch, $chp[$i] ) }
            }
            $xc=$#ch;
            if ($opt_v || $opt_V )
            { print "SELECTION GRANTED:   @ch\n" }
if ( $opt_N ){ print "@ch\n"; exit }
	  }
        }
#
# here channel prep dep. on argv
# 
        if ( $_=~ /\[GRAVITY_CAL_1_UGAL_V\]\W+(.*)./ )
        { if (!exists($conv{"Grav-1"}))
	  { %conv = ( %conv, "Grav-1" => "factor,".$1 );
            if ( $opt_v ){ print "Grav-1 factor: $1 from file\n" }
	  }
        }
        if ( $_=~ /\[PRESSURE_CAL_MBAR_V]\W+(.*)/ )
        { if (!exists($conv{"Baro-1"}))
	  { %conv = ( %conv, "Baro-1" => "factor,".$1 ); 
            if ( $opt_v ){ print "Baro-1 factor: $1\n" }
	  }}

        if ( $_ =~ /\[CHANNELS\]/ )
        { $q=1; 
          #print "Start at $_\n" 
        }    
	if ( $_ =~ /\[COMMENT\]/ )
        { $co=1; 
          if ( $opt_v ){ print "COMMENTS>\n" }
        }
	if ( $_ =~ /\[UNDETVAL\]\W+(.*)/ )
        { $uv=$1;
          if ( $opt_v || $opt_c ){ print "MRS: $uv\n" }
        }
  } 
}
}
# Here's the right brace of the speed-up 'else {'
#
if ( $opt_v )
{ @cs = ( keys %conv );
  foreach $key ( @cs ){ print "$key :: $conv{$key}\n" }
}
if ( $opt_V || $opt_S || $opt_I )
{ 
    for ( $j=0; $j<=$xc; $j++)
    {
        $units{$name[$ch[$j]-1]} =~ m/([%.,\/\w]+)/; $un=$1;
        $conv{$name[$ch[$j]-1]}  =~ m/([%.,\/\w]+)/; $cv=$1;
        
        print "CONVERSION $ch[$j] $name[$ch[$j]-1] $cv $un\n"
    }
}
exit
