function ocr3000_plot(configfile)
% function ocr3000_plot(configfile)
% Reads ocr3000 data from netcdf files, 
% plots last 1 day of averaged data, combining datasets
% plots last 7 days of data, one data set at a time
% plots full series of data, one data set at a time
% and saves the processed results in several netCDF files.
%
% Created 2004-06-11 Michael Godin, godin@mbari.org
%
%-------------------------------------------------------------------------
% Requires:
%   readProperties.m
%   getProperty.m
%-------------------------------------------------------------------------

% Add netcdf to path if necessary
if(isempty (strfind(path, 'netcdf'))),
    path(path,[pwd '/netcdf']);
    path(path,[pwd '/netcdf/ncfiles']);    
    path(path,[pwd '/netcdf/nctype']);
    path(path,[pwd '/netcdf/ncutility']);
end

% set up some useful constants

wavelengths=[412 440 490 510 555 610 650 670 700];
suffixes    = {'_mean', '_median'};
suffixTitles    = {'Mean', 'Median'};

Eds  = 1;
Ed5  = 2;
Lu5  = 3;
Ed10 = 4;
Lu10 = 5;
% use similar mapping as ocr3000_derive
Kd_0_5 =  1 + 5;
Kd_0_10 = 2 + 5;
Kd_5_10 = 3 + 5;
Ku_5_10 = 4 + 5;
Ku_0_5  = 5 + 5;
Lw_Kd_0_5   = 6 + 5;
Lw_Ku_5_10  = 7 + 5;
Lw_Ku_0_5   = 8 + 5;
RSR_Kd_0_5  = 9 + 5;
RSR_Ku_5_10 = 10 + 5;
RSR_Ku_0_5  = 11 + 5;
bb670_Kd_0_5      = 12 + 5;
bb670_Ku_5_10     = 13 + 5;
bb670_Ku_0_5      = 14 + 5;
chl_490_Kd_0_5    = 15 + 5;
chl_490_Ku_5_10   = 16 + 5;
chl_490_Ku_0_5    = 17 + 5;
chl_oc2v4_Kd_0_5  = 18 + 5;
chl_oc2v4_Ku_5_10 = 19 + 5;
chl_oc2v4_Ku_0_5  = 20 + 5;

files     = {    'Eds','Ed5','Lu5','Ed10','Lu10',...
                 'Kd_0_5','Kd_0_10','Kd_5_10','Ku_5_10',...
                 'Ku_0_5',...
                 'Lw_Kd_0_5','Lw_Ku_5_10','Lw_Ku_0_5',...
                 'RSR_Kd_0_5','RSR_Ku_5_10','RSR_Ku_0_5',...
                 'bb670_Kd_0_5','bb670_Ku_5_10','bb670_Ku_0_5',...
                 'chl_490_Kd_0_5','chl_490_Ku_5_10','chl_490_Ku_0_5',...
                 'chl_oc2v4_Kd_0_5','chl_oc2v4_Ku_5_10','chl_oc2v4_Ku_0_5'
                 };
           
% Read in common properties
basePath     = pwd;
prop         = readProperties(configfile);
mooring.name = getProperty(prop, 'mooring');
radiometers = getProperty(prop, [mooring.name '.radiometers']);

% Parse out radiometers into cell array radiometer{1:radiometer_count}
radiometer_count = 0;
while(size(radiometers,2) > 0)
    radiometer_count = radiometer_count + 1;
    [radiometer{radiometer_count} radiometers] = strtok(radiometers,', ');
    % Ignore invalid entries
    if(size(radiometer{radiometer_count},2) < 1) 
        radiometer_count = radiometer_count - 1;
        continue
    end
    names{radiometer_count} = getProperty(prop, [radiometer{radiometer_count} '.name']);
end

netCDF_path = MBARI_path(getProperty(prop, ['netCDF_path']));
plot_path = MBARI_path(getProperty(prop, ['plot_path']));

