%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Written by Joel Brito - 17/Aug/2014
%Modifications by Milena Ponczek - jun 2019 %
% %
%This routine reads the PTR tdms files and plots daily aerosol conc. %
%values for a quick look. It saves a .mat with read data. %
% %
%Latest changes: %
%To improve processing time now the routine reads the previous mat %
%file (if existing) and imports only updated ascii files. The same is%
%valid for making new plots %
% %
%Milena:
% %
% Bug at BG routine: We need to remove the PI and impurities from BG %
% matrix %
% %
% Added
% Primary , H2O cluster and impurities analysis and plots %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%

clear
close all
clc

%Choose the folder we want
[tmp,mydir]=fileparts(uigetdir(pwd));
%%%[tmp,mydir]=fileparts(pwd);

%Choose the title of your series of experiment
station=['RioBranco-noBG'];
fig_dir = 'fig/';
mkdir(fig_dir);

is_BG=0;

if is_BG==1
fl_dir='PTRMS_BG/';
else
fl_dir = strcat(mydir,'/');
end

disp(['looking for files: ' fl_dir])

%{
if exist('mat-files/PTRMS.mat')>0
load mat-files/PTRMS.mat
count_rem=0;
clear rem_idx
fl=dir([fl_dir,'*.tdms']);

for i=1:max(size(fl_old))
for j=1:max(size(fl))
if strcmp(fl(j).name,fl_old(i).name)
if fl(j).bytes==fl_old(i).bytes
count_rem=count_rem+1;
rem_idx(count_rem)=j;
end;
end;
end;
end;

if count_rem>0
fl(rem_idx)=[];
fl_old=[fl_old;fl];
end;

else
%}
%fl=dir([fl_dir,'*.tdms']);

fl=dir(fl_dir);
dirFlags = [fl.isdir];
dirFlags(1:2)=0; %força os dois pimeiros diretorios . e .. a terem a dirFlag = 0
fl = fl(dirFlags);
count=0;
count_old=1;
cps=[];
DO=[]; %DO means Digital Output
fl_old=fl;
time_cps=[];
time_BG=[];
time_Cal=[];
Cal=[];
BG=[];
%end;
%%
%--------------------------------------------------------------------
% Read new files
%--------------------------------------------------------------------

error_flag=1;
%Check if there are files in the choosen folder
if size(fl,1)==0 && count==0
fl
disp('No files found')
return

elseif min(size(fl))>0
%look for .tdms files
for fl_number=1:max(size(fl))

filenames=dir([fl_dir,fl(fl_number).name,'/*.tdms']);

for filenames_number=1:max(size(filenames))

%tests if the filename.mat exist already (the fc exist returns 2)
if exist([fl_dir,fl(fl_number).name,'/',filenames(filenames_number).name(1:end-4),'mat'])==2
disp(['The mat file exist. Loading file: ',fl_dir,fl(fl_number).name,'/',filenames(filenames_number).name])
load([fl_dir,fl(fl_number).name,'/',filenames(filenames_number).name(1:end-4),'mat']);
else
%disp(['converting file: ',fl_dir,fl(fl_number).name,'/',filenames(filenames_number).name])
tic
[ConvertedData,ConvertVer,ChanNames]=convertTDMS(false,[fl_dir,fl(fl_number).name,'/',filenames(filenames_number).name]);
toc
save([fl_dir,fl(fl_number).name,'/',filenames(filenames_number).name(1:end-4),'mat'],'ConvertedData','ConvertVer','ChanNames');
end

%@Milena comments:
% ChanNames faz parte da estrutura formada pelo 'convertTDMS', é
% uma celula (cell array) que tem os canais usados pelo PTRMS
% (sempre 166)
% check if (???) file is valid and has data 'ChanNames' should
% be greater than 115
if count==0 && size(ChanNames{1,1},1)>115
header_cps=ConvertedData.Data.MeasuredData(1,3).Data';
cps=ones(1,max(size(header_cps)));
DO=ones(1,8);
error_flag=0; %Changes to zero, first file OK
end

clear dummy
clear dummy_DO

if max(size(ConvertedData.Data.MeasuredData(1,3).Data))==max(size(header_cps))
check_DO=0;

