#!/usr/bin/perl -w

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

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

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

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

  -r: Unapply NVSS flux registration factors in $nvssfrff.
      Correction factors < 0.1 (i.e. 0 = not ready yet) aren't unapplied.
      If file1 is blah.ext the corresponding output wiil be blah_JY.ext.
};

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

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

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

  if($mdef{UNIT} =~ /Jy/){
    print "Skipping $filn ($mind)...it is already in $mdef{UNIT}.\n";
    next;
  }
  if($mdef{UNIT} !~ /K/){
    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.2221e9 / ($f2beam * $mdef{SCALFACT});

  if($mdef{UNIT} =~ /mK/){
    $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(/$surv ([0-9.]+)/){
	print "$surv $1\n";
	if($1 > 0.1){            # 0 if not ready.
	  $convfac /= $1;        # The "correction" factors are given as
	                         # $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/\.(K_)?/.YJ_/;
  $mdef{FILENAME} = reverse $gnirts;

  $mdef{UNIT}     = 'Jy/beam';
  $mdef{SCALFACT} = 0.001;

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

  manip("f$j_ind = f$mind / $convfac");
}