Purewater_SpecBuilder.m

% Purewater_SpecBuilder
% August 29, 2016
% Jesse Bausell 
%
%
% This matlab script reads in an ac-s pure-water calibration file and
% out of it creates a single spectrum. The file is typically generated by
% running deionized water through ac-s flow chamber(s) while the instrument
% is recording. 
%
% This code was updated and modified from existing code written by Jennifer
% Schullien.

 %Clear all variables, close any files/figures, clear command window
clear all; close all; clc;

%% 1. First section of code allows user to interactively select a pure-water file

worD = {' ABSORPTION (a) ';' ATTENUATION (c) '};
% Create cell variable with words that can be used interchangeably in
% subsequent interactive strings.

while 1
    % First while loop prompts user to select whether he/she will 'build'
    % an absorption or an attenuation pure-water calibration spectrum
    
    % Ask user to specify an absorption or an attenuation pure-water file.
    mdquestion = input('Press "a" for absorption and "c" for attenuation:','s');
    if strcmpi(mdquestion,'a') %User selects absorption (a)
       hh = 1; clc; worD1 = 'a'; %structure element 1, clear command window, variable 'a'
       break %break the while loop
   elseif strcmpi(mdquestion,'c') %User selects attenuation (c)
       hh = 2; clc; worD1 = 'c'; %structure element 2, clear command window, variable 'c'
       break %break the while loop
   else %User doesn't follow instructions and selects something nonsensical      
       disp('Invalid entry. Try again.') % Tells user to try again. Re-runs while loop       
    end    
end

