function ocr3000(configfile, startDateTime)
% function ocr3000(configfile, startDateTime)
%
% Reads ocr3000 data from SSDS, processes the data, and 
% saves the processed results in several netCDF files.
% Inputs: startDateTime (GMT)
% If the startDateTime string is specified, then data is 
% read from then until to the present (for initialization).
% Otherwise (recommended), the existing netcdf files are queried.
%
% Created 2004-06-11 Michael Godin, godin@mbari.org
%-------------------------------------------------------------------------
% Requires:
%   readProperties.m
%   getProperty.m
%   string2escape.m
%   urlread.m
%   stringtokenizer.m
%   hex2uint8.m
%-------------------------------------------------------------------------

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

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
    light_ids{radiometer_count} = getProperty(prop, [radiometer{radiometer_count} '.light_id']);
    dark_ids{radiometer_count}  = getProperty(prop, [radiometer{radiometer_count} '.dark_id']);
end


% Determine when to start reading in raw data
GMT_offset = str2num(getProperty(prop, 'GMT_offset'));
netCDF_path = MBARI_path(getProperty(prop, ['netCDF_path']));
if(nargin < 2),
    for radiometer_index = 1:radiometer_count,
        ncfilename = [netCDF_path '/' radiometer{radiometer_index} '_mean.nc'];
        ncfile = netcdf(ncfilename, 'nowrite');
        if(size(ncfile,2)==0),
            lastTime(radiometer_index) = 0;
        else
            lastTime(radiometer_index) = ncfile{'time'}(end)/86400 + datenum(1970,1,1); % was + GMT_offset/24;
        end
        close(ncfile);
    end
    startFrom = min(lastTime)+1/86400;
else
    startFrom = datenum(startDateTime);
end

% Invalid start time not handled for now
if startFrom == 0
    return
end

% Build up data access URL
URL_base                  = getProperty(prop, 'URL_base');
URL_identification        = strrep(getProperty(prop, 'URL_identification'),'%3D','=');
URL_options               = strrep(getProperty(prop, 'URL_options'),'%3D','=');
URL_start_label           = getProperty(prop, 'URL_start_label');
URL_end_label             = getProperty(prop, 'URL_end_label');
URL_datetime_form         = getProperty(prop, 'URL_datetime_form');
URL_datetime_convert_from = getProperty(prop, 'URL_datetime_convert_from');
URL_datetime_convert_to   = getProperty(prop, 'URL_datetime_convert_to');

% Read in calibraton files
cal_path = MBARI_path(getProperty(prop, 'calibration_path'));
cal_data = read_ocr3000_cal(prop, cal_path, light_ids, dark_ids);

hoursTillNow = ceil((now- GMT_offset/24-startFrom)*24);
hoursPerDownload = 24;
while hoursTillNow >= 1,
    disp (['hoursTillNow = ' num2str(hoursTillNow)]);
    start_string = strrep(datestr(startFrom,URL_datetime_form),URL_datetime_convert_from,URL_datetime_convert_to);
    end_time = startFrom + hoursPerDownload/24.0-1/86400;
    end_string = strrep(datestr(end_time,URL_datetime_form),URL_datetime_convert_from,URL_datetime_convert_to);
    
    URL = [URL_base '?' URL_identification '&' URL_options '&' ...
            URL_start_label '=' start_string '&' URL_end_label '=' end_string]  % was int2str(hoursTillNow)
	
    startFrom = startFrom + hoursPerDownload/24.0;

    % Read in raw data for each radiometer
	raw_data = read_ocr3000_raw(URL, light_ids, dark_ids);
	
    if isempty(raw_data), 
        hoursTillNow = hoursTillNow - hoursPerDownload;
    	continue; 
    end
    
    % Calibrate the raw data into engineering units
	eng_data = ocr3000_calc_eng(raw_data, cal_data);

	% Subtract dark data, and reduce data samples to one per time stamp
	[c_time mean_data median_data std_data n_samples] = ocr3000_dark_sub(eng_data);
	
	% Grid data to a uniform wavelength set
	% 2005_01_05 Changed wave_grid from an output to an input
	wavelength.min = str2num(getProperty(prop, 'wavelength.min'));
	wavelength.max = str2num(getProperty(prop, 'wavelength.max'));
	wavelength.inc = str2num(getProperty(prop, 'wavelength.inc'));
	wave_grid = [wavelength.min:wavelength.inc:wavelength.max];
        [mean_grid median_grid std_grid] = ...
            ocr3000_grid(mean_data, median_data, std_data, cal_data, wave_grid);
        %[mean_grid median_grid std_grid wave_grid] = ...
        %    ocr3000_grid(mean_data, median_data, std_data, cal_data);

        
	% Modify engineering units
	[mean_mod median_mod std_mod units] = ...
        ocr3000_mod_units(mean_grid, median_grid, std_grid, cal_data, 1/100, 'uW/cm^2', 'W/m^2');
	
	% Store data in netcdf files
    disp('Storing data in netCDF files...');
	ocr3000_netCDF_store(c_time, prop, radiometer, light_ids, netCDF_path, wave_grid, units, n_samples, mean_mod, '_mean');
	ocr3000_netCDF_store(c_time, prop, radiometer, light_ids, netCDF_path, wave_grid, units, n_samples, mean_mod, '_median');
	ocr3000_netCDF_store(c_time, prop, radiometer, light_ids, netCDF_path, wave_grid, units, n_samples, std_mod, '_std');    

    hoursTillNow = hoursTillNow - hoursPerDownload;

