function ocr3000_derive(configfile,Ed_cutoff,Lu_cutoff)
% function ocr3000_derive
%
% Reads ocr3000 data from netcdf files, calcualtes Kds and Kus, 
% 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

if(nargin < 2)
    Ed_cutoff = 0.001 % limit for usable derived data
end
if(nargin < 3)
    Lu_cutoff = 0.00005 % limit for usable derived data
end

cutoff = [Ed_cutoff, Ed_cutoff, Lu_cutoff,  Ed_cutoff, Lu_cutoff];
    
% set up some useful constants
Eds  = 1;
Ed5  = 2;
Lu5  = 3;
Ed10 = 4;
Lu10 = 5;
Kd_0_5 =  1;
Kd_0_10 = 2;
Kd_5_10 = 3;
Ku_5_10 = 4;
Ku_0_5  = 5;
Lw_Kd_0_5   = 6;
Lw_Ku_5_10  = 7;
Lw_Ku_0_5   = 8;
RSR_Kd_0_5  = 9;
RSR_Ku_5_10 = 10;
RSR_Ku_0_5  = 11;
bb670_Kd_0_5      = 12;
bb670_Ku_5_10     = 13;
bb670_Ku_0_5      = 14;
chl_490_Kd_0_5    = 15;
chl_490_Ku_5_10   = 16;
chl_490_Ku_0_5    = 17;
chl_oc2v4_Kd_0_5  = 18;
chl_oc2v4_Ku_5_10 = 19;
chl_oc2v4_Ku_0_5  = 20;

out_files     = {'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'
                 };
out_variables = {'attenuation_downwelling_0m_5m',...
                 'attenuation_downwelling_0m_10m',...
                 'attenuation_downwelling_5m_10m',...
                 'attenuation_upwelling_5m_10m',...
                 'attenuation_upwelling_derived_0m_5m',...
                 'water_leaving_radiance_downwelling_0m_5m',...
                 'water_leaving_radiance_upwelling_5m_10m',...
                 'water_leaving_radiance_derived_0m_5m',...
                 'surface_remote_sensing_reflectance_downwelling_0m_5m',...
                 'surface_remote_sensing_reflectance_upwelling_5m_10m',...
                 'surface_remote_sensing_reflectance_derived_0m_5m',...
                 'backscatter_670_downwelling_0m_5m',...
                 'backscatter_670_upwelling_5m_10m',...
                 'backscatter_670_derived_0m_5m',...
                 'chlorophyll_490_downwelling_0m_5m',...
                 'chlorophyll_490_upwelling_5m_10m',...
                 'chlorophyll_490_derived_0m_5m',...
                 'chlorophyll_oc2v4_downwelling_0m_5m',...
                 'chlorophyll_oc2v4_upwelling_5m_10m',...
                 'chlorophyll_oc2v4_derived_0m_5m'
                 };
out_long_names = {'Downwelling irradiance attenuation from 0m to 5m',...
                 'Downwelling irradiance attenuation from 0m to 10m',...
                 'Downwelling irradiance attenuation from 5m to 10m',...
                 'Upwelling radiance attenuation from 5m to 10m',...
                 'Upwelling attenuation derived from 5m to 10m',...
                 'Water leaving radiance from downwelling 0m to 5m',...
                 'Water leaving radiance from upwelling 5m to 10m',...
                 'Water leaving radiance from derived 0m to 5m',...
                 'Surface Remote sensing reflectance from downwelling 0m to 5m',...
                 'Surface Remote sensing reflectance from upwelling 5m to 10m',...
                 'Surface Remote sensing reflectance from derived 0m to 5m',...
                 'Backscatter at 670nm from downwelling 0m to 5m',...
                 'Backscatter at 670nm from upwelling 5m to 10m',...
                 'Backscatter at 670nm from derived 0m to 5m',...
                 'Chlorophyll 490 from downwelling 0m to 5m',...
                 'Chlorophyll 490 from upwelling 5m to 10m',...
                 'Chlorophyll 490 from derived 0m to 5m',...
                 'Chlorophyll oc2v4 from downwelling 0m to 5m',...
                 'Chlorophyll oc2v4 from upwelling 5m to 10m',...
                 'Chlorophyll oc2v4 from derived 0m to 5m'
                 };