for suffix_index = 1:length(suffixes)
    %if(0==1) % for debugging
    plot_lastday(netCDF_path, plot_path,...
            [suffixTitles{suffix_index} ' Downwelling Irradiances'],...
            { [files{Eds}  suffixes{suffix_index}],...
              [files{Ed5}  suffixes{suffix_index}],...
              [files{Ed10} suffixes{suffix_index}]},...
              ['Ed' suffixes{suffix_index}]);
       
    plot_lastday(netCDF_path, plot_path,...
            [suffixTitles{suffix_index} ' Upwelling Radiances'],...
            { [files{Lu5}  suffixes{suffix_index}],...
              [files{Lu10} suffixes{suffix_index}]},...
              ['Lu' suffixes{suffix_index}]);
      
    plot_lastday(netCDF_path, plot_path,...
            [suffixTitles{suffix_index} ' Attenuation Coefficients'],...
            { [files{Kd_0_5}  suffixes{suffix_index}],...
              [files{Kd_0_10}  suffixes{suffix_index}],...
              [files{Kd_5_10}  suffixes{suffix_index}],...
              [files{Ku_5_10}  suffixes{suffix_index}],...
              [files{Ku_0_5}  suffixes{suffix_index}]},...
              ['K' suffixes{suffix_index}]);
      
    plot_lastday(netCDF_path, plot_path,...
            [suffixTitles{suffix_index} ' Water Leaving Radiances'],...
            { [files{Lw_Kd_0_5}  suffixes{suffix_index}],...
              [files{Lw_Ku_5_10}  suffixes{suffix_index}],...
              [files{Lw_Ku_0_5}  suffixes{suffix_index}]},...
              ['Lw' suffixes{suffix_index}]);
      
    plot_lastday(netCDF_path, plot_path,...
            [suffixTitles{suffix_index} ' Remote Sensing Reflectances'],...
            { [files{RSR_Kd_0_5}  suffixes{suffix_index}],...
              [files{RSR_Ku_5_10}  suffixes{suffix_index}],...
              [files{RSR_Ku_0_5}  suffixes{suffix_index}]},...
              ['RSR' suffixes{suffix_index}]);
    %end % for debugging
    for(file_index = 1:length(files)),
        plot_timeseries(netCDF_path, plot_path,...
                [files{file_index}  suffixes{suffix_index}],...
                wavelengths, [suffixTitles{suffix_index} ' ']);
    end
end
disp ('ocr3000_plot done');
%quit
return 

function plot_lastday(netCDF_path, plot_path, plot_title, fileroots, png_root)
disp(['Plotting ' plot_title]);
figure('visible','off'); 
rect = [242 349 479 365]; 
set(gcf,'Position',rect);
set(gcf,'PaperPosition',[0 0 3.75 2.75]);
c = cool(length(fileroots));
for(file_index = 1: length(fileroots)),
    fileroot = fileroots{file_index};
    nc_file = netcdf([netCDF_path '/' fileroot '.nc']);
    max_valid = nc_file.max_valid(:);
    if(isempty(max_valid)), max_valid = inf; end
    min_valid = nc_file.min_valid(:);
    if(isempty(min_valid)), min_valid = -inf; end
    % find the multi-dimensional variable
    for vars = length(ncnames(var(nc_file))) : -1 : 1,
        var_dims = length(dim(nc_file{vars}));
        if(var_dims > 1)
            var_name = char(ncnames(nc_file{vars}));
            break;
        end
    end
    all_times = nc_file{'time'}(:);
    first_time = max(all_times) - 86400;
    start_index = min(find(all_times > first_time));
    spectra = squeeze(nc_file{var_name}(start_index:end,:));
    spectra(find(spectra<min_valid)) = NaN;
    spectra(find(spectra>max_valid)) = NaN;
    spectrum = nanmean(spectra);
    waves = nc_file{'wavelength'}(:);
    %c(file_index,:) = c(file_index,:) / (floor(file_index/length(fileroots))/2+1);
    plot(waves', spectrum, 'Color',c(file_index,:),'LineWidth',1)
    if(file_index == 1)
           xlabel('\lambda (nm)');
           ylabel(nc_file{var_name}.units(:));
    end
    labels{file_index} = nc_file{var_name}.long_name(:);
    sorty = spectrum;
    %sorty(find(isnan(sorty)))=max(sorty);
    sorty = sort(sorty);
    if(length(sorty) > sum(isnan(sorty))),
        minY(file_index) = nanmedian(sorty(1:ceil((length(sorty)-sum(isnan(sorty)))/5)));
        maxY(file_index) = nanmedian(sorty(end-sum(isnan(sorty))-floor((length(sorty)-sum(isnan(sorty)))/5):end-sum(isnan(sorty))));
    else
        minY(file_index) = 0;
        maxY(file_index) = 1;
    end
    hold on;
    close(nc_file);
end
minnY = min(minY);
maxxY = max(maxY);
deltaY = maxxY - minnY;
minnY = minnY-deltaY/5;
maxxY = maxxY+deltaY/5;
%disp([plot_title ' ' num2str(minnY) ' ' num2str(maxxY)]);
if(isnan(minnY)) minnY = 0; end
if(isnan(maxxY)) maxxY = 1; end
set(gca,'YLim',[minnY, maxxY]);
hold off;
title([plot_title ' - 24 hour mean']) %  until '...
          %  datestr(all_times(start_index)/86400+datenum(1970,1,1))...
          %          ' GMT']);
