#!/usr/bin/perl


# load packages

use PGPLOT;
use POSIX;


# set parameters

$file = "table.txt";


# read data

open (FILE, $file) or die "Cannot open $file: $!";
while (<FILE>) {
   next if /^$/;
   next if /^#/;
   ($name[$i], $caid[$i], $band[$i], $f_pred[$i], $s_pred[$i], $f_meas,
      $s_meas, $dce_time[$i], $sp_type[$i], $cal_factor) = split;
   ($f_meas[$i], $f_err[$i]) = split ('&plusmn;', $f_meas);
   ($s_meas[$i], $s_err[$i]) = split ('&plusmn;', $s_meas);
   ($cal_factor[$i], $cal_err[$i]) = split ('&plusmn;', $cal_factor);
   $f_pred[$i] = log10 ($f_pred[$i]);
   $s_pred[$i] = log10 ($s_pred[$i]);
   $temp[$i] = log10 (&temp ($sp_type[$i]));
   $i++;
}
close (FILE);


# compute range of data

$min70 = 1e9;  $max70 = 0;
$min160 = 1e9;  $max160 = 0;
$fmin = 1e9;  $fmax = -1e9;
$smin = 1e9;  $smax = -1e9;
$tmin = 1e9;  $tmax = -1e9;

for ($i = 0; $i <= $#cal_factor; $i++) {
   $low = $cal_factor[$i] - $cal_err[$i];
   $high = $cal_factor[$i] + $cal_err[$i];
   if ($band[$i] == 70) {
      if ($low < $min70) {$min70 = $low;}
      if ($high > $max70) {$max70 = $high;}
   }
   if ($band[$i] == 160) {
      if ($low < $min160) {$min160 = $low;}
      if ($high > $max160) {$max160 = $high;}
   }
   if ($f_pred[$i] < $fmin) {$fmin = $f_pred[$i];}
   if ($f_pred[$i] > $fmax) {$fmax = $f_pred[$i];}
   if ($s_pred[$i] < $smin) {$smin = $s_pred[$i];}
   if ($s_pred[$i] > $smax) {$smax = $s_pred[$i];}
   if ($temp[$i] < $tmin) {$tmin = $temp[$i];}
   if ($temp[$i] > $tmax) {$tmax = $temp[$i];}
}

$fmin -= 0.05 * ($fmax - $fmin);  $fmax += 0.05 * ($fmax - $fmin);
$smin -= 0.05 * ($smax - $smin);  $smax += 0.05 * ($smax - $smin);
$tmin -= 0.05 * ($tmax - $tmin);  $tmax += 0.05 * ($tmax - $tmin);
$min70 -= 0.05 * ($max70 - $min70);  $max70 += 0.05 * ($max70 - $min70);
$min160 -= 0.05 * ($max160 - $min160);  $max160 += 0.05 * ($max160 - $min160);


# initialize plot

pgopen ('?');
pgsubp (2, 2);
pgscf (1);
pgsch (2);
pgslw (3);
$lab70 = "\\gmJy/arcsec\\u2\\d/MIPS\\d70\\uUnit";
$lab160 = "\\gmJy/arcsec\\u2\\d/MIPS\\d160\\uUnit";


# plot calibration factor vs. flux (linearity)

pgsci (1);
pgenv ($fmin, $fmax, $min70, $max70, 0, -2);
pgbox (BCLNTS, 0, 0, "", 0, 0);
pglab ("Flux Density (Jy)", "", "");
pgsci (4);
pgbox ("", 0, 0, BNTS, 0, 0);
pgmtxt (L, 2.5, 0.5, 0.5, $lab70);
pgmtxt (T, 1, 0.4, 1, chr(4) . " = 70\\gmm");
for ($i = 0; $i <= $#name; $i++) {
   if ($band[$i] == 70) {
      pgpt (1, $f_pred[$i], $cal_factor[$i], 4);
      pgerrb (6, 1, $f_pred[$i], $cal_factor[$i], $cal_err[$i], 1);
   }
}
pgswin ($fmin, $fmax, $min160, $max160);
pgsci (2);
pgbox ("", 0, 0, CMTS, 0, 0);
pgmtxt (R, 2.5, 0.5, 0.5, $lab160);
pgmtxt (T, 1, 0.6, 0, chr(13) . " = 160\\gmm");
for ($i = 0; $i <= $#name; $i++) {
   if ($band[$i] == 160) {
      pgpt (1, $f_pred[$i], $cal_factor[$i], 13);
      pgerrb (6, 1, $f_pred[$i], $cal_factor[$i], $cal_err[$i], 1);
   }
}


# plot calibration factor vs. predicted sky (background effects)

