function drawFoveatedRFs(RFscalingRate, maxEccentricity, ...
                         minEccentricity, gap, col)

% drawFoveatedRFs(RFscalingRate, maxEccentricity, minEccentricity, gap, col)
%
% In the current figure, draw a cartoon diagram of circular Receptive
% Fields that grow linearly eccentricity.
%
% - RFscalingRate: expresses the ratio of RF radius to  eccentricity
% - maxEccentricity: maximum distance from origin (units are arbitrary)
% - minEccentricity: minimum distance at which RFs are drawn.
%      Optional, defaulting to maEccentricty/20 
% - gap: space between RFs, as a fraction of RF radius.  Must be in
%        the interval [0, 1].  Optional, defaulting to 0.15
% - col: color specification (single-character string, or rgb triple)
%
% Examples:
%   figure, drawFoveatedRFs(.5, 15, .2, .2, [1 0 0])
%   figure, drawFoveatedRFs(.3, 20, .2, .2, [0 1 0])
%
%  human V1 based on Kay et al 2013
%   figure(1), set(gcf, 'Position', get(0, 'ScreenSize')); drawFoveatedRFs(2/20, 40, 1, 0.15, [1 0 0]); % V1
%   figure(2), drawFoveatedRFs(5/20, 40, 4, 0.15, [0.8 0.8 0]);  % V3
%   figure(3), drawFoveatedRFs(10/20, 40, 1, 0.15, [0 0.8 0.8]); % hV4
%   hgexport(1, '~/Desktop/prf_V1.eps')
%   hgexport(2, '~/Desktop/prf_V3.eps')
%   hgexport(3, '~/Desktop/prf_hV4.eps')

% EPS, 8/2009.  Revised 9/2014 by HH and JW.

if (nargin < 4.5)
    col = [0 0 0];
end

if (nargin < 3.5)
    gap = 0.15;
elseif ((gap >= 1) | (gap <= 0))
    warning('Gap must be between 0 and 1.  Setting to 0.15');
    gap = 0.15;
end

if (nargin < 2.5)
    minEccentricity = maxEccentricity/20;  %*** Magic number
end

clf

numPerRing = floor(pi/RFscalingRate); 
angles = 2*pi*([1:numPerRing]-1)/numPerRing;
Xctrs = cos(angles);  Yctrs = sin(angles);

% geometric calculation of spacing between rings
cf = cos((angles(2)-angles(1))/2);
cfPlusR2 = cf + RFscalingRate^2;
eFactor = (cfPlusR2 - sqrt( cfPlusR2^2 + cf^2 * (RFscalingRate^2 - 1))) / cf^2;

cNumPoints = 60; %** Magic number
cAngles = 2*pi*([1:cNumPoints]'-1)/(cNumPoints-1);
% coordinates for one circle, horizontally displayed, at eccentricity=1
X1 = (1-gap)*RFscalingRate*cos(cAngles);  Y1 = (1-gap)*RFscalingRate*sin(cAngles);

Xall = repmat(Xctrs, cNumPoints, 1) + repmat(X1,1,numPerRing);
Yall = repmat(Yctrs, cNumPoints, 1) + repmat(Y1,1,numPerRing);

rotAngle = pi/numPerRing;
XallR = cos(rotAngle)*Xall + sin(rotAngle)*Yall;
YallR = -sin(rotAngle)*Xall + cos(rotAngle)*Yall;

ecc  = (maxEccentricity/(1+RFscalingRate));
hold on
while (ecc > minEccentricity)
    h = fill(ecc*Xall, ecc*Yall, col);
    for n=1:length(h)
        set(h(n), 'LineStyle', 'none');
    end
    ecc = eFactor * ecc;
    tmp = XallR;  XallR = Xall; Xall = tmp;  % toggle to rotated
    tmp = YallR;  YallR = Yall; Yall = tmp;  % circles
end
hold off
axis equal
axis(maxEccentricity*[-1 1 -1 1]);
%axis off

return

%% ----------------
% Examples:
drawFoveatedRFs(0.1, 80)
drawFoveatedRFs(0.22, 80)
figure(99); clf, drawFoveatedRFs(23/32/2, 40, 4, 0.1, [0.9 0.7 0.5])
axis off