setfontsize(7);
if ~isempty(labels)
    lgnd = legend(char(labels),0);
    set(lgnd, 'Fontsize',5);
end 
axis tight;
ylim = get(gca,'YLim');
dy = ylim(2)-ylim(1);
txt = text(min(waves), ylim(1)-dy/10, ['Last data: ' datestr(all_times(end)/86400+datenum(1970,1,1))], 'Fontsize', 5);
print('-dpng', '-painters', [plot_path '/' png_root '_day.png']);
close(gcf);
return

function plot_timeseries(netCDF_path, plot_path, fileroot, wavelengths, title_prefix)
disp(['Plotting ' fileroot]);
nc_file = netcdf([netCDF_path '/' fileroot '.nc']);
% First, look for the multidimensional variable
var_name = '';
for vars = length(ncnames(var(nc_file))) : -1 : 1,
    var_dims = length(dim(nc_file{vars}));
    if(var_dims > 1)
        var_name = char(ncnames(nc_file{vars}));
        break;
    end
end
if(isempty (var_name))
    % No multi-dimensional variable - find the last defined single
    % dimensional variable
	for vars = length(ncnames(var(nc_file))) : -1 : 1,
        var_dims = length(dim(nc_file{vars}));
        if(var_dims == 1)
            var_name = char(ncnames(nc_file{vars}));
            break;
        end
	end
end   
all_times = nc_file{'time'}(:);
% plot 7 day data
first_time = max(all_times) - 7*86400;
start_index = min(find(all_times > first_time));
generate_timeseries(fileroot, nc_file, start_index, var_name, var_dims, wavelengths, plot_path, '_7day', title_prefix, ' - 7 days')
% plot all data
generate_timeseries(fileroot, nc_file, 1, var_name, var_dims, wavelengths, plot_path, '_all', title_prefix, ' - all data')
close(nc_file);
return

