% Analysis of MID strip from LLO MY 10/16/2014
% MEZ 
%
% open and parse Pfeiffer Quadera ascii export file
RGAfile = '20141012 134518 Zucker_MID_beamtube_10_Oct_14_final.asc';
fid = fopen(RGAfile);
RGAscanstart = textscan(fid, '%*s %*s %s %s %s', 1,  'Headerlines', 3, 'CollectOutput', 1);
% disp(['start time ='  , data{1}(1,:)])
frewind(fid); % restart from top of file
RGAdata = textscan(fid, ...
    '%s %s %s %f %f %s %s %s %f %f %s %s %s %f %f %s %s %s %f %f %s %s %s %f %f %s %s %s %f %f %s %s %s %f %f', ...
    'HeaderLines', 8, 'CollectOutput', 1);
fclose(fid);
RGAstopstr = '16-Oct-2014 09:12:16';  % manual date in matlab date format
RGAstopnum = datenum(RGAstopstr); 
disp(['Scan ending ' , RGAstopstr] );

% RGA scan sequence: AMU 2 14 18 28 32 40 5
AMU2 = RGAdata{2};
AMU14 = RGAdata{4};
AMU18 = RGAdata{6};
AMU28 = RGAdata{8};
AMU32 = RGAdata{10};
AMU40 = RGAdata{12};
AMU5 = RGAdata{14};
%
lenamu = length(AMU5)-1;
AMUSUM=AMU5(1:lenamu,:); % establish array size
AMUSUM(:,2)= AMU2(1:lenamu,2)+ AMU28(1:lenamu,2)+ AMU40(1:lenamu,2)+ AMU18(1:lenamu,2); % omit AMU 14 & 32
%
volume = 1e3*3947*pi*0.62^2;  % tube volume in liters
filterSize = 300; % length of smoothing filter
trendLength = 40; % length of history to fit for trend, hours
RGAscale = 4.16e+04; % torr/ampere estimated from gauge mean and total ion current
speriod = mean(diff(AMU5(:,1))); % measure average sample period
nsamp = round(trendLength * 3600/speriod); % no. samples to go back
maxsamp = length(AMU5(:,1))-1 ; % no. samples in original series
% filt1 = gausswin(filterSize);
% filt1 = hann(filterSize);
filt1 = ones(filterSize); 
filt1 = filt1/sum(filt1); % normalize filter
AMU5(:,2) = filter(filt1,1,AMU5(:,2)); % lowpass null channel
%
tAMU5 = AMU5(maxsamp-nsamp:maxsamp,:); % snip section, size trendLength
tAMU2 = AMU2(maxsamp-nsamp:maxsamp,:);
tAMU14 = AMU14(maxsamp-nsamp:maxsamp,:);
tAMU18 = AMU18(maxsamp-nsamp:maxsamp,:);
tAMU28 = AMU28(maxsamp-nsamp:maxsamp,:);
tAMU32 = AMU32(maxsamp-nsamp:maxsamp,:);
tAMU40 = AMU40(maxsamp-nsamp:maxsamp,:);
%
tAMU2(:,2)=abs(tAMU2(:,2)-tAMU5(:,2)); % remove preamp offset by subtracting smoothed null channel 
tAMU14(:,2)=abs(tAMU14(:,2)-tAMU5(:,2));
tAMU18(:,2)=abs(tAMU18(:,2)-tAMU5(:,2));
tAMU28(:,2)=abs(tAMU28(:,2)-tAMU5(:,2));
tAMU32(:,2)=abs(tAMU32(:,2)-tAMU5(:,2));
tAMU40(:,2)=abs(tAMU40(:,2)-tAMU5(:,2));
%
tAMUSUM = tAMU5; % establish array size
tAMUSUM(:,2)=tAMU2(:,2)+tAMU28(:,2)+tAMU40(:,2)+tAMU18(:,2); % omit AMU 14 & 32
% =======================================
%
%
% Simple linear fit
f=fittype('a*x+b');
[model28,gof28]=fit(tAMU28(:,1),tAMU28(:,2),f,'Startpoint',[1  1]);
[model2,gof2]=fit(tAMU2(:,1),tAMU2(:,2),f,'Startpoint',[1  1]);
[models,gofs]=fit(tAMUSUM(:,1),tAMUSUM(:,2),f,'Startpoint',[1  1]);
% 
LR28 = volume * RGAscale * model28.a ; % 
LR2 = volume * RGAscale * model2.a ; %
LRs = volume * RGAscale * models.a;
disp(' ');
disp('Simple regression: ');
disp('--------------------------------------------------------------');
disp(['AMU 28 slope = ', num2str(LR28), ' Tl/s', '  (R squared = ', num2str(gof28.adjrsquare), ')' ]);
disp(['AMU 2 slope = ', num2str(LR2), ' Tl/s', '  (R squared = ', num2str(gof2.adjrsquare), ')']);
disp(['SUM slope = ', num2str(LRs), ' Tl/s', '  (R squared = ', num2str(gofs.adjrsquare), ')']);
areacm = pi*124*3947*100;
ionxsfactor=2.4; % correction for H2 current to N2 equivalent
JH2 = ionxsfactor*LR2/areacm;
Tavg = mean(TempMean(length(TempMean)-3600*trendLength:length(TempMean)));
disp(['Implied H2 desorption = ', num2str(JH2), ' Tl/s/cm^2 at <T> = ', num2str(Tavg), ' C']); 
disp(' ');
%    
    figure(3);
    plot(tAMU28(:,1),tAMU28(:,2),tAMU28(:,1),model28.a*tAMU28(:,1)+model28.b);
    hold all; 
    grid on; 
    plot(tAMU2(:,1),tAMU2(:,2),tAMU2(:,1),model2.a*tAMU2(:,1)+model2.b);
    %ylim([0,4e-7]);
    xlabel('t ');
    ylabel('P ');
    legend('AMU28','28 Fit',' AMU 2 ','2 Fit', 'Location', 'SouthEast');
    %ylim([1e-9 1e-6]);
    ylim(7e-12*[0 1]);
    % xlim([0 8.7e4]);
    xlabel('time (sec)');
    ylabel(['ion current (A)'] );
    title({'MY RGA'  ['mez ' , datestr(now)]});
    hold off;
