function flightrack(varargin)
%
% function flighttrack(args)
%   Allows you to colour any variable from the FAAM  core netcdf data file
%   along the flighttrack, and/or vector any two variables along the flighttrack. 
%   In addition AVHRR data can be coloured if available.
%
% =======================================================
% usage:
%      flighttrack(property1,value1,property2,value2,...)
% =======================================================
%
% Valid properties:
%
%
%   'level'     -   Can take values 'lower' or 'upper'. Lower will plot
%                   lower altitudes above higher ones, and vice-versa.
%                   If this option is omitted, the default will plot newer
%                   data above older data.
%
%   'corefile' -   String pointing to the netcdf file containing core data
%                   for the flight.
%
%   'var'       -   Variable to be coloured along the flighttrack. This can
%                   be a string naming the variable in the netcdf file, or
%                   a variable in the current matlab workspace. If 0 is
%                   passed as the value to this property, no variable will
%                   be coloured along the flighttrack.
%
%   'domain'    -   Define the domain in the form
%                   ...'domain',[lat_min lat_max lon_min lon_max],...
%                   If no domain is specified the domain will be based on
%                   1) The domain of AVHRR data if 'satfile' is set, or
%                   2) The domain of the flighttrack, with some padding.
%
%   'vector'    -   Plot vectors along the flighttrack. The value should be
%                   a two element column vector containing strings matching
%                   the names of variables in the netcdf datafile.
%                   E.g. ...,'vector',['Uvar';'Vvar'],...
%
%   'vec_zrange'-   Two element row vector defining minimum and maximum
%                   altitudes to plot vectors over (metres)
%
%   'range'     -   A row vector containing the min and max values of var
%                   to be plotted.
%
%   'satfile'   -   String pointing to a .mat file containing AVHRR data.
%
%   'skip'      -   A value to reduce the resolution of the AVHRR data. For
%                   example, a value of 5 will only use every 5th data
%                   point.
%
%   'dotsize'   -   defines the size of dots to use in generating a 
%                   coloured flight path. This generally wants to be larger
%                   the smaller the domain.
%
%   'coast'     -   Can be set to 'low' to plot with m_coast, otherwise
%                   m_gshhs_i is used.
%
%   'vecscale'  -   Allows you to vary the scale of vectors when producing
%                   a vector plot. Default is 0.001.
%
%   'vec_skip'  -   Plot vector data desampled by this factor (below 0.1
%                   Hz)
%
%   'cloud_map' -   Colormap for colouring the satellite data.
%
%   'data_map'  -   Colormap for data along the flight track
%
%   Note: this function requires:   arrow.m
%                                    > http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=278&objectType=file
%
%                                   m_map toolbox (optional high-res
%                                   coastlines recommended)
%
%                                   Flight coredata netCDF file
%                                   AVHRR workspace (Optional: for sat imagery)


if(length(varargin) == 0)
    help flighttrack;
    return;
end


%===========================
% DEFINE PARAMETER DEFAULTS
%===========================
flight_level = 0;           % Plot data in flight-time order
core_file = 'data.nc';      % If no core file is given, assume 'data.nc'
range = [0 0];              % data range
nc = 1;                     % Reading from netcdf(1) or workspace(0)
sat = 0;                    % Colour AVHRR imagery
skip = 1;                   % Sampling interval (km) for AVHRR
track = 0;                  % Plot a flight track?
vecs = 0;                   % plot vectors
var_set = 0;                % Has a variable to plot been defined?
coast = 1;                  % High res(1) or low res(0) coastlines
domain = [0 0 0 0];         % domain to plot
scatter_size = 4;           % Size of dots for the scatter plot
vecscale = 0.001;           % Scaling parameter for vectors
vec_height_range = [0 0];   % Sub range of altitudes to plot vectors over
vec_skip = 10;              % Vector desampling factor
cbar = 0;                   % Add colorbars? NOT CURRENTLY WORKING
data_map = jet;
cloud_map = 0;

