function [handle,Zi,grid,Xi,Yi] = cm_eeg_topoplot_20140226(Values,loc_file,MAPLIMITS)

%%  functions needed:
%   - THG_eeg_readlocs
%   - THG_eeg_finputcheck

%%%%%%%%%%%%%%%%%%%%%%%% Set defaults %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% whitebk = 'on';  % by default, make gridplot background color = EEGLAB screen background color

rmax = 0.5;                 % actual head radius - Don't change this!

GRID_SCALE = 57;            % plot map on a n x n grid
CIRCGRID   = 201;           % number of angles to use in drawing circles
AXHEADFAC = 1.3;            % head to axes scaling factor
CONTOURNUM = 6;             % number of contour levels to plot
STYLE = 'both';             % default 'style': both,straight,fill,contour,blank

HEADCOLOR = [0 0 0];        % default head color (black = [0 0 0])
BACKCOLOR = [1 1 1];

ELECTRODES = 'on';          % default 'electrodes': on|off|label - set below

EMARKER = '.';              % mark electrode locations with small disks
EMARKERSIZE = 6;

ECOLOR = [0 0 0];           % default electrode color = black

EMARKERLINEWIDTH = 1;       % default edge linewidth for emarkers

HLINEWIDTH = 1;             % default linewidth for head, nose, ears
BLANKINGRINGWIDTH = .05;    % width of the blanking ring 
HEADRINGWIDTH    = .01;     % width of the cartoon head ring
% SHADING = 'flat';         % default 'shading': flat|interp




cmap = colormap;
cmaplen = size(cmap,1);




%% Read the channel location information
[tmpeloc labels Th Rd indices] = cm_eeg_readlocs_20140226( loc_file );
Th = pi/180*Th;                              
allchansind = 1:length(Th);

%%  Channels to plot    
    plotchans = indices;




[x,y]       = pol2cart(Th,Rd);  % transform electrode locations from polar to cartesian coordinates
plotchans   = abs(plotchans);   % reverse indicated channel polarities
allchansind = allchansind(plotchans);
Th          = Th(plotchans);
Rd          = Rd(plotchans);
x           = x(plotchans);
y           = y(plotchans);
labels      = labels(plotchans); % remove labels for electrodes without locations
labels      = strvcat(labels); % make a label string matrix
Values      = Values(plotchans);

%%  Read plotting radius from chanlocs  
    plotrad = min(1.0,max(Rd)*1.02);            % default: just outside the outermost electrode location
    plotrad = max(plotrad,0.5);                 % default: plot out to the 0.5 head boundary
    default_intrad = 1;                         % indicator for (no) specified intrad
    intrad = min(1.0,max(Rd)*1.02);             % default: just outside the outermost electrode location

%%  Set radius of head cartoon
    headrad = rmax;  % (anatomically correct)

%%  Find plotting channels 
    pltchans = find(Rd <= plotrad);             % plot channels inside plotting circle
    intchans = find(x <= intrad & y <= intrad); % interpolate and plot channels inside interpolation square

%%  Eliminate channels not plotted 
    allx      = x;
    ally      = y;
    intchans; % interpolate using only the 'intchans' channels
    pltchans; % plot using only indicated 'plotchans' channels

    intValues = Values(intchans);
    Values = Values(pltchans);
    % now channel parameters and values all refer to plotting channels only

    allchansind = allchansind(pltchans);
    intTh = Th(intchans);           % eliminate channels outside the interpolation area
    intRd = Rd(intchans);
    intx  = x(intchans);
    inty  = y(intchans);
    Th    = Th(pltchans);           % eliminate channels outside the plotting area
    Rd    = Rd(pltchans);
    x     = x(pltchans);
    y     = y(pltchans);

    labels = labels(pltchans,:);

%%  Squeeze channel locations to <= rmax
    squeezefac  = rmax/plotrad;
    intRd       = intRd*squeezefac; % squeeze electrode arc_lengths towards the vertex
    Rd          = Rd*squeezefac;       % squeeze electrode arc_lengths towards the vertex
                              % to plot all inside the head cartoon
    intx        = intx*squeezefac;   
    inty        = inty*squeezefac;  
    x           = x*squeezefac;    
    y           = y*squeezefac;   
    allx        = allx*squeezefac;    
    ally        = ally*squeezefac;   
    % Note: Now outermost channel will be plotted just inside rmax

    % rotate channels based on chaninfo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     rotate = 0;