for j=1:max(size(ConvertedData.Data.MeasuredData))
%verifica onde dentro de 'ConvertedData.Data.MeasuredData'
%está a string 'Time / Cycle/Abs Time [sec]', no caso
%pode-se ver que sempre quando j=52
if strcmp(ConvertedData.Data.MeasuredData(1,j).Name,'Time / Cycle/Abs Time [sec]')
j;
%quando j = 52, temos 'Time / Cycle/Abs Time
%[sec]'. Os valores são convertidos na operação:
time_dummy = (ConvertedData.Data.MeasuredData(1,j).Data./(24*60*60)+ datenum(1904,1,1,0,0,0));
%porque datenum(1904,1,1,0,0,0) 1904???
end

% Misc/X1 DO tem informacao do digital output (DO)
% usado para controlar as valvulas que fazem a
% aquisicao do BG
if max(size(ConvertedData.Data.MeasuredData(1,j).Name))==12
dummy_name=ConvertedData.Data.MeasuredData(1,j).Name;
if strcmp(dummy_name(1:11),'Misc/X1 DO ')
dummy_size=ConvertedData.Data.MeasuredData(1,j).Total_Samples;
%dummy_channel recebe o numero apos 'Misc/X1 DO '
dummy_channel=str2num(dummy_name(12));
dummy_DO(:,dummy_channel)=ConvertedData.Data.MeasuredData(1,j).Data(1:dummy_size);
check_DO=1;
end
end
if max(size(ConvertedData.Data.MeasuredData(1,j).Name))>24
if strcmp(ConvertedData.Data.MeasuredData(1,j).Name(1:24),'Raw signal intensities/m')
for k=1:max(size(header_cps))
dummy_header=num2str(header_cps(k));
if strcmp(ConvertedData.Data.MeasuredData(1,j).Name(25:25+max(size(dummy_header))-1),dummy_header)
dummy(:,k)=ConvertedData.Data.MeasuredData(1,j).Data;
end
end
end
end
end

count=count+size(cps,1);
cps = [cps ; dummy];
DO = [DO ; dummy_DO];
time_cps = [time_cps ; time_dummy];

end
end
end
end

%--------------------------------------------------------------------

disp('Removing duplicate data...')
[time_cps,idx]=unique(time_cps);
cps=cps(idx,:);
DO=DO(idx,:);

cps(time_cps==0,:)=[];
DO(time_cps==0,:)=[];
time_cps(time_cps==0)=[];

%remove first 5 lines
cps(1:5,:)=[];
DO(1:5,:)=[];
time_cps(1:5,:)=[];

if exist('mat-files', 'dir') ~= 7
% tests whether the 'mat-files' directory already exists,
% if not, create the folder
disp(strcat('mat-files',' was created'))
mkdir mat-files/
end

if is_BG==1
save mat-files/PTRMS_BG.mat
else

filename_mat_file = strcat('mat-files/',station, '.mat');
save (filename_mat_file);

end

%%
%--------------------------------------------------------------------
% Processing files
%--------------------------------------------------------------------

%@Milena: why this?????
Readme_PTRMS='Data from PTR-MS. Data filtered for OFR cycles.';

label_21='Primary ion (cps)';
label_28='m/z 28 ';
label_30='NO+ (cps)';
label_31='Formaldehyde ';
label_32='O2+ (cps)';
label_33='Methanol ';
label_37='H3O+(H2O) (cps)';
label_42='Acetonitrile ';
label_45='Acetaldehyde ';
label_46='m/z 46 ';
label_47='Ethanol/Formic Acid';
label_57='m/z 57 ';
label_59='Acetone ';
label_60='m/z 60 ';
label_61='Acetic Acid ';
label_63='DMS ';
label_69='Isoprene ';
label_71='MVK + MACR ';
label_73='MEK ';
label_77='PAN ';
label_79='Benzene ';
label_81='Monoterpenes ';
label_83='eighty-three ';
label_89='m/z 89 ';
label_93='Toluene ';
label_95='m/z 95 ';
label_105='m/z 105 ';
label_107='C8 Aromatics ';
label_121='m/z 121 ';
label_125='one-two-five ';
label_129='C9 Aromatics ';
label_137='Monoterpenes ';
label_205='Sesquiterpenes ';

vocs=cps;
time_PTRMS=time_cps;

