function [ Alt Az ] = xy2AltAz(x,y,x0,y0)
%xy2AltAz computes alt and azimuth from pixel x,y with x0 y0 at zenith
% this is for a Sigma 15mm lens and a Canon 5D Mark III camera
% Christopher STubbs Jan 16 2014. 

% Azoffset is amount y pixel axis is to the E of N, in degrees
% changed Jan 17 2014 to use empirical values

Azoffset=-4.4;

% default values, can override with additional argument if desired
if nargin == 2
        x0=1408;
        y0=924;
end

% For no distortion, the radius in pixels is proportional to sin(zenith/2).

% Let's assume the lens goes to 90 deg of zenith angle at the diagonal 
% corner, so that is a distance of (monochromatic pixels= 2888 x 1924)
% sqrt(1444^2+962^2)=1735 pixels. So rpix=alpha*sind(zenith/2), 
% and 1735=alpha*sind(90/2) so alpha=1735/.707=2454. 

% empirical value is 2474.4

% radial distance in pixels, away from zenith
rpix=sqrt((x-x0).^2+(y-y0).^2);

% convert to zenith angle, in degrees
zenith=abs(2*asind(rpix./2474.4));

% convert from zenith to altitude, in degrees
Alt=90-zenith;

% now figure out the azimuth angle, in radians
PixAz1=atan2(-(x-x0),(y-y0));

% go to degrees
PixAz=180*PixAz1/pi;

% account for the offset
Az=PixAz-Azoffset;
Az=mod(Az,360);


end
