#!/usr/local/bin/perl -w
#
# Usage: cirsi_grid2.pl telescope catalogfn outputfn
#
# you may need to change the perl location. type where perl to find
# where it is.
#
# where telescope is either INT or DUPONT, catalogfn is an input file of 
# field centres with format as described below, and outputfn is the name of 
# the output coordinate file (for use by the TCS).
#
# Catalogue file format:
#
# RA  DEC  FIELDNAME
#
# ie, the catalogue consists of the J2000 RA, DEC of each mosaic image field 
# centre (30'x30' on INT, and 13'x13' on duPont), and the fieldname of that 
# mosaic image, eg:
#
# 10:52:13 57:28:48 LockmanHole
#
# note: the colons are required, and FIELDNAME should have no spaces. 
#
# Output files:
#
# The TCS coordinates of the 4 pointings to make a filled mosaic image are
# written to outputfn, and GAIA format catalogue files (.gaia extension) are
# created for use as finder charts.
#
# Modification history:
#
# 22/5/00 MS  -- created original version
# 11/8/00 CNS -- added tel option, fixed -ve DEC bug, minor changes 
#

die "Usage: cirsi_grid2.pl telescope catalogfn outputfn\n"
    unless ($#ARGV == 2);

if ($ARGV[0] eq "INT") {
    $secpix = 0.457;
} elsif ($ARGV[0] eq "DUPONT") {
    $secpix = 0.197;
} else {
    die "telescope must be INT or DUPONT\n"
}

print "Telescope = $ARGV[0], scale = $secpix [arcsec/pixel]\n";

open(CENTRES, $ARGV[1]) or 
    die "Unable to open coordinate file: $ARGV[1]\n";

open(OUTPUT,"> $ARGV[2]") or
    die "Unable to open TCS output file: $ARGV[2]\n";
    
print OUTPUT "! TCS coordinates - centre of each of 4 pointings\n!\n";

$n = 0;

while (<CENTRES>) {                 # read RA DEC FIELDNAME list
    chomp;

    if (m/^\s*\d/) {
        ($ras[$n], $decs[$n], $fields[$n]) = split(/\s+/, $_);
        $n++;
    }
}

for ($k = 0; $k < $n; $k++) {            # loop over each field centre

    $ra  = $ras[$k];
    $dec = $decs[$k];

    open(GAIA, "> $fields[$k].gaia")
        or die "Unable to open: $fields[$k].gaia\n";

    &convert_ra;
    &convert_dec;                                # convert to decimal

    # Declination step with 5% overlap

    $stepdec = (($secpix * 1024.) - (0.05 * $secpix * 1024.)) / 60.0;

    $width=($secpix *1024)/(2*3600);    # half width of plotted square in gaia
    $stepdecdeg=$stepdec/60; # convert to degrees
    $stepra=($stepdec)/(cos($decdeg*3.141592/180)); # declination correction
    $stepradeg=$stepra/60; # convert to degrees

    # The next bit has the correct tab format for a gaia catalogue
    # do not change or the catalogue won't be read in    

    print GAIA "Query Result\n\n";
    print GAIA "# Config entry for catalogue:\n";
    print GAIA  'symbol: id {square {} {} {} {$id} {}} { ' . $width .
                ' {deg 2000}}' . "\n";
    print GAIA "# End config entry\n\n";
    print GAIA "id\tra\tdec\t\n";
    print GAIA "--\t--\t---\t\n";

    # Now loop round and calculate the 16 chip centres
    # This bit is used to produce the files for the finding chart 

    for ($i=0;$i<4;$i++) {
	for ($j=0;$j<4;$j++) {
	    $id=($i+1).($j+1);
	    $dec=($decdeg-2*$stepdecdeg)+($i*$stepdecdeg+$stepdecdeg/2);
	    $ra=($radeg-2*$stepradeg)+($j*$stepradeg+$stepradeg/2);
	    $ra=$ra/15; # convert ra
	    if ($ra<0) {$ra=24+$ra;} # Trap RA wrapraround
	    &decimalra_to_formatra;
	    &decimaldec_to_formatdec;
	    print GAIA "$id\t$rah $ram $ras\t$decd $decm $decs\n";
	}
    }

    close(GAIA);

# Now loop round and produce the TCS coordinates for each pointing
# There are 4 pointings per mosaic
# the outputted names look like: (16 chips in 4-pattern mosaic)
#
#     N
#    ____ 
#   |    |
# E | DA | W 
#   | CB | 
#   |____|
#   
#     S
#
    for ($i=1;$i>-1;$i--) {
	for ($j=1;$j>-1;$j--) {
	    #Get each id of each pointing correct by appending
	    # a letter according to the pattern above
	    if ($i==1 && $j==1) { $id=$fields[$k].'D';}
	    if ($i==0 && $j==1) { $id=$fields[$k].'C';}
	    if ($i==1 && $j==0) { $id=$fields[$k].'A';}
	    if ($i==0 && $j==0) { $id=$fields[$k].'B';}
	    $dec=($decdeg-$stepdecdeg/2)+($i*$stepdecdeg);
	    $ra=($radeg-$stepradeg/2)+($j*$stepradeg);
	    $ra=$ra/15;
	    if ($ra<0) {$ra=24+$ra;}
	    &decimalra_to_formatra;
	    &decimaldec_to_formatdec;
	    print OUTPUT "$id $rah $ram $ras  $decd $decm $decs  J2000.0\n";
	}
    }

    $ra  = $ras[$k];
    $dec = $decs[$k];
    $id  = $fields[$k].'X';               # write mosaic center position

    &convert_ra;
    &convert_dec;
 
    $ra = $radeg;
    $dec = $decdeg;
    $ra=$ra/15; # convert ra
    if ($ra<0) {$ra=24+$ra;} # Trap RA wrapraround

    &decimalra_to_formatra;
    &decimaldec_to_formatdec;

    print OUTPUT "$id $rah $ram $ras  $decd $decm $decs  J2000.0\n";
    print OUTPUT "!\n";
}

close(OUTPUT);

# Subroutines below for RA-DEC conversion
    
sub convert_ra {
    ($rah,$ram,$ras)=split(':',$ra,3);
    $radeg=$rah*15+($ram/60)*15+($ras/3600)*15;
}

sub convert_dec {
    ($decd,$decm,$decs)=split(':',$dec,3);

    if (substr($decd,0,1)eq'-') {
	
	$decdeg=$decd-($decm/60)-($decs/3600);

    }

    else {

    $decdeg=$decd+($decm/60)+($decs/3600);

    }
}

sub decimalra_to_formatra {
    $rah=(sprintf('%2.2d',(int($ra))));
    $ram=(sprintf('%2.2d',(int(($ra-$rah)*60.01))));
    $ras=(sprintf('%4.2f',(abs(($ra-$rah-($ram/60))*3600))));
    if (substr($ras,1,1) eq '.')  {$ras='0'.$ras;}
}

sub decimaldec_to_formatdec {

    if ($dec >=0) {
        $decd='+'.(sprintf('%2.2d',(int($dec))));
        $decm=(sprintf('%2.2d',(int(($dec-$decd)*60.01))));
        $decs=(sprintf('%3.1f',(abs(($dec-$decd-($decm/60))*3600))));

        if (substr($decs,1,1) eq '.')  {
            $decs='0'.$decs;
        }
    } else {
        $decd=(sprintf('%2.2d',(int($dec))));
        $decm=(sprintf('%2.2d',(abs(int(($dec-$decd)*60.010)))));
        $decs=(sprintf('%3.1f',(abs(($dec-$decd+($decm/60))*3600))));

        if (substr($decs,1,1) eq '.')  {
            $decs='0'.$decs;
        }

        if ($decd eq "00") {
            $decd = '-' . $decd;
        }
    }
}