function ocr3000_derive(deployment,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
% Modified 2006-09-27 Seth Bushinsky, sethb@mbari.org
% 
% Attenuation coefficient equation can be found in: Morel and Maritorena,
% "Bio-Optical properties of oceanic waters: A reappraisal", JGR - Oceans,
% Vol 106. 2001.
%    Kd(wavelength) =  Z^-1 ln[Eds(wavelength)/EdZ (wavelength)]
%
%-------------------------------------------------------------------------
% Requires:
%   readProperties.m
%   getProperty.m
%-------------------------------------------------------------------------

% Add netcdf to path if necessary
toolbox_dir = '/u/ssdsadmin/dev/DPforSSDS/';
if(isempty (strfind(path, 'netcdf_toolbox'))),
    addpath('/usr/local/DODS_ml/bin');
    addpath([toolbox_dir 'mexnc']);
    addpath([toolbox_dir 'mexcdf/snctools']);
    addpath([pwd '/netcdf_toolbox/netcdf']);
    addpath([pwd '/netcdf_toolbox/netcdf/nctype']);
    addpath([pwd '/netcdf_toolbox/netcdf/ncutility']);
end

if(nargin < 3)
    Ed_cutoff = 0.0005; % limit for usable derived data
    disp(['Ed_cutoff' Ed_cutoff])
end
if(nargin < 4)
    Lu_cutoff = 0.00001; % limit for usable derived data
    disp(['Lu_cutoff' Lu_cutoff])
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;
chl_490_Kd_0_5    = 21;
chl_490_Kd_0_10    = 22;
Kd_0_5_95 = 23;
Kd_0_10_95 = 24;

% To turn on other derived variables, change the outfile_used parameter
outfile_used = [1:5,21,22];



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',...
                 'chl_490_Kd_0_5','chl_490_Kd_0_10','Kd_0_5_95','Kd_0_10_95'
                 };
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',...
                 'chlorophyll_490_downwelling_0m_5m_M88',...
                 'chlorophyll_490_downwelling_0m_10m_M88',...
                 'attenuation_downwelling_0m_5m95',...
                 'attenuation_downwelling_0m_10m95'
                 };
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',...
                 'chlorophyll 490 downwelling 0m 5m Morel 1988',...
                 'chlorophyll 490 downwelling 0m 10m Morel 1988',...
                 'Downwelling irradiance attenuation from 0m to 5m - 95',...
                 'Downwelling irradiance attenuation from 0m to 10m - 95'
                 };
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',...
                 'mg/m^3',...
                 'mg/m^3',...
                 'm^-1',...
                 'm^-1'
                 };
upper_indices = [Eds, Eds,  Ed5,  Lu5,  Eds, Eds, Lu5,  Eds, Eds, Lu5,  Eds, Eds, Lu5,  Eds, Eds, Lu5,  Eds, Eds, Lu5,  Eds, Eds, Eds,  Eds, Eds];
lower_indices = [Ed5, Ed10, Ed10, Lu10, Lu5, Ed5, Lu10, Lu5, Ed5, Lu10, Lu5, Ed5, Lu10, Lu5, Ed5, Lu10, Lu5, Ed5, Lu10, Lu5, Ed5, Ed10, Ed5, Ed10];

