% dsspos: Convert DSS (Digital Sky Survey) pixel positions into RA and DEC.
%
% [RA,Dec]=dsspos(xpix,ypix,header);
%
% ARGUMENTS
%  xpix		X pixel position of target
%  ypix		Y pixel position of target
%  header	Header of FITS file, as returned by fitsload().
%
% RETURNS
%  [RA, Dec]	The RA and Declination of the target.
%
% SEE ALSO
%  fitsload

% AUTHOR: Eric Tittley
%  Routine adapted from the WCSLib file, dsspos.c
%
%  File saoimage/wcslib/dsspos.c
%  October 21, 1999
%  By Doug Mink, Harvard-Smithsonian Center for Astrophysics
%
%  Module:     dsspos.c (Plate solution WCS conversion)
%  Purpose:    Compute WCS from Digital Sky Survey plate fit
%  Subroutine: dsspos() converts from pixel location to RA,Dec 
%  Subroutine: dsspix() converts from RA,Dec to pixel location
%
%   These functions are based on the astrmcal.c portion of GETIMAGE by
%   J. Doggett and the documentation distributed with the Digital Sky Survey.
%
% HISTORY
%  00 04 28: Conversion done
%  06 03 01 Standardized header comments.

function [RA,dec]=dsspos(xpix,ypix,header)

  cond2r = 1.745329252e-2;
  cons2r = 206264.8062470964;
  twopi = 6.28318530717959;

% Convert from image pixels to plate pixels
  x_pixel_offset=fitsfindkeyvalue(header,'CNPIX1');
  y_pixel_offset=fitsfindkeyvalue(header,'CNPIX2');
  x = xpix + x_pixel_offset - 1.0 + 0.5;
  y = ypix + y_pixel_offset - 1.0 + 0.5;

% Convert from pixels to millimeters
  ppo_coeff=zeros(1,6);
  for i=1:6
   ppo_coeff(i)=fitsfindkeyvalue(header,['PPO',int2str(i)]);
  end
  x_pixel_size=fitsfindkeyvalue(header,'XPIXELSZ');
  y_pixel_size=fitsfindkeyvalue(header,'YPIXELSZ');
  xmm = (ppo_coeff(3) - x * x_pixel_size) / 1000.0;
  ymm = (y * y_pixel_size - ppo_coeff(6)) / 1000.0;
  xmm2 = xmm .* xmm;
  ymm2 = ymm .* ymm;
  xmm3 = xmm .* xmm2;
  ymm3 = ymm .* ymm2;
  x2y2 = xmm2 + ymm2;

%  Compute coordinates from x,y and plate model */
  x_coeff=zeros(1,12);
  y_coeff=x_coeff;
  for i=1:20
   x_coeff(i)=fitsfindkeyvalue(header,['AMDX',int2str(i)]);
   y_coeff(i)=fitsfindkeyvalue(header,['AMDY',int2str(i)]);
  end

  xi =  x_coeff(1) *xmm       + x_coeff(2) *ymm + ...
	x_coeff(3)            + x_coeff(4) *xmm2 + ...
	x_coeff(5) *xmm.*ymm  + x_coeff(6) *ymm2 + ...
	x_coeff(7) *(x2y2)    + x_coeff(8) *xmm3 + ...
	x_coeff(9) *xmm2.*ymm + x_coeff(10)*xmm.*ymm2 + ...
	x_coeff(11)*ymm3      + x_coeff(12)*xmm.*x2y2 + ...
	x_coeff(13)*xmm.*x2y2.*x2y2;

  eta =	y_coeff(1) *ymm       + y_coeff(2) *xmm + ...
	y_coeff(3)            + y_coeff(4) *ymm2 + ...
	y_coeff(5) *xmm.*ymm  + y_coeff(6) *xmm2 + ...
	y_coeff(7) *x2y2      + y_coeff(8) *ymm3 + ...
	y_coeff(9) *ymm2.*xmm + y_coeff(10)*ymm.*xmm2 + ...
	y_coeff(11)*xmm3      + y_coeff(12)*ymm.*x2y2 + ...
	y_coeff(13)*ymm.*x2y2.*x2y2;

% Convert to radians

  xir = xi / cons2r;
  etar = eta / cons2r;

% Convert to RA and Dec
  PlateDecSign=fitsfindkeystring(header,'PLTDECSN');
  if(PlateDecSign(1)=='-')
   PlateDecSign=-1;
  else
   PlateDecSign=1;
  end
  PlateDecD=fitsfindkeyvalue(header,'PLTDECD');
  PlateDecM=fitsfindkeyvalue(header,'PLTDECM');
  PlateDecS=fitsfindkeyvalue(header,'PLTDECS');
  plate_dec=PlateDecSign*dms2deg(PlateDecD,PlateDecM,PlateDecS)/180*pi;

  PlateRAH=fitsfindkeyvalue(header,'PLTRAH');
  PlateRAM=fitsfindkeyvalue(header,'PLTRAM');
  PlateRAS=fitsfindkeyvalue(header,'PLTRAS');
  plate_ra=dms2deg(PlateRAH,PlateRAM,PlateRAS)/12*pi;
 
  ctan = tan(plate_dec);
  ccos = cos(plate_dec);

  raoff = atan2(xir / ccos, 1.0 - etar * ctan);
  ra = raoff + plate_ra;
  if (ra < 0.0), ra = ra + twopi; end
  RA = ra / cond2r;

  dec = atan (cos (raoff) .* ((etar + ctan) ./ (1.0 - (etar * ctan))));
  Dec = dec / cond2r;