%
% =======================================% %
% Following assumes pressure & temp data in
% workspace already, w/ same value for trendLength
% 
toffset = round(daterun - RGAstopnum) * 3600 * 24; % subtract matlab dates, conv to seconds
Ptime = max(tAMU5(:,1)) - 3600 * trendLength +1 : max(tAMU5(:,1));
Ptime = Ptime + toffset;
% =======================================% %
    figure(4);
    semilogy(Ptime,trendData);
    % plot(Ptime,trendData);
    hold all;
    plot(tAMUSUM(:,1),(tAMUSUM(:,2)*RGAscale), '-.');
    plot(tAMU2(:,1),tAMU2(:,2)*RGAscale);
    % plot(tAMU14(:,1),tAMU14(:,2)*RGAscale);
    plot(tAMU18(:,1),tAMU18(:,2)*RGAscale);
    plot(tAMU28(:,1),tAMU28(:,2)*RGAscale);
    % plot(tAMU32(:,1),tAMU32(:,2)*RGAscale);
    plot(tAMU40(:,1),tAMU40(:,2)*RGAscale);
    legend('CC Gauge','RGA sum','2 ','18', '28','40', 'Location', 'SouthEast');
    ylim([1e-9 1e-6]);
    % ylim(1e-7*[0 4]);
    % xlim([0 8.7e4]);
    xlabel('time (sec)');
    ylabel(['P (torr) or ion current (A) x ', num2str(RGAscale) , ' Torr/A'] );
    title({'MY RGA'  ['mez ' , datestr(now)]});
    grid;
    hold off;
% 
%
disp('Linear model with tempco: ');
disp('--------------------------------------------------------------');
% snip a temperature trace corresponding to the RGA measurement period
% 
tempDelay2 = round(3.1*3600); % apparent delay between T and effect on AMU28?
tslong = timeseries(TempMean, max(tAMU28(:,1))+toffset-length(TempMean)+tempDelay2:max(tAMU28(:,1))+toffset-1+tempDelay2);
tsshort = resample(tslong, tAMU28(1:length(tAMU28(:,1))-1 , 1) ) ; % temperature resampled at RGA sample times
tAMU28new = tAMU28(1:length(tAMU28(:,1))-1,:);
PAMU28 = tAMU28new;
PAMU28(:,2) = PAMU28(:,2)*RGAscale;
t=PAMU28(:,1);
tempIn = tsshort.data;
meantemp = mean(tempIn);
pptemp = max(tempIn) - min(tempIn);
tempIn = tempIn-meantemp; % remove mean
X = [t, tempIn];
mdl1 = fitlm(X, PAMU28(:,2), 'linear');
u = mdl1.Coefficients.Estimate;
% disp(['mean = ', num2str(u(1)), ' Torr']);
disp(['AMU28 slope = ', num2str(u(2)*volume), ' Tl/s']);
disp(['tempco = ', num2str(u(3)), ' T/K']);
disp(['R squared = ', num2str(mdl1.Rsquared.Adjusted)] );
disp(['at mean temp = ', num2str(meantemp), ' C, excursions ', num2str(pptemp), ' C p-p']);
disp(['using delay = ', num2str(tempDelay2/3600), ' h']);
%
    figure(5);
    plot(t/3600,PAMU28(:,2) ...
    , t/3600, RGAscale * (model28.b + model28.a*t) ...
    , t/3600, u(1) + u(2)*t+u(3)*tempIn  ...
    , t/3600, u(1) + u(2)*t ...
    );      % , tAMUSUM(:,1)/3600, tAMUSUM(:,2)*RGAscale ...
            % , 1/3600 * ( t(1):t(1)+length(trendData)-1 ), trendData );
    hold all;
    title({'Fit to AMU28 dP/dt'; ['mez', datestr(now)]})
    xlabel('time (h)');
    ylabel('P (torr)');
    grid on; 
    legend('28 AMU', 'Linear', 'Linear + Tempco', 'Revised' ...
        );   % , 'RGA sum', 'CC Gauge' );
    hold off;
%
%
    figure(6);
    Ptimelong = max(AMUSUM(:,1)) - length(PressureMean) + 1 : max(AMUSUM(:,1));
    Ptimelong = Ptimelong + toffset;
        plot(AMUSUM(:,1)/3600, RGAscale*AMUSUM(:,2)...
        , AMU28(:,1)/3600, RGAscale*AMU28(:,2)...
        , AMU2(:,1)/3600, RGAscale*AMU2(:,2)...
        , Ptimelong/3600,PressureMean, '--' );
        hold all;
        grid on;
        xlim([10 75])
        ylim(1e-7*[0 4])
        xlabel('time (h)');
        ylabel('P (torr)');
        legend('RGA Sum', '28 AMU', '2 AMU', 'Cold Cathode', 'Location', 'SouthEast');
        title({'CC vs. RGA over 70 hours'; ['mez', datestr(now)]})
        %
        hold off;
        
        
 
    