% 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 '.' deployment]);
mooring.name = getProperty(prop, 'mooring');
radiometers = getProperty(prop, [mooring.name '.radiometers']);
latitude = str2double(getProperty(prop, [mooring.name '.latitude']));
longitude = str2double(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'));
netCDF_path = [netCDF_path '/' deployment];
for out_file_index = 1:length(out_files),
    ncfilename = [netCDF_path '/' out_files{out_file_index} '_mean.nc'];
    ncid = fopen(ncfilename);
    if (ncid~=-1),
        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);
        close(ncfile);
    else
        last_time_stamp(out_file_index) = 0;
    end
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 start_index==0,
            start_index=1;
        end
        
        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,:));
        else
            disp ('Missing Radiometer Data: No Data in file');
            time{radiometer_index} = 0; %ncfile_prev{'time'}(start_index_prev : end);
            spectra{radiometer_index} = 0; %squeeze(ncfile_prev{names{radiometer_index}}(start_index_prev : end));
        end
        disp (['Read ' num2str(length(time{radiometer_index})) ' samples.']);
        depth(radiometer_index) = ncfile.depth(:);
        if(radiometer_index == 1), wave_grid = ncfile{'wavelength'}(:) ; end
        close(ncfile);
    end
   
    % Apply cutoff
    disp('Applying cutoffs'); 
    for radiometer_index = 1:radiometer_count,
        spectra{radiometer_index}(find(spectra{radiometer_index} < cutoff(radiometer_index))) = NaN;
    end
        
    % Calculate index of refraction
    disp('Calculating 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
    % Do so for each different radiometer (different times...)
    for rad = 1:radiometer_count,
        %mTime is matlab time, time is in epoch seconds
        mTime{rad} = time{rad}/86400+datenum(1970,1,1);
        if (length(mTime{rad}) == 0), return, end
        disp('Calculating solar zenith angles');
        solZen{rad} = solarzenith(mTime{rad}(:), latitude, longitude);  % 8 hours difference between time stamps and solar time
        disp('Calculating transmissions');
%         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
%         for(time_index = 1:length(mTime)),
        index2{rad} = zeros(length(mTime{rad}(:)), length(index_ref));
        Transmission{rad} = index2{rad}(:,:);
        for wave_index = 1:length(index_ref),
            Transmission{rad}(:,wave_index) = 1-fresnel_(solZen{rad}(:), index_ref(wave_index));
            index2{rad}(:,wave_index) = index_ref(wave_index)^2; % index2 is the squared index of refraction
        end

        Transmission{rad}(find(Transmission{rad} <= 0)) = NaN;
        Transmission{rad}(find(Transmission{rad} >= 1)) = NaN;

    end
    %end


    %creates arrays Ed5_match and Eds_match that contain the indices for
    %each time match between the two radiometers.  This section is
    %necessary because the derive program no longer inserts blanks to match
    %each time record to the other radiometers.  As a result, for each
    %derived variable, a new time record that only contains matching times
    %must be created.
    for tt=1:length(time{Eds}),
        if tt==1, Ed5_match_Eds=[];Eds_match_Ed5=[]; end
        temp= find(time{Ed5}==time{Eds}(tt));
        if isempty(temp),
            continue
        else
            Ed5_match_Eds(end+1)=temp;
            Eds_match_Ed5(end+1)=tt;
        end
        clear temp;
    end

    %time matches between Eds and Ed10
    for tt=1:length(time{Eds}),
        if tt==1, Ed10_match_Eds=[];Eds_match_Ed10=[]; end
        temp= find(time{Ed10}==time{Eds}(tt));
        if isempty(temp),
            continue
        else
            Ed10_match_Eds(end+1)=temp;
            Eds_match_Ed10(end+1)=tt;
        end
        clear temp;
    end

    %time matches between Ed5 and Ed10
    for tt=1:length(time{Ed5}),
        if tt==1, Ed10_match_Ed5=[];Ed5_match_Ed10=[]; end
        temp = find(time{Ed10}==time{Ed5}(tt));
        if isempty(temp),
            continue
        else
            Ed10_match_Ed5(end+1)=temp;
            Ed5_match_Ed10(end+1)=tt;
        end
        clear temp;
    end

    %time matches between Lu5 and Lu10
    for tt=1:length(time{Lu5}),
        if tt==1, Lu10_match_Lu5=[];Lu5_match_Lu10=[]; end
        temp = find(time{Lu10}==time{Lu5}(tt));
        if isempty(temp),
            continue
        else
            Lu10_match_Lu5(end+1)=temp;
            Lu5_match_Lu10(end+1)=tt;
        end
        clear temp;
    end


    time_case = cell(1,length(out_files));
    %makes a time log
    time_case{Kd_0_5} = time{Eds}(Eds_match_Ed5);
    time_case{Kd_0_10} = time{Eds}(Eds_match_Ed10);
    time_case{Kd_5_10} = time{Ed5}(Ed5_match_Ed10);
    time_case{Ku_5_10} = time{Lu5}(Lu5_match_Lu10);

    %matches Kd_0_5, Ku_5_10, Kd_5_10
    for tt=1:length(time_case{Ku_5_10}),
        if tt==1, Kd_0_5_match_both=[];Ku_5_10_match_both=[];Kd_5_10_match_both=[]; end
        temp1 = find(time_case{Kd_0_5}==time_case{Ku_5_10}(tt));
        temp2 = find(time_case{Kd_5_10}==time_case{Ku_5_10}(tt));
        if (~isempty(temp1)&&~isempty(temp2)),
            Kd_0_5_match_both (end+1) = temp1;
            Kd_5_10_match_both (end+1)= temp2;
            Ku_5_10_match_both (end+1)= tt;
            time_case{Ku_0_5} (end+1) = time_case{Ku_5_10}(tt);
        end
        clear temp1;
        clear temp2;
    end

    %matches Lu5 and Kd_5_10
    for tt=1:length(time_case{Kd_5_10}),
        if tt==1, Lu5_match_Kd_5_10=[];Kd_5_10_match_Lu5=[]; end
        temp1 = find(time{Lu5}==time_case{Kd_5_10}(tt));
        if (~isempty(temp1)),
            Lu5_match_Kd_5_10 (end+1) = temp1;
            Kd_5_10_match_Lu5 (end+1)= tt;
            time_case{Lw_Kd_0_5} (end+1) = time_case{Kd_5_10}(tt);
        end
        clear temp1;
    end

    %matches Lu5 and Ku_5_10
    for tt=1:length(time_case{Ku_5_10}),
        if tt==1, Lu5_match_Ku_5_10=[];Ku_5_10_match_Lu5=[]; end
        temp1 = find(time{Lu5}==time_case{Ku_5_10}(tt));
        if (~isempty(temp1)),
            Lu5_match_Ku_5_10 (end+1) = temp1;
            Ku_5_10_match_Lu5 (end+1)= tt;
            time_case{Lw_Ku_5_10}(end+1) = time_case{Ku_5_10}(tt);
        end
        clear temp1;
    end
    
    %matches Lu5 and Ku_0_10
    for tt=1:length(time_case{Ku_0_5}),
        if tt==1, Lu5_match_Ku_0_5=[];Ku_0_5_match_Lu5=[]; end
        temp1 = find(time{Lu5}==time_case{Ku_0_5}(tt));
        if (~isempty(temp1)),
            Lu5_match_Ku_0_5 (end+1) = temp1;
            Ku_0_5_match_Lu5 (end+1)= tt;
            time_case{Lw_Ku_0_5}(end+1) = time_case{Ku_0_5}(tt);
        end
        clear temp1;
    end

    %matches Eds and Lw_Kd_0_5
    for tt=1:length(time_case{Lw_Kd_0_5}),
        if tt==1, Eds_match_Lw_Kd_0_5=[];Lw_Kd_0_5_match_Eds=[]; end
        temp1 = find(time{Eds}==time_case{Lw_Kd_0_5}(tt));
        if (~isempty(temp1)),
            Eds_match_Lw_Kd_0_5 (end+1) = temp1;
            Lw_Kd_0_5_match_Eds (end+1)= tt;
            time_case{RSR_Kd_0_5} (end+1) = time_case{Lw_Kd_0_5}(tt);
        end
        clear temp1;
    end
    %matches Eds and Lw_Ku_5_10
    for tt=1:length(time_case{Lw_Ku_5_10}),
        if tt==1, Eds_match_Lw_Ku_5_10=[];Lw_Ku_5_10_match_Eds=[]; end
        temp1 = find(time{Eds}==time_case{Lw_Ku_5_10}(tt));
        if (~isempty(temp1)),
            Eds_match_Lw_Ku_5_10 (end+1) = temp1;
            Lw_Ku_5_10_match_Eds (end+1)= tt;
            time_case{RSR_Ku_5_10} (end+1) = time_case{Lw_Ku_5_10}(tt);
        end
        clear temp1;
    end
    %matches Eds and Lw_Ku_0_5
    for tt=1:length(time_case{Lw_Ku_0_5}),
        if tt==1, Eds_match_Lw_Ku_0_5=[];Lw_Ku_0_5_match_Eds=[]; end
        temp1 = find(time{Eds}==time_case{Lw_Ku_0_5}(tt));
        if (~isempty(temp1)),
            Eds_match_Lw_Ku_0_5 (end+1) = temp1;
            Lw_Ku_0_5_match_Eds (end+1)= tt;
            time_case{RSR_Ku_0_5} (end+1) = time_case{Lw_Ku_0_5}(tt);
        end
        clear temp1;
    end

    % Create output files
    for j = 1:length(outfile_used),
        out_file_index = outfile_used(j);
   % for out_file_index = 1:length(out_files),
        channel = out_variables{out_file_index};
        disp(['Calculating ' channel '...']);
        ncfilename = [netCDF_path '/' out_files{out_file_index} suffixes{suffix_index} '.nc'];
        ncid = fopen(ncfilename);
        if (ncid~=-1)
            ncfile = netcdf(ncfilename, 'write');
        else
            %        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) || (out_file_index >= Kd_0_5_95))
                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}(Eds_match_Ed5,:).* Transmission{Eds}(Eds_match_Ed5,:)) - ...
                    log(spectra{Ed5}(Ed5_match_Eds,:)) )  / depth(Ed5) ;
                max_valid = 10;
                min_valid = -1;
                temp = find(F{Kd_0_5}<0);
                F{Kd_0_5}(temp) = NaN;
            case{Kd_0_5_95}
                F{Kd_0_5_95} =   (log(spectra{Eds}(Eds_match_Ed5,:).* 0.95) - ...
                    log(spectra{Ed5}(Ed5_match_Eds,:)) )  / depth(Ed5) ;
                max_valid = 10;
                min_valid = -1;
                time_case{Kd_0_5_95} = time_case{Kd_0_5};
            case {Kd_0_10} %changed from old formula to one based on Morel 1988
                F{Kd_0_10} = (log(spectra{Eds}(Eds_match_Ed10,:).* Transmission{Eds}(Eds_match_Ed10,:)) - log(spectra{Ed10}(Ed10_match_Eds,:)) ) / ...
                    depth(Ed10) ;
                max_valid = 10;
                min_valid = -1;
                temp = find(F{Kd_0_10}<0);
                F{Kd_0_10}(temp) = NaN;
            case {Kd_0_10_95}
                F{Kd_0_10_95} = (log(spectra{Eds}(Eds_match_Ed10,:).* 0.95) - log(spectra{Ed10}(Ed10_match_Eds,:)) ) / ...
                    depth(Ed10) ;
                max_valid = 10;
                min_valid = -1;
                time_case{Kd_0_10_95} = time_case{Kd_0_10};
            case {Kd_5_10}
                F{Kd_5_10} = (log(spectra{Ed5}(Ed5_match_Ed10,:)) - log(spectra{Ed10}(Ed10_match_Ed5,:)))  / ...
                    (depth(Ed10) - depth(Ed5)) ;
                max_valid = 10;
                min_valid = -1;
            case{chl_490_Kd_0_5}
                F{chl_490_Kd_0_5}  = ((interp1(wave_grid,F{Kd_0_5}',490))./0.121).^(1/.428);
                time_case{chl_490_Kd_0_5} = time_case{Kd_0_5};
                max_valid = 1000;
                min_valid = 0;
            case{chl_490_Kd_0_10} %only chl product used - based on Morel 1988
                F{chl_490_Kd_0_10}  = ((interp1(wave_grid,F{Kd_0_10}',490))./0.121).^(1/.428);
                time_case{chl_490_Kd_0_10} = time_case{Kd_0_10};
                max_valid = 1000;
                min_valid = 0;
            case {Ku_5_10}
                F{Ku_5_10} = (log(spectra{Lu5}(Lu5_match_Lu10,:)) - log(spectra{Lu10}(Lu10_match_Lu5,:)))  / ...
                    (depth(Lu10) - depth(Lu5)) ;
                max_valid = 10;
                min_valid = -1;
            case {Ku_0_5}
                F{Ku_0_5}  = F{Kd_0_5}(Kd_0_5_match_both,:) + F{Ku_5_10}(Ku_5_10_match_both,:) - F{Kd_5_10}(Kd_5_10_match_both,:); % was 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}(Lu5_match_Kd_5_10,:) ./ index2{Lu5}(Lu5_match_Kd_5_10,:) .* Transmission{Lu5}(Lu5_match_Kd_5_10,:) .* exp(F{Kd_5_10}(Kd_5_10_match_Lu5,:) * depth(Lu5));%check if Kd_5_10 is the correct one to use here
                max_valid = 1;
                min_valid = 0;
            case {Lw_Ku_5_10}
                F{Lw_Ku_5_10} = spectra{Lu5}(Lu5_match_Ku_5_10,:) ./ index2{Lu5}(Lu5_match_Ku_5_10,:) .* Transmission{Lu5}(Lu5_match_Ku_5_10,:) .* exp(F{Ku_5_10}(Ku_5_10_match_Lu5,:) * depth(Lu5));
                max_valid = 1;
                min_valid = 0;
            case {Lw_Ku_0_5}
                F{Lw_Ku_0_5}  = spectra{Lu5}(Lu5_match_Ku_0_5,:) ./ index2{Lu5}(Lu5_match_Ku_0_5,:) .* Transmission{Lu5}(Lu5_match_Ku_0_5,:) .* exp(F{Ku_0_5}(Ku_0_5_match_Lu5,:)  * depth(Lu5));
                max_valid = 1;
                min_valid = 0;
            case {RSR_Kd_0_5}
                F{RSR_Kd_0_5}  = F{Lw_Kd_0_5}(Lw_Kd_0_5_match_Eds,:)  ./ spectra{Eds}(Eds_match_Lw_Kd_0_5,:);
                max_valid = 10;
                min_valid = 0;
            case {RSR_Ku_5_10}
                F{RSR_Ku_5_10} = F{Lw_Ku_5_10}(Lw_Ku_5_10_match_Eds,:) ./ spectra{Eds}(Eds_match_Lw_Ku_5_10,:);
                max_valid = 10;
                min_valid = 0;
            case {RSR_Ku_0_5}
                F{RSR_Ku_0_5}  = F{Lw_Ku_0_5}(Lw_Ku_0_5_match_Eds,:)   ./ spectra{Eds}(Eds_match_Lw_Ku_0_5,:);
                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;
                time_case{bb670_Kd_0_5} = time_case{RSR_Kd_0_5};
                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;
                time_case{bb670_Ku_5_10} = time_case{RSR_Ku_5_10};
                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;
                time_case{bb670_Ku_0_5} = time_case{RSR_Ku_0_5};
                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;
                time_case{chl_490_Kd_0_5} = time_case{Lw_Kd_0_5};
                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;
                time_case{chl_490_Ku_5_10} = time_case{Lw_Ku_5_10};
                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;
                time_case{chl_490_Ku_0_5} = time_case{Lw_Ku_0_5};
                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);
                time_case{chl_oc2v4_Kd_0_5} = time_case{RSR_Kd_0_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);
                time_case{chl_oc2v4_Ku_5_10} = time_case{RSR_Ku_5_10};
                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);
                time_case{chl_oc2v4_Ku_0_5} = time_case{RSR_Ku_0_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 = find(ncfile{'time'}==last_time_stamp(out_file_index))+1;
        %         if isempty(start_index),
        start_index=1;
        %         end

        max_index = length(time_case{out_file_index});
    
        while(~isempty(nc_end_time) && start_index < max_index && time_case{out_file_index}(start_index) <= nc_end_time ),
           start_index = start_index + 1;
        end

        if isempty(nc_end_time),
            nc_end_time = 0;
        end
        
        if((max_index~=0) && (nc_end_time<time_case{out_file_index}(max_index))),
            disp(['Writing data to to file.']);
      
            ncfile{'time'}(nc_index + 1: nc_index + max_index - start_index +1) = ...
                time_case{out_file_index}(start_index : max_index);
            if(out_file_index < bb670_Kd_0_5) || (out_file_index >= Kd_0_5_95)
                ncfile{channel}(nc_index + 1: nc_index + max_index - start_index + 1,:) = ...
                    F{out_file_index}(start_index : max_index,:);
            else
                ncfile{channel}(nc_index + 1: nc_index + max_index - start_index + 1)   = ...
                    F{out_file_index}(start_index : max_index);
            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