end

disp('ocr3000 done.');
quit

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cal_data = read_ocr3000_cal(prop, cal_path, light_ids, dark_ids);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read in OCR3000 calibraton files
    
disp('Reading calibration data...');
all_ids = [light_ids, dark_ids];
id_size = size(light_ids,2);

for ld_index = 1:2,
    for index = 1:id_size,
        cal_filename = [cal_path '/' getProperty(prop, [all_ids{(ld_index-1)*id_size + index} '.calibration'])];
        [cal_data_input gotdata]= urlread(['file://' cal_filename]);
        if(gotdata)
            linefeeds = findstr(char(10),cal_data_input);
        else
            linefeeds = 0;
        end
        wl_index = 1;
        dark_sum = 0;
        for line = 1:size(linefeeds,2)-1,
            data_line = cal_data_input(linefeeds(line)+1:linefeeds(line+1)-1);
            %parse out the data lines into cal file records
            spaces = findstr(' ',data_line);
            %test to make sure this is a valid record 
            if size(spaces,2) == 6 && data_line(1) ~= '#' && isempty(strmatch(data_line(1:3),'AUX')) ...
                    && data_line(spaces(1)+1) >= '1' && data_line(spaces(1)+1) <= '9',
                % Read in the wavelength and units
		        cal_data{ld_index, index}.wavelength(wl_index) = str2num(data_line(spaces(1)+1:spaces(2)-1));
                if(wl_index == 1),
                    cal_data{ld_index, index}.units = data_line(spaces(2)+2:spaces(3)-2);
                end
                % Now read in the calibration coefficients
                line = line + 1;
                data_line = cal_data_input(linefeeds(line)+1:linefeeds(line+1)-1);
                tabs = [findstr(char(9),data_line) size(data_line,2)];
                if(size(tabs,2) == 4),
                    if data_line(tabs(end)) == char(13), tabs(end) = tabs(end) -1; end
                    cal_data{ld_index, index}.a0(wl_index) = str2num(data_line(1:tabs(1)-1));
                    cal_data{ld_index, index}.a1(wl_index) = str2num(data_line(tabs(1)+1:tabs(2)-1));
                    cal_data{ld_index, index}.Im(wl_index) = str2num(data_line(tabs(2)+1:tabs(3)-1));
                    cal_data{ld_index, index}.CInt(wl_index) = str2num(data_line(tabs(3)+1:tabs(4)-1));
                    dark_sum = dark_sum + cal_data{ld_index, index}.a0(wl_index);
                end
                cal_data{ld_index, index}.dark_average = dark_sum/(wl_index);
                wl_index = wl_index +1;
            end
        end
    end
end
disp('Done.');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [raw_data] = read_ocr3000_raw(URL, light_ids, dark_ids)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reads in an OCR3000 raw data file, 
% prefixed with a SSDS timestamp and recond#,
% with the data in hex
% Organizes data based on values in 
% cell arrays light_ids and dark_ids
%
% Created 2004-06-11 Michael Godin, godin@mbari.org
%

disp('Reading data from internet...')
raw_data = [];
ssds_data = urlread(URL);
%ssds_data = urlread('file:///C:/Monterey_Moorings/CIMT/200406/SSDSData.htm');
disp(['Read ' num2str(length(ssds_data)) ' bytes.']);
fid = fopen('urlread.txt','wb');
fwrite(fid,ssds_data,'uint8');
fclose(fid);