out_units =      {'m^-1',...
                 'm^-1',...
                 'm^-1',...
                 'm^-1',...
                 'm^-1',...
                 'W/m^2/nm/sr',...
                 'W/m^2/nm/sr',...
                 'W/m^2/nm/sr',...
                 '1/sr',...
                 '1/sr',...
                 '1/sr',...
                 'm^-1',...
                 'm^-1',...
                 'm^-1',...
                 'ug/l',...
                 'ug/l',...
                 'ug/l',...
                 'ug/l',...
                 'ug/l',...
                 'ug/l'                 
                 };
upper_indices = [Eds, Eds,  Ed5,  Lu5,  Eds, Eds, Lu5,  Eds, Eds, Lu5,  Eds, Eds, Lu5,  Eds, Eds, Lu5,  Eds, Eds, Lu5,  Eds];
lower_indices = [Ed5, Ed10, Ed10, Lu10, Lu5, Ed5, Lu10, Lu5, Ed5, Lu10, Lu5, Ed5, Lu10, Lu5, Ed5, Lu10, Lu5, Ed5, Lu10, Lu5];

% Scientific constants
aw_670 = 0.430; % for calculating bb_670
a_oc2v4 = [0.319 -2.366 0.879 -0.135 -0.071];  % for oc2v4 chl algorythm
                 
% Read in common properties
basePath     = pwd;
prop         = readProperties(configfile);
mooring.name = getProperty(prop, 'mooring');
radiometers = getProperty(prop, [mooring.name '.radiometers']);
latitude = str2num(getProperty(prop, [mooring.name '.latitude']));
longitude = str2num(getProperty(prop, [mooring.name '.longitude']));

% 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

% Determine when to start reading in raw data
netCDF_path = MBARI_path(getProperty(prop, ['netCDF_path']));
for out_file_index = 1:length(out_files),
    ncfilename = [netCDF_path '/' out_files{out_file_index} '_mean.nc'];
    ncfile = netcdf(ncfilename, 'nowrite');
    if(size(ncfile,2)==0 || isempty(ncfile{'time'}(end))),
        last_time_stamp(out_file_index) = 0;
    else
        last_time_stamp(out_file_index) = ncfile{'time'}(end);
    end
    close(ncfile);
end
start_time_stamp = min(last_time_stamp);

