%# csq=csquare2(beta,f,g) returns cos(pi*beta*(f^2+g^2))
%# beta is equal to lambda*D/l^2, i.e. the Fresnelnumber
%# for the complete spectrum (from 0 to fsample)
%# f_isall = 0
%# only f and g from 0 to fsample/2 should be given, the other part is 
%# calculated using the periodicity in frequency-space
%# f_isall = 1
%# complete frequency range is given
%# allows for astigmatism (5th argument)

function csq=csquare2(betah,f,g,varargin)

f_isall = 0;
switch nargin
case 3
    betav = betah;
case 4
    f_isall = varargin{1};
    betav = betah;
case 5
    f_isall = varargin{1};
    betav = varargin{2};
end

csq=cos((pi*betah)*f.^2+(pi*betav)*g.^2);

if (f_isall~=1)
    [n m]=size(f);
    n=n-1;
    m=m-1;
    csq=cat(2,cat(1,csq,flipdim(csq(2:n,:,:),1)),cat(1,flipdim(csq(:,2:m,:),2),flipdim(flipdim(csq(2:n,2:m,:),1),2)));
end

end
