#!/usr/bin/perl -w

# use: finishmos moscode

# Grab Q, U, and weight mosaics from /igps/delrizzo/mos/p21/,

# make coverage (sqrt(wt)) and mask maps from weight map,

# tmapconvrt SD.MPA21 and SD.MCPI21 to local frames,
#     Can SD.MCPI21 easily account for different noises of 26m and Eff.?

# make corrected polarized intensity (.MCPI21) using madr's manipulator.
# makepolmap does not correct for Ricean bias.

# Fill in edges of ST Q, U, and CPI mosaics with single dish data,
# clean up with rm and neg everything except filled in mosaics,
# ktranslate filled in mosaics to miriad format,
# regrid them to q_a1 frame, options=offset,
# rm -rf original miriad dirs,

# The big mosaics should be made from CPI, Q, and U, not angle.  Unfortunately
# impol, miriad's task for producing polarization angle and intensity maps, can
# only properly correct for Ricean bias when the noise is uniform.  It can of
# course make the angle map out of Q and U, though.


use DRAORed;
use strict;

# Redefine $mind to not have a no-data value (it starts at 0.0000).
sub nonodata {
  my $mind = $_[0];

  my %mdef = altreaddef($mind);
  delete($mdef{NODATA});
  writedef($mind, $mind, %mdef);

  return %mdef;
}

# Returns the theoretical rms noise in K at the center of a mosaic that has
# been multiplied by the square root of its weight mosaic, assuming an rms
# noise of 0.04K for MIJ2 (empirically determined).
sub getmosnoise {
  my $letter = uc($_[0]);
  my $mind   = fileindex("${letter}_Q");
  my %mdef   = readdef($mind);
  my $dec    = 60.0;

  # Get dec.
  open(GC, "|galcon > galcon.log") or die "error opening |galcon > galcon.log";
  print GC "\n";                 # Select epoch of equatorial coordinates [J]:
  print GC "\n";                 # Do you wish to convert galactic to
                                 # equatorial? [Y]:
  print GC "$mdef{REFCX}\n";     # Longitude (give units) [exit]:
  print GC "$mdef{REFCY}\n";     # Latitude (give units) [exit]:

  # writes answer

  print GC "\n";                 # Longitude (give units) [exit]:
  close(GC) or die "error closing |galcon > galcon.log";
  open(GC, "galcon.log") or die "error opening galcon.log";
  while(<GC>){
    my $line = $_;
    if($line =~ /RA = [^,]+, DEC = (\d+)D (\d+)\' ([0-9.]+)\"/){
      $dec = $1 + $2 / 60.0 + $3 / 3600.0;
      last;
    }
  }
  close(GC) or die "error closing galcon.log";

  # MIJ2's dec is 55.2294166666667.
  return 0.04869484928 * sin(0.01745329252 * $dec);
}

# Grab Q, U, and weight mosaics from /igps/delrizzo/mos/p21/,
#my $davedir = "/igps/delrizzo/mos/p21/";  # May have to ln -s this to something
#                                          # shorter, if not already done.
my $davedir = "dm";