header=header_cps;
header_BG=header_cps;
header_Cal=header_cps;

%DO indicates the valves control

%for this, I understand if
%DO1 and DO4 and DO5 = 1, it represents a BG measurement
time_BG=time_PTRMS(DO(:,1)==1 & DO(:,4)==1 & DO(:,6)~=1);
BG=vocs(DO(:,1)==1 & DO(:,4)==1 & DO(:,6)~=1,:);

%for this, I understand if
%DO1 and DO2 and DO3 and DO5 = 1, it represents a Calibration measurement
% time_Cal=time_PTRMS(DO(:,1)==1 & DO(:,2)==1 & DO(:,3)==1 & DO(:,5)==1);
% Cal=vocs(DO(:,1)==1 & DO(:,2)==1 & DO(:,3)==1 & DO(:,5)==1,:);

%remove data when DO(1) = 1
%if comented, no BG removing!
time_PTRMS(DO(:,1)==1)=[];
vocs(DO(:,1)==1,:)=[];

k=1;
%Identifying the Primary Ion mz=21 -> Prim
for i=1:size(header,2);
if header(i)==21
Prim=vocs(:,i);
if mean(Prim)<1e6
Prim=Prim.*500;
end;
idx_remove(k)=i;
k=k+1;
%Identifying the mz=25
elseif header(i)==25
idx_remove(k)=i;
k=k+1;
%Identifying the mz=30 -> NO+ (impurity)
elseif header(i)==30
m30=vocs(:,i);
idx_remove(k)=i;
k=k+1;
%Identifying the mz=34 -> O2+ (impurity)
elseif header(i)==32
m32=vocs(:,i);
idx_remove(k)=i;
k=k+1;
%Identifying the mz=37 -> H30+(H20) (water cluster)
elseif header(i)==37
m37=vocs(:,i);
idx_remove(k)=i;
k=k+1;
end;

% idx_remove represents the header of removed ions (21, 25, 30, 34, 37)
end;

%Remove the m/z of the primary and impurities ions that do not need to have
%mixing ratios calculated.
if k>1
vocs(:,idx_remove)=[];
header(idx_remove)=[];
% Bug at BG routine: We need to remove the PI and impurities from BG
% matrix as well, otherwise the BG routine does not work properly.
% If not running the BG, comment the 2 lines below
BG(:,idx_remove)=[];
header_BG(idx_remove)=[];
end;
%%
%----------------------------------------------------------------
%Primary , H2O cluster and impurities analysis
%Added by Milena
%----------------------------------------------------------------

ratioCluster = (m37./Prim) *100;
ratioO2 = ((m32)./Prim) *100;
ratioNO = (m30./Prim) *100;
ratioImpurities = ratioO2 + ratioNO;

meanCluster = mean(ratioCluster);
meanImpurities = mean(ratioImpurities);

% plot(time_PTRMS,Prim,'r')
% % hold on
% % plot(time_PTRMS,m37),'b')

%@Milena: a pensar: colocar a datavalidation aqui!!!

%Correct for fluctuations in primary
%@Milena: para corrigir, ele multiplica o valor atual de VOC_cps pela media (PI) e
%divide pelo valor de PI atual
for i=1:size(vocs,1)
vocs(i,1:end)=vocs(i,1:end).*mean(Prim)./Prim(i);
end;

%%
%Retrieve VOC index
for i=1:max(size(header))
eval(['Idx',num2str(header(i)),'=i;']);
eval(['title_',num2str(header(i)),'=label_',num2str(header(i)),';']);
end;

title_PTRMS='VOCs ';
label_PTRMS='Mixing ratio (ppbv)';

if is_BG==1
return
end
%%

%--------------------------------------------------------
% % %Here it loads and applies the BG corrections
% BG_PTRMS
%Ignore BG for now

% % load PTRMS\BG_dry.mat
% % load PTRMS\BG_wet.mat
% %
% %
% % idx=find(time_PTRMS<datenum(2018,09,26,14,0,0));
% % vocs(idx,:)=vocs(idx,:)-BG_wet;
% %
% % idx=find(time_PTRMS>datenum(2018,09,26,14,0,0));
% % vocs(idx,:)=vocs(idx,:)-BG_dry;
% % idx=1;
%%
%--------------------------------------------------------
%Here it loads and applies the Cal corrections
Cal_PTRMS