while 1
    % Second while loop prompts the user to select an
    % absorption/attenuation file (predicated on his/her choice of which
    % type of spectrum to build.
    disp(['choose a pure water ' worD{hh} ' file']); %Instructs user to select a pure-water file
    [wcal_a_infile, a_DIR] = uigetfile('.dat'); %User selects a file.
    disp(wcal_a_infile); %Display the file that the user selected
    % Asks the user whether or not he/she selected the correct file
    y = input(['Is ' wcal_a_infile ' the correct' worD{hh} ' file? (y/n): '],'s'); 
    % Allows the user to accept the file or select another one if
    % unsatisfied
    if isequal(strcmpi(y,'y'),1)
        break
    end
end

%% 2. Section of code will extract the data from the selected pure-water file

fid = fopen([a_DIR wcal_a_infile]); % Open pure-water file
temp = textscan(fid, '%s','Delimiter',','); % Scan the entire file into matlab
fclose(fid); clear ans fid % close the file and clear variables 'ans' and 'fid'
stra = temp{:}; % organize file lines into a cell array (line by line)

for ii = 1:length(stra)
    % This for-loop goes through pure-water data (extracted from the file)
    % line by line. It separates the relevant portion of the file (a/c time
    % series) from the irrelevant parts.
    linE = stra{ii}; % Variable grabs one line at a time (in sequential order)
    if regexpi(linE,'output wavelengths') 
        % Locates the line with number of ac-s channel wavelengths
        inD = regexpi(linE,';'); % Index the semicolon
        wvl_num = str2num(linE(1:inD-1)); % Locates the number of channels and makes it an integer
    end

    if regexpi(linE,'acquisition binsize') 
        % Locates the line that says 'aqcuisition binsize'
        break % breaks the for-loop
    end
end

header_lines = stra(1:ii-4); % Puts the entire header onto one variable
wcal_a_temp = stra(ii+2:end); % Puts all of the a/c data into one variable
wcal_a_waveLENGTHS = stra{ii+1}; %added to Jen's code: this takes wave lengths line from the file
a_waveIND = regexpi(wcal_a_waveLENGTHS,worD1); %this finds the locations of 'a' or 'c' in wavelength line
% These bottom two lines convert a or c data from string to matrix form
a = cellfun(@str2num,wcal_a_temp,'UniformOutput',0); 
wcal_a = cell2mat(a); 

% The if statement below indexes absorption (a) or attenuation (c) channels 
% based on user's previous selection (see section 1).
if isequal(hh,1) % If user selected a for absorpption
    cLR = {'b';'r'}; % Color code the plots in the next section
    a_channels = wcal_a(:,length(a_waveIND)+2:(2*length(a_waveIND))+2); %index a values
elseif isequal(hh,2) % If user selected c for attenuation
    cLR = {'r';'b'}; % Color code the plots in the next section
    a_channels = wcal_a(:,2:length(a_waveIND)+1); % index c values 
end

%% 3. The final section of code allows user to interactively build a pure-water spectrum

% Empty arrays for averaging indices. The two variables below will be used
% to index the ranges at which a/c will be averaged at each channel.
col_starT = nan(1,length(a_waveIND)); % empty nan array for starting averaging index
col_stoP = nan(1,length(a_waveIND)); % empty nan array for ending averaging index
% Create empty structure to place final product:
ACS(hh).lambda = nan(1,length(a_waveIND)); % empty nan array for wavelength
ACS(hh).spectra = nan(1,length(a_waveIND)); % empty nan array for a or c
ACS(hh).T_cal = nan(1,length(a_waveIND)); % empty nan array for temperature
% Indexing variables:
ii = 1; % index variable that goes from wavelength to wavelength
keY =0; % index variable that facilitates different parts of the while loop (see below)

while ii <= length(a_waveIND)
    % This while loop builds the average a/c spectrum. The pure-water data
    % measured by ac-s is a time series. Pure-water runs through the flow
    % chamber(s) and the ac-s consistently measures pure-water absorption
    % and attenuation. The loop allows the user to visualize the time
    % series one channel at a time, and then select a portion that is
    % reliable (e.g. values are realistic and reasonable, not artificially
    % high or low). This loop will then average the a/c for the
    % aforementioned channel, allowing the user to continue to the next
    % one.
    close; clc; % Close all plots, clear the command window.
    % Add channel wavelength to the matlab structure (final product)
    ACS(hh).lambda(ii) = str2num(wcal_a_waveLENGTHS(a_waveIND(ii)+1:a_waveIND(ii)+5));
    % Plot the a/c timeline. Print the wavelength (as lambda) on the figure
    h = plot(a_channels(:,ii)); text(0.9*length(a_channels(:,ii)),0.9*max(a_channels(:,ii)),['\lambda = ' num2str(ACS(hh).lambda(ii))],'FontSize',20);
    hold on % hold the figure to allow for future user input.

    % Following 'if' statement selects minimum index for each
    % channel (wavelength) that will be used in creating an average
    % spectrum (one channel/wavelength at a time)
    if isequal(keY,0) % No minumum or maximum index value is selected
        starT = input('Enter minimum limit: '); clc; %Prompts user to select minimum index
        % Plots minimum index onto the time series graph so user can
        % visualize it (line below)
        plot([starT starT],[min(a_channels(:,ii)) max(a_channels(:,ii))],'r','LineWidth',2);
        % Asks user to confirm selection of miminum index (line below)
        quE = input(['Is ' num2str(starT) ' an acceptable minimum? (y/n) '],'s');
        
        if strcmpi(quE,'y') % User confirms or rejects selection using 'y' key
            % keY variable changes to shift into next section of while loop
            keY = keY + 1; 
        else 
            close; % 
        end
    end

    % Following 'if' statement selects maximum index for each
    % channel (wavelength) that will be used in creating an average
    % spectrum (one channel/wavelength at a time)
    if isequal(keY,1) % No maximum index is selected
        stoP = input('Enter maximum limit: '); % Prompts user to select maximum index
        % Plots maximum index onto the time series graph so user can
        % visualize it (line below)
        plot([stoP stoP],[min(a_channels(:,ii)) max(a_channels(:,ii))],'r','LineWidth',2);
        % Asks user to confirm selection of miminum index (line below)
        quE = input(['Is ' num2str(stoP) ' an acceptable maximum? (y/n) '],'s');

        if strcmpi(quE,'y') % User confirms or rejects selection using 'y' key
            % keY variable changes to shift into next section of while loop
             keY = keY + 1;
        end
    end

    if isequal(keY,2) % Minimum and maximum indices selected and approved by user
        % Plots minimum and maximum indices as vertical lines on time
        % series. Gives user a visualization of what he/she selected. keY
        % is kept equal to 2 in the while loop continuously unless user
        % says otherwise. This is to save user time of redundantly
        % selecting the same min/max indices over and over again.
        plot([starT starT],[min(a_channels(:,ii)) max(a_channels(:,ii))],'r','LineWidth',2);
        plot([stoP stoP],[min(a_channels(:,ii)) max(a_channels(:,ii))],'r','LineWidth',2);
        % Asks user to confirm his/her min/max index selection (line below)
        quE = input(['Are ' num2str(starT) ' and ' num2str(stoP) ' acceptable limits? (y/n) '],'s');

        if strcmpi(quE,'y') % User confirms selection
            ACS(hh).spectra(ii) = nanmean(a_channels(starT:stoP,ii)); %a/c is averaged for given channel
            ACS(hh).T_cal(ii) = nanmean(wcal_a(starT:stoP,183)); %temperature is averaged for given channel
            col_starT(ii) = starT; % minimum index is put into an array
            col_stoP(ii) = stoP; % maximum index is put into an array
            hold off; % Plot hold is removed allowing for a creation of a fresh new plot
            ii = ii + 1; % channel index increased by 1, allowing while loop to move to next channel
        elseif strcmpi(quE,'n') % User rejects selection
            keY = 0; % keY brought from 2 to zero. Index selection process is repeated
        end
    end
end

% Newly created average a/c pure-water spectrum is put into a .mat file
fileIND = regexp(wcal_a_infile,'.dat') - 1; % Isolate name of pure-water calibration file
save([a_DIR wcal_a_infile(1:fileIND) '_pwavg'],'ACS'); % Save pure-water calibration spectrum as .mat file

% Choose one randomly-selected time series of ac-s channel and save as a
% time series plot.
rr = round(1 + (length(ACS(hh).spectra)-1).*rand(1)); %Select a wavelength at random
% Plot a or c time series associated with selected channel (a/c vs. index)
h = plot(a_channels(:,rr),cLR{1}); text(0.9*length(a_channels(:,rr)),0.8*max(a_channels(:,rr)),['\lambda = ' num2str(ACS(hh).lambda(rr))],'FontSize',20);
hold on; plot([col_starT(rr) col_starT(rr)],[min(a_channels(:,rr)) max(a_channels(:,rr))],cLR{2},'LineWidth',2);
plot([col_stoP(rr) col_stoP(rr)],[min(a_channels(:,rr)) max(a_channels(:,rr))],cLR{2},'LineWidth',2);
ylabel(worD{hh});
saveas(gcf,[a_DIR wcal_a_infile(1:fileIND) '_gateSAMPLE'],'jpeg'); % Save plot as jpeg
close; %close aforementioned plot

% Last plot is a subplot. The top panel is average  the a/c spectrum. 
% Bottom plot is the average water temperature of each channel. This may
% vary slightly if different indices are used for different channels.
subplot(2,1,1) % Top panel
plot(ACS(hh).lambda,ACS(hh).spectra,cLR{2},'LineWidth',2) % plot a/c spectrum
xlabel('\lambda  (nm)'); ylabel([worD{hh}(1:end-4) '(m ^-^1)']); % label the plot

subplot(2,1,2) % Bottom panel
plot(ACS(hh).lambda,ACS(hh).T_cal,'k') % plot the temperature panel
xlabel('\lambda  (nm)'); ylabel('Temperature (^oC)'); % label the plot
saveas(gcf,[a_DIR wcal_a_infile(1:fileIND) '_spectrum'],'jpeg'); % save subplots as jpeg
close; % close the subplot