%parse out the data file, one line at a time
linefeeds = findstr(char(10),ssds_data);
record_num = zeros(2,size(light_ids,2));
for index = 1:length(linefeeds)-1,
    data_line = ssds_data(linefeeds(index)+1:linefeeds(index+1)-1);
    
    %parse out the data lines into timestamp, record#, and datasets
    comma_spaces = findstr(', ',data_line);
    %test to make sure this is a valid record 
    if length(comma_spaces) == 3 && data_line(1) >= '1' && data_line(1) <= '9' ...
            && data_line(comma_spaces(2)-3) == 'G', 
 %       some_data = zeros(2,length(light_ids));
        time_stamp = str2num (data_line(comma_spaces(2)+2:comma_spaces(3)-4));  % This omits the last 3 sig figs of the time stamp
        dataset = data_line(comma_spaces(3)+2:end);
        %now we parse out the data line itself into measuremnts
        packets_SATH = findstr('53415448',dataset); % look for encoded 'SATH'
        packets_SATP = findstr('53415450',dataset); % look for encoded 'SATP'
        packets = sort([packets_SATH packets_SATP size(dataset,2)+1]);
        for index_data = 1:size(packets,2)-1,
            hexdata = dataset(packets(index_data):packets(index_data+1)-1);
            if(size(hexdata,2) ~= 622 && size(hexdata,2) ~= 623)
                disp('Bad Data !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!');
                continue
            end
            id = char(hex2uint8(hexdata(1:20)));
            int_time = hex2dec(hexdata(21:24))/1000;
            sample_delay = hex2dec(hexdata(25:28));
            if(sample_delay) > 32767, 
                sample_delay = sample_delay - 65536 ;
            end
                        
            % Pixel Read
            spectra = hex2uint16(hexdata(29:576));

            %dark_samples = hex2dec(hexdata(577:578);
            dark_average = hex2dec(hexdata(579:582));
            %aux1 = hex2uint8(hexdata(583:590));
            %aux2 = hex2uint8(hexdata(591:594));
            frame_counter =  hex2dec(hexdata(595:596));
            timer = str2num(char(hex2uint8(hexdata(597:616))));
            %checksum = hex2dec(hexdata(617:618);
            %terminator = char(hex2uint8(hexdata(619:622)))

            light_match = strmatch(id, light_ids);
            dark_match = strmatch(id, dark_ids);
            
            ld_index = 0;
            if (light_match),
                ld_index = 1;
            elseif (dark_match),
                ld_index = 2;
            end
            if ld_index > 0,
                sensor = max([light_match dark_match]);
                record_num(ld_index, sensor) = record_num(ld_index, sensor) + 1;
                raw_data{ld_index, sensor}(record_num(ld_index, sensor)).time_stamp = time_stamp;
                raw_data{ld_index, sensor}(record_num(ld_index, sensor)).int_time = int_time;
                raw_data{ld_index, sensor}(record_num(ld_index, sensor)).sample_delay = sample_delay;
                raw_data{ld_index, sensor}(record_num(ld_index, sensor)).spectra = spectra;
                raw_data{ld_index, sensor}(record_num(ld_index, sensor)).dark_average = dark_average;
                raw_data{ld_index, sensor}(record_num(ld_index, sensor)).frame_counter = frame_counter;
                raw_data{ld_index, sensor}(record_num(ld_index, sensor)).timer = timer;
%                some_data(ld_index, sensor)=1;
            end               
        end
%        if(sum(sum(some_data)) > 0),
%            for sensor = 1 : length(light_ids),
%                if(some_data(1,sensor) == 0),
%                    record_num(1, sensor) = record_num(1, sensor) + 1;
%                    raw_data{1, sensor}(record_num(1, sensor)).time_stamp = time_stamp;
%                    raw_data{1, sensor}(record_num(1, sensor)).int_time = NaN;
%                    raw_data{1, sensor}(record_num(1, sensor)).sample_delay = NaN;
%                    raw_data{1, sensor}(record_num(1, sensor)).spectra = NaN;
%                    raw_data{1, sensor}(record_num(1, sensor)).dark_average = NaN;
%                    raw_data{1, sensor}(record_num(1, sensor)).frame_counter = NaN;
%                    raw_data{1, sensor}(record_num(1, sensor)).timer = NaN;
%                end
%            end
%        end
        disp(['Processed raw data from ' datestr(time_stamp/86400 + datenum(1970,1,1))]);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function eng_data = ocr3000_calc_eng(raw_data, cal_data);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calibrate the raw data into engineering units
%
disp('Calculating Engineering Units');
for i = 1: size(raw_data,1),
    for j = 1 : size(raw_data,2),
        for k = 1: size(raw_data{i,j},2),
            eng_data{i,j}(k).time_stamp = raw_data{i,j}(k).time_stamp;
            eng_data{i,j}(k).timer = raw_data{i,j}(k).timer;
            eng_data{i,j}(k).int_time = raw_data{i,j}(k).int_time;
            eng_data{i,j}(k).spectra    = cal_data{i,j}.Im .* cal_data{i,j}.a1 ...
                .* ( raw_data{i,j}(k).spectra - cal_data{i,j}.a0 ...
                    + cal_data{i,j}.dark_average - raw_data{i,j}(k).dark_average) ...
                .* ( cal_data{i,j}.CInt ./ raw_data{i,j}(k).int_time);
        end
    end
end    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [c_time,mean_data,median_data,std_data,n_samples] = ocr3000_dark_sub(eng_data);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subtract dark data, and reduce data samples to one per time stamp
% Keep track of the int_time for each instrument...  
%   as long as it's lower than the previous, discard the previous data.

disp('Subtracting Dark Data');
for j = 1 : size(eng_data,2),   % this steps from Es through Lu10
    time_stamp = 0;
    record_num = 1;
    dark_index = size(eng_data{2,j},2);
    for k = size(eng_data{1,j},2) : -1 : 1,
        if eng_data{1,j}(k).time_stamp ~= time_stamp,
            time_stamp = eng_data{1,j}(k).time_stamp;
            int_time = eng_data{1,j}(k).int_time;
            max_index = k;
            clear array_data;
            n_samples(j,record_num) = 0;
            done = 0;
        end
        if done, continue; end
        
        % Find the next dark data to subtract...
        for dark_index = dark_index : -1  : 1,
            if eng_data{2,j}(dark_index).time_stamp <= time_stamp,
                break;
            end
        end
        
        % Subtract the dark data 
        if(dark_index > 0) 
            eng_data{1,j}(k).spectra = eng_data{1,j}(k).spectra ...
                                     - eng_data{2,j}(dark_index).spectra;
        end
        
        % Add to array data
        n_samples(j,record_num) = n_samples(j,record_num) + 1;
        array_data(n_samples(j,record_num),:) = eng_data{1,j}(k).spectra;
        
        % Check if this is the last record of the set
        if k==1 || eng_data{1,j}(k-1).time_stamp ~= time_stamp ...
                    ||  eng_data{1,j}(k-1).int_time ~= int_time,
            done = 1;
            c_time(j,record_num) = time_stamp;
            if(n_samples(j,record_num) > 1)
                mean_data(j,:,record_num) = mean(array_data);
                median_data(j,:,record_num) = median(array_data);
                std_data(j,:,record_num) = std(array_data);
            else
                mean_data(j,:,record_num) = array_data;
                median_data(j,:,record_num) = array_data;
                std_data(j,:,record_num) = array_data * 0;
            end
            record_num = record_num + 1;
        end
    end
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [mean_grid,median_grid,std_grid] = ...
    ocr3000_grid(mean_data, median_data, std_data, cal_data, wave_grid);
% 2005_01_05 Changed wave_grid from an output to an input
%function [mean_grid,median_grid,std_grid,wave_grid] = ...
%    ocr3000_grid(mean_data, median_data, std_data, cal_data);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grid data
%
wave_grid_size = 10/3;

% Find limits
%for i = 1 : size(cal_data,2),
%    min_array(i) = cal_data{1,i}.wavelength(1);
%    max_array(i) = cal_data{1,i}.wavelength(end);
%end
%min_wave = ceil(max(min_array)/wave_grid_size)*wave_grid_size;
%max_wave = floor(min(max_array)/wave_grid_size)*wave_grid_size;
%wave_grid = [min_wave : wave_grid_size : max_wave];

for i = 1 : size(mean_data,1),
    for k = 1: size(mean_data,3),
        mean_grid(i,:,k) = interp1( ...
            cal_data{1,i}.wavelength, mean_data(i,:,k), wave_grid(:),'cubic');
        median_grid(i,:,k) = interp1( ...
            cal_data{1,i}.wavelength, median_data(i,:,k), wave_grid(:),'cubic');
        std_grid(i,:,k) = interp1( ...
            cal_data{1,i}.wavelength, std_data(i,:,k), wave_grid(:),'cubic');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [mean_mod,median_mod,std_mod,units] = ...
    ocr3000_mod_units(mean_data, median_data, std_data, cal_data, factor, old_unit, new_unit);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Modify engineering units

mean_mod = mean_data * factor;
median_mod = median_data * factor;
std_mod = std_data * factor;

% Get new unit names
for i = 1 : size(cal_data,2), 
    units{i} = strrep(cal_data{1,i}.units, old_unit, new_unit);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ocr3000_netCDF_store(c_time, prop, radiometer, light_ids, netCDF_path, ...
            wave_grid, units, n_samples, dataset, suffix)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Store data in a netCDF file

mooring.name = getProperty(prop, 'mooring');

for radiometer_index = 1:length(radiometer),
    channel = radiometer{radiometer_index};
    channel_name = getProperty(prop, [channel '.name']);
    ncfilename = [netCDF_path '/' channel suffix '.nc'];
    ncfile = netcdf(ncfilename, 'write');
    if(size(ncfile,2)==0 || isempty(ncfile{'time'}(end))),
        lastTime(radiometer_index) = 0;
    else
        lastTime(radiometer_index) = ncfile{'time'}(end);
    end
    start_index = size(c_time,2);
     
    while start_index >= 1 ...
            && lastTime(radiometer_index) >= c_time(radiometer_index,start_index),
        start_index = start_index - 1;
    end
    if start_index >= 1,
        if(lastTime(radiometer_index) == 0),
                  
            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'}(:)          = str2num(getProperty(prop, [mooring.name '.latitude']));
            ncfile.latitude                = str2num(getProperty(prop, [mooring.name '.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'}(:)         = str2num(getProperty(prop, [mooring.name '.longitude']));
            ncfile.latitude                = str2num(getProperty(prop, [mooring.name '.longitude']));

            ncfile('depth')                = 1;
            ncfile{'depth'}                = ncdouble('depth'); 
            ncfile{'depth'}.long_name      = getProperty(prop, 'depth.long_name'); 
            ncfile{'depth'}.units          = getProperty(prop, 'depth.units');
            ncfile{'depth'}.AXIS           = getProperty(prop, 'depth.AXIS');
            ncfile{'depth'}.positive       = getProperty(prop, 'depth.positive');
            ncfile{'depth'}(:)             = str2num(getProperty(prop, [channel '.depth']));
            ncfile.depth                   = str2num(getProperty(prop, [channel '.depth']));

            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      
            ncfile{channel_name}                = ncfloat('time','depth','latitude','longitude','wavelength');
            %ncfile{channel_name}                = ncfloat('time','wavelength');
            ncfile{channel_name}.long_name      = getProperty(prop, [channel '.long_name']);
            ncfile{channel_name}.units          = units{radiometer_index};
            ncfile{channel_name}.immersed       = getProperty(prop, [channel '.immersed']);
            ncfile{channel_name}.light_id       = getProperty(prop, [channel '.light_id']);
            ncfile{channel_name}.dark_id        = getProperty(prop, [channel '.dark_id']);
            ncfile{channel_name}.FillValue_     = -9999;
            ncfile{channel_name}.missing_value  = -9999;

            ncfile{'sample_count'}              = ncshort('time');
            ncfile{'sample_count'}.long_name    = 'samples per measurement';
        end            
        
        [Year Month Day]                        = datevec(now);
        dateS                                   = sprintf('%02i/%02i/%04i', Month, Day, Year);
        ncfile.lastmodified                     = dateS;
        nc_index = length(ncfile{'time'}(:));
        %for time_index = start_index : -1 : 1,
        %    ncfile{'time'}(nc_index)                 = c_time(radiometer_index,time_index);
        %    ncfile{channel_name}(nc_index,:)         = squeeze(dataset(radiometer_index,:,time_index));
        %    ncfile{'sample_count'}(nc_index)         = n_samples(radiometer_index,time_index);
        %    nc_index = nc_index + 1;
        %end
        ncfile{'time'}(nc_index + 1 : nc_index + start_index)= ...
            flipdim(c_time(radiometer_index, 1 : start_index),2);
        ncfile{channel_name}(nc_index + 1 : nc_index + start_index,:) = ... 
            flipdim(permute(dataset(radiometer_index, : ,1: start_index),[1 3 2]),2);
        ncfile{'sample_count'}(nc_index + 1 : nc_index + start_index) = ...
            flipdim(n_samples(radiometer_index, 1 : start_index),2);
    end
    close(ncfile);
end
return