%%
filename_mat_file = strcat('mat-files/',station, '.mat');
save (filename_mat_file);

%%
% % %Remove data above 50ppbv
% % [max_val,idx_time]=max(vocs);
% % [max_max_val,idx_voc]=max(max_val);
% %
% % % [M,I] = max(___) M represent the max values of each colum of A and I find the indices of the maximum values of A and returns
% % % them in output vector I.
% %
% % remove_50ppb_counter=0;
% % while max_max_val>50
% % time_PTRMS(idx_time)=[];
% % vocs(idx_time,:)=[];
% % remove_50ppb_counter=remove_50ppb_counter+1;
% %
% % [max_val,idx_time]=max(vocs);
% % [max_max_val,idx_voc]=max(max_val);
% %
% % end;

%%
%--------------------------------------------------------------------
% Plot section below
%--------------------------------------------------------------------

%%
%--------------------------------------------------------------------
% Plot full time series
%--------------------------------------------------------------------

disp('Plot: time series')

if isunix
fig1 = figure('visible','off');
else
fig1=figure;
end;

set(fig1,'InvertHardcopy','on');
set(gca, 'FontSize', 12, 'LineWidth', 2);

plot(time_PTRMS,vocs(:,Idx45),'Color',[0.8 0 0.2],'LineWidth', 1.2)
hold on
%plot(time_PTRMS,vocs(:,Idx42),'Color', [0.3216 0.1882 0.1882],'LineWidth', 1.2)
plot(time_PTRMS,vocs(:,Idx59),'Color', [0.929, 0.694, 0.125],'LineWidth', 1.2)
plot(time_PTRMS,vocs(:,Idx69),'Color', [0 0.4980 0],'LineWidth', 1.2)
plot(time_PTRMS,vocs(:,Idx71),'Color',[0.75, 0.75, 0],'LineWidth', 1.2)
plot(time_PTRMS,vocs(:,Idx73),'Color',[1 0.3216 0.4392],'LineWidth', 1.2)
plot(time_PTRMS,vocs(:,Idx79),'Color',[0.4118 0.3255 0.2392],'LineWidth', 1.2)
plot(time_PTRMS,vocs(:,Idx81),'Color', [0.2000 0.4196 0.6667],'LineWidth', 1.2)
plot(time_PTRMS,vocs(:,Idx93),'Color',[0.25, 0.25, 0.25],'LineWidth', 1.2)
title([title_PTRMS,station])
legend('Acetaldehyde','Acetone','Isoprene','MVK+MACR','MEK','Benzene', 'Monoterpenes', 'Toluene')
xlabel('Date')
ylabel(label_PTRMS)

%ylim(axes1,[0 nanmax(nanmax([vocs(idx_st:idx_end,Idx45) vocs(idx_st:idx_end,Idx69) vocs(idx_st:idx_end,Idx59)]))])
ylim([0 15])
box on
dynamicDateTicks([], [], 'dd/mm');

set(gcf,'units','inch','position',[0,0,8,4]); %
set(gcf,'PaperPositionMode','auto') %salva a figura na dimensao que estamos vendo na tela

%nome=['../RioBranco/fig/PTRMS_' mydir '_Time_series'];
nome=[fig_dir, 'Time_series'];

print(fig1,'-dpng',[nome,'.png']);
hold off

%%
%--------------------------------------------------
%plot mz21 corrected (Prim) and mz37
%--------------------------------------------------
fig1=figure;
plot(time_cps,cps(:,1),'m') %mz21 raw
hold on
plot(time_cps,(cps(:,1)*500),'r') %mz21 x 500 -> primary ion
plot(time_cps,cps(:,4),'b') %mz37 -> water cluster
title(['Prim and Water Cluster'])
legend('mz21','Primary',' H_2O Cluster')
xlabel('Date')
ylabel('Intensity (cps)')

%ylim(axes1,[0 nanmax(nanmax([vocs(idx_st:idx_end,Idx45) vocs(idx_st:idx_end,Idx43) vocs(idx_st:idx_end,Idx59)]))])
box on
dynamicDateTicks([], [], 'dd/mm');

set(gcf,'units','inch','position',[0,0,8,4]); %
set(gcf,'PaperPositionMode','auto')