% Deal with user input
k=1;
while(k<=length(varargin)),
    switch(lower(varargin{k})),
        case('level'),
            switch(lower(varargin{k+1})),
                case('lower'),
                    flight_level = 2;
                case('upper'),
                    flight_level = 1;
            end
        case('corefile'),
            track = 1;
            core_file = varargin{k+1};
        case('range'),
            range = varargin{k+1};
            if(size(range) ~= size([0 0]))
                error('range must be a two element vector');
            end
        case('coast'),
            temp = varargin{k+1};
            switch(temp),
                case('low'),
                    coast = 0;
            end
        case('var'),
            x = varargin{k+1};
            var_set = 1;
            if(length(x) > 4)
                if(x(1:4) ~= 'PARA')
                    nc = 0;
                end
            else
                nc = 1;
            end
            if(x == 0)
                var_set = 0;
            end
        case('domain'),
            if(size(varargin{k+1}) ~= size(domain))
                warning('Domain incorrectly specified, using a default...');
            else
                domain = varargin{k+1};
            end
        case('vector'),
            vec_vars = varargin{k+1};
            vecs = 1;
        case('vecscale'),
            vecscale = varargin{k+1};
            if(vecs == 0)
                warning('''vecscale'' given but vector plotting is off');
            end
        case('vec_zrange'),
            vec_height_range = varargin{k+1};
            if(size(vec_height_range) ~= size([0 0]))
                error('''vec_zrange'' is invalid');
            end
        case('satfile'),
            sat = 1;
            sat_file = varargin{k+1};
        case('skip'),
            skip = varargin{k+1};
            skip = floor(skip);
            if(skip < 1)
                error('skip must be an integer >= 1');
            end
        case('dotsize'),
            scatter_size = varargin{k+1};
            if(scatter_size < 0)
                error(['''dotsize'' is invalid']);
            end
            if(scatter_size == 0)
                warning('''dotsize'' is set to zero');
            end
        case('vec_skip'),
            vec_skip = varargin{k+1};
            if(vec_skip ~= floor(vec_skip))
                error('''vec_skip'' must be an integer');
            end
            if(vec_skip <= 0)
                error('''vec_skip'' must be >= 1');
            end
        case('cloud_map'),
            cloud_map = varargin{k+1};
        case('data_map'),
            data_map = varargin{k+1};
        case('colorbar'),
            cbar = varargin{k+1};
    end
    k=k+2;
    if(mod(length(varargin),2)~=0)
        error('mod(nargin, 2) ~= 0. Check your inputs!');
    end
end

if(sat ~= 0)
    if(isempty(sat_file))
        error('Satellite imagery requested, but no file given');
    end

    % Load AVHRR workspace
    eval(['load ' sat_file])

    % Reduce Imagery resolution...
    if(domain ~= [0 0 0 0])
        Lat = Lat(:,1);
        Lon = Lon(1,:);

        a = find(Lat <= domain(2));
        b = find(Lat >= domain(1));
        c = find(Lon <= domain(4));
        d = find(Lon >= domain(3));
        Lata = Lat(a);
        Latb = Lat(b);
        Lat = intersect(Lata,Latb);
        Lona = Lon(c);
        Lonb = Lon(d);
        Lon = intersect(Lona, Lonb);

        [Lon Lat] = meshgrid(Lon, Lat);
        Image = Image(intersect(a,b),intersect(c,d));

        clear a b c d Lata Latb Lona Lonb
    end

    Lat = Lat(1:skip:end,1:skip:end);
    Lon = Lon(1:skip:end,1:skip:end);
    Image = Image(1:skip:end,1:skip:end);
   
    % Load colourmap for satellite data
    if(cloud_map == 0)
        map = colmap;
        map = 1.-map;   % make clouds light rather than dark
    else
        map = cloud_map;
    end
    
    % Nomalise image
    Image = Image./max(max(Image)); 
end

if(vecs == 1)
    % We're using vectors - read in the data and the quality flags
    vec_var1 = readnc(core_file, vec_vars(1,:));
    vec_var2 = readnc(core_file, vec_vars(2,:));
    vec_var1_flag = readnc(core_file, strcat(vec_vars(1,:),'FLAG'));
    vec_var2_flag = readnc(core_file, strcat(vec_vars(2,:),'FLAG'));
end

if(track == 1)
    % Read lat/lon coordinates from .nc file
    north = readnc(core_file, 'PARA0580');
    east = readnc(core_file, 'PARA0581');
    north_flag = readnc(core_file, 'PARA0580FLAG');
    east_flag = readnc(core_file, 'PARA0581FLAG');

    % Read flight profile from .nc file
    z = readnc(core_file, 'PARA0582');
    z_flag = readnc(core_file, 'PARA0582FLAG');
end

% Are we reading a variable from .nc file for variable as an arguement?
if(track == 1 && var_set ~= 0)
    if(nc)
        ft_var = readnc(core_file, x);
        var_flag = readnc(core_file, strcat(x, 'FLAG'));
    else
        ft_var = x;
    end
end

% Resample data higher than 1 Hz down to 1 Hz
if(track == 1 && var_set ~= 0)
    if(min(size(ft_var)) > 1)
        ft_var = ft_var(1,:);
    end
end
if(vecs == 1)
    if(min(size(vec_var1) > 1))
        vec_var1 = vec_var1(1,:);
        vec_var1_flag = vec_var1_flag(1,:);
    end
    if(min(size(vec_var2) > 1))
        vec_var2 = vec_var2(1,:);
        vec_var2_flag = vec_var2_flag(1,:);
    end
end
%% N.B. If data is lower than 1 Hz will need to resample lat,lon,z before
%% this. Also var no doubt. Prolly vecs as well. Dammit!

if(track == 1)
    % NaN out dodgy data
    z(find(z_flag > 0)) = NaN;
    east(find(east_flag > 0)) = NaN;
    north(find(north_flag > 0)) = NaN;
    
    if(var_set ~= 0)
        ft_var(find(var_flag > 0)) = NaN;
    end
    
    %clear *flag
end

if(track == 1)
    % Sort arrays to overlay/underlay flight levels
    if(flight_level == 1)
        [z i] = sort(z, 'descend');
        ft_var = ft_var(i);
        north = north(i);
        east = east(i);
    elseif(flight_level == 2)
        [z i] = sort(z, 'ascend');
        ft_var = ft_var(i);
        north = north(i);
        east = east(i);
    end
end

%% Create Image bounds...
if(domain ~= [0 0 0 0])         % Domain has been specified
    lat_range = [domain(1) domain(2)];
    lon_range = [domain(3) domain(4)];
elseif(sat ~= 0)                % No domain specified, but have sat imagery
    lat_range = [min(min(Lat)) max(max(Lat))];
    lon_range = [min(min(Lon)) max(max(Lon))];
elseif(track == 1)              % No domain or sat data specified, base
    lat_range = [min(north)-3 max(north)+3];        % domain on flighttrack
    lon_range = [min(east)-3 max(east)+3];
end

% Create a projection
m_proj('lambert', 'lats', lat_range, 'lons', lon_range);

%% SATELLITE IMAGERY STUFF %%
if(sat == 1)
    [x y] = m_ll2xy(Lon, Lat);
    clear Lon Lat
    colormap(map);

    pcolor(x,y,Image),shading interp;
end
%% DONE WITH IMAGERY %%

hold on          %err... turn on hold ;)

% Plot coastlines 
if(coast == 0)
    if(sat ~= 0)
        m_coast('color','y');
    else
        m_coast('patch',[.5 .5 .5]);
    end 
else
    if(sat ~= 0)
        m_gshhs_i('color', 'y');
    else
        m_gshhs_i('patch', [.5 .5 .5]);
    end
end

%Add a colorbar if required
if(cbar ~= 0)
    colorbar;
end

% name the current axis
ax1 = gca;

% Make a funky grid
m_grid('linestyle', ':', 'box', 'fancy');

if(track == 1)

    % Create [x y] variables, so don't have to use m_[plot function]
    [lon lat] = m_ll2xy(east, north);

    % desample geographic data to 0.1 Hz
    lon = lon(1:10:end);
    lat = lat(1:10:end);
    if(var_set ~= 0)
        ft_var = ft_var(1:10:end);
    end
    z = z(1:10:end);
    clear east north        % tidy up

    % Define a second layer axis.
    ax2 = axes('position', get(ax1, 'position'));

    
    % If no data range has been specified, use full range of the data
    if(range == [0 0])
        if(var_set ~= 0)
            range = [min(ft_var) max(ft_var)];
        end
    end
    % NaN out data outside range.
    if(var_set ~= 0)
        a = find(ft_var < range(1));
        b = find(ft_var > range(2));
        ft_var(a) = nan;
        ft_var(b) = nan;
    end
    
    % If a variable has been specified, colour it along the flighttrack,
    % otherwise just plot a red line
    if(var_set ~= 0)
        scatter(lon,lat,scatter_size,ft_var,'filled');
    else
        plot(lon,lat,'-r');
    end

    % Draw a standard-ish grid
    m_grid('linestyle', ':', 'box', 'fancy');

    % Double check axes are aligned
    set(ax1,'position',get(ax1,'position'));        
    
    if(sat == 1)
        %% Do con. colourmap stuff - i.e. paste two colourmaps together so
        %% that we can plot both clouds and colour along the flight track
        colormap([map; data_map]);
        CmLength   = length(colormap);   % Colormap length
        BeginSlot1 = 1;                  % Beginning slot
        EndSlot1   = length(map);        % Ending slot
        BeginSlot2 = EndSlot1 + 1;
        EndSlot2   = CmLength;
        CLim1      = get(ax1,'CLim');  % CLim values for each axis
        CLim2      = get(ax2,'CLim');
        set(ax1,'CLim',newclim(BeginSlot1,EndSlot1,CLim1(1),CLim1(2),CmLength))
        set(ax2,'CLim',newclim(BeginSlot2,EndSlot2,CLim2(1),CLim2(2),CmLength))
        set(ax2,'visible','off');
    else
        colormap(jet);
    end
end

if(vecs == 1) 
    % Quality control the vector data
    vec_var1(find(vec_var1_flag>0)) = nan;
    vec_var2(find(vec_var2_flag>0)) = nan;
    
    % Desample the vector data to 0.1 Hz
    vec_var1 = vec_var1(1:10:end);
    vec_var2 = vec_var2(1:10:end);
    
    % If a height range has been given, NaN out any vectors that are not in
    % this range.
    if(vec_height_range(2) > 0)
        a = find(z <= vec_height_range(1));
        b = find(z >= vec_height_range(2));
        vec_var1(a) = nan;
        vec_var1(b) = nan;
        vec_var2(a) = nan;
        vec_var2(b) = nan;
    end

    for(i = 1:vec_skip:length(lon))
        do_plot = 1;
        if(isnan(lon(i)) || isnan(lat(i)) || isnan(vec_var1(i)) || isnan(vec_var2(i)))
            do_plot = 0;
        end

        if(do_plot == 1)
           warning off
           arrow([lon(i),lat(i)], [lon(i)+vecscale*vec_var1(i), lat(i)+vecscale*vec_var2(i)],'tipangle',5,'baseangle',12,'length',10); 
        end
    end
end

% Make the top plot transparent
mgc = findobj('tag','m_grid_color');
set(mgc(1),'visible','off')

%Add second colorbar
% if(cbar ~= 0)
%     colorbar;
% end

end



function out = colmap

% NOAA colourmap for AVHRR data. Note the inverse of this is actually
% coloured, to use this exact colourmap, you will have to pass it as an
% argument to flighttrack

out = [
    1.0000    1.0000    1.0000;
    0.9961    0.9961    0.9961;
    0.9922    0.9922    0.9922;
    0.9882    0.9882    0.9882;
    0.9843    0.9843    0.9843;
    0.9804    0.9804    0.9804;
    0.9765    0.9765    0.9765;
    0.9725    0.9725    0.9725;
    0.9686    0.9686    0.9686;
    0.9647    0.9647    0.9647;
    0.9608    0.9608    0.9608;
    0.9569    0.9569    0.9569;
    0.9529    0.9529    0.9529;
    0.9490    0.9490    0.9490;
    0.9451    0.9451    0.9451;
    0.9412    0.9412    0.9412;
    0.9373    0.9373    0.9373;
    0.9333    0.9333    0.9333;
    0.9294    0.9294    0.9294;
    0.9255    0.9255    0.9255;
    0.9216    0.9216    0.9216;
    0.9176    0.9176    0.9176;
    0.9137    0.9137    0.9137;
    0.9098    0.9098    0.9098;
    0.9059    0.9059    0.9059;
    0.9020    0.9020    0.9020;
    0.8980    0.8980    0.8980;
    0.8941    0.8941    0.8941;
    0.8902    0.8902    0.8902;
    0.8863    0.8863    0.8863;
    0.8824    0.8824    0.8824;
    0.8784    0.8784    0.8784;
    0.8745    0.8745    0.8745;
    0.8706    0.8706    0.8706;
    0.8667    0.8667    0.8667;
    0.8627    0.8627    0.8627;
    0.8588    0.8588    0.8588;
    0.8549    0.8549    0.8549;
    0.8510    0.8510    0.8510;
    0.8471    0.8471    0.8471;
    0.8431    0.8431    0.8431;
    0.8392    0.8392    0.8392;
    0.8353    0.8353    0.8353;
    0.8314    0.8314    0.8314;
    0.8275    0.8275    0.8275;
    0.8235    0.8235    0.8235;
    0.8196    0.8196    0.8196;
    0.8157    0.8157    0.8157;
    0.8118    0.8118    0.8118;
    0.8078    0.8078    0.8078;
    0.8039    0.8039    0.8039;
    0.8000    0.8000    0.8000;
    0.7961    0.7961    0.7961;
    0.7922    0.7922    0.7922;
    0.7882    0.7882    0.7882;
    0.7843    0.7843    0.7843;
    0.7804    0.7804    0.7804;
    0.7765    0.7765    0.7765;
    0.7725    0.7725    0.7725;
    0.7686    0.7686    0.7686;
    0.7647    0.7647    0.7647;
    0.7608    0.7608    0.7608;
    0.7569    0.7569    0.7569;
    0.7529    0.7529    0.7529;
    0.7490    0.7490    0.7490;
    0.7451    0.7451    0.7451;
    0.7412    0.7412    0.7412;
    0.7373    0.7373    0.7373;
    0.7333    0.7333    0.7333;
    0.7294    0.7294    0.7294;
    0.7255    0.7255    0.7255;
    0.7216    0.7216    0.7216;
    0.7176    0.7176    0.7176;
    0.7137    0.7137    0.7137;
    0.7098    0.7098    0.7098;
    0.7059    0.7059    0.7059;
    0.7020    0.7020    0.7020;
    0.6980    0.6980    0.6980;
    0.6941    0.6941    0.6941;
    0.6902    0.6902    0.6902;
    0.6863    0.6863    0.6863;
    0.6824    0.6824    0.6824;
    0.6784    0.6784    0.6784;
    0.6745    0.6745    0.6745;
    0.6706    0.6706    0.6706;
    0.6667    0.6667    0.6667;
    0.6627    0.6627    0.6627;
    0.6588    0.6588    0.6588;
    0.6549    0.6549    0.6549;
    0.6510    0.6510    0.6510;
    0.6471    0.6471    0.6471;
    0.6431    0.6431    0.6431;
    0.6392    0.6392    0.6392;
    0.6353    0.6353    0.6353;
    0.6314    0.6314    0.6314;
    0.6275    0.6275    0.6275;
    0.6235    0.6235    0.6235;
    0.6196    0.6196    0.6196;
    0.6157    0.6157    0.6157;
    0.6118    0.6118    0.6118;
    0.6078    0.6078    0.6078;
    0.6039    0.6039    0.6039;
    0.6000    0.6000    0.6000;
    0.5961    0.5961    0.5961;
    0.5922    0.5922    0.5922;
    0.5882    0.5882    0.5882;
    0.5843    0.5843    0.5843;
    0.5804    0.5804    0.5804;
    0.5765    0.5765    0.5765;
    0.5725    0.5725    0.5725;
    0.5686    0.5686    0.5686;
    0.5647    0.5647    0.5647;
    0.5608    0.5608    0.5608;
    0.5569    0.5569    0.5569;
    0.5529    0.5529    0.5529;
    0.5490    0.5490    0.5490;
    0.5451    0.5451    0.5451;
    0.5412    0.5412    0.5412;
    0.5373    0.5373    0.5373;
    0.5333    0.5333    0.5333;
    0.5294    0.5294    0.5294;
    0.5255    0.5255    0.5255;
    0.5216    0.5216    0.5216;
    0.5176    0.5176    0.5176;
    0.5137    0.5137    0.5137;
    0.5098    0.5098    0.5098;
    0.5059    0.5059    0.5059;
    0.5020    0.5020    0.5020;
    0.4980    0.4980    0.4980;
    0.4941    0.4941    0.4941;
    0.4902    0.4902    0.4902;
    0.4863    0.4863    0.4863;
    0.4824    0.4824    0.4824;
    0.4784    0.4784    0.4784;
    0.4745    0.4745    0.4745;
    0.4706    0.4706    0.4706;
    0.4667    0.4667    0.4667;
    0.4627    0.4627    0.4627;
    0.4588    0.4588    0.4588;
    0.4549    0.4549    0.4549;
    0.4510    0.4510    0.4510;
    0.4471    0.4471    0.4471;
    0.4431    0.4431    0.4431;
    0.4392    0.4392    0.4392;
    0.4353    0.4353    0.4353;
    0.4314    0.4314    0.4314;
    0.4275    0.4275    0.4275;
    0.4235    0.4235    0.4235;
    0.4196    0.4196    0.4196;
    0.4157    0.4157    0.4157;
    0.4118    0.4118    0.4118;
    0.4078    0.4078    0.4078;
    0.4039    0.4039    0.4039;
    0.4000    0.4000    0.4000;
    0.3961    0.3961    0.3961;
    0.3922    0.3922    0.3922;
    0.3882    0.3882    0.3882;
    0.3843    0.3843    0.3843;
    0.3804    0.3804    0.3804;
    0.3765    0.3765    0.3765;
    0.3725    0.3725    0.3725;
    0.3686    0.3686    0.3686;
    0.3647    0.3647    0.3647;
    0.3608    0.3608    0.3608;
    0.3569    0.3569    0.3569;
    0.3529    0.3529    0.3529;
    0.3490    0.3490    0.3490;
    0.3451    0.3451    0.3451;
    0.3412    0.3412    0.3412;
    0.3373    0.3373    0.3373;
    0.3333    0.3333    0.3333;
    0.3294    0.3294    0.3294;
    0.3255    0.3255    0.3255;
    0.3216    0.3216    0.3216;
    0.3176    0.3176    0.3176;
    0.3137    0.3137    0.3137;
    0.3098    0.3098    0.3098;
    0.3059    0.3059    0.3059;
    0.3020    0.3020    0.3020;
    0.2980    0.2980    0.2980;
    0.2941    0.2941    0.2941;
    0.2902    0.2902    0.2902;
    0.2863    0.2863    0.2863;
    0.2824    0.2824    0.2824;
    0.2784    0.2784    0.2784;
    0.2745    0.2745    0.2745;
    0.2706    0.2706    0.2706;
    0.2667    0.2667    0.2667;
    0.2627    0.2627    0.2627;
    0.2588    0.2588    0.2588;
    0.2549    0.2549    0.2549;
    0.2510    0.2510    0.2510;
    0.2471    0.2471    0.2471;
    0.2431    0.2431    0.2431;
    0.2392    0.2392    0.2392;
    0.2353    0.2353    0.2353;
    0.2314    0.2314    0.2314;
    0.2275    0.2275    0.2275;
    0.2235    0.2235    0.2235;
    0.2196    0.2196    0.2196;
    0.2157    0.2157    0.2157;
    0.2118    0.2118    0.2118;
    0.2078    0.2078    0.2078;
    0.2039    0.2039    0.2039;
    0.2000    0.2000    0.2000;
    0.1961    0.1961    0.1961;
    0.1922    0.1922    0.1922;
    0.1882    0.1882    0.1882;
    0.1843    0.1843    0.1843;
    0.1804    0.1804    0.1804;
    0.1765    0.1765    0.1765;
    0.1725    0.1725    0.1725;
    0.1686    0.1686    0.1686;
    0.1647    0.1647    0.1647;
    0.1608    0.1608    0.1608;
    0.1569    0.1569    0.1569;
    0.1529    0.1529    0.1529;
    0.1490    0.1490    0.1490;
    0.1451    0.1451    0.1451;
    0.1412    0.1412    0.1412;
    0.1373    0.1373    0.1373;
    0.1333    0.1333    0.1333;
    0.1294    0.1294    0.1294;
    0.1255    0.1255    0.1255;
    0.1216    0.1216    0.1216;
    0.1176    0.1176    0.1176;
    0.1137    0.1137    0.1137;
    0.1098    0.1098    0.1098;
    0.1059    0.1059    0.1059;
    0.1020    0.1020    0.1020;
    0.0980    0.0980    0.0980;
    0.0941    0.0941    0.0941;
    0.0902    0.0902    0.0902;
    0.0863    0.0863    0.0863;
    0.0824    0.0824    0.0824;
    0.0784    0.0784    0.0784;
    0.0745    0.0745    0.0745;
    0.0706    0.0706    0.0706;
    0.0667    0.0667    0.0667;
    0.0627    0.0627    0.0627;
    0.0588    0.0588    0.0588;
    0.0549    0.0549    0.0549;
    0.0510    0.0510    0.0510;
    0.0471    0.0471    0.0471;
    0.0431    0.0431    0.0431;
    0.0392    0.0392    0.0392;
    0.0353    0.0353    0.0353;
    0.0314    0.0314    0.0314;
    0.0275    0.0275    0.0275;
    0.0235    0.0235    0.0235;
    0.0196    0.0196    0.0196;
    0.0157    0.0157    0.0157;
    0.0118    0.0118    0.0118;
    0.0078    0.0078    0.0078;
    0.0039    0.0039    0.0039;
    0         0         0];
end



function CLim = newclim(BeginSlot,EndSlot,CDmin,CDmax,CmLength)
% Hashes together a couple of different colourspaces

   % 				Convert slot number and range
   % 				to percent of colormap
   PBeginSlot    = (BeginSlot - 1) / (CmLength - 1);
   PEndSlot      = (EndSlot - 1) / (CmLength - 1);
   PCmRange      = PEndSlot - PBeginSlot;
   % 				Determine range and min and max 
   % 				of new CLim values
   DataRange     = CDmax - CDmin;
   ClimRange     = DataRange / PCmRange;
   NewCmin       = CDmin - (PBeginSlot * ClimRange);
   NewCmax       = CDmax + (1 - PEndSlot) * ClimRange;
   CLim          = [NewCmin,NewCmax];
end % of function

function output = readnc(filename, varname)
% Returns a named variable from a netcdf file

[ncid, status] = mexnc('open', filename);
[varid, status] = mexnc('inq_varid', ncid, varname);
[output, status] = mexnc('get_var_double', ncid, varid);

mexnc('close', ncid);
end %of function