diff --git a/lopap/Lopap20210311_0951.m b/lopap/Lopap20210311_0951.m index e63c282..1e47835 100644 --- a/lopap/Lopap20210311_0951.m +++ b/lopap/Lopap20210311_0951.m @@ -1,656 +1,656 @@ -%% Lopap Template -% Uses data from 11.03.2021 09h51 -% -%% Info: How to use this -% -% # Make a copy of this file for each Lopap measurement. Name it -% "Lopap"_date_time.m and store it in the folder with the Lopap data. % -% Remove all the general text and comments in your file and add notes and -% observations about the experiment. -% -% # You can run this script, publish it, or evaluate it section wise. For -% new data run it section wise to carefully update the settings and to -% observe the results. -% # Once the data is loaded, you can also run individual sections only. -% # Use the LabBook :: Chrono Notes section as lab book. -% # First, you need to give the date and time of the Lopap measurement to -% load and process the data. -% # All sections that can be modified are marked with the word ":". -% # You also need to give time regions such as for the background fit, as -% well as to calculate means. There are several options how to do this -% -% * While running the experiments and observing the signal in the Lopap -% software, you can enter the date and times directly in the LabBook : -% Select data ranges section. -% * You can manually select the regions in a plot. A window will pop up to -% select the region to modify. To modify several regions, re-run this -% section of the script. This overwrites the current settings. -% * The times of the regions are stored in the "TimeRanges.mat" file. If it -% exists, this is read; the LabBook : Select data ranges section is then -% ignored. -%% ToDO -% -% * Analyse time delay of Lopap and and to exported variables. We need this -% when looking at flow tube data. - -%% Info: Output & Variables -% -% * *lopapReferenceWavelength*: This _string variable_ holds text with -% stating the reference wavelength for channel 1 and 2 used to process the -% data by the Lopap software. -% * *LopapDataCh1 & LopapDataCh2*: These _time tables_ hold the data that -% have been processed by the Lopap Software during the data acquisition. It -% contains the absorption for a number -% of selected wavelength. The absorption is the log of the intensity at the -% reference wavelength relative to the intensity at each wavelength of -% interest. The timetables hold the time in each row and the wavelength in -% the columns. The content is thus identical to that shown in the -% "measurement window" of the Lopap software. -% * *LopapSpecCh1,LopapSpecCh2*: The _time tables_ hold the full spectrum -% for each time. Time is in the rows, wavelength in the columns. The -% wavelength are read from the *.prj file, the spectra data from the *.spec -% file. The information is identical to that shown in the "spectrum window" -% of the Lopap software. -% * *LopapTimeT* This is the main results table: -% -% * time -% * raw signal channel 1 -% * raw signal channel 2 -% * background subtracted channel 1 -% * background subtracted channel 2 -% * calibrated signal channel 1 in ppt -% * calibrated signal channel 2 in ppt -% * HONO in ppt (channel 2 -channel 1) -% -% -% * *TimeRangesCalib.mat*: Time ranges after they have been manually modified. If -% this file exists, the time range definitions in this script are ignored. -% * *LopapMatlabDataDATE_TIMECalib.mat*: Raw data and experimental settings. -% Background subtraction and calibration are not saved in this file. When -% this file exists, data will not be read from the windows server. -% * *LopapDATE_TIMECalib.mat*: Processed data - absorption with time (column 1 -% and 2: 'Ch #1' and 'Ch #2') after background (columns 3 and 4: 'Ch #1 -% Bkgrd' and 'Ch #2 Bkgrd') subtraction and calibration (columns 5 and 6: -% 'Ch #1 Calib' and 'Ch #2 Calib') for one wavelength - that can be used -% for further comparison with other experiments and discussion. - - -%% Matlab: Initialise -% -% This starts Matlab and defines some common settings such as the path to -% the Lopap data. Any modifications can be done at the end of this file in -% where the local function _localInitializeMatlab_ is stored. - - -% clear variables in workspace, clear text in command window -clearvars -clc -% set the path to the data and define some custom colours for plotting -[fileLocation, myColorMap] = localInitializeMatlab; -% : Date when Lopap recording was started. Defines names of Lopap -% files and folders. Seconds are irrelevant, but format is important. -LopapParameters.MeasurementDate = datenum('11-March-2021 09:51:00'); - -%% - -% Print info to command window - -disp('The data is read from:'); -disp(fullfile(fileLocation.server,fileLocation.dataPath,datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM'))); -disp(' ') -disp('Date of experiment:') -disp(datestr(LopapParameters.MeasurementDate)); -disp(' ') -disp('Date of analyis:') -disp(datestr(today)); - -%% Process: Import data -% -% Reads and imports data from Lopap *.dat, *.spec, and *.prj files into: -% -% * *lopapReferenceWavelength*: This _string variable_ holds text with -% stating the reference wavelength for channel 1 and 2 used to process the -% data by the Lopap software. -% * *LopapDataCh1 & LopapDataCh2*: These _time tables_ hold the data that -% have been processed by the Lopap Software during the data acquisition. It -% contains the absorption for a number -% of selected wavelength. The absorption is the log of the intensity at the -% reference wavelength relative to the intensity at each wavelength of -% interest. The timetables hold the time in each row and the wavelength in -% the columns. The content is thus identical to that shown in the -% "measurement window" of the Lopap software. -% * *LopapSpecCh1,LopapSpecCh2*: The _time tables_ hold the full spectrum -% for each time. Time is in the rows, wavelength in the columns. The -% wavelength are read from the *.prj file, the spectra data from the *.spec -% file. The information is identical to that shown in the "spectrum window" -% of the Lopap software. - - -if isfile(['LopapMatlabData' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat']) - load(['LopapMatlabData' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat']) -else -[LopapParameters.instrumentReferenceWavelength, LopapDataCh1,LopapDataCh2,LopapSpecCh1,LopapSpecCh2] ... - = importLopap(LopapParameters.MeasurementDate,fileLocation); -end - -%% Matlab: Run mode -% -% : -% -% * runs without user interaction and without ability to modify data -% ranges, when set to 0; use this to publish. -% * if set to 1: stops for user interaction, -% modify data ranges is enabled.Use this to run. -LopapParameters.Interactive = 1; - -%% LabBook : Measurement Settings % from 20.12.2020 template -% -% : list the Lopap settings here: -LopapParameters.concStandard = 0.01; % mg/l concentration of nitrite in standard used for calibration -LopapParameters.flowReagents = 5./[22+0/60 25+26/60]; % ml/min, enter time for 5 ml in minutes, -% channel 1 and channel 2, first number is minutes, second is seconds. -LopapParameters.flowGas = 396; % ml/min, at room temp and pressure, take average of 10. -LopapParameters.temperature = 298; % temp of gas flow in K -LopapParameters.pressure = 1013; % pressure of gas flow in mbar -LopapParameters.flowGas = LopapParameters.flowGas.*LopapParameters.pressure ./LopapParameters.temperature... - .*298.15./1013.25; % STP ml/min (as used in original LOPAP documents) - -% : Choose a reference wavelength. Default is 650 nm. -LopapParameters.wavelengthR = 650.3; %nm, reference intensity -% : Choose a signal wavelength to calculate the absorption with time -% from the spectra -LopapParameters.wavelengthS = 530; %nm, signal intensity - -%% LabBook : Define data ranges -% -% These settings are ignored if the TimesRanges.mat file exist. If you do -% not want to use the data stored in the TimesRanges.mat file, you need to -% delete it before running this section. -% -% These settings can be modified by selecting data ranges in a plot later. -% So, you do not necessarily need to type the dates and times here. - -if LopapParameters.Interactive == 1 - answer = questdlg('How would you like to set the data ranges of the experiment?', ... - 'Data range selection', ... - 'Read from script code','Read from mat file','Use from workspace','Read from script code'); - % Handle response - switch answer - case 'Read from script code' - disp([answer 'OK. Data range will set as coded in the script.']) - tmpDataRangeSource = 1; - case 'Read from mat file' - disp([answer 'OK. Data range will loaded from mat file in current directory.']) - tmpDataRangeSource = 2; - case 'Use from workspace' - disp([answer 'OK. Data range in workspace variable will be used.']) - tmpDataRangeSource = 3; - case 'Cancel' - disp('OK. Do nothing') - tmpDataRangeSource = 0; - end -else - if isfile(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') '.mat']) - tmpDataRangeSource = 2; - disp([answer ' OK. Data range will loaded from mat file in current directory.']) - elseif ~exist('TimeRange','var') - tmpDataRangeSource = 1; - disp([answer 'OK. Data range will set as coded in the script.']) - elseif exist('TimeRange','var') - tmpDataRangeSource = 3; - disp(' OK. Data range in workspace variable will be used.') - end -end - -if tmpDataRangeSource == 2 - load(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') '.mat']) - disp('TimeRange in TimeRange*.mat file in current directory used') - disp('You need to delete the TimeRange variable in the workspace AND the TimeRange*.mat file in the current directory') - disp('to use the dates entered in the script') - disp('You keep to delete the TimeRange variable in the workspace AND delete TimeRange*.mat file in the current directory') - disp('to use the dates saved in the workspace') -elseif tmpDataRangeSource == 1 - % (Optional - comment if no needed) state the start and end - % time of the time that you want to look at (to exclude very noisy - % signal in the beginning and at the end of the data) - TimeRange.calibrationRange = [datenum('15-Nov-2021 16:30:00') datenum('15-Nov-2021 17:20:00')]; - % : State the time for the first spectrum that you like to - % compare (to confirm nitrite) - TimeRange.specCompareOne = datenum('15-Nov-2021 16:49:30 '); - % : State the time for one spectrum of the second spectrum that - % you like to compare (to confirm nitrite) - TimeRange.specCompareTwo = datenum('15-Nov-2021 17:02:00 '); - % : Choose a start and end time for one as many intervals that can - % be used to fit the background. Background is the signal when R1&R2 was - % attached to the Lopap with zero air entering the sampling coil. - TimeRange.backgroundFit = [datenum('15-Nov-2021 15:10:00') datenum('15-Nov-2021 15:50:00');... - datenum('0-Dec-2021 16:26:00') datenum('0-Dec-2021 16:26:00')]; - % : State the start and end time of data that can be used to - % calculate the mean of the zero air signal without nitrite - TimeRange.MeanZeroAir = [datenum('04-Feb-2022 11:49:07') datenum('04-Feb-2022 11:59:14 ')]; - % : State the start and end time of data that can be used to - % calculate the mean of the zero air signal with nitrite - TimeRange.MeanNitriteStandard = [datenum('04-Feb-2022 12:26:28') datenum('04-Feb-2022 12:39:46')]; - - save(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'TimeRange'); - disp('TimeRange as coded in script used') - -elseif tmpDataRangeSource == 3 - if exist('TimeRange','var') - disp('TimeRange in Workspace used') - disp('You need to delete the TimeRange variable in the workspace AND the TimeRange*.mat file in the current directory') - disp('to use the dates entered in the script') - elseif ~exist('TimeRange','var') - warning('TimeRange in the Workspace was selected, but there is no TimeRange variable') - end -end - -% confirm that all dates are reasonable -tmp = fieldnames(TimeRange); -for ii = 1:length(tmp) - mustBeLopapReasonableDate(TimeRange.(tmp{ii})) -end -clearvars tmp* - -%% Process: Modify data ranges -% -% Select the range to modify in the pop-up window and the number of fit -% ranges for the background fit in the 2nd pop-up window. Then, a figure -% will show, select the data ranges in the figure: -% -% * a line for calibration range, mean zero air and mean -% nitrite standard. Select the range by keeping the mouse pushed while -% moving from start to end time. -% * a line for each background region, default is 2. Select the range by -% keeping the mouse pushed while moving from start to end time. -% * a point for spectrum with reagents and spectrum with nitrite. Click -% once to select a specific time. -% -% * LopapParameters.Interactive needs to be 1 for this to run. -% -% # Run this section several times to modify several time regions -% # Press cancel in the pop-up window if you do not want to modify any -% region. - -TimeRange = modifyLopapTimeRangesCalib(TimeRange, LopapParameters, LopapDataCh1, LopapDataCh2); -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* -save(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') '.mat'],'TimeRange'); - -%% Process: Find closest match of selected dates in data set -% -% The data is sampled at specific times with a certain interval. The -% selected or chosen dates might not necessarily exactly match the data. -% Here, closest match to the given dates are selected. - -tmp.d1 = TimeRange.specCompareOne; -tmp.d2 = LopapSpecCh1.Time; -[~,ind1] = min(abs(bsxfun(@minus,datenum(tmp.d2),datenum(tmp.d1)'))); -TimeRange.specCompareOne = tmp.d2(ind1,:); - -tmp.d1 = TimeRange.specCompareTwo ; -tmp.d2 = LopapSpecCh1.Time; -[~,ind1] = min(abs(bsxfun(@minus,datenum(tmp.d2),datenum(tmp.d1)'))); -TimeRange.specCompareTwo = tmp.d2(ind1,:); - -% save current state of -clearvars -except Lopap* fileLocation myColorMap TimeRange* -save(['LopapMatlabData' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'LopapDataCh1','LopapDataCh2','LopapParameters','LopapSpecCh1','LopapSpecCh2') -save(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'TimeRange'); - -%% LabBook: Data ranges used in this work -% -clc -fprintf('Calibration range: %s to %s \n',datestr(TimeRange.calibrationRange(1)), datestr(TimeRange.calibrationRange(2))) -for ii=1:size(TimeRange.backgroundFit,1) - fprintf('Background range: %s to %s \n',datestr(TimeRange.backgroundFit(ii,1)), datestr(TimeRange.backgroundFit(ii,2))) -end -fprintf('Zero mean: %s to %s \n',datestr(TimeRange.MeanZeroAir(1)), datestr(TimeRange.MeanZeroAir(2))) -fprintf('Nitrite mean: %s to %s \n',datestr(TimeRange.MeanNitriteStandard(1)), datestr(TimeRange.MeanNitriteStandard(2))) -fprintf('Comparison Spectrum #1: %s \n',datestr(TimeRange.specCompareOne(1))) -fprintf('Comparison Spectrum #2: %s \n',datestr(TimeRange.specCompareTwo(1))) - - -%% Plot all data -% -% Here, we plot the data as shown in the "Measurement window" of the Lopap -% software. That is, the processed $I^{*}_{i}$ data for all wavelength as -% defined in the settings of the Lopap software. -% -% -% For every measurement wavelength I_i the program calculates the -% measurement (absorption) value I_i^* from the spectral intensity on the -% measurement wavelength I_i and a reference wavelength I_r with the -% following operation: -% -% -% $$I^{*}_{i} = log_{10} (\frac{I_r}{I_i})$$ -% -% -% -% The plots show time on the x-axis and the ABS intensity on the y axis for -% each wavelength in separate graphs. One Figure for each Channel. - -figure( 'Name', 'LopapChannel1','NumberTitle','off','Position',[10 100 550 700]); -ch1s=stackedplot(LopapDataCh1); -title(LopapParameters.instrumentReferenceWavelength(1)) -for ii = 1:width(LopapDataCh1) - % : The numbers in the brackets give the lower and upper limit - % of the y-axis. Adapt these to fit your data and remove the % sign. - % ch1s.AxesProperties(ii).YLimits = [-0.4 -0.1]; -end -grid ON -% Add title and axis labels -xlabel('Lab Time') - -figure( 'Name', 'LopapChannel2','NumberTitle','off','Position',[600 100 550 700]); -ch2s=stackedplot(LopapDataCh2); -title(LopapParameters.instrumentReferenceWavelength(2)) -for ii = 1:width(LopapDataCh2) - % : The numbers in the brackets give the lower and upper limit - % of the y-axis. Adapt these to fit your data and remove the % sign. - % ch2s.AxesProperties(ii).YLimits = [-0.2 0.3]; -end -grid ON -% Add title and axis labels -xlabel('Lab Time') -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* -if LopapParameters.Interactive == 1 - distFig - fprintf('press any key if you want to proceed\n') - pause - %% - % - close all -end - -%% Plot all data :: selected time range -% - -TR = timerange(datestr(TimeRange.calibrationRange(1)),datestr(TimeRange.calibrationRange(2))); -figure( 'Name', 'LopapChannel1','NumberTitle','off','Position',[10 100 550 700]); -ch1s=stackedplot(LopapDataCh1(TR,:)); -title(LopapParameters.instrumentReferenceWavelength(1)) -for ii = 1:width(LopapDataCh1) - % : The numbers in the brackets give the lower and upper limit - % of the y-axis. Adapt these to fit your data and remove the % sign. - % ch1s.AxesProperties(ii).YLimits = [-0.4 -0.1]; -end -grid ON -% Add title and axis labels -xlabel('Lab Time') - -figure( 'Name', 'LopapChannel2','NumberTitle','off','Position',[600 100 550 700]); -ch2s=stackedplot(LopapDataCh2(TR,:)); -title(LopapParameters.instrumentReferenceWavelength(2)) -for ii = 1:width(LopapDataCh2) - % : The numbers in the brackets give the lower and upper limit - % of the y-axis. Adapt these to fit your data and remove the % sign. - % ch2s.AxesProperties(ii).YLimits = [-0.2 0.3]; -end -grid ON -% Add title and axis labels -xlabel('Lab Time') -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* -if LopapParameters.Interactive == 1 - distFig - fprintf('press any key if you want to proceed\n') - pause - %% - % - close all -end - - -%% Plot selected data in selected time range -% -% Here we plot the $I^{*}_{i}$ data for one wavelength only. - -TR = timerange(datestr(TimeRange.calibrationRange(1)),datestr(TimeRange.calibrationRange(2))); -figure( 'Name', 'LopapChannel1 550','NumberTitle','off','Position',[10 100 550 800]); - % : The numbers in the curled brackets give wavelength for which you'd - % like to display the data. The numbers need to be exactly match the - % labels in the previous figure. - ch1s=stackedplot(LopapDataCh1(TR,:),{'550.7'}); - title(LopapParameters.instrumentReferenceWavelength(1)) - % : The numbers in the brackets give the lower and upper limit - % of the y-axis. Adapt these to fit your data and remove the % sign. - %ch1s.AxesProperties.YLimits = [-0.4 -0.1]; - grid ON - % Add title and axis labels - xlabel('Lab Time') - - -figure( 'Name', 'LopapChannel2 550','NumberTitle','off','Position',[600 100 550 800]); - % : The numbers in the curled brackets give wavelength for which you'd - % like to display the data. The numbers need to be exactly match the - % labels in the previous figure - ch2s=stackedplot(LopapDataCh2(TR,:),{'550.5'}); - title(LopapParameters.instrumentReferenceWavelength(2)) - % : The numbers in the brackets give the lower and upper limit - % of the y-axis. Adapt these to fit your data and remove the % sign. - % ch2s.AxesProperties.YLimits = [-0.2 0.3]; - grid ON - % Add title and axis labels - xlabel('Lab Time') - -clearvars -except Lopap* fileLocation myColorMap lopap* TimeRange* -if LopapParameters.Interactive == 1 - distFig - fprintf('press any key if you want to proceed\n') - pause - %% - % - close all -end - -%% Plot selected spectra -% -% Here, we plot a spectrum without ('R1&R2') and with nitrite -% ('R1+nitrite&R2'') for each channel or which ever combination of spectra -% that we have selected. -clc -close all -figure('Name', 'LopapSelectedSpectra') -plot(str2double(LopapSpecCh1.Properties.VariableNames),LopapSpecCh1{datestr(TimeRange.specCompareOne),:},... - 'DisplayName',['#1 ' datestr(TimeRange.specCompareOne) ]) % a time where only R1&R2 was connected to the Lopap -hold on -plot(str2double(LopapSpecCh2.Properties.VariableNames),LopapSpecCh2{datestr(TimeRange.specCompareOne),:},... - 'DisplayName',['#2 ' datestr(TimeRange.specCompareOne)],'LineStyle','--') % a time where only R1&R2 was connected to the Lopap -plot(str2double(LopapSpecCh1.Properties.VariableNames),LopapSpecCh1{datestr(TimeRange.specCompareTwo),:},... - 'DisplayName',['#1 ' datestr(TimeRange.specCompareTwo)]) % a time where nitrite in R1 was connected to the Lopap -plot(str2double(LopapSpecCh2.Properties.VariableNames),LopapSpecCh2{datestr(TimeRange.specCompareTwo),:},... - 'DisplayName',['#2 ' datestr(TimeRange.specCompareTwo)],'LineStyle','--') % a time where nitrite in R1 was connected to the Lopap -legend -box ON -grid ON -% Add title and axis labels -xlabel('wavelength [nm]') -ylabel('Intensity') - -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* -if LopapParameters.Interactive == 1 - distFig - fprintf('press any key if you want to proceed\n') - pause - %% - % - close all -end -%% Process: Absorption spectrum -% -% Here, we reprocess the data to derive an absorption spectrum. This should -% match the spectra shown for nitrite in the Lopap manual. - - -clc -close all -indexR(1) = find(str2double(LopapSpecCh1.Properties.VariableNames) >= LopapParameters.wavelengthR,1); -indexR(2) = find(str2double(LopapSpecCh2.Properties.VariableNames) >= LopapParameters.wavelengthR,1); - -figure('Name', 'LopapAbsorbtionSpectra') -plot(str2double(LopapSpecCh1.Properties.VariableNames),log10(LopapSpecCh1{datestr(TimeRange.specCompareTwo),... - indexR(1)}./LopapSpecCh1{datestr(TimeRange.specCompareTwo),:})-log10(LopapSpecCh1{datestr(TimeRange.specCompareOne),... - indexR(1)}./LopapSpecCh1{datestr(TimeRange.specCompareOne),:}),'DisplayName','#1') -hold on -plot(str2double(LopapSpecCh2.Properties.VariableNames),log10(LopapSpecCh2{datestr(TimeRange.specCompareTwo),... - indexR(2)}./LopapSpecCh2{datestr(TimeRange.specCompareTwo),:})-log10(LopapSpecCh2{datestr(TimeRange.specCompareOne),... - indexR(2)}./LopapSpecCh2{datestr(TimeRange.specCompareOne),:}),'DisplayName','#2','LineStyle','--') -legend -box ON -grid ON -% Add title and axis labels -xlabel('wavelength [nm]') -ylabel('ABS (log I_{r}/I_{i}') -title([datestr(TimeRange.specCompareOne) ' & ' datestr(TimeRange.specCompareTwo)]) -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* -if LopapParameters.Interactive == 1 - distFig - fprintf('press any key if you want to proceed\n') - pause - %% - % - close all -end -%% Process :: Absorption with time -% -% Here, we reprocess the data to get the intensity with time for 1 chosen -% signal wavelength and 1 reference wavelength. The data is identical to that -% shown in the "measurement window" of the Lopap software. -% -% The data is written to the _LopapTimeT_ time table. Rows hold the time, -% columns hold the processed signal for channel 1 and channel 2. -% -% It is important that the same signal wavelength and reference wavelength -% are taken for the calibration and for the experimental data. - -clc -close all -indexS(1) = find(str2double(LopapSpecCh1.Properties.VariableNames) >= LopapParameters.wavelengthS,1); -indexS(2) = find(str2double(LopapSpecCh2.Properties.VariableNames) >= LopapParameters.wavelengthS,1); -indexR(1) = find(str2double(LopapSpecCh1.Properties.VariableNames) >= LopapParameters.wavelengthR,1); -indexR(2) = find(str2double(LopapSpecCh2.Properties.VariableNames) >= LopapParameters.wavelengthR,1); -LopapTimeT = timetable(LopapSpecCh1.Time,log10(LopapSpecCh1{:,indexR(1)}./LopapSpecCh1{:,indexS(1)}),... - log10(LopapSpecCh2{:,indexR(2)}./LopapSpecCh2{:,indexS(2)})); -LopapTimeT.Properties.VariableNames ={'Ch #1','Ch #2'}; -figure('Name', 'Lopap TimeTable') -stackedplot(LopapTimeT,{{'Ch #1','Ch #2'}}); -title(['Lopap' datestr(LopapParameters.MeasurementDate)]) -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* -if LopapParameters.Interactive == 1 - distFig - fprintf('press any key if you want to proceed\n') - pause - %% - % - close all -end -%% Process :: Background -% -% Currently a linear fit forced through 0. -clc -close all -LopapTimeT = fitLopapBkgrd(LopapTimeT,'BkgrdPeriod',TimeRange.backgroundFit); -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* -if LopapParameters.Interactive == 1 - distFig - fprintf('press any key if you want to proceed\n') - pause - %% - % - close all -end -%% -% -clc -close all -% Plot -figure('Name', 'Lopap TimeTable Bkgrd') -s=stackedplot(LopapTimeT,{{'Ch #1 Bkgrd','Ch #2 Bkgrd'}}); -title(['Lopap' datestr(LopapParameters.MeasurementDate)]) -% s.AxesProperties.YLimits = [-0.01 0.6]; -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* - -%% Process :: Calibrate -% -clc -close all - -S = timerange(datestr(TimeRange.MeanZeroAir(1)), datestr(TimeRange.MeanZeroAir(2))); -TZ = LopapTimeT(S,:); -S = timerange(datestr(TimeRange.MeanNitriteStandard(1)), datestr(TimeRange.MeanNitriteStandard(2))); -TN = LopapTimeT(S,:); -calib.mean = mean(TN{:,:},1) - mean(TZ{:,:},1); -calib.std = std(TN{:,:},1); -calib.mixingRatio = LopapParameters.concStandard./1000.*LopapParameters.flowReagents./1000./46.*6.022E23./... - (246E17.*LopapParameters.flowGas).*1E12; - -[LopapTimeT, calib.fitresult] = fitLopapCalibration(LopapTimeT,calib); -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* fitresult calib - -fprintf('Detection limit channel #1 is %.0f ppt \n',calib.fitresult.calibChannel1(3.*calib.std(3))) -fprintf('Detection limit channel #2 is %.0f ppt \n',calib.fitresult.calibChannel2(3.*calib.std(4))) - -%% Plot calibrated data -% -clc -close all -% Plot -figure('Name', 'Lopap TimeTable Calibrated') -s=stackedplot(LopapTimeT,{{'Ch #1 Calib','Ch #2 Calib'}}); -title(['HONO in ppt, Lopap' datestr(LopapParameters.MeasurementDate)]) -% s.AxesProperties.YLimits = [-100 2500]; -clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* fitresult calib -if LopapParameters.Interactive == 1 - distFig - fprintf('press any key if you want to proceed\n') - pause - %% - % - close all -end -%% Save the data -% -save(['LopapMatlabData' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'LopapDataCh1','LopapDataCh2','LopapParameters','calib','LopapSpecCh1','LopapSpecCh2') -save(['Lopap' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'LopapParameters','LopapTimeT','TimeRange') -%% Local functions -function [fileLocation, myColorMap] = localInitializeMatlab -%INITIALIZEMATLAB Sets Lopap specific settings -% Adds processing functions to path and defines file locations for mac, windows, -% and Matlab online (which runs on unix) - -% change working directory -path = matlab.desktop.editor.getActiveFilename; -[filepath,~,~] = fileparts(path); -cd(filepath); - - -% Data Server Address -if ismac - fileLocation.server = '/Volumes/'; - fileLocation.dataPath = 'data_2411/Lopap/'; - % add Lopap functions and general add-ons to path - addpath('/Users/bartels-rausch/SwitchDrive/PSI_Matlab/Add-Ons'); - addpath('/Users/bartels-rausch/SwitchDrive/PSI_Matlab/Lopap/1.LOPAP_Restliche Daten'); -elseif ispc - fileLocation.server = '\\fs00\'; %UNC name for file server that hosts Napp shared drive - fileLocation.dataPath = 'data_2411\Lopap\'; - % add Lopap functions and general add-ons to path - - addpath('C:\Users\vieill_l\switchdrive\Add-Ons'); - - - addpath('C:\Users\vieill_l\switchdrive\Lopap\1.LOPAP_Restliche Daten'); - - -elseif isunix - fileLocation.server = '/'; - fileLocation.dataPath = 'Lopap/'; -end - -% Color maps -myColorMap.red = colormap(autumn(7)); -myColorMap.blue = colormap(winter(7)); - -close all -end - +%% Lopap Template +% Uses data from 11.03.2021 09h51 +% +%% Info: How to use this +% +% # Make a copy of this file for each Lopap measurement. Name it +% "Lopap"_date_time.m and store it in the folder with the Lopap data. % +% Remove all the general text and comments in your file and add notes and +% observations about the experiment. +% +% # You can run this script, publish it, or evaluate it section wise. For +% new data run it section wise to carefully update the settings and to +% observe the results. +% # Once the data is loaded, you can also run individual sections only. +% # Use the LabBook :: Chrono Notes section as lab book. +% # First, you need to give the date and time of the Lopap measurement to +% load and process the data. +% # All sections that can be modified are marked with the word ":". +% # You also need to give time regions such as for the background fit, as +% well as to calculate means. There are several options how to do this +% +% * While running the experiments and observing the signal in the Lopap +% software, you can enter the date and times directly in the LabBook : +% Select data ranges section. +% * You can manually select the regions in a plot. A window will pop up to +% select the region to modify. To modify several regions, re-run this +% section of the script. This overwrites the current settings. +% * The times of the regions are stored in the "TimeRanges.mat" file. If it +% exists, this is read; the LabBook : Select data ranges section is then +% ignored. +%% ToDO +% +% * Analyse time delay of Lopap and and to exported variables. We need this +% when looking at flow tube data. + +%% Info: Output & Variables +% +% * *lopapReferenceWavelength*: This _string variable_ holds text with +% stating the reference wavelength for channel 1 and 2 used to process the +% data by the Lopap software. +% * *LopapDataCh1 & LopapDataCh2*: These _time tables_ hold the data that +% have been processed by the Lopap Software during the data acquisition. It +% contains the absorption for a number +% of selected wavelength. The absorption is the log of the intensity at the +% reference wavelength relative to the intensity at each wavelength of +% interest. The timetables hold the time in each row and the wavelength in +% the columns. The content is thus identical to that shown in the +% "measurement window" of the Lopap software. +% * *LopapSpecCh1,LopapSpecCh2*: The _time tables_ hold the full spectrum +% for each time. Time is in the rows, wavelength in the columns. The +% wavelength are read from the *.prj file, the spectra data from the *.spec +% file. The information is identical to that shown in the "spectrum window" +% of the Lopap software. +% * *LopapTimeT* This is the main results table: +% +% * time +% * raw signal channel 1 +% * raw signal channel 2 +% * background subtracted channel 1 +% * background subtracted channel 2 +% * calibrated signal channel 1 in ppt +% * calibrated signal channel 2 in ppt +% * HONO in ppt (channel 2 -channel 1) +% +% +% * *TimeRangesCalib.mat*: Time ranges after they have been manually modified. If +% this file exists, the time range definitions in this script are ignored. +% * *LopapMatlabDataDATE_TIMECalib.mat*: Raw data and experimental settings. +% Background subtraction and calibration are not saved in this file. When +% this file exists, data will not be read from the windows server. +% * *LopapDATE_TIMECalib.mat*: Processed data - absorption with time (column 1 +% and 2: 'Ch #1' and 'Ch #2') after background (columns 3 and 4: 'Ch #1 +% Bkgrd' and 'Ch #2 Bkgrd') subtraction and calibration (columns 5 and 6: +% 'Ch #1 Calib' and 'Ch #2 Calib') for one wavelength - that can be used +% for further comparison with other experiments and discussion. + + +%% Matlab: Initialise +% +% This starts Matlab and defines some common settings such as the path to +% the Lopap data. Any modifications can be done at the end of this file in +% where the local function _localInitializeMatlab_ is stored. + + +% clear variables in workspace, clear text in command window +clearvars +clc +% set the path to the data and define some custom colours for plotting +[fileLocation, myColorMap] = localInitializeMatlab; +% : Date when Lopap recording was started. Defines names of Lopap +% files and folders. Seconds are irrelevant, but format is important. +LopapParameters.MeasurementDate = datenum('11-March-2021 09:51:00'); + +%% + +% Print info to command window + +disp('The data is read from:'); +disp(fullfile(fileLocation.server,fileLocation.dataPath,datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM'))); +disp(' ') +disp('Date of experiment:') +disp(datestr(LopapParameters.MeasurementDate)); +disp(' ') +disp('Date of analyis:') +disp(datestr(today)); + +%% Process: Import data +% +% Reads and imports data from Lopap *.dat, *.spec, and *.prj files into: +% +% * *lopapReferenceWavelength*: This _string variable_ holds text with +% stating the reference wavelength for channel 1 and 2 used to process the +% data by the Lopap software. +% * *LopapDataCh1 & LopapDataCh2*: These _time tables_ hold the data that +% have been processed by the Lopap Software during the data acquisition. It +% contains the absorption for a number +% of selected wavelength. The absorption is the log of the intensity at the +% reference wavelength relative to the intensity at each wavelength of +% interest. The timetables hold the time in each row and the wavelength in +% the columns. The content is thus identical to that shown in the +% "measurement window" of the Lopap software. +% * *LopapSpecCh1,LopapSpecCh2*: The _time tables_ hold the full spectrum +% for each time. Time is in the rows, wavelength in the columns. The +% wavelength are read from the *.prj file, the spectra data from the *.spec +% file. The information is identical to that shown in the "spectrum window" +% of the Lopap software. + + +if isfile(['LopapMatlabData' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat']) + load(['LopapMatlabData' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat']) +else +[LopapParameters.instrumentReferenceWavelength, LopapDataCh1,LopapDataCh2,LopapSpecCh1,LopapSpecCh2] ... + = importLopap(LopapParameters.MeasurementDate,fileLocation); +end + +%% Matlab: Run mode +% +% : +% +% * runs without user interaction and without ability to modify data +% ranges, when set to 0; use this to publish. +% * if set to 1: stops for user interaction, +% modify data ranges is enabled.Use this to run. +LopapParameters.Interactive = 1; + +%% LabBook : Measurement Settings % from 20.12.2020 template +% +% : list the Lopap settings here: +LopapParameters.concStandard = 0.01; % mg/l concentration of nitrite in standard used for calibration +LopapParameters.flowReagents = 5./[22+0/60 25+26/60]; % ml/min, enter time for 5 ml in minutes, +% channel 1 and channel 2, first number is minutes, second is seconds. +LopapParameters.flowGas = 396; % ml/min, at room temp and pressure, take average of 10. +LopapParameters.temperature = 298; % temp of gas flow in K +LopapParameters.pressure = 1013; % pressure of gas flow in mbar +LopapParameters.flowGas = LopapParameters.flowGas.*LopapParameters.pressure ./LopapParameters.temperature... + .*298.15./1013.25; % STP ml/min (as used in original LOPAP documents) + +% : Choose a reference wavelength. Default is 650 nm. +LopapParameters.wavelengthR = 650.3; %nm, reference intensity +% : Choose a signal wavelength to calculate the absorption with time +% from the spectra +LopapParameters.wavelengthS = 530; %nm, signal intensity + +%% LabBook : Define data ranges +% +% These settings are ignored if the TimesRanges.mat file exist. If you do +% not want to use the data stored in the TimesRanges.mat file, you need to +% delete it before running this section. +% +% These settings can be modified by selecting data ranges in a plot later. +% So, you do not necessarily need to type the dates and times here. + +if LopapParameters.Interactive == 1 + answer = questdlg('How would you like to set the data ranges of the experiment?', ... + 'Data range selection', ... + 'Read from script code','Read from mat file','Use from workspace','Read from script code'); + % Handle response + switch answer + case 'Read from script code' + disp([answer 'OK. Data range will set as coded in the script.']) + tmpDataRangeSource = 1; + case 'Read from mat file' + disp([answer 'OK. Data range will loaded from mat file in current directory.']) + tmpDataRangeSource = 2; + case 'Use from workspace' + disp([answer 'OK. Data range in workspace variable will be used.']) + tmpDataRangeSource = 3; + case 'Cancel' + disp('OK. Do nothing') + tmpDataRangeSource = 0; + end +else + if isfile(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') '.mat']) + tmpDataRangeSource = 2; + disp([answer ' OK. Data range will loaded from mat file in current directory.']) + elseif ~exist('TimeRange','var') + tmpDataRangeSource = 1; + disp([answer 'OK. Data range will set as coded in the script.']) + elseif exist('TimeRange','var') + tmpDataRangeSource = 3; + disp(' OK. Data range in workspace variable will be used.') + end +end + +if tmpDataRangeSource == 2 + load(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') '.mat']) + disp('TimeRange in TimeRange*.mat file in current directory used') + disp('You need to delete the TimeRange variable in the workspace AND the TimeRange*.mat file in the current directory') + disp('to use the dates entered in the script') + disp('You keep to delete the TimeRange variable in the workspace AND delete TimeRange*.mat file in the current directory') + disp('to use the dates saved in the workspace') +elseif tmpDataRangeSource == 1 + % (Optional - comment if no needed) state the start and end + % time of the time that you want to look at (to exclude very noisy + % signal in the beginning and at the end of the data) + TimeRange.calibrationRange = [datenum('15-Nov-2021 16:30:00') datenum('15-Nov-2021 17:20:00')]; + % : State the time for the first spectrum that you like to + % compare (to confirm nitrite) + TimeRange.specCompareOne = datenum('15-Nov-2021 16:49:30 '); + % : State the time for one spectrum of the second spectrum that + % you like to compare (to confirm nitrite) + TimeRange.specCompareTwo = datenum('15-Nov-2021 17:02:00 '); + % : Choose a start and end time for one as many intervals that can + % be used to fit the background. Background is the signal when R1&R2 was + % attached to the Lopap with zero air entering the sampling coil. + TimeRange.backgroundFit = [datenum('15-Nov-2021 15:10:00') datenum('15-Nov-2021 15:50:00');... + datenum('0-Dec-2021 16:26:00') datenum('0-Dec-2021 16:26:00')]; + % : State the start and end time of data that can be used to + % calculate the mean of the zero air signal without nitrite + TimeRange.MeanZeroAir = [datenum('04-Feb-2022 11:49:07') datenum('04-Feb-2022 11:59:14 ')]; + % : State the start and end time of data that can be used to + % calculate the mean of the zero air signal with nitrite + TimeRange.MeanNitriteStandard = [datenum('04-Feb-2022 12:26:28') datenum('04-Feb-2022 12:39:46')]; + + save(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'TimeRange'); + disp('TimeRange as coded in script used') + +elseif tmpDataRangeSource == 3 + if exist('TimeRange','var') + disp('TimeRange in Workspace used') + disp('You need to delete the TimeRange variable in the workspace AND the TimeRange*.mat file in the current directory') + disp('to use the dates entered in the script') + elseif ~exist('TimeRange','var') + warning('TimeRange in the Workspace was selected, but there is no TimeRange variable') + end +end + +% confirm that all dates are reasonable +tmp = fieldnames(TimeRange); +for ii = 1:length(tmp) + mustBeLopapReasonableDate(TimeRange.(tmp{ii})) +end +clearvars tmp* + +%% Process: Modify data ranges +% +% Select the range to modify in the pop-up window and the number of fit +% ranges for the background fit in the 2nd pop-up window. Then, a figure +% will show, select the data ranges in the figure: +% +% * a line for calibration range, mean zero air and mean +% nitrite standard. Select the range by keeping the mouse pushed while +% moving from start to end time. +% * a line for each background region, default is 2. Select the range by +% keeping the mouse pushed while moving from start to end time. +% * a point for spectrum with reagents and spectrum with nitrite. Click +% once to select a specific time. +% +% * LopapParameters.Interactive needs to be 1 for this to run. +% +% # Run this section several times to modify several time regions +% # Press cancel in the pop-up window if you do not want to modify any +% region. + +TimeRange = modifyLopapTimeRangesCalib(TimeRange, LopapParameters, LopapDataCh1, LopapDataCh2); +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* +save(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') '.mat'],'TimeRange'); + +%% Process: Find closest match of selected dates in data set +% +% The data is sampled at specific times with a certain interval. The +% selected or chosen dates might not necessarily exactly match the data. +% Here, closest match to the given dates are selected. + +tmp.d1 = TimeRange.specCompareOne; +tmp.d2 = LopapSpecCh1.Time; +[~,ind1] = min(abs(bsxfun(@minus,datenum(tmp.d2),datenum(tmp.d1)'))); +TimeRange.specCompareOne = tmp.d2(ind1,:); + +tmp.d1 = TimeRange.specCompareTwo ; +tmp.d2 = LopapSpecCh1.Time; +[~,ind1] = min(abs(bsxfun(@minus,datenum(tmp.d2),datenum(tmp.d1)'))); +TimeRange.specCompareTwo = tmp.d2(ind1,:); + +% save current state of +clearvars -except Lopap* fileLocation myColorMap TimeRange* +save(['LopapMatlabData' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'LopapDataCh1','LopapDataCh2','LopapParameters','LopapSpecCh1','LopapSpecCh2') +save(['TimeRanges' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'TimeRange'); + +%% LabBook: Data ranges used in this work +% +clc +fprintf('Calibration range: %s to %s \n',datestr(TimeRange.calibrationRange(1)), datestr(TimeRange.calibrationRange(2))) +for ii=1:size(TimeRange.backgroundFit,1) + fprintf('Background range: %s to %s \n',datestr(TimeRange.backgroundFit(ii,1)), datestr(TimeRange.backgroundFit(ii,2))) +end +fprintf('Zero mean: %s to %s \n',datestr(TimeRange.MeanZeroAir(1)), datestr(TimeRange.MeanZeroAir(2))) +fprintf('Nitrite mean: %s to %s \n',datestr(TimeRange.MeanNitriteStandard(1)), datestr(TimeRange.MeanNitriteStandard(2))) +fprintf('Comparison Spectrum #1: %s \n',datestr(TimeRange.specCompareOne(1))) +fprintf('Comparison Spectrum #2: %s \n',datestr(TimeRange.specCompareTwo(1))) + + +%% Plot all data +% +% Here, we plot the data as shown in the "Measurement window" of the Lopap +% software. That is, the processed $I^{*}_{i}$ data for all wavelength as +% defined in the settings of the Lopap software. +% +% +% For every measurement wavelength I_i the program calculates the +% measurement (absorption) value I_i^* from the spectral intensity on the +% measurement wavelength I_i and a reference wavelength I_r with the +% following operation: +% +% +% $$I^{*}_{i} = log_{10} (\frac{I_r}{I_i})$$ +% +% +% +% The plots show time on the x-axis and the ABS intensity on the y axis for +% each wavelength in separate graphs. One Figure for each Channel. + +figure( 'Name', 'LopapChannel1','NumberTitle','off','Position',[10 100 550 700]); +ch1s=stackedplot(LopapDataCh1); +title(LopapParameters.instrumentReferenceWavelength(1)) +for ii = 1:width(LopapDataCh1) + % : The numbers in the brackets give the lower and upper limit + % of the y-axis. Adapt these to fit your data and remove the % sign. + % ch1s.AxesProperties(ii).YLimits = [-0.4 -0.1]; +end +grid ON +% Add title and axis labels +xlabel('Lab Time') + +figure( 'Name', 'LopapChannel2','NumberTitle','off','Position',[600 100 550 700]); +ch2s=stackedplot(LopapDataCh2); +title(LopapParameters.instrumentReferenceWavelength(2)) +for ii = 1:width(LopapDataCh2) + % : The numbers in the brackets give the lower and upper limit + % of the y-axis. Adapt these to fit your data and remove the % sign. + % ch2s.AxesProperties(ii).YLimits = [-0.2 0.3]; +end +grid ON +% Add title and axis labels +xlabel('Lab Time') +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* +if LopapParameters.Interactive == 1 + distFig + fprintf('press any key if you want to proceed\n') + pause + %% + % + close all +end + +%% Plot all data :: selected time range +% + +TR = timerange(datestr(TimeRange.calibrationRange(1)),datestr(TimeRange.calibrationRange(2))); +figure( 'Name', 'LopapChannel1','NumberTitle','off','Position',[10 100 550 700]); +ch1s=stackedplot(LopapDataCh1(TR,:)); +title(LopapParameters.instrumentReferenceWavelength(1)) +for ii = 1:width(LopapDataCh1) + % : The numbers in the brackets give the lower and upper limit + % of the y-axis. Adapt these to fit your data and remove the % sign. + % ch1s.AxesProperties(ii).YLimits = [-0.4 -0.1]; +end +grid ON +% Add title and axis labels +xlabel('Lab Time') + +figure( 'Name', 'LopapChannel2','NumberTitle','off','Position',[600 100 550 700]); +ch2s=stackedplot(LopapDataCh2(TR,:)); +title(LopapParameters.instrumentReferenceWavelength(2)) +for ii = 1:width(LopapDataCh2) + % : The numbers in the brackets give the lower and upper limit + % of the y-axis. Adapt these to fit your data and remove the % sign. + % ch2s.AxesProperties(ii).YLimits = [-0.2 0.3]; +end +grid ON +% Add title and axis labels +xlabel('Lab Time') +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* +if LopapParameters.Interactive == 1 + distFig + fprintf('press any key if you want to proceed\n') + pause + %% + % + close all +end + + +%% Plot selected data in selected time range +% +% Here we plot the $I^{*}_{i}$ data for one wavelength only. + +TR = timerange(datestr(TimeRange.calibrationRange(1)),datestr(TimeRange.calibrationRange(2))); +figure( 'Name', 'LopapChannel1 550','NumberTitle','off','Position',[10 100 550 800]); + % : The numbers in the curled brackets give wavelength for which you'd + % like to display the data. The numbers need to be exactly match the + % labels in the previous figure. + ch1s=stackedplot(LopapDataCh1(TR,:),{'550.7'}); + title(LopapParameters.instrumentReferenceWavelength(1)) + % : The numbers in the brackets give the lower and upper limit + % of the y-axis. Adapt these to fit your data and remove the % sign. + %ch1s.AxesProperties.YLimits = [-0.4 -0.1]; + grid ON + % Add title and axis labels + xlabel('Lab Time') + + +figure( 'Name', 'LopapChannel2 550','NumberTitle','off','Position',[600 100 550 800]); + % : The numbers in the curled brackets give wavelength for which you'd + % like to display the data. The numbers need to be exactly match the + % labels in the previous figure + ch2s=stackedplot(LopapDataCh2(TR,:),{'550.5'}); + title(LopapParameters.instrumentReferenceWavelength(2)) + % : The numbers in the brackets give the lower and upper limit + % of the y-axis. Adapt these to fit your data and remove the % sign. + % ch2s.AxesProperties.YLimits = [-0.2 0.3]; + grid ON + % Add title and axis labels + xlabel('Lab Time') + +clearvars -except Lopap* fileLocation myColorMap lopap* TimeRange* +if LopapParameters.Interactive == 1 + distFig + fprintf('press any key if you want to proceed\n') + pause + %% + % + close all +end + +%% Plot selected spectra +% +% Here, we plot a spectrum without ('R1&R2') and with nitrite +% ('R1+nitrite&R2'') for each channel or which ever combination of spectra +% that we have selected. +clc +close all +figure('Name', 'LopapSelectedSpectra') +plot(str2double(LopapSpecCh1.Properties.VariableNames),LopapSpecCh1{datestr(TimeRange.specCompareOne),:},... + 'DisplayName',['#1 ' datestr(TimeRange.specCompareOne) ]) % a time where only R1&R2 was connected to the Lopap +hold on +plot(str2double(LopapSpecCh2.Properties.VariableNames),LopapSpecCh2{datestr(TimeRange.specCompareOne),:},... + 'DisplayName',['#2 ' datestr(TimeRange.specCompareOne)],'LineStyle','--') % a time where only R1&R2 was connected to the Lopap +plot(str2double(LopapSpecCh1.Properties.VariableNames),LopapSpecCh1{datestr(TimeRange.specCompareTwo),:},... + 'DisplayName',['#1 ' datestr(TimeRange.specCompareTwo)]) % a time where nitrite in R1 was connected to the Lopap +plot(str2double(LopapSpecCh2.Properties.VariableNames),LopapSpecCh2{datestr(TimeRange.specCompareTwo),:},... + 'DisplayName',['#2 ' datestr(TimeRange.specCompareTwo)],'LineStyle','--') % a time where nitrite in R1 was connected to the Lopap +legend +box ON +grid ON +% Add title and axis labels +xlabel('wavelength [nm]') +ylabel('Intensity') + +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* +if LopapParameters.Interactive == 1 + distFig + fprintf('press any key if you want to proceed\n') + pause + %% + % + close all +end +%% Process: Absorption spectrum +% +% Here, we reprocess the data to derive an absorption spectrum. This should +% match the spectra shown for nitrite in the Lopap manual. + + +clc +close all +indexR(1) = find(str2double(LopapSpecCh1.Properties.VariableNames) >= LopapParameters.wavelengthR,1); +indexR(2) = find(str2double(LopapSpecCh2.Properties.VariableNames) >= LopapParameters.wavelengthR,1); + +figure('Name', 'LopapAbsorbtionSpectra') +plot(str2double(LopapSpecCh1.Properties.VariableNames),log10(LopapSpecCh1{datestr(TimeRange.specCompareTwo),... + indexR(1)}./LopapSpecCh1{datestr(TimeRange.specCompareTwo),:})-log10(LopapSpecCh1{datestr(TimeRange.specCompareOne),... + indexR(1)}./LopapSpecCh1{datestr(TimeRange.specCompareOne),:}),'DisplayName','#1') +hold on +plot(str2double(LopapSpecCh2.Properties.VariableNames),log10(LopapSpecCh2{datestr(TimeRange.specCompareTwo),... + indexR(2)}./LopapSpecCh2{datestr(TimeRange.specCompareTwo),:})-log10(LopapSpecCh2{datestr(TimeRange.specCompareOne),... + indexR(2)}./LopapSpecCh2{datestr(TimeRange.specCompareOne),:}),'DisplayName','#2','LineStyle','--') +legend +box ON +grid ON +% Add title and axis labels +xlabel('wavelength [nm]') +ylabel('ABS (log I_{r}/I_{i}') +title([datestr(TimeRange.specCompareOne) ' & ' datestr(TimeRange.specCompareTwo)]) +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* +if LopapParameters.Interactive == 1 + distFig + fprintf('press any key if you want to proceed\n') + pause + %% + % + close all +end +%% Process :: Absorption with time +% +% Here, we reprocess the data to get the intensity with time for 1 chosen +% signal wavelength and 1 reference wavelength. The data is identical to that +% shown in the "measurement window" of the Lopap software. +% +% The data is written to the _LopapTimeT_ time table. Rows hold the time, +% columns hold the processed signal for channel 1 and channel 2. +% +% It is important that the same signal wavelength and reference wavelength +% are taken for the calibration and for the experimental data. + +clc +close all +indexS(1) = find(str2double(LopapSpecCh1.Properties.VariableNames) >= LopapParameters.wavelengthS,1); +indexS(2) = find(str2double(LopapSpecCh2.Properties.VariableNames) >= LopapParameters.wavelengthS,1); +indexR(1) = find(str2double(LopapSpecCh1.Properties.VariableNames) >= LopapParameters.wavelengthR,1); +indexR(2) = find(str2double(LopapSpecCh2.Properties.VariableNames) >= LopapParameters.wavelengthR,1); +LopapTimeT = timetable(LopapSpecCh1.Time,log10(LopapSpecCh1{:,indexR(1)}./LopapSpecCh1{:,indexS(1)}),... + log10(LopapSpecCh2{:,indexR(2)}./LopapSpecCh2{:,indexS(2)})); +LopapTimeT.Properties.VariableNames ={'Ch #1','Ch #2'}; +figure('Name', 'Lopap TimeTable') +stackedplot(LopapTimeT,{{'Ch #1','Ch #2'}}); +title(['Lopap' datestr(LopapParameters.MeasurementDate)]) +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* +if LopapParameters.Interactive == 1 + distFig + fprintf('press any key if you want to proceed\n') + pause + %% + % + close all +end +%% Process :: Background +% +% Currently a linear fit forced through 0. +clc +close all +LopapTimeT = fitLopapBkgrd(LopapTimeT,'BkgrdPeriod',TimeRange.backgroundFit); +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* +if LopapParameters.Interactive == 1 + distFig + fprintf('press any key if you want to proceed\n') + pause + %% + % + close all +end +%% +% +clc +close all +% Plot +figure('Name', 'Lopap TimeTable Bkgrd') +s=stackedplot(LopapTimeT,{{'Ch #1 Bkgrd','Ch #2 Bkgrd'}}); +title(['Lopap' datestr(LopapParameters.MeasurementDate)]) +% s.AxesProperties.YLimits = [-0.01 0.6]; +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* + +%% Process :: Calibrate +% +clc +close all + +S = timerange(datestr(TimeRange.MeanZeroAir(1)), datestr(TimeRange.MeanZeroAir(2))); +TZ = LopapTimeT(S,:); +S = timerange(datestr(TimeRange.MeanNitriteStandard(1)), datestr(TimeRange.MeanNitriteStandard(2))); +TN = LopapTimeT(S,:); +calib.mean = mean(TN{:,:},1) - mean(TZ{:,:},1); +calib.std = std(TN{:,:},1); +calib.mixingRatio = LopapParameters.concStandard./1000.*LopapParameters.flowReagents./1000./46.*6.022E23./... + (246E17.*LopapParameters.flowGas).*1E12; + +[LopapTimeT, calib.fitresult] = fitLopapCalibration(LopapTimeT,calib); +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* fitresult calib + +fprintf('Detection limit channel #1 is %.0f ppt \n',calib.fitresult.calibChannel1(3.*calib.std(3))) +fprintf('Detection limit channel #2 is %.0f ppt \n',calib.fitresult.calibChannel2(3.*calib.std(4))) + +%% Plot calibrated data +% +clc +close all +% Plot +figure('Name', 'Lopap TimeTable Calibrated') +s=stackedplot(LopapTimeT,{{'Ch #1 Calib','Ch #2 Calib'}}); +title(['HONO in ppt, Lopap' datestr(LopapParameters.MeasurementDate)]) +% s.AxesProperties.YLimits = [-100 2500]; +clearvars -except Lopap* fileLocation myColorMap LopapParameters TimeRange* fitresult calib +if LopapParameters.Interactive == 1 + distFig + fprintf('press any key if you want to proceed\n') + pause + %% + % + close all +end +%% Save the data +% +save(['LopapMatlabData' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'LopapDataCh1','LopapDataCh2','LopapParameters','calib','LopapSpecCh1','LopapSpecCh2') +save(['Lopap' datestr(LopapParameters.MeasurementDate,'yyyymmdd_HHMM') 'Calib.mat'],'LopapParameters','LopapTimeT','TimeRange') +%% Local functions +function [fileLocation, myColorMap] = localInitializeMatlab +%INITIALIZEMATLAB Sets Lopap specific settings +% Adds processing functions to path and defines file locations for mac, windows, +% and Matlab online (which runs on unix) + +% change working directory +path = matlab.desktop.editor.getActiveFilename; +[filepath,~,~] = fileparts(path); +cd(filepath); + + +% Data Server Address +if ismac + fileLocation.server = '/Volumes/'; + fileLocation.dataPath = 'data_2411/Lopap/'; + % add Lopap functions and general add-ons to path + addpath('/Users/bartels-rausch/SwitchDrive/PSI_Matlab/Add-Ons'); + addpath('/Users/bartels-rausch/SwitchDrive/PSI_Matlab/Lopap/1.LOPAP_Restliche Daten'); +elseif ispc + fileLocation.server = '\\fs00\'; %UNC name for file server that hosts Napp shared drive + fileLocation.dataPath = 'data_2411\Lopap\'; + % add Lopap functions and general add-ons to path + + addpath('C:\Users\vieill_l\switchdrive\Add-Ons'); + + + addpath('C:\Users\vieill_l\switchdrive\Lopap\1.LOPAP_Restliche Daten'); + + +elseif isunix + fileLocation.server = '/'; + fileLocation.dataPath = 'Lopap/'; +end + +% Color maps +myColorMap.red = colormap(autumn(7)); +myColorMap.blue = colormap(winter(7)); + +close all +end + diff --git a/utils/outerjoinRowToTable.m b/utils/outerjoinRowToTable.m index 0f08ad7..1d37b4c 100644 --- a/utils/outerjoinRowToTable.m +++ b/utils/outerjoinRowToTable.m @@ -6,7 +6,7 @@ function Table = outerjoinRowToTable(Table,TableRowToBeAdded) % % * Strings in the table are not yet supported. % * The row that is added needs to be formated as table -% * Adds only one row at a time% +% * Adds only one row at a time. Update 2024: Adds several rows, that is tables % % T = outerjoinRowToTable(TableToAddTo,RowToBeAdded) % @@ -21,11 +21,14 @@ if ~all(ismember(Table.Properties.VariableNames,TableRowToBeAdded.Properties.Var tmpMissingVariable = ~ismember(Table.Properties.VariableNames,TableRowToBeAdded.Properties.VariableNames); tmpVariablesToBeAdded = Table.Properties.VariableNames(tmpMissingVariable); for ii = 1:size(tmpVariablesToBeAdded,2) - if matches(class(Table.(tmpVariablesToBeAdded{ii})),'double') - tmpContentToBeAdded = NaN(1,size(Table.(tmpVariablesToBeAdded{ii}),2)); + if matches(class(Table.(tmpVariablesToBeAdded{ii})),'double') || matches(class(Table.(tmpVariablesToBeAdded{ii})),'int32') + tmpContentToBeAdded = NaN(size(TableRowToBeAdded,1),size(Table.(tmpVariablesToBeAdded{ii}),2)); TableRowToBeAdded = addvars(TableRowToBeAdded,tmpContentToBeAdded,'NewVariableNames',tmpVariablesToBeAdded{ii}); elseif matches(class(Table.(tmpVariablesToBeAdded{ii})),'cell') - tmpContentToBeAdded = {' '}; + tmpContentToBeAdded = repmat({' '},size(TableRowToBeAdded,1),size(Table.(tmpVariablesToBeAdded{ii}),2)); + TableRowToBeAdded = addvars(TableRowToBeAdded,tmpContentToBeAdded,'NewVariableNames',tmpVariablesToBeAdded{ii}); + elseif matches(class(Table.(tmpVariablesToBeAdded{ii})),'string') + tmpContentToBeAdded = repmat({" "},size(TableRowToBeAdded,1),size(Table.(tmpVariablesToBeAdded{ii}),2)); TableRowToBeAdded = addvars(TableRowToBeAdded,tmpContentToBeAdded,'NewVariableNames',tmpVariablesToBeAdded{ii}); end end @@ -35,11 +38,14 @@ if ~all(ismember(TableRowToBeAdded.Properties.VariableNames,Table.Properties.Var tmpMissingVariable = ~ismember(TableRowToBeAdded.Properties.VariableNames,Table.Properties.VariableNames); tmpVariablesToBeAdded = TableRowToBeAdded.Properties.VariableNames(tmpMissingVariable); for ii = 1:size(tmpVariablesToBeAdded,2) - if matches(class(TableRowToBeAdded.(tmpVariablesToBeAdded{ii})),'double') + if matches(class(TableRowToBeAdded.(tmpVariablesToBeAdded{ii})),'double') || matches(class(TableRowToBeAdded.(tmpVariablesToBeAdded{ii})),'int32') tmpContentToBeAdded = NaN(size(Table,1),size(TableRowToBeAdded.(tmpVariablesToBeAdded{ii}),2)); Table = addvars(Table,tmpContentToBeAdded,'NewVariableNames',tmpVariablesToBeAdded{ii}); elseif matches(class(TableRowToBeAdded.(tmpVariablesToBeAdded{ii})),'cell') - tmpContentToBeAdded = repmat({' '},size(Table,1),1); + tmpContentToBeAdded = repmat({' '},size(Table,1),size(TableRowToBeAdded.(tmpVariablesToBeAdded{ii}),2)); + Table = addvars(Table,tmpContentToBeAdded,'NewVariableNames',tmpVariablesToBeAdded{ii}); + elseif matches(class(TableRowToBeAdded.(tmpVariablesToBeAdded{ii})),'string') + tmpContentToBeAdded = repmat({" "},size(Table,1),size(TableRowToBeAdded.(tmpVariablesToBeAdded{ii}),2)); Table = addvars(Table,tmpContentToBeAdded,'NewVariableNames',tmpVariablesToBeAdded{ii}); end end