pgsci (1);
pgenv ($smin, $smax, $min70, $max70, 0, -2);
pgbox (BCLNTS, 0, 0, "", 0, 0);
pglab ("Sky Background (MJy/sr)", "", "");
pgsci (4);
pgbox ("", 0, 0, BNTS, 0, 0);
pgmtxt (L, 2.5, 0.5, 0.5, $lab70);
pgmtxt (T, 1, 0.4, 1, chr(4) . " = 70\\gmm");
for ($i = 0; $i <= $#name; $i++) {
   if ($band[$i] == 70) {
      pgpt (1, $s_pred[$i], $cal_factor[$i], 4);
      pgerrb (6, 1, $s_pred[$i], $cal_factor[$i], $cal_err[$i], 1);
   }
}
pgswin ($smin, $smax, $min160, $max160);
pgsci (2);
pgbox ("", 0, 0, CMTS, 0, 0);
pgmtxt (R, 2.5, 0.5, 0.5, $lab160);
pgmtxt (T, 1, 0.6, 0, chr(13) . " = 160\\gmm");
for ($i = 0; $i <= $#name; $i++) {
   if ($band[$i] == 160) {
      pgpt (1, $s_pred[$i], $cal_factor[$i], 13);
      pgerrb (6, 1, $s_pred[$i], $cal_factor[$i], $cal_err[$i], 1);
   }
}


# plot calibration factor vs. DCE time (integration time effects)

pgsci (1);
pgenv (2, 11, $min70, $max70, 0, -2);
pgbox (BCNTS, 0, 0, "", 0, 0);
pglab ("DCE time (MIPS sec.)", "", "");
pgsci (4);
pgbox ("", 0, 0, BNTS, 0, 0);
pgmtxt (L, 2.5, 0.5, 0.5, $lab70);
pgmtxt (T, 1, 0.4, 1, chr(4) . " = 70\\gmm");
for ($i = 0; $i <= $#name; $i++) {
   if ($band[$i] == 70) {
      pgpt (1, $dce_time[$i], $cal_factor[$i], 4);
      pgerrb (6, 1, $dce_time[$i], $cal_factor[$i], $cal_err[$i], 1);
   }
}
pgswin (2, 11, $min160, $max160);
pgsci (2);
pgbox ("", 0, 0, CMTS, 0, 0);
pgmtxt (R, 2.5, 0.5, 0.5, $lab160);
pgmtxt (T, 1, 0.6, 0, chr(13) . " = 160\\gmm");
for ($i = 0; $i <= $#name; $i++) {
   if ($band[$i] == 160) {
      pgpt (1, $dce_time[$i], $cal_factor[$i], 13);
      pgerrb (6, 1, $dce_time[$i], $cal_factor[$i], $cal_err[$i], 1);
   }
}


# plot calibration factor vs. temperature (spectral type effects)

pgsci (1);
pgenv ($tmin, $tmax, $min70, $max70, 0, -2);
pgbox (BCLNTS, 0, 0, "", 0, 0);
pglab ("T\\deff\\u (K)", "", "");
pgsci (4);
pgbox ("", 0, 0, BNTS, 0, 0);
pgmtxt (L, 2.5, 0.5, 0.5, $lab70);
pgmtxt (T, 1, 0.4, 1, chr(4) . " = 70\\gmm");
for ($i = 0; $i <= $#name; $i++) {
   if ($band[$i] == 70) {
      pgpt (1, $temp[$i], $cal_factor[$i], 4);
      pgerrb (6, 1, $temp[$i], $cal_factor[$i], $cal_err[$i], 1);
   }
}
pgswin ($tmin, $tmax, $min160, $max160);
pgsci (2);
pgbox ("", 0, 0, CMTS, 0, 0);
pgmtxt (R, 2.5, 0.5, 0.5, $lab160);
pgmtxt (T, 1, 0.6, 0, chr(13) . " = 160\\gmm");
for ($i = 0; $i <= $#name; $i++) {
   if ($band[$i] == 160) {
      pgpt (1, $temp[$i], $cal_factor[$i], 13);
      pgerrb (6, 1, $temp[$i], $cal_factor[$i], $cal_err[$i], 1);
   }
}


# close plot

pgend;




# SUBROUTINES

# convert spectral type to effective temperature using tables from Lang's
# "Astrophysical Data"