%nome=['../RioBranco/fig/PTRMS_' mydir '_Time_series'];
nome2=[fig_dir, 'Time_series_Primary'];

%print(fig1,'-dpng',[nome2,'.png']);

%%
%--------------------------------------------------
%plot Prim and impurities analysis
%--------------------------------------------------
fig1=figure;

plot(time_PTRMS,ratioCluster,'c')
hold on
plot(time_PTRMS,ratioO2,'b')
plot(time_PTRMS,ratioNO,'r')
plot(time_PTRMS,ratioImpurities,'k')

title(['H_2O Cluster and Impurities'])
legend('H_2O Cluster','O_2^+','NO^+', 'Total Impurity')
xlabel('Date')
ylabel('Ratio')

%ylim(axes1,[0 nanmax(nanmax([vocs(idx_st:idx_end,Idx45) vocs(idx_st:idx_end,Idx43) vocs(idx_st:idx_end,Idx59)]))])
box on
dynamicDateTicks([], [], 'dd/mm');

set(gcf,'units','inch','position',[0,0,8,4]); %
set(gcf,'PaperPositionMode','auto')

nome3 = [fig_dir, 'Time_series_Impurities_ratio'];

print(fig1,'-dpng',[nome3 '.png']);

%%
%--------------------------------------------------------------------
% Plot 1-day, for all days
%--------------------------------------------------------------------

days_PTRMS=unique(floor(time_PTRMS));
% floor(x) = rounds each element of x to the nearest integer less than or equal to that element.

%if you don't want to consider all days, stablish here:
%days_PTRMS_OK =

count_days=0;
if exist('days_PTRMS_OK','var')==1
if min(size(days_PTRMS_OK))>0
days_PTRMS_OK=unique(days_PTRMS_OK);
% for i=1:max(size(days_PTRMS))
% if max(days_PTRMS(i)==days_PTRMS_OK)>0
% count_days=count_days+1;
% rem_days(count_days)=i;
% end;
% end;
days_PTRMS(rem_days)=[];
end;
else
days_PTRMS_OK=[];
end;
%%
if min(size(days_PTRMS))>0
for i=1:max(size(days_PTRMS))

fig_name = [fig_dir 'PTRMS_',station,'_',datestr(days_PTRMS(i),29)];

clear fig1;

quick_time_PTRMS=(time_PTRMS-days_PTRMS(i)).*24;

[diff_st,idx_st]=min(abs(quick_time_PTRMS));
[diff_end,idx_end]=min(abs(quick_time_PTRMS - 24));

if diff_st<1/24 && diff_end<1/24
days_PTRMS_OK=[days_PTRMS_OK days_PTRMS];
end;

%Calculate sunrise/sunset
%hmjb - we need a way to automatize this
% coordinates São Paulo - USP:
% Latitude: -23.561815 | Longitude: -46.734703

%coordinates Rio Branco:
% Latitude:-9.97 | Longitude:-67.81

[rs,t,d,z,b,r]=suncycle(-9.97,-67.81,i);

%hmjb
% bug from Joel?? he was adding +8h... maybe because he
% passed longitute positive (manaus shoould be -60 instead
% of +60?)
%rs=rs+8;

if idx_end>idx_st+30

if nanmax(nanmax([vocs(idx_st:idx_end,Idx45) vocs(idx_st:idx_end,Idx59) vocs(idx_st:idx_end,Idx71) vocs(idx_st:idx_end,Idx79)]))>1
if isunix
fig1 = figure('visible','off');
else
fig1=figure;
end
set(fig1,'InvertHardcopy','on');
set(gca,'FontSize', 12, 'LineWidth', 2);

axis('off');
axes1 = axes('Parent',fig1,...
'XTickLabel',{'0','2','4','6','8','10','12','14','16','18','20','22','24'},...
'XTick',[0 2 4 6 8 10 12 14 16 18 20 22 24]);

% axes1 = axes('Parent',fig1)

hold on
rectangle('Parent',axes1,'Position',[rs(1),-1000,rs(2)-rs(1),200000],'FaceColor',[1 1 0.7],'EdgeColor',[1 1 0.7])
%
% plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx45),'m')
% %plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx42),'b')
% plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx59),'c')
% plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx69),'g')
% plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx71),'r')
% plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx79),'k')
% plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx81),'b')
% plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx93),'y')
% legend('Acetaldehyde','Acetone','Isoprene','MVK+MACR','Benzene', 'Monoterpenes', 'Toluene')

plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx45),'Color',[0.8 0 0.2],'LineWidth', 1.2)
%plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx42),'b')
plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx59),'Color', [0.929, 0.694, 0.125],'LineWidth', 1.2)
plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx69),'Color', [0 0.4980 0],'LineWidth', 1.2)
plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx71),'Color',[0.75, 0.75, 0],'LineWidth', 1.2)
plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx73),'Color', [1 0.3216 0.4392],'LineWidth', 1.2)
plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx79),'Color',[0.4118 0.3255 0.2392],'LineWidth', 1.2)
plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx81),'Color',[0.2 0.4196 0.6667],'LineWidth', 1.2)
plot(quick_time_PTRMS(idx_st:idx_end),vocs(idx_st:idx_end,Idx93),'Color',[0.25, 0.25, 0.25],'LineWidth', 1.2)
legend('Acetaldehyde','Acetone','Isoprene','MVK+MACR','MEK','Benzene', 'Monoterpenes', 'Toluene')

%ylim(axes1,[0 nanmax(nanmax([vocs(idx_st:idx_end,Idx45) vocs(idx_st:idx_end,Idx71) vocs(idx_st:idx_end,Idx79)]))])
%colocando o ylim fixo (em 10 ppb)
ylim([0 12])

xlim(axes1,[0 24])

title([title_PTRMS,station,' - ',datestr(days_PTRMS(i),1)])
xlabel('Time (UTC)')
ylabel(label_PTRMS)
box on

set(gcf,'units','inch','position',[0,0,8,4]); %
set(gcf,'PaperPositionMode','auto')

print(fig1,'-dpng',[fig_name,'.png']);

end
end;
end;
end;

%%
% Data validation filters (for further processing)
%

%%
%Creating a mat file with the processed data
% %
% % label = {label_31; label_33; label_42; label_45; label_47; label_59; label_69; label_71; label_73; label_79; label_93; label_107; label_81; label_137; label_205;label_129};
% % file_name_mat = strcat(station,'.mat');
% % save (file_name_mat, 'time_PTRMS', 'vocs' , 'label' );

%%
%CHECK WHAT TO WRITE IN THE ASCII FILE!!

% % Results_time = time_PTRMS;
% % Results_date = datestr(Results_time,0); %converts the time from number to date format dd-mmm-yyyy HH:MM:SS
% % %choose which fragments we want to save in the ASCII file
% % mz_69 = vocs(:,Idx69);
% % mz_71 = vocs(:,Idx71);
% % mz_79 = vocs(:,Idx79);
% % mz_93 = vocs(:,Idx93);
% % Results_VOCs = vocs(:,:);
% %
% % Results_Table = table(Results_time,Results_date,mz_69, mz_71, mz_79, mz_93);
% % %writetable(Results_Table,'myData.txt','Delimiter','\t')
% % file_name = strcat(station,'.txt');
% % writetable(Results_Table,file_name,'Delimiter','\t');
% % %teste
% % Results_Table = table(Results_time,Results_date,Results_VOCs);
% % file_name2 = strcat(station,'_2','.txt');
% % writetable(Results_Table,file_name2,'Delimiter','\t');
% %
% % if isunix
% %
% % %CHECK WHAT TO WRITE IN THE ASCII FILE!!
% % %{
% % fid = fopen('0_Ascii-files/PTRMS.csv','wt');
% % fprintf(fid,'Date(UTC), Number concentration (ug m-3)\n');
% % for i=1:max(size(time_PTRMS_avg))
% % fprintf(fid,'%s, %2.2f\n',datestr(time_PTRMS_avg(i)),CPC_avg(i));
% % end;
% % fclose(fid);
% % %}
% % %!chmod 777 ../RioBranco/0_Ascii-files/PTRMS.csv
% % !rm -f ../RioBranco/mat-files/PTRMS.mat
% % end;
% %
% %
% % close all
% % % save ../RioBranco/mat-files/PTRMS.mat
% %
% % if isunix
% % !chmod 777 ../RioBranco/mat-files/PTRMS.mat
% % end