% Returns colormaps needed for various moment plots
function [x] = cmaps(mapname, arg)
    bgcolor = [0.89 0.87 0.83];
    switch mapname
        case 'zmap'
            % Standard reflectivity colormap
            if nargin<2
                zlim = [0 75];
            else
                zlim = arg;
            end
            if zlim(1)<0,zlim(1) = 0; end
            if zlim(2)>75,zlim(2) = 75; end
            tmp = zlim(1):5:zlim(2);
            tmp = floor(tmp/5)+1;
            std_map = [bgcolor; ...
                       0.20 1.00 1.00; ...
                       0.20 0.60 1.00; ...
                       0.00 0.00 1.00; ...
                       0.30 1.00 0.00; ...
                       0.10 0.80 0.00; ...
                       0.00 0.60 0.00; ...
                       1.00 1.00 0.00; ...
                       1.00 0.75 0.00; ...
                       1.00 0.50 0.00; ...
            %           1.00 0.30 0.30; ...
                       1.00 0.00 0.00; ...
                       0.75 0.00 0.00; ...
                       0.50 0.00 0.00; ...
                       1.00 0.00 0.80; ...
                       0.60 0.30 1.00; ...
                       1.00 1.00 1.00];
            x = std_map(tmp,:);
            
        case 'rgmap'
            % Red-green velocity colormap
            if nargin<2
                num = size(colormap,1);
            else
                num = arg;
            end
            % Fixed color points
            pts = [0
                0.30
                0.50
                0.70
                1.00];
            % Fixed colors at each point
            rgb = [0.00 0.20 0.00;...
                0.00 0.80 0.00;...
                0.85 0.85 0.85;...
                0.80 0.00 0.00;...
                0.20 0.00 0.00];
            if any(rgb>1),fprintf('PT Matrix cannot have value > 1.0.\n');
                return;
            end
            % Interpolation points
            base = linspace(0,1,num);
            base = base(:);
            x = interp1(pts, rgb, base, 'linear');
            x(x<0) = 0;
            x(x>1) = 1;
            
        case 'rapmap'
            % Intensity-scaling colormap; often used for spectrum width
            if nargin < 2
                npts = size(colormap,1);
            else
                npts = arg;
            end
            cromIQ = [0.00, 0.00; ...
                     -0.07, 0.13; ...
                     -0.11, 0.20; ...
                     -0.03, 0.21; ...
                      0.23, 0.19; ...
                      0.39, 0.11; ...
                      0.41,-0.01; ...
                      0.39,-0.12; ...
                      0.32,-0.19; ...
                      0.20,-0.15; ...
                      0.00, 0.00];
            oldx = 1:11;
            newx = linspace(1,11,npts).';
            yiqmap(1:npts,1) = linspace(0,1,npts).';
            yiqmap(1:npts,2) = interp1(oldx,cromIQ(:,1),newx,'pchip');
            yiqmap(1:npts,3) = interp1(oldx,cromIQ(:,2),newx,'pchip');
            x=yiqmap*[1.0000, 1.0000, 1.0000; ...
                      0.9562,-0.2727,-1.1037; ...
                      0.6214,-0.6468, 1.7006];
            x(x<0)=0; x(x>1)=1;
            
        case 'dmap'
            % Differential reflectivity colormap
            % Similar to zmap, but with rgmap-style interpolation
            % Fixed color points (forming color regions)
            pts = [0, 9, ...
                10, 39, ... 
                40, 69, ...
                70, 99, ...
                100, 129, ...
                130, 159, ...
                160, 189, ...
                190, 219, ...
                220, 249, ...
                250, 254] ./ 254.0;
            % Fixed colors at each point
            rgb = [0.30, 0.45, 0.50;0.60, 0.90, 1.00    
                0.45, 0.20, 0.80;   0.70, 0.40, 1.00    
                0.50, 0.20, 0.35;   1.00, 0.50, 0.85    
                0.70, 0.50, 0.15;   1.00, 1.00, 0.85    
                1.00, 1.00, 1.00;   0.00, 0.35, 1.00    
                0.10, 1.00, 0.50;   0.00, 0.50, 0.00    
                1.00, 1.00, 0.00;   1.00, 0.50, 0.00    
                1.00, 0.00, 0.00;   0.50, 0.00, 0.00    
                1.00, 0.00, 1.00;   0.50, 0.00, 0.50    
                1.00, 1.00, 1.00;   0.60, 1.00, 1.00];
            if any(rgb>1),fprintf('PT Matrix cannot have value > 1.0.\n');
                return;
            end
            % Interpolation points (51 values)
            base = linspace(0,1,51);
            base = base(:);
            x = interp1(pts, rgb, base, 'linear');
            x(x<0) = 0;
            x(x>1) = 1;
            
        case 'rmap'
            % Correlation coefficient colormap
            % Fixed color points (forming color regions)
            pts = [0, 6, ...
                7,  11, ...
                12, 16, ... 
                17, 21, ...
                22, 26, ...
                27, 31, ...
                32, 36, ...
                37, 41] ./ 41.0;
            % Fixed colors at each point
            rgb = [0.00, 0.00, 0.00; 0.50, 0.60, 0.70;
                0.00, 1.00, 1.00; 0.00, 0.00, 0.85;
                0.00, 1.00, 0.00; 0.00, 0.50, 0.00;
                1.00, 1.00, 0.00; 1.00, 0.50, 0.00;
                1.00, 0.00, 0.00; 0.50, 0.00, 0.00;
                1.00, 0.00, 1.00; 0.50, 0.00, 0.50;
                0.60, 0.22, 1.00; 0.35, 0.11, 0.55;
                0.40, 0.45, 1.00; 0.20, 0.22, 0.60];
            if any(rgb>1),fprintf('PT Matrix cannot have value > 1.0.\n');
                return;
            end
            % Interpolation points (42 values)
            base = linspace(0,1,42);
            base = base(:);
            x = interp1(pts, rgb, base, 'linear');
            x(x<0) = 0;
            x(x>1) = 1;
    end
end