Merge pull request #26 from tsipkens/developer

Minor, mostly cosmetic updates, including io and Grid improvements

+invert/gen_Gd.m

@@ -0,0 +1,99 @@
+
+% GEN_GD  Generate the Gd matrix from a series of name-value pairs.
+% Author: Timothy Sipkens, 2020-07-17
+% 
+% Inputs: 
+%   OPTION 1: One can enter a series of name-value pairs.
+%   At least, three of the following arguments must be supplied:
+%       'l1' or 's1' - Correlation length for variable 1 (e.g. mass).
+%       'l2' or 's2' - Correlation length for variable 2 (e.g. mobility diameter).
+%       'R12' - Correlation between variables.
+%       'Dm' - Mass-mobility exponent or slope between variables
+%         e.g., varargin = {'l1',0.2,'Dm',3,'R12',0.99}
+%   
+%   OPTION 2: Alternatively, one can supply Gd, which is used to fill out the 
+%   sg structure, containing the above properties. Remaining name-value
+%   pairs are used to output a new Gd by modifying l2.
+%         e.g., varargin = {Gd,'R12',0.99}
+%=========================================================================%
+
+function [Gd, sg] = gen_Gd(varargin)
+
+sg = struct('l1',[],'l2',[],...
+    'R12',[],'Dm',[]); % structure of Gd properties
+
+
+%== OPTION 2: Gd was supplied, return correlation lengths, corr., etc. ===%
+%   Update parameters if additional arguments are supplied.
+if isa(varargin{1}, 'double')
+    Gd = varargin{1};
+
+    % evaluate parameter set implied by Gd, before modifiction
+    sg.l1 = sqrt(Gd(1,1));
+    sg.l2 = sqrt(Gd(2,2));
+    sg.R12 = Gd(1,2) / (sg.l1 * sg.l2);
+    sg.Dm = Gd(1,2) / Gd(2,2);
+
+    % loop through name-values pairs
+    for ii=2:2:length(varargin)
+        if strcmp(varargin{ii},'s1'); varargin{ii} = 'l1'; end % allow for 's1' in the place of 'l1'
+        if or(strcmp(varargin{ii},'s2'),strcmp(varargin{ii},'12')) % cannot use l2
+            warning
+            continue;
+        end
+        sg.(varargin{ii}) = varargin{ii+1}; % copy Gd properties to structure
+    end
+    sg.l2 = sg.R12 / sg.Dm * sg.l1;
+
+
+%== OPTION 1: Otherwise, name-value pairs supplied =======================%
+else
+    % loop through name-values pairs
+    for ii=1:2:length(varargin) 
+        % allow for 's1' or 's2' in the place of 'l1' and 'l2'
+        if strcmp(varargin{ii},'s1'); varargin{ii} = 'l1'; end
+        if strcmp(varargin{ii},'s2'); varargin{ii} = 'l2'; end
+        sg.(varargin{ii}) = varargin{ii+1}; % copy Gd properties to structure
+    end
+
+    % update correlation lengths, if necessary
+    if isempty(sg.l1) % use Dm, R12, and l2 if l1 is missing
+        sg.l1 = sg.Dm / sg.R12 * sg.l2;
+    elseif isempty(sg.l2) % use Dm, R12, and l1 if l2 is missing
+        sg.l2 = sg.R12 / sg.Dm * sg.l1;
+    end
+end
+
+
+% initiate Gd using correlation lengths
+% re-evaluates to Gd if only Gd was supplied as input to method
+Gd = [sg.l1^2, 0; ...
+    0, sg.l2^2];
+
+
+% incorporate correlation lengths
+if ~isempty(sg.R12) % OPTION 1: By default, use supplied value of R12
+    Gd(1,2) = sg.R12 * sg.l1 * sg.l2;
+
+else % OPTION 1: Alternatively, use mass-mobility exponent / slope
+    Gd(1,2) = sg.Dm * sg.l2 .^ 2;
+
+end
+Gd(2,1) = Gd(1,2);
+
+
+%-- Update sg structure --------------------------------------------------%
+sg.l1 = sqrt(Gd(1,1));
+sg.l2 = sqrt(Gd(2,2));
+sg.R12 = Gd(1,2) / (sg.l1 * sg.l2);
+sg.Dm = Gd(1,2) / Gd(2,2);
+%-------------------------------------------------------------------------%
+
+
+% return error if correlation is unphysical
+if sg.R12>=1
+    error('Correlation implied by Gd matrix exceeds unity.');
+end
+
+end
+