foreach my $davelet ( <$davedir/md*> ){
  my $letter = $davelet;                   # Also includes number.
  my $srcd   = "$davelet/abcd";
  my %mosind = ();
  my %sdind  = ();

  $letter =~ s,^$davedir/,,;

  foreach my $stokes ( "q", "u" ) {
    system("lnrget $srcd/$stokes \U${letter}.P21_${stokes}ABCD.MOSAIC ${letter}_$stokes\E");
    # I have to copy so I can modify this file.
    unlink("\U${letter}_$stokes\E") or die "error unlinking \U${letter}_$stokes\E";
    system("cp $srcd/$stokes/\U${letter}.P21_${stokes}ABCD.MOSAIC ${letter}_$stokes\E");
    $mosind{$stokes} = fileindex("\U${letter}_$stokes\E");
    my %mdef = nonodata($mosind{$stokes});

    system("tmapconvrt \USD.M${stokes}G21 ${letter}_$stokes ${letter}_SD_$stokes\E");
    $sdind{$stokes} = fileindex("\U${letter}_SD_$stokes\E");
    %mdef = nonodata($sdind{$stokes});
  }

  my $wtfilna = "\U$letter\E.WEIGHT.ARRAY";
  if(-f "$srcd/u/${wtfilna}.2"){
    $wtfilna = "${wtfilna}.2";
  }
  system("lnrget $srcd/u $wtfilna \U$letter\E_WT");
  # I have to copy so I can modify this file.
  unlink("\U${letter}_WT\E") or die "error unlinking \U${letter}_WT\E";
  system("cp $srcd/u/$wtfilna \U${letter}_WT\E");

# make coverage (sqrt(wt)) and mask maps from weight map,
  my $wtind = fileindex("\U$letter\E_WT");
  my %wtdef = nonodata($wtind);

  my $maskind = get_first_slot();
  $wtdef{FILENAME} = "\U$letter\E_MASK";
  writedef($wtind, $maskind, %wtdef);

  my $clipind = get_first_slot();
  $wtdef{FILENAME} = "\U$letter\E_CPICLIP";
  writedef($wtind, $clipind, %wtdef);

  my $cpiind = get_first_slot();
  %wtdef = altreaddef($mosind{q});
  $wtdef{FILENAME}     = "\U$letter\E_CPI";
  $wtdef{POLARIZATION} = "Corr. Pol. Int."; # Maximum of 16 characters.
  writedef($mosind{q}, $cpiind, %wtdef);

  my $squaredbias = getmosnoise($letter);
  $squaredbias *= 1.2 * $squaredbias;

  open(MADR, "|madr > /dev/null 2>&1") or die "Could not |madr > /dev/null 2>&1";
  print MADR "m\n";
  print MADR "f$maskind = clip(f$wtind, 0.001)\n";
  print MADR "f${wtind}m = f$wtind + 1.0 - f$maskind\n";

  # Fill in edges of ST Q and U mosaics with single dish data,
  print MADR "f$mosind{q}m = f$mosind{q} * f$maskind + f$sdind{q} - f$sdind{q} * f$maskind\n";
  print MADR "f$mosind{u}m = f$mosind{u} * f$maskind + f$sdind{u} - f$sdind{u} * f$maskind\n";

  #     Can SD.MCPI21 easily account for different noises of 26m and Eff.?  I
  #     don't think it even matters at this point, when you consider systematic
  #     errors.  They all (after combination of their respective spatial
  #     frequencies) probably have "noise" of about 0.04K.
#  #        So just make the noise uniform so miriad's impol can handle it, and
#  #        convert the corrected polarization intensity back using maths.
#  print MADR "f$mosind{q}m = f$mosind{q} * f$covind\n";
#  print MADR "f$mosind{u}m = f$mosind{u} * f$covind\n";

  print MADR "f$cpiind = f$mosind{q} * f$mosind{q} + f$mosind{u} * f$mosind{u} - $squaredbias / f$wtind\n";
  print MADR "f$clipind = clip(f$cpiind, 0.000001)\n";
  print MADR "f${cpiind}m = power(f$cpiind * f$clipind, 0.5)\n";

  print MADR "exit\n";  # Exit manip.

  # clean up with rm and neg everything except filled in mosaics,
  # The local SD and mask files are not needed anymore.
  foreach my $mind ( $sdind{q}, $sdind{u}, $maskind ){
    print MADR "neg f$mind\n";
    print MADR "y\n";       # Do you wish to delete the file? [N]: y
    print MADR "y\n";       # Are you sure you wish to delete file(s)? [N]: y
  }

  print MADR "exit\n";
  close(MADR) or die "Error closing |madr > /dev/null 2>&1";



#   foreach my $t ( "q", "u", "cov" ){
#     my $draof = "\U${letter}_$t\E";
#     my $mirf  = "${t}_$letter";

#     system("ln -s $draof ${draof}.drao");
#     system("/usr/local/karma/bin/ktranslate -outtype miriad - ${draof}.drao $mirf");

#     system("/usr/local/miriad/bin/linux/puthd in=$mirf/epoch value=2000.0");
#     system("/usr/local/miriad/bin/linux/regrid in=$mirf out=${mirf}.rg axes=1,2 tin=q_a1 options=offset");

#     system("/bin/rm -rf $mirf ${draof}.drao");
#   }
}


# The big mosaics should be made from CPI, Q, and U, not angle.  Unfortunately
# impol, miriad's task for producing polarization angle and intensity maps, can
# only properly correct for Ricean bias when the noise is uniform.  It can of
# course make the angle map out of Q and U, though.