```%BIPOLAR returns an M-by-3 matrix containing a blue-red colormap, in
%	in which red stands for positive, blue stands for negative,
%	and white stands for 0.
%
%  Usage: cmap = bipolar(M, lo, hi, contrast);  or  cmap = bipolar;
%
%  cmap:  output M-by-3 matrix for BIPOLAR colormap.
%  M:	  number of shades in the colormap. By default, it is the
%	  same length as the current colormap.
%  lo:	  the lowest value to represent.
%  hi:	  the highest value to represent.
%
%  Inspired from the LORETA PASCAL program:
%	http://www.unizh.ch/keyinst/NewLORETA
%
%  jimmy@rotman-baycrest.on.ca
%
%----------------------------------------------------------------
function cmap = bipolar(M, lo, hi, contrast)

if ~exist('contrast','var')
contrast = 128;
end

if ~exist('lo','var')
lo = -1;
end

if ~exist('hi','var')
hi = 1;
end

if ~exist('M','var')
cmap = colormap;
M = size(cmap,1);
end

steepness = 10 ^ (1 - (contrast-1)/127);
pos_infs = 1e-99;
neg_infs = -1e-99;

doubleredc = [];
doublebluec = [];

if lo >= 0		% all positive

if lo == 0
lo = pos_infs;
end

for i=linspace(hi/M, hi, M)
t = exp(log(i/hi)*steepness);
doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
end

cmap = doubleredc;

elseif hi <= 0	% all negative

if hi == 0
hi = neg_infs;
end

for i=linspace(abs(lo)/M, abs(lo), M)
t = exp(log(i/abs(lo))*steepness);
doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
end

cmap = flipud(doublebluec);

else

if hi > abs(lo)
maxc = hi;
else
maxc = abs(lo);
end

for i=linspace(maxc/M, hi, round(M*hi/(hi-lo)))
t = exp(log(i/maxc)*steepness);
doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
end

for i=linspace(maxc/M, abs(lo), round(M*abs(lo)/(hi-lo)))
t = exp(log(i/maxc)*steepness);
doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
end

cmap = [flipud(doublebluec); doubleredc];

end

return;					% bipolar

```