Index: trunk/ippScripts/scripts/ipp_mops_translate.pl =================================================================== --- trunk/ippScripts/scripts/ipp_mops_translate.pl (revision 21201) +++ trunk/ippScripts/scripts/ipp_mops_translate.pl (revision 24190) @@ -15,7 +15,5 @@ use constant EXTNAME => 'MOPS_TRANSIENT_DETECTIONS'; # Extension name for output table -use constant POSITION_ERROR => 0.2 / 3600; # Error in positions, degrees use constant ZERO_POINT => 25; # Magnitude zero point -use constant MAG_ERROR => 0.1; # Error in magnitudes use constant OBSERVATORY_CODE => 566; # IAU Observatory Code use constant FAKE_LIMITING_MAG => 23.0; # Fake limiting magnitude to report @@ -25,5 +23,4 @@ $extname, # Name of extension containing photometry $skycell, # Skycell file with WCS - $diff_image_id, # difference image id $output, # Name of output file $save_temps, # Save temporary files? @@ -33,17 +30,13 @@ 'input=s' => \$input, 'extname=s' => \$extname, - 'skycell=s' => \$skycell, 'output=s' => \$output, - 'diff_image_id=s' => \$diff_image_id, 'save-temps' => \$save_temps, ) or pod2usage( 2 ); pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV; -pod2usage( -msg => "Required options: --input --extname --skycell --diff_image_id --output", +pod2usage( -msg => "Required options: --input --extname --output", -exitval => 3) unless defined $input and defined $extname - and defined $skycell - and defined $diff_image_id, and defined $output; @@ -64,7 +57,7 @@ # Header translation table my $headers = { + # IPP name => MOPS name, type, comment 'FPA.RA' => { name => 'RA', type => TSTRING, comment => 'Right ascension of boresight'}, 'FPA.DEC' => { name => 'DEC', type => TSTRING, comment => 'Declination of boresight' }, - 'FPA.OBS' => { name => 'FPA_ID', type => TSTRING, comment => 'Exposure identifier' }, 'FPA.FILTER' => { name => 'FILTER', type => TSTRING, comment => 'Filter name' }, 'EXPTIME' => { name => 'EXPTIME', type => TDOUBLE, comment => 'Exposure time' }, @@ -72,87 +65,48 @@ 'FPA.ALT' => { name => 'TEL_ALT', type => TDOUBLE, comment => 'Telescope altitude' }, 'FPA.AZ' => { name => 'TEL_AZ', type => TDOUBLE, comment => 'Telescope azimuth' }, - }; - - -# Read the skycell header for the WCS -my $skycellResolved = $ipprc->file_resolve( $skycell ); -my $status = 0; -my $skyfile_fits = Astro::FITS::CFITSIO::open_file( $skycellResolved, READONLY, $status ); - die("failed to open skycell file: $skycellResolved: $status") if $status; - -my $wcsheader; -($wcsheader, $status) = Astro::FITS::CFITSIO::fits_read_header( $skyfile_fits ); -die("Unable to read skycell header: $status") if $status; - -# Get the useful header keywords -my $naxis1 = $$wcsheader{'NAXIS1'}; -my $naxis2; -if ($naxis1) { - $naxis2 = $$wcsheader{'NAXIS2'} or die("Can't find NAXIS2"); -} else { - # if the skyfile is compressed then the WCS won't be in the primary header - my $hdutype; - $skyfile_fits->movrel_hdu(1, $hdutype, $status); - die("Unable to movrel_hdu: $status") if $status; - - ($wcsheader, $status) = Astro::FITS::CFITSIO::fits_read_header( $skyfile_fits ); - die("Unable to read extension header: $status") if $status; - my $xtension = $$wcsheader{'XTENSION'} or die("Can't find XTENSION"); - die("XTENSION found: $xtension") if $xtension ne "'BINTABLE'"; - $naxis1 = $$wcsheader{'ZNAXIS1'} or die("Can't find ZNAXIS1"); - $naxis2 = $$wcsheader{'ZNAXIS2'} or die("Can't find ZNAXIS2"); -} - -my $ctype1 = $$wcsheader{'CTYPE1'} or die("Can't find CTYPE1"); -my $ctype2 = $$wcsheader{'CTYPE2'} or die("Can't find CTYPE2"); -my $cdelt1 = $$wcsheader{'CDELT1'} or die("Can't find CDELT1"); -my $cdelt2 = $$wcsheader{'CDELT2'} or die("Can't find CDELT2"); -my $crval1 = deg2rad($$wcsheader{'CRVAL1'}) or die("Can't find CRVAL1"); -my $crval2 = deg2rad($$wcsheader{'CRVAL2'}) or die("Can't find CRVAL2"); -my $crpix1 = $$wcsheader{'CRPIX1'} or die("Can't find CRPIX1"); -my $crpix2 = $$wcsheader{'CRPIX2'} or die("Can't find CRPIX2"); -my $pc11 = $$wcsheader{'PC001001'} or die("Can't find PC001001"); -my $pc12 = $$wcsheader{'PC001002'} or die("Can't find PC001002"); -my $pc21 = $$wcsheader{'PC002001'} or die("Can't find PC002001"); -my $pc22 = $$wcsheader{'PC002002'} or die("Can't find PC002002"); -my $crota1 = $$wcsheader{'CROTA1'}; -my $crota2 = $$wcsheader{'CROTA2'}; - -die("Unexpected projection: $ctype1 and $ctype2.") unless $ctype1 =~ /^\'RA---TAN\s*\'$/ and - $ctype2 =~ /^\'DEC--TAN\s*\'$/; -die("Can't determine size of skycell ($naxis1,$naxis2)") unless $naxis1 > 0 and $naxis2 > 0; -die("Can't determine scale of skycell ($cdelt1,$cdelt2)") if not defined $cdelt1 or $cdelt1 == 0 or - not defined $cdelt2 or $cdelt2 == 0; -die("We don't know how to handle rotations ($crota1,$crota2)") if defined $crota1 or defined $crota2; - - -# Read the input table -my $inputResolved = $ipprc->file_resolve($input); + 'IMAGEID' => { name => 'DIFFIMID', type => TINT, comment => 'Difference image identifier'}, + 'FPA.OBS' => { name => 'FPA_ID', type => TSTRING, comment => 'Exposure identifier' }, +}; + + +# Read header +my $inputResolved = $ipprc->file_resolve($input); # Resolved filename +my $status = 0; # CFITSIO status my $inFits = Astro::FITS::CFITSIO::open_file( $inputResolved, READONLY, $status ); # FITS file handle +die("failed to open input file: $inputResolved: $status") if $status; my $inHeader = $inFits->read_header(); # Header for input check_fitsio($status); + +# Plate scales +my $cdelt1 = $$inHeader{'CDELT1'} or die("Can't find CDELT1"); +my $cdelt2 = $$inHeader{'CDELT2'} or die("Can't find CDELT2"); + +# Read table data $inFits->movnam_hdu(BINARY_TBL, $extname, 0, $status) and check_fitsio($status); my $numRows; # Number of rows in table $inFits->get_num_rows($numRows, $status) and check_fitsio($status); -my ($xCol, $yCol, $magCol, $ensCol); # Column numbers for x and y -$inFits->get_colnum(0, 'X_PSF', $xCol, $status) and check_fitsio($status); -$inFits->get_colnum(0, 'Y_PSF', $yCol, $status) and check_fitsio($status); -$inFits->get_colnum(0, 'PSF_INST_MAG', $magCol, $status) and check_fitsio($status); -$inFits->get_colnum(0, 'EXT_NSIGMA', $ensCol, $status) and check_fitsio($status); - -my ($xType, $yType, $magType, $ensType); # Column types -$inFits->get_coltype($xCol, $xType, undef, undef, $status) and check_fitsio($status); -$inFits->get_coltype($yCol, $yType, undef, undef, $status) and check_fitsio($status); -$inFits->get_coltype($magCol, $magType, undef, undef, $status) and check_fitsio($status); -$inFits->get_coltype($ensCol, $ensType, undef, undef, $status) and check_fitsio($status); - -my ($x, $y, $mag, $ens); # Sources arrays -$inFits->read_col($xType, $xCol, 1, 1, $numRows, 0, $x, undef, $status) and check_fitsio($status); -$inFits->read_col($yType, $yCol, 1, 1, $numRows, 0, $y, undef, $status) and check_fitsio($status); -$inFits->read_col($magType, $magCol, 1, 1, $numRows, 0, $mag, undef, $status) and check_fitsio($status); -$inFits->read_col($ensType, $ensCol, 1, 1, $numRows, 0, $ens, undef, $status) and check_fitsio($status); +my $ra = column($inFits, 'RA_PSF', $numRows); # Right Ascension, degrees +my $dec = column($inFits, 'DEC_PSF', $numRows); # Declination, degrees +my $mag = column($inFits, 'PSF_INST_MAG', $numRows); # Magnitude +my $magErr = column($inFits, 'PSF_INST_MAG_SIG', $numRows); # Magnitude error +my $ens = column($inFits, 'EXT_NSIGMA', $numRows); # Significance of extension +my $xErr = column($inFits, 'X_PSF_SIG', $numRows); # Error in x position, pixels +my $yErr = column($inFits, 'Y_PSF_SIG', $numRows); # Error in y position, pixels +my $scale = column($inFits, 'PLTSCALE', $numRows); # Plate scale, arcsec/pixel +my $angle = column($inFits, 'POSANGLE', $numRows); # Position angle, degrees $inFits->close_file( $status ); + +my ($raErr, $decErr); # Error in ra, dec +for (my $i = 0; $i < $numRows; $i++) { + my $cosPA = cos($$angle[$i]); + my $sinPA = sin($$angle[$i]); + + # XXX Not sure about the transformation here --- check + $$raErr[$i] = $$scale[$i] * ($cosPA * $xErr + $sinPA * $yErr) / 3600; + $$decErr[$i] = $$scale[$i] * ($sinPA * $xErr - $cosPA * $yErr) / 3600; +} + # Parse the list of columns @@ -169,12 +123,10 @@ # Convert the input data into the output formats for (my $i = 0; $i < $numRows; $i++) { - my ($ra, $dec) = pixels_to_sky($$x[$i], $$y[$i]); # Sky coordinates - push @{$colData{'RA_DEG'}}, $ra; push @{$colData{'DEC_DEG'}}, $dec; - push @{$colData{'RA_SIG'}}, POSITION_ERROR; - push @{$colData{'DEC_SIG'}}, POSITION_ERROR; + push @{$colData{'RA_SIG'}}, $raErr; + push @{$colData{'DEC_SIG'}}, $decErr; push @{$colData{'FLUX'}}, $$mag[$i] + ZERO_POINT; - push @{$colData{'FLUX_SIG'}}, MAG_ERROR; + push @{$colData{'FLUX_SIG'}}, $magErr; push @{$colData{'ANG'}}, 0.0; push @{$colData{'ANG_SIG'}}, 0.0; @@ -211,9 +163,4 @@ } -# set DIFFIMID -# XXX TODO: we could put the diff_image_id in the .cmf file and get it from there -$outFits->write_key( TLONGLONG, 'DIFFIMID', $diff_image_id, 'Difference image identifier', $status); -check_fitsio( $status ); - # Adjust the time from start to mid-point { @@ -260,4 +207,24 @@ +# Read a column specified by name +sub column +{ + my $fits = shift; # FITS file + my $name = shift; # Name of column + my $rows = shift; # Number of rows + my $status = 0; # Status of CFITSIO + + my $num; # Column number + $fits->get_colnum(0, $name, $num, $status) and check_fitsio($status); + my $type; # Type of column + $inFits->get_coltype($num, $type, undef, undef, $status) and check_fitsio($status); + my $data; # Array with data + $inFits->read_col($type, $num, 1, 1, $rows, 0, $data, undef, $status) and check_fitsio($status); + + return $data; +} + + + # From Astro::FITS::CFITSIO demo @@ -273,40 +240,4 @@ } -# Convert pixel coordinates to sky coordinates -sub pixels_to_sky -{ - my ($x, $y) = @_; # Coordinates - - # Pixel coordinate relative to reference - $x -= $crpix1; - $y -= $crpix2; - - # Coordinates on tangent plane - my $xi = $pc11 * ($x) + $pc12 * ($y); - my $eta = $pc21 * ($x) + $pc22 * ($y); - $xi *= $cdelt1; - $eta *= $cdelt2; - - # Coordinates on rotated celestial sphere - my $phi = atan2($eta,$xi) + pi/2; - my $theta = atan(180 / pi / sqrt($xi**2 + $eta**2)); - - # Coordinates on celestial sphere - my $ra = $crval1 + atan2(cos($theta) * sin($phi), - sin($theta) * cos($crval2) + cos($theta) * sin($crval2) * cos($phi)); - my $dec = asin(sin($theta) * sin($crval2) - cos($theta) * cos($crval2) * cos($phi)); - - $ra = rad2deg($ra); - $dec = rad2deg($dec); - - $ra += 360 if $ra < 0; - $dec += 360 if $dec < 0; - $ra -= 360 if $ra >= 360; - $dec -= 360 if $dec >= 360; - - return ($ra, $dec); -} - - __END__