function s_out = polar(s, delth)
%SAXS/POLAR Convert a SAXS image object to polar coordinates.
%   POLAR(S) unrolls an image from rectangular coordinates to a polar
%   transformation, where x is the distance from the center of the image
%   and y is the angle of rotation, in degrees.  The positive x axis is
%   zero degrees, with angles increasing in a counterclockwise direction.
%
%   POLAR(S, DELTA) specifies the spacing in the y axis (angle of
%   rotation).  The default is 0.5 degrees.
%
%   The resulting image includes data points which do not correspond to the
%   original image.  Those pixels have an intensity value of NaN, a special
%   numeric code that indicates "not a number".
%
%   Since we are interpolating between the original image's pixels we have
%   mostly fractional intensities even if the original image was made up of
%   integer intensities (counts).
%
%   See also SAXS/CENTER, SAXS/IMAGE.
% 

%   
%delr=2;
if isequal(s.type, 'rth')  %image already in polar form.
    s_out = s;
    %But what if a different delth is specified?  Ah, well.
    return
end

if isequal(pixelcenter(s), [0 0])  % image probably not centered
    s = center(s, 'nostatus');
end

if nargin > 1  % optional delta theta argument is here
    if ~ (isnumeric(delth) && isequal(length(delth), 1))
        error('Invalid argument.')
    end
else
    delth =2;  % default delta theta = 2 degree by Karsten March 2005
end

% If calibration data are present, prepare for non-equal x and y calibrations.
if ~ isempty(s.koffset) && ~ isempty(s.wavelength)
    kcaltemp=tan(2*asin(s.kcal*s.wavelength/4/pi))/2*4*pi/s.wavelength;
    p2k = (kcaltemp - s.koffset)./ s.pixelcal;
    p2kaverage=(p2k(2)+p2k(1))/2;
else
    p2kaverage=1;
    p2k=[1 1];
end

% Get image size in pixels.
[ysize xsize] = size(s.pixels);

% Get maximum radius from center.
imgcenter = pixel(s, [0, 0])';
corners = [1 1 xsize xsize; 1 ysize 1 ysize];
rmax = max(sqrt(sum(((corners - repmat(imgcenter, 1, 4)).*repmat([p2k(1)/p2kaverage;p2k(2)/p2kaverage],1,4)).^2)));
if ~ isempty(s.koffset) && ~ isempty(s.wavelength)
    twok=4*pi/s.wavelength;
    rmax=sin(0.5*atan(rmax*p2kaverage/twok*2))*twok/p2kaverage;
    rmax=min(rmax,1/p2kaverage*twok);
end
% Get image center in fractional pixel numbers.
pixc = pixelfrac(s, [0, 0]);
xc = pixc(1);
yc = pixc(2);

% Create arrays containing the x and y pixel coordinates.
[x y] = meshgrid(1:xsize, 1:ysize);

% Create a set of polar coordinates to compute.

delr=rmax/200;
[xip yip] = meshgrid([delr/10:200]*delr,(0:delth:360)*2*pi/360);
%[xip yip] = meshgrid(delr:delr:rmax, (0:delth:360-delth)*2*pi/360);
%[xip yip] = meshgrid(.5:.5:rmax, (0:delth:360-delth)*2*pi/360);

% Get the cartesian equivalents of our polar grid.
if ~ isempty(s.koffset) && ~ isempty(s.wavelength)
    % this takes into account wide angle non-linearity
    Q=xip*p2kaverage;
    phi=yip;
    twok=4*pi/s.wavelength;
    onek=twok/2;
    A=Q.*(Q/twok<1);
    B=(A./A).*A;
    Q_star=onek*tan(2*asin(B/twok))/p2kaverage;
    
    xi=floor(Q_star.*cos(phi));
    yi=floor(Q_star.*sin(phi));
    [xxi yyi] = pol2cart(yip, xip);
% for ii=1:size(xi)(1)
%     for jj=1:size(xi)(2)
%         A(ii,jj)=s(xi,yi)
%     end
% end
%     [xi yi]= 
else
    [xi yi] = pol2cart(yip, xip);
end

% Interpolate our polar coordinate values from the image array.
imgrth = interp2((x - xc)*p2k(1)/p2kaverage, (y - yc)*p2k(2)/p2kaverage, s.pixels, xi, yi,'*nearest');
s_out = s;
%s_out.raw_pixels=s_out.raw_pixels ; %saving original data in s.raw_pixels
s_out.raw_center=[xc yc];
s_out.pixels = imgrth;
s_out.type = 'rth';
s_out.xaxis = [delr rmax];

% If calibration data are present, convert xaxis to K.
if ~ isempty(s_out.koffset) && ~ isempty(s_out.wavelength)
    p2k = (s_out.kcal - s_out.koffset)./ s_out.pixelcal;
    s_out.xaxis = s_out.xaxis .* p2kaverage + s_out.koffset;
end

s_out.yaxis = [0 360];
record = 'unroll to polar coordinates (r/theta)';
s_out.history = strvcat(s_out.history, record);