+invert/tikhonov.m

@@ -1,7 +1,7 @@
 
 % TIKHONOV  Performs inversion using various order Tikhonov regularization in 2D.
 % Author:   Timothy Sipkens, 2018-11-21
-%-------------------------------------------------------------------------%
+% 
 % Inputs:
 %   A           Model matrix
 %   b           Data

+io/import_2.m

@@ -0,0 +1,86 @@
+
+% IMPORT_2  Read tandem CPMA-SMPS data using SMPS data export and a CPMA log file.
+% Author: Timothy Sipkens, 2019-12-17
+%
+% Note: fn_cpma can be either a single file name or a cell array of file
+% names. The times are combined / stacked vertically before a search for
+% the times corresponding to the DMA time stamps is performed. CPMA data
+% files should be given in chronological order.
+%=========================================================================%
+
+function [data,d_star,sp,prop_dma,prop_pma] = import_2(fn_smps,fn_cpma,prop_pma)
+
+
+%== Read SMPS/CPC file =========================================%
+[data,d_star,prop_dma,time_dma] = io.import_dma(fn_smps);
+%===============================================================%
+
+
+
+%== Read CPMA file =============================================%
+addpath('tfer_pma');
+
+if ~iscell(fn_cpma); fn_cpma = {fn_cpma}; end
+
+opts = detectImportOptions(fn_cpma{1}, ...
+    'FileType', 'text'); % default options
+for ff=1:length(fn_cpma)
+    if ~exist('time_cpma','var'); idx0 = 1;
+    else; idx0 = [idx0, length(time_cpma)]; end
+
+    t0 = readtable(fn_cpma{ff}, opts); % read CPMA file
+    time_cpma0 = datetime(t0.Date_Time); % get times to compare against DMA data
+
+    if ~exist('t1','var'); t1 = t0; time_cpma = time_cpma0;
+    else; t1 = [t1;t0]; time_cpma = [time_cpma;time_cpma0]; end
+end
+
+
+[gr1,gr2] = ndgrid(time_dma, time_cpma); % grid of CPMA/SMPS times
+[~,idx_t] = max(gr1<gr2,[],2); % get first CPMA time corresponding to DMA scan
+
+if any(idx_t==idx0)
+    warning(['CPMA data files do not contain correct times. ' ...
+        'Only DMA data output.']);
+    sp = get_setpoint();
+    return;
+end
+
+
+% voltage and speed used to define CPMA setpoint
+omega = t1.ClassSpFB_rad_s__(idx_t); % speed of CPMA (no averaging)
+V = t1.AdjVoltageFB_V__(idx_t); % voltage of CPMA (no averaging)
+
+
+% alternate voltages / speeds (unused)
+V_raw = t1.RawVoltageFB_V__(idx_t);
+omega1 = t1.InnerSpFB_rad_s__(idx_t);
+omega2 = t1.OuterSpFB_rad_s__(idx_t);
+
+
+% if prop_pma is not specified
+if ~exist('prop_pma','var'); prop_pma = []; end
+if isempty(prop_pma); prop_pma = kernel.prop_pma; end
+
+
+prop_pma.T = mean(t1.Temperature_C__(idx_t)) + 273.15; % average temperature
+prop_pma.p = mean(t1.ClassPress_Pa__(idx_t) ./ t1.RefPress_Pa__(idx_t)); % average pressure
+sp = get_setpoint(prop_pma,...
+    'omega', omega, 'V', V); % get CPMA setpoint information
+%===============================================================%
+
+
+
+%== Sort based on mass-setpoints ===============================%
+m_star = [sp.m_star];
+[~,idx_sort] = sort(m_star);
+
+sp = sp(idx_sort);
+data = data(:,idx_sort)';
+%===============================================================%
+
+end
+
+
+
+

+io/import_c.m

@@ -10,78 +10,77 @@
 function [data0,d_star0,sp0,prop_dma,prop_pma,out] = ...
     import_c(fnames,prop_pma)
 
+%-- Parse inputs ---------------------------------------------------------%
+addpath tfer_pma; % add tfer_pma dependency
+
+% check if input is a filename structure, or convert to structure
 if ~isstruct(fnames)
     t0 = fnames;
     fnames = struct();
     [fnames.folder,fnames.name,ext] = fileparts(t0);
     fnames.name = [fnames.name,ext];
 end
 
+% if not supplied get properties structure
+if ~exist('prop_pma','var'); prop_pma = []; end
+if isempty(prop_pma); prop_pma = kernel.prop_pma(' CPMA'); end
+%-------------------------------------------------------------------------%
+
+
+% initiate diameters
 data0 = [];
 d_star0 = [];
 