%%  Make the plot 
    xmin = min(-rmax,min(intx)); xmax = max(rmax,max(intx));
    ymin = min(-rmax,min(inty)); ymax = max(rmax,max(inty));

    % Interpolate scalp map data
    xi         = linspace(xmin,xmax,GRID_SCALE);   % x-axis description (row vector)
    yi         = linspace(ymin,ymax,GRID_SCALE);   % y-axis description (row vector)
    [Xi,Yi,Zi] = griddata(inty,intx,intValues,yi',xi,'cubic'); % interpolate data ('invdist' replaced by)
    
    % Mask out data outside the head
    mask      = (sqrt(Xi.^2 + Yi.^2) <= rmax); % mask outside the plotting circle
    ii        = find(mask == 0);
    Zi(ii)    = NaN;                         % mask non-plotting voxels with NaNs
    grid      = plotrad;                     % unless 'noplot', then 3rd output arg is plotrad

    % Colormap limits
    amin    = MAPLIMITS(1);
    amax    = MAPLIMITS(2);
    delta   = xi(2)-xi(1); % length of grid entry

    % Scale the axes
    cla         % clear current axis
    hold on
    h = gca;    % uses current axes


  set(gca,'Xlim',[-rmax rmax]*AXHEADFAC,'Ylim',[-rmax rmax]*AXHEADFAC,'color','w','xcolor','w','ycolor','w','zcolor','w','xtick',[],'ytick',[]);

%%  plot
    tmph = surface(Xi-delta/2,Yi-delta/2,zeros(size(Zi)),Zi,'EdgeColor','none','FaceColor','flat'); % 'FaceColor','interp' or 'flat'
    colormap jet
    
%%  Set color axis
    caxis([amin amax]) % set coloraxis



    handle = gca;


%%  Plot filled ring to mask jagged grid boundary
    hwidth  = HEADRINGWIDTH;                        % width of head ring 
    hin     = squeezefac*headrad*(1- hwidth/2);     % inner head ring radius

    rwidth  = BLANKINGRINGWIDTH;                    % width of blanking outer ring
    rin     =  rmax*(1-rwidth/2);                   % inner ring radius
    if hin>rin
        rin = hin;                                  % don't blank inside the head ring
    end



%%  mask the jagged border around rmax
    circ    = linspace(0,2*pi,CIRCGRID);
    rx      = sin(circ); 
    ry      = cos(circ); 
    ringx   = [[rx(:)' rx(1) ]*(rin+rwidth)  [rx(:)' rx(1)]*rin];
    ringy   = [[ry(:)' ry(1) ]*(rin+rwidth)  [ry(:)' ry(1)]*rin];

    ringh   = patch(ringx,ringy,0.01*ones(size(ringx)),BACKCOLOR,'edgecolor','none'); hold on 


%%  Plot head outline
    headx = [[rx(:)' rx(1) ]*(hin+hwidth)  [rx(:)' rx(1)]*hin];
    heady = [[ry(:)' ry(1) ]*(hin+hwidth)  [ry(:)' ry(1)]*hin];
    ringh  = patch(headx,heady,ones(size(headx)),'k','edgecolor','none'); hold on

%%  Plot ears and nose
    base  = rmax-.0046;
    basex = 0.18*rmax;          % nose width
    tip   = 1.15*rmax; 
    tiphw = .04*rmax;           % nose tip half width
    tipr  = .01*rmax;           % nose tip rounding
    q     = .04;                % ear lengthening
    EarX  = [.497-.005  .510  .518  .5299 .5419  .54    .547   .532   .510   .489-.005]; % rmax = 0.5
    EarY  = [q+.0555 q+.0775 q+.0783 q+.0746 q+.0555 -.0055 -.0932 -.1313 -.1384 -.1199];
    sf    = headrad/plotrad;                                          % squeeze the model ears and nose 
                                                                      % by this factor
    plot3([basex;tiphw;0;-tiphw;-basex]*sf,[base;tip-tipr;tip;tip-tipr;base]*sf,...
         2*ones(size([basex;tiphw;0;-tiphw;-basex])),...
         'Color',HEADCOLOR,'LineWidth',HLINEWIDTH);                                       % plot nose
    plot3(EarX*sf,EarY*sf,2*ones(size(EarX)),'color',HEADCOLOR,'LineWidth',HLINEWIDTH)    % plot left ear
    plot3(-EarX*sf,EarY*sf,2*ones(size(EarY)),'color',HEADCOLOR,'LineWidth',HLINEWIDTH)   % plot right ear


%%  Show electrode information
     plotax = gca;
     axis square                                      
     axis off

     pos = get(gca,'position');
     xlm = get(gca,'xlim');
     ylm = get(gca,'ylim');
     axis square % make textax square

     pos = get(gca,'position');
     set(plotax,'position',pos);

     xlm = get(gca,'xlim');
     set(plotax,'xlim',xlm);

     ylm = get(gca,'ylim');
     set(plotax,'ylim',ylm);                               % copy position and axis limits again


%%  Mark electrode locations only
    ELECTRODE_HEIGHT = 2.1;  % z value for plotting electrode information (above the surf)
    hp2 = plot3(y*.95,x*.95,ones(size(x))*ELECTRODE_HEIGHT,...
            EMARKER,'Color',ECOLOR,'markersize',EMARKERSIZE,'linewidth',EMARKERLINEWIDTH);


%%%%%%%%%%%%% Set EEGLAB background color to match head border %%%%%%%%%%%%%%%%%%%%%%%%

hold off
axis off
return