sub temp {

   my $sp_type = $_[0];
   my ($i, $temp, $type, $num, $lum);

   my @l1 = (
      ["O", 3, 47300],
      ["O", 4, 44100],
      ["O", 5, 40300],
      ["O", 6, 39000],
      ["O", 7, 35700],
      ["O", 8, 34200],
      ["O", 9, 32600],
      ["B", 0, 26000],
      ["B", 1, 20800],
      ["B", 2, 18500],
      ["B", 3, 16200],
      ["B", 5, 13600],
      ["B", 6, 13000],
      ["B", 7, 12200],
      ["B", 8, 11200],
      ["B", 9, 10300],
      ["A", 0, 9730],
      ["A", 1, 9230],
      ["A", 2, 9080],
      ["A", 3, 8770],
      ["A", 5, 8510],
      ["A", 7, 8150],
      ["A", 8, 7950],
      ["F", 0, 7700],
      ["F", 2, 7350],
      ["F", 5, 6900],
      ["F", 8, 6100],
      ["G", 0, 5550],
      ["G", 2, 5200],
      ["G", 5, 4850],
      ["G", 8, 4600],
      ["K", 0, 4420],
      ["K", 1, 4330],
      ["K", 2, 4250],
      ["K", 3, 4080],
      ["K", 4, 3950],
      ["K", 5, 3850],
      ["K", 7, 3700],
      ["M", 0, 3650],
      ["M", 1, 3550],
      ["M", 2, 3450],
      ["M", 3, 3200],
      ["M", 4, 2980],
      ["M", 5, 2800],
      ["M", 6, 2600]
   );

   my @l3 = (
      ["O", 3, 50000],
      ["O", 4, 45500],
      ["O", 5, 42500],
      ["O", 6, 39500],
      ["O", 7, 37000],
      ["O", 8, 34700],
      ["O", 9, 32000],
      ["B", 0, 29000],
      ["B", 1, 24000],
      ["B", 2, 20300],
      ["B", 3, 17100],
      ["B", 5, 15000],
      ["B", 6, 14100],
      ["B", 7, 13200],
      ["B", 8, 12400],
      ["B", 9, 11000],
      ["A", 0, 10100],
      ["A", 1, 9480],
      ["A", 2, 9000],
      ["A", 3, 8600],
      ["A", 5, 8100],
      ["A", 7, 7650],
      ["A", 8, 7450],
      ["F", 0, 7150],
      ["F", 2, 6870],
      ["F", 5, 6470],
      ["F", 8, 6150],
      ["G", 0, 5850],
      ["G", 2, 5450],
      ["G", 5, 5150],
      ["G", 8, 4900],
      ["K", 0, 4750],
      ["K", 1, 4600],
      ["K", 2, 4420],
      ["K", 3, 4200],
      ["K", 4, 4000],
      ["K", 5, 3950],
      ["K", 7, 3850],
      ["M", 0, 3800],
      ["M", 1, 3720],
      ["M", 2, 3620],
      ["M", 3, 3530],
      ["M", 4, 3430],
      ["M", 5, 3330],
      ["M", 6, 3240]
   );

   my @l5 = (
      ["O", 3, 52500],
      ["O", 4, 48000],
      ["O", 5, 44500],
      ["O", 6, 41000],
      ["O", 7, 38000],
      ["O", 8, 35800],
      ["O", 9, 33000],
      ["B", 0, 30000],
      ["B", 1, 25400],
      ["B", 2, 22000],
      ["B", 3, 18700],
      ["B", 5, 15400],
      ["B", 6, 14000],
      ["B", 7, 13000],
      ["B", 8, 11900],
      ["B", 9, 10500],
      ["A", 0, 9520],
      ["A", 1, 9230],
      ["A", 2, 8970],
      ["A", 3, 8720],
      ["A", 5, 8200],
      ["A", 7, 7850],
      ["A", 8, 7580],
      ["F", 0, 7200],
      ["F", 2, 6890],
      ["F", 5, 6440],
      ["F", 8, 6200],
      ["G", 0, 6030],
      ["G", 2, 5860],
      ["G", 5, 5770],
      ["G", 8, 5570],
      ["K", 0, 5250],
      ["K", 1, 5080],
      ["K", 2, 4900],
      ["K", 3, 4730],
      ["K", 4, 4590],
      ["K", 5, 4350],
      ["K", 7, 4060],
      ["M", 0, 3850],
      ["M", 1, 3720],
      ["M", 2, 3580],
      ["M", 3, 3470],
      ["M", 4, 3370],
      ["M", 5, 3240],
      ["M", 6, 3050],
      ["M", 7, 2940],
      ["M", 8, 2640]
   );

   $type = substr ($sp_type, 0, 1);
   $num = substr ($sp_type, 1, 1);

   if ($sp_type =~ "I") {$lum = \@l1;}
   if ($sp_type =~ "III") {$lum = \@l3;}
   if ($sp_type =~ "V") {$lum = \@l5;}
   if (!defined ($lum)) {
      print STDERR "Could not find luminosity class for $sp_type.\n";
   }

   for ($i = 0; $i <= $#{$lum}; $i++) {
      if (${$lum}[$i][0] eq $type) {
         if (${$lum}[$i][1] >= $num) {
            $temp = ${$lum}[$i-1][2] + (${$lum}[$i][2] - ${$lum}[$i-1][2]) *
               ($num - ${$lum}[$i-1][1]) / (${$lum}[$i][1] - ${$lum}[$i-1][1]);
            last;
         }
         if ((${$lum}[$i][1] < $num) && (${$lum}[$i][0] ne ${$lum}[$i+1][0])) {
            $temp = (${$lum}[$i][2] + ${$lum}[$i+1][2]) / 2;
            last;
         }
      }
   }

   if (!defined ($temp)) {
      print STDERR "Could not compute temperature for $sp_type.\n";
   }

   return ($temp);

}
