function fRc = UWiAtRt(a, atRT, thetas, rotMethod, filter, d)
%IATRT inversion of the attenuated Radon transform
% 
% INPUT ARGS:  a         : attenuation map
%              atRT      : attenuated Radon Transform data
%              thetas    : the angles corresponding to each column of atRT
%              rotMethod : rotation method used by 'rotateImage'
%              filter    : filter type used by 'hilbTrans'
%              d         : cropping frequency used by 'hilbTrans'
%
% OUTPUT ARGS: fRc    : reconstructed map
%
% See also ROTATEIMAGE, HILBTRANS

n = size(a,1); h = 2/n;

% define the reconstruction grid
x = -1+1/n:h:1-1/n;
x = repmat(x, n, 1);

y = (1-1/n:-2/n:-1+1/n).';
y = repmat(y, 1, n);

nr = size(atRT,1);

% definition of h1 and h2 as seen in class

h1 = %%% FILL OUT
h2 = %%% FILL OUT

% construction of g_a slice by slice
ga = zeros(size(h1));
for j=1:size(ga,2),
    ga(:,j) = %%% FILL OUT
end

% weighted backprojection
tempx = zeros(n); 
tempy = zeros(n);

for th=1:length(thetas),
    % compute exp(Da(x,\theta^\perp))
    Da = %%% FILL OUT
    
    % contribution to the integral over S1
    proj = ga(:,th);
    t = x.*cos(thetas(th)) + y.*sin(thetas(th));
    t = (t-(-nr/n+1/n))*n/2+1;
    ft = floor(t);
    
    tempx = tempx + cos(thetas(th))*Da.*((t-ft).*proj(ft+1) + (ft+1-t).*proj(ft));
    tempy = tempy + sin(thetas(th))*Da.*((t-ft).*proj(ft+1) + (ft+1-t).*proj(ft));
end

% rescaling constants
tempx = tempx*2*pi/length(thetas)/(4*pi);
tempy = tempy*2*pi/length(thetas)/(4*pi);

% figure; 
% subplot(1,2,1); dispmat(tempx, 'tempx');
% subplot(1,2,2); dispmat(tempy, 'tempy');

% take the divergence using forward derivatives
dxtempx = %%% FILL OUT
dytempy = %%% FILL OUT

fRc = dxtempx + dytempy;

%--- take out the part outside the unit disk
fRc = fRc.*double(x.^2+y.^2<=0.95);

end