function generate_timeseries(fileroot, nc_file, start_index, var_name, var_dims, wavelengths, plot_path, plot_suffix, title_prefix, title_suffix)
figure('visible','off'); 
rect = [242 349 479 365]; 
set(gcf,'Position',rect);
set(gcf,'PaperPosition',[0 0 3.75 2.75]);
max_valid = nc_file.max_valid(:);
if(isempty(max_valid)), max_valid = inf; end
min_valid = nc_file.min_valid(:);
if(isempty(min_valid)), min_valid = -inf; end
if(size(squeeze(nc_file{var_name}(:)),2) > 1)
    waves = nc_file{'wavelength'}(:);
    wave_indices = 1: length(waves);
    for wave_index = 1: length(wavelengths),
        interp_index = interp1(waves, wave_indices, wavelengths(wave_index));
        lower_index = floor(interp_index);
        upper_index = ceil(interp_index);
        if(lower_index ~= upper_index)
            lower_frac = upper_index - interp_index;
            upper_frac = interp_index - lower_index;
            switch(var_dims)
                case 2,
                    lower_series = squeeze(nc_file{var_name}(start_index:end,lower_index));
                    upper_series = squeeze(nc_file{var_name}(start_index:end,upper_index));
                case 3,
                    lower_series = squeeze(nc_file{var_name}(start_index:end,:,lower_index));
                    upper_series = squeeze(nc_file{var_name}(start_index:end,:,upper_index));
                case 4,
                    lower_series = squeeze(nc_file{var_name}(start_index:end,:,:,lower_index));
                    upper_series = squeeze(nc_file{var_name}(start_index:end,:,:,upper_index));
                case 5,
                    lower_series = squeeze(nc_file{var_name}(start_index:end,:,:,:,lower_index));
                    upper_series = squeeze(nc_file{var_name}(start_index:end,:,:,:,upper_index));
                case 6,
                    lower_series = squeeze(nc_file{var_name}(start_index:end,:,:,:,:,lower_index));
                    upper_series = squeeze(nc_file{var_name}(start_index:end,:,:,:,:,upper_index));
            end
            series = lower_frac * lower_series + upper_frac * upper_series;
        else
            switch(var_dims)
                case 2,
                    series = squeeze(nc_file{var_name}(start_index:end,lower_index));
                case 3,
                    series = squeeze(nc_file{var_name}(start_index:end,:,lower_index));
                case 4,
                    series = squeeze(nc_file{var_name}(start_index:end,:,:,lower_index));
                case 5,
                    series = squeeze(nc_file{var_name}(start_index:end,:,:,:,lower_index));
                case 6,
                    series = squeeze(nc_file{var_name}(start_index:end,:,:,:,:,lower_index));
            end
        end
        c = wave2rgb(wavelengths(wave_index),1);
        series(find(series>max_valid))=nan;
        series(find(series<min_valid))=nan;
        plot(nc_file{'time'}(start_index:end)/86400+datenum(1970,1,1), series, 'Color',c);
        labels{wave_index} = [num2str(wavelengths(wave_index)) ' nm'];
        sorty = series;
        %sorty(find(isnan(sorty)))=max(sorty);
        sorty = sort(sorty);
        if(length(sorty) > sum(isnan(sorty))),
            minY(wave_index) = nanmedian(sorty(1:ceil((length(sorty)-sum(isnan(sorty)))/5)));
            maxY(wave_index) = nanmedian(sorty(end-sum(isnan(sorty))-floor((length(sorty)-sum(isnan(sorty)))/5):end-sum(isnan(sorty))));
        else
            minY(wave_index) = 0;
            maxY(wave_index) = 1;
        end
        hold on;
    end % for wave_index = 1: length(wavelengths),
else % if(var_dims > 1)
    series = squeeze(nc_file{var_name}(start_index:end));
    series(find(series>max_valid))=nan;
    series(find(series<min_valid))=nan;
    plot(nc_file{'time'}(start_index:end)/86400+datenum(1970,1,1), series);
    sorty = series;
    %sorty(find(isnan(sorty)))=max(sorty);
    sorty = sort(sorty);
    if(length(sorty) > sum(isnan(sorty))),
        minY(1) = nanmedian(sorty(1:ceil((length(sorty)-sum(isnan(sorty)))/5)));
        maxY(1) = nanmedian(sorty(end-sum(isnan(sorty))-floor((length(sorty)-sum(isnan(sorty)))/5):end-sum(isnan(sorty))));
    else
        minY(1) = 0;
        maxY(1) = 1;
    end
    labels = {};
end
ylabel(nc_file{var_name}.units(:));
minnY = min(minY);
maxxY = max(maxY);
deltaY = maxxY - minnY;
minnY = minnY-deltaY/5;
maxxY = maxxY+deltaY/5;
%disp([plot_title ' ' num2str(minnY) ' ' num2str(maxxY)]);
if (~isnan(minnY) && ~isnan(maxxY) && (minnY < maxxY)),
    set(gca,'YLim',[minnY, maxxY]);
end
hold off;
title([title_prefix nc_file{var_name}.long_name(:) title_suffix]);
setfontsize(7);
if ~isempty(labels)
    lgnd = legend(char(labels),2);
    set(lgnd, 'Fontsize',5);
end 
%axis tight;
minTime = nc_file{'time'}(start_index)/86400+datenum(1970,1,1);
maxTime = nc_file{'time'}(end)/86400+datenum(1970,1,1);
dTime = maxTime-minTime;
timeStep = ceil(dTime/7);
datetick('x',25);
xlabel('date (GMT)','FontWeight','bold');
set(gca,'xtick',[floor(minTime):timeStep:ceil(maxTime)]);
ylim = get(gca,'YLim');
dy = ylim(2)-ylim(1);
xlim = get(gca,'XLim');
txt = text(xlim(1), ylim(1)-dy/10, ['Last data: ' datestr(maxTime)], 'Fontsize', 5);
print('-dpng', '-painters', [plot_path '/' fileroot plot_suffix '.png']);
close(gcf);

return