-disp('Reading files...');
 
+disp('Reading files...');
 N = length(fnames); % number of files to read
-if N>1; tools.textbar(0); end
-for ff=1:N
+if N>1; tools.textbar(0); end % initiate commandline textbar
+for ff=1:N % loop through files
     fn = [fnames(ff).folder,'\',fnames(ff).name];
+        % ffth file name
 
 
-    %== Open file and read data =========%
-    fid = fopen(fn);
+    opts = detectImportOptions(fn,'NumHeaderLines',1);
+        % default options
 
-    ii = 0; % automatically detect size of header
-    line = {''};
-    while ~strcmp(line{1},'Point')
-        line = split(fgets(fid),',');
-        ii = ii+1;
-    end
+    %== Proceed with reading data ================%
+    tab0 = readtable(fn, opts);
 
-    [~] = fclose(fid);
+    %-- PMA setpoints -----------------------------------------------%
+    m_star = tab0.m_fg_;
+    V = tab0.V_V__1;
+    omega = tab0.w_rad_s_; % centerline radial speed
+    Rm = []; % CPMA implied resolution, currently N/A
+    p = tab0.P_Pa_; % PMA pressure
+    T = tab0.T_C_; % PMA tempreature
 
 
-    %== Proceed with reading data ========%
-    t = readtable(fn, 'HeaderLines',ii-1);
 
-    %-- PMA setpoints -----------------------------------------------%
-    m_star = table2array(t(:,3));
-    V = table2array(t(:,4));
-    omega = table2array(t(:,5)); % centerline radial speed
-    Rm = table2array(t(:,6)); % CPMA implied resolution
-    p = table2array(t(:,7)); % PMA pressure
-    T = table2array(t(:,8)); % PMA tempreature
-
-    if ~exist('prop_pma','var'); prop_pma = []; end
-    if isempty(prop_pma); prop_pma = kernel.prop_pma(' CPMA'); end
-
-    prop_pma.Q = mean(table2array(t(:,9)))/1000/60; % current averages over all setpoint
+    prop_pma.Q = mean(tab0.Qa_LPM__1) / 1000 / 60; % current assumed constant over all setpoint
             % and does not support changing Q during measurements
     prop_pma.v_bar = prop_pma.Q/prop_pma.A; % average flow velocity
 
     clear sp; % reset PMA setpoints
-    for ii=1:length(V)
-        prop_pma.T = T(ii)+273.15; % pressure converted to Kelvin
-        prop_pma.p = p(ii)/101325; % pressure converted to atm
-        % sp(ii,1) = tfer_pma.get_setpoint(prop_pma,...
-        %     'm_star',m_star(ii).*1e-18,'Rm',Rm(ii));
-        sp(ii,1) = tfer_pma.get_setpoint(prop_pma,...
-            'V',V(ii),'omega',omega(ii));
-    end
+    prop_pma.T = mean(T) + 273.15; % pressure converted to Kelvin
+    prop_pma.p = mean(p) / 101325; % pressure converted to atm
+    % sp(ii,1) = tfer_pma.get_setpoint(prop_pma,...
+    %     'm_star',m_star(ii).*1e-18,'Rm',Rm(ii));
+    sp = get_setpoint(prop_pma,...
+        'V', V, 'omega', omega);
     %----------------------------------------------------------------%
 
 
     %-- DMA setpoints -----------------------------------------------%
-    d_star = table2array(t(:,12)); % DMA setpoints
+    d_star = tab0.dm_nm_; % DMA setpoints
     prop_dma = kernel.prop_dma(' Electrostatic Classifier Model 3080');
-    prop_dma.Q_a = mean(table2array(t(:,9)))/1000/60; % aerosol flow (same as PMA)
+    prop_dma.Q_a = prop_pma.Q; % aerosol flow (same as PMA)
     prop_dma.Q_s = prop_dma.Q_a;
-    prop_dma.Q_c = mean(table2array(t(:,14)))/1000/60; % shield flow
+    prop_dma.Q_c = mean(tab0.Qsh_LPM_) / 1000 / 60; % shield flow
     prop_dma.Q_m = prop_dma.Q_c;
     %----------------------------------------------------------------%
 
 
-    data = table2array(t(:,20)); % data as a vector
+    data = tab0.N_p_cc_; % data as a vector
 
 
     %-- Append to arrays --------------------------------------------%