#!/usr/bin/perl -w

# Given input in mJy/beam, creates a copy in K.
# Optionally applies the NVSS flux registration.  If the flux correction factor
# is < 0.1, (i.e. 0 = not ready yet), it is not applied.

use Astangles;
use DRAORed;
use Getopt::Std;
use strict;

my $nvssfrff = "/home/rreid/f_fluxreg_c21.text";

my $use = q {
  Use: kelvinize [-r] file1 [ file2 file3 ... ]
  Makes copies in K of file1 [ file2 file3 ... ] (in mJy/beam).
  If file1 is blah.ext the corresponding output will be blah_K.ext.
  Logs to manip.madrout and the standard writedef logs.

  -r: Apply NVSS flux registration factors in $nvssfrff to the K copies.
      Correction factors < 0.1 (i.e. 0 = not ready yet) aren't applied.
      If file1 is LHblah.ext (it is strongly recommended that it start with
      the survey code!), the corresponding output wiil be LHblah_FRK.ext.
};

our($opt_r);
getopts('r');

unless(@ARGV > 0) {
  die $use;
}

foreach my $filn (@ARGV) {
  my $mind  = fileindex($filn);
  my %mdef  = readdef($mind);
  my $k_ind = get_first_slot();

  if($mdef{UNIT} =~ /K/){
    print "Skipping $filn ($mind)...it is already in $mdef{UNIT}.\n";
    next;
  }
  if($mdef{UNIT} !~ /Jy\/beam/){
    print "Cannot handle $mdef{UNIT} in $filn ($mind)...skipping it.\n";
    next;
  }

  unless($mdef{CHF}){
    print "Central Hel. Freq. not found in $filn ($mind)...skipping it.\n";
    next;
  }
  unless($mdef{BMAJ} && $mdef{BMIN}){
    print "Beam not found in $filn ($mind)...skipping it.\n";
    next;
  }

  my $f2beam = $mdef{CHF} * dms2s($mdef{BMAJ});
  $f2beam *= $mdef{CHF} * dms2s($mdef{BMIN});
  my $convfac = 1.2221e12 * $mdef{SCALFACT} / $f2beam;

  if($mdef{UNIT} =~ /mJy/){
    $convfac *= 0.001;
  }

  my $indicator = '.K_';

  if($opt_r){
    my $surv = survey($filn);

    open(FR, $nvssfrff)
       or die "Could not open $nvssfrff to get flux registrations";
    while(<FR>){
      if(/\b$surv ([0-9.]+)/){
	print "$surv $1\n";
	if($1 > 0.1){            # 0 if not ready.
	  $convfac /= $1;        # The "correction" factors are given as
	  $indicator = '.KRF_';  # $f_{DRAO} / f_{NVSS}$.
	}
	last;
      }
    }
    close(FR) or warn "Error closing $nvssfrff";
  }

  # <a href="http://www.perlmonks.org/index.pl?node=japhy">japhy</a> claims
  # this is the fastest way of replacing only the final occurrence.
  my $gnirts = reverse $mdef{FILENAME};
  $gnirts =~ s/\./$indicator/;
  $mdef{FILENAME} = reverse $gnirts;

  $mdef{UNIT}     = 'K (Tb)';
  $mdef{SCALFACT} = 1.0;

  writedef($mind, $k_ind, %mdef);

  manip("f$k_ind = $convfac * f$mind");
}