suffixes    = {'_mean', '_median'};
for suffix_index = 1:length(suffixes)
    
    % Read in new data
    for radiometer_index = 1:radiometer_count,
        ncfile = netcdf([netCDF_path '/' radiometer{radiometer_index} suffixes{suffix_index} '.nc']);
        disp(['Reading ' radiometer{radiometer_index} suffixes{suffix_index} '...']);
        alltime = ncfile{'time'}(:);
        start_index = min(find(alltime >= start_time_stamp))-1;
        if ( isempty(start_index)), start_index= 0; end
        if(start_index > 0)
            time{radiometer_index} = ncfile{'time'}(start_index : end);
            %spectra{radiometer_index} = squeeze(ncfile{names{radiometer_index}}(start_index : end,1,1,1,:));
            spectra{radiometer_index} = squeeze(ncfile{names{radiometer_index}}(start_index : end,:));
        else
            time{radiometer_index} = ncfile{'time'}(:);
            %spectra{radiometer_index} = squeeze(ncfile{names{radiometer_index}}(:,1,1,1,:));
            spectra{radiometer_index} = squeeze(ncfile{names{radiometer_index}}(:,:));
        end
        %depth(radiometer_index) = ncfile{'depth'}(:);
        depth(radiometer_index) = ncfile.depth(:);
        if(radiometer_index == 1), wave_grid = ncfile{'wavelength'}(:) ; end
        close(ncfile);
    end
    
    % Apply cutoff
    for radiometer_index = 1:radiometer_count,
        spectra{radiometer_index}(find(spectra{radiometer_index} < cutoff(radiometer_index))) = NaN;
    end
        
    % Calculate index of refraction
    index_ref = wave_grid * 0 + 1.34; % water_index(wave_grid);  % index of refraction
    
    % Calculate time-dependent transmission, and index of refraction squared
    mTime = time{1}/86400+datenum(1970,1,1);
    if (length(mTime) == 0) return, end
    solZen = solarzenith(mTime, latitude, longitude); % 8 hours difference between time stamps and solar time
    for(time_index = 1:length(mTime)),
        for(wave_index = 1:length(index_ref)),
             Transmission(time_index,wave_index) = 1-fresnel_(solZen(time_index), index_ref(wave_index));
             index2(time_index,wave_index) = index_ref(wave_index)^2;
        end
    end     
    Transmission(find(Transmission <= 0)) = NaN;
    Transmission(find(Transmission >= 1)) = NaN;
    
    % Create output files
    for out_file_index = 1:length(out_files),
        ncfilename = [netCDF_path '/' out_files{out_file_index} suffixes{suffix_index} '.nc'];
        ncfile = netcdf(ncfilename, 'write');
        channel = out_variables{out_file_index};
        disp(['Calculating ' channel '...']);
        if(size(ncfile,2)==0 || isempty(ncfile{'time'}(end))),
            close(ncfile);    
            ncfile = netcdf(ncfilename, 'clobber');
            %Initialize Attributes
            ncfile.waterdepth                               = getProperty(prop, [mooring.name '.waterdepth']);
            ncfile.waterdepth_units                         = getProperty(prop, [mooring.name '.waterdepth_units']);
			[Year Month Day]                                = datevec(now);
			dateS                                           = sprintf('%02i/%02i/%04i', Month, Day, Year);
			ncfile.creationdate            = dateS; %datestr(now)
            ncfile.mooring                 = mooring.name;

            % Initialize Dimensions        
            ncfile('time')                 = 0; % Dimension 'time' is the Record Dimension
            ncfile{'time'}                 = nclong('time'); % Varaible 'time' stored as double
            ncfile{'time'}.long_name       = getProperty(prop, 'time.long_name'); 
            ncfile{'time'}.units           = getProperty(prop, 'time.units');

            %ncfile('latitude')             = 1;
            %ncfile{'latitude'}             = ncdouble('latitude'); 
            %ncfile{'latitude'}.long_name   = getProperty(prop, 'latitude.long_name'); 
            %ncfile{'latitude'}.units       = getProperty(prop, 'latitude.units');
            %ncfile{'latitude'}.AXIS        = getProperty(prop, 'latitude.AXIS');
            %ncfile{'latitude'}(:)          = latitude;
            ncfile.latitude                 = latitude;

            %ncfile('longitude')            = 1;
            %ncfile{'longitude'}            = ncdouble('longitude'); 
            %ncfile{'longitude'}.long_name  = getProperty(prop, 'longitude.long_name'); 
            %ncfile{'longitude'}.units      = getProperty(prop, 'longitude.units');
            %ncfile{'longitude'}.AXIS       = getProperty(prop, 'longitude.AXIS');
            %ncfile{'longitude'}(:)         = longitude;
            ncfile.longitude               = longitude;

            %ncfile('depth_upper')          = 1;
            %ncfile{'depth_upper'}          = ncdouble('depth_upper'); 
            %ncfile{'depth_upper'}.long_name= 'water depth of upper measuremnt'; 
            %ncfile{'depth_upper'}.units    = getProperty(prop, 'depth.units');
            %ncfile{'depth_upper'}.AXIS     = getProperty(prop, 'depth.AXIS');
            %ncfile{'depth_upper'}.positive = getProperty(prop, 'depth.positive');
            %ncfile{'depth_upper'}(:)       = depth(upper_indices(out_file_index));
            ncfile.depth_upper             = depth(upper_indices(out_file_index));

            %ncfile('depth_lower')          = 1;
            %ncfile{'depth_lower'}          = ncdouble('depth_lower'); 
            %ncfile{'depth_lower'}.long_name= 'water depth of lower measuremnt'; 
            %ncfile{'depth_lower'}.units    = getProperty(prop, 'depth.units');
            %ncfile{'depth_lower'}.AXIS     = getProperty(prop, 'depth.AXIS');
            %ncfile{'depth_lower'}.positive = getProperty(prop, 'depth.positive');
            %ncfile{'depth_lower'}(:)       = depth(lower_indices(out_file_index));
            ncfile.depth_lower             = depth(lower_indices(out_file_index));

            ncfile('wavelength')           = length(wave_grid);
            ncfile{'wavelength'}           = ncdouble('wavelength'); 
            ncfile{'wavelength'}.long_name = getProperty(prop, 'wavelength.long_name'); 
            ncfile{'wavelength'}.units     = getProperty(prop, 'wavelength.units');
            ncfile{'wavelength'}(:)        = wave_grid;
            
            % Initialize Variables      
            if(out_file_index < bb670_Kd_0_5)
                %ncfile{channel}                = ncfloat('time','latitude','longitude','depth_upper','depth_lower','wavelength');
                ncfile{channel}                = ncfloat('time','wavelength');
            else
                %ncfile{channel}                = ncfloat('time','latitude','longitude','depth_upper','depth_lower');
                ncfile{channel}                = ncfloat('time');
            end
            ncfile{channel}.units          = out_units{out_file_index};
            ncfile{channel}.long_name      = out_long_names{out_file_index};
            ncfile{channel}.FillValue_     = -9999;
            ncfile{channel}.missing_value  = -9999;    
        end            
        
        % do the calculations....
        switch(out_file_index)
            case {Kd_0_5}
                F{Kd_0_5}  = (log(spectra{Eds}.* Transmission) - log(spectra{Ed5}) )  / ...
                                depth(Ed5) ;
                max_valid = 10;
                min_valid = -1;
            case {Kd_0_10}
                F{Kd_0_10} = (log(spectra{Eds}.* Transmission) - log(spectra{Ed10}) )  / ...
                                depth(Ed10) ;
                max_valid = 10;
                min_valid = -1;
            case {Kd_5_10}
                F{Kd_5_10} = (log(spectra{Ed5}) - log(spectra{Ed10}))  / ...
                                (depth(Ed10) - depth(Ed5)) ;
                max_valid = 10;
                min_valid = -1;
            case {Ku_5_10}
                F{Ku_5_10} = (log(spectra{Lu5}) - log(spectra{Lu10}))  / ...
                                (depth(Lu10) - depth(Lu5)) ;
                max_valid = 10;
                min_valid = -1;
            case {Ku_0_5}
                F{Ku_0_5}  = F{Kd_0_5} .* F{Ku_5_10} ./ F{Kd_5_10};
                max_valid = 10;
                min_valid = -1;
            case {Lw_Kd_0_5}
                F{Lw_Kd_0_5}  = spectra{Lu5} ./ index2 .* Transmission .* exp(F{Kd_5_10} * depth(Lu5));
                max_valid = 1;
                min_valid = 0;
            case {Lw_Ku_5_10}
                F{Lw_Ku_5_10} = spectra{Lu5} ./ index2 .* Transmission .* exp(F{Ku_5_10} * depth(Lu5));
                max_valid = 1;
                min_valid = 0;
            case {Lw_Ku_0_5}
                F{Lw_Ku_0_5}  = spectra{Lu5} ./ index2 .* Transmission .* exp(F{Ku_0_5}  * depth(Lu5));
                max_valid = 1;
                min_valid = 0;
            case {RSR_Kd_0_5}
                F{RSR_Kd_0_5}  = F{Lw_Kd_0_5}  ./ spectra{Eds}; 
                max_valid = 10;
                min_valid = 0;
            case {RSR_Ku_5_10}
                F{RSR_Ku_5_10} = F{Lw_Ku_5_10} ./ spectra{Eds}; 
                max_valid = 10;
                min_valid = 0;
            case {RSR_Ku_0_5}
                F{RSR_Ku_0_5}  = F{Lw_Ku_0_5}  ./ spectra{Eds}; 
                max_valid = 10;
                min_valid = 0;
            case {bb670_Kd_0_5}
                F{bb670_Kd_0_5}  = interp1(wave_grid,F{RSR_Kd_0_5}',670) * aw_670 / 0.051;
                max_valid = 10;
                min_valid = 0;
            case {bb670_Ku_5_10}
                F{bb670_Ku_5_10} = interp1(wave_grid,F{RSR_Ku_5_10}',670)* aw_670 / 0.051;
                max_valid = 10;
                min_valid = 0;
            case {bb670_Ku_0_5}
                F{bb670_Ku_0_5}  = interp1(wave_grid,F{RSR_Ku_0_5}',670) * aw_670 / 0.051;
                max_valid = 10;
                min_valid = 0;
            case {chl_490_Kd_0_5}
                F{chl_490_Kd_0_5}  = 0.2818 * (( interp1(wave_grid,F{Lw_Kd_0_5}',412) + ...
                                               interp1(wave_grid,F{Lw_Kd_0_5}',565) ) ./ ...
                                                interp1(wave_grid,F{Lw_Kd_0_5}',490)) .^ 3.497;
                max_valid = 1000;
                min_valid = 0;
            case {chl_490_Ku_5_10}
                F{chl_490_Ku_5_10} = 0.2818 * (( interp1(wave_grid,F{Lw_Ku_5_10}',412) + ...
                                               interp1(wave_grid,F{Lw_Ku_5_10}',565) ) ./ ...
                                                interp1(wave_grid,F{Lw_Ku_5_10}',490)) .^ 3.497;
                max_valid = 1000;
                min_valid = 0;
            case {chl_490_Ku_0_5}
                F{chl_490_Ku_0_5}  = 0.2818 * (( interp1(wave_grid,F{Lw_Ku_0_5}',412) + ...
                                               interp1(wave_grid,F{Lw_Ku_0_5}',565) ) ./ ...
                                                interp1(wave_grid,F{Lw_Ku_0_5}',490)) .^ 3.497;
                max_valid = 1000;
                min_valid = 0;
            case {chl_oc2v4_Kd_0_5}
                R = log10(interp1(wave_grid,F{RSR_Kd_0_5}',490) ./ interp1(wave_grid,F{RSR_Kd_0_5}',555));
                F{chl_oc2v4_Kd_0_5}  = 10.^(a_oc2v4(1) + a_oc2v4(2)*R + a_oc2v4(3)*R.^2 + a_oc2v4(4)*R.^3) + a_oc2v4(5);
                max_valid = 1000;
                min_valid = 0;
            case {chl_oc2v4_Ku_5_10}
                R = log10(interp1(wave_grid,F{RSR_Ku_5_10}',490) ./ interp1(wave_grid,F{RSR_Ku_5_10}',555));
                F{chl_oc2v4_Ku_5_10}  = 10.^(a_oc2v4(1) + a_oc2v4(2)*R + a_oc2v4(3)*R.^2 + a_oc2v4(4)*R.^3) + a_oc2v4(5);
                max_valid = 1000;
                min_valid = 0;
            case {chl_oc2v4_Ku_0_5}
                R = log10(interp1(wave_grid,F{RSR_Ku_0_5}',490) ./ interp1(wave_grid,F{RSR_Ku_0_5}',555));
                F{chl_oc2v4_Ku_0_5}  = 10.^(a_oc2v4(1) + a_oc2v4(2)*R + a_oc2v4(3)*R.^2 + a_oc2v4(4)*R.^3) + a_oc2v4(5);
                max_valid = 1000;
                min_valid = 0;
        end
        % We only care about real values here.
        F{out_file_index}(find(~isreal(F{out_file_index}))) = NaN;
        % Write the calculated variables
        nc_index = length(ncfile{'time'}(:));
        nc_end_time = ncfile{'time'}(end);
        start_index = 1;
        max_index = length(time{1});
        while(~isempty(nc_end_time) && start_index < max_index && time{1}(start_index) <= nc_end_time ),
           start_index = start_index + 1;
        end
        
        if(start_index <= max_index),
            for time_index = start_index : max_index,
                ncfile{'time'}(nc_index)              = time{1}(time_index);
                if(out_file_index < bb670_Kd_0_5)
                    %ncfile{channel}(nc_index,1,1,1,:) = F{out_file_index}(time_index,:);
                    ncfile{channel}(nc_index,:) = F{out_file_index}(time_index,:);
                else
                    %ncfile{channel}(nc_index,1,1,1)   = F{out_file_index}(time_index);
                    ncfile{channel}(nc_index)   = F{out_file_index}(time_index);
                end
                nc_index = nc_index + 1;         
            end
        else
            disp(['No new ' out_files{out_file_index} ' data to add to files.']);
        end
        ncfile.min_valid = min_valid;
        ncfile.max_valid = max_valid;
        % Close output file
        close(ncfile);
    end
end
disp ('ocr3000_derive done');
%quit