--- /dev/null
+%# ssq=ssquare2(beta,f,g) returns sin(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 ssq=ssquare1(betah,f,varargin)
+
+f_isall = 0;
+switch nargin
+case 3
+ f_isall = varargin{1};
+end
+
+ssq=sin((pi*betah)*f.^2);
+
+if (f_isall~=1)
+ [n m]=size(f);
+ n=n-1;
+ m=m-1;
+ ssq=cat(2,cat(1,ssq,flipdim(ssq(2:n,:,:),1)),cat(1,flipdim(ssq(:,2:m,:),2),flipdim(flipdim(ssq(2:n,2:m,:),1),2)));
+end
+
+end