
% leakmodel.m
% MEZ 10/14/2014
% takes in vacuum data from .mat file generated by get_pressures_lite.m
% using gps start time defined in workspace 
% (e.g., from recent execution of get_pressures_lite)
% 
% filename=['/data/zucker/pressure',num2str(gpsstart),'.mat'];
filename=['./pressure',num2str(gpsstart),'.mat'];
load(filename);
disp(['time series from GPS  ' , num2str(gpsstart), ' to ', num2str(gpsstart + duration)]);
disp(['thus ending ', datestr(daterun)]);
disp(['(span of ', num2str(duration/3600), ' hours)']);
%
volume = 4.73e6;  % tube volume in liters
filterSize = 1000; % length of smoothing filter
trendLength = 95; % length of history to fit for trend, hours
trendEnd = 7.1778 ; % hours ago to stop 
tempDelay = 3600*0 ; % empirical delay time for thermal effect
% tempDelay = 0; 
    %
    PT723 = press(1).data; % EY
    PT990 = press(2).data; % Y2-6
    PT984 = press(3).data; % Y2-4
    PT753 = press(4).data; % MY
    PT970 = press(5).data; % Y1-6
    PT964 = press(6).data; % Y1-4
    PT958 = press(7).data; % Y1-2
    PT624 = press(8).data; % LY
    %
    IP758 = capdat(1).data; % cap ion pump
    FL755 = capdat(4).data; % cap turbo foreline
    CP756 = capdat(5).data; % cap Inficon
    IP757 = capdat(6).data; % ion pump Inficon
    %
    TY27 = tempdat(1).data; % Y2-7 degrees C
    TY24 = tempdat(2).data; % Y2-4 degrees C 
    TY21 = tempdat(3).data; % Y2-1 degrees C
    TY17 = tempdat(4).data; % Y1-7 degrees C 
    TY14 = tempdat(5).data; % Y1-4 degrees C 
    TY11 = tempdat(6).data; % Y1-1 degrees C 
    %
%PressureMean = (PT723+PT990+PT984+PT753+PT970+PT964+PT958+PT624)/8;
%PressureMean = (PT723+PT990+PT753+PT970+PT964+PT958+PT624)/7; % PT984 flaky
PressureMean = (PT990+PT753+PT984+PT624)/4; % other gauges flaky or have data dropouts
TempMean = (TY27+TY24+TY21+TY17)/4; % average tube temp. Y11 and Y14 died
% TempMean = TY21;
%      figure(1);
 PTime = 1:size(PT753);
 PTime = PTime - duration;
%      % plot(PTime/3600,PressureMean*1e7,PTime/3600, TempMean);
%      %plot(PTime/3600,PressureMean);
plot(PTime/3600,TempMean);
%     % %
%
% tt=PTime/3600;
% plot(tt,PT990,tt,PT753,tt,PT984,tt,PT624,tt,PressureMean);
% figure;
%
%
% filt1 = gausswin(filterSize);
% filt1 = hann(filterSize);
filt1 = ones(filterSize); 
filt1 = filt1/sum(filt1); % normalize filter
filt2 = filt1; 
%
PressureMean = filter(filt1,1,PressureMean); % lowpass pressure
TempMean = filter(filt1,1,TempMean); % lowpass temperature
%
trendData= PressureMean(size(PressureMean)-round(trendEnd*3600)-trendLength*3600:size(PressureMean)-round(trendEnd*3600)-1); % snip tail
% trendData= PT753(size(PT753)-trendLength*3600:size(PT753)-1); % snip tail
%
disp(' ');
disp(['Over ', num2str(max(size(trendData))/3600), ' hours:']);
disp(' ');
disp('Linear Fit');
disp('================================');
%
t=-size(trendData)-3600*trendEnd:-3600*trendEnd-1;
f=fittype('a*x+b');
[model,gof]=fit(t',trendData,f,'Startpoint',[1  1]);
    figure(2);
    plot(t/3600,trendData,t/3600,model.a*t+model.b);
    grid on; 
    hold all;
    ylim([0,5e-8]);
    xlabel('t (h)');
    ylabel('P (Torr)');
    LR = model.a * volume ; % torr/sec * liters
    disp(['slope = ', num2str(LR), ' Tl/s']);
    disp(['R squared = ', num2str(gof.adjrsquare)] );
    disp(' ');
    disp(' ');
%
disp('Linear Fit with Temperature Regression');
disp('================================');
%
    tempInput = TempMean(size(TempMean)-trendLength*3600-tempDelay:size(TempMean)-1-tempDelay); % snip tail
    tempInput = tempInput - mean(tempInput); % remove mean temperature so no DC term gets fitted (??)
    X = [t' tempInput];
    mdl1 = LinearModel.fit(X, trendData, 'linear');
    u = mdl1.Coefficients.Estimate;
    % disp(['mean = ', num2str(u(1)), ' Torr']);
    disp(['slope = ', num2str(u(2)*volume), ' Tl/s']);
    disp(['tempco = ', num2str(u(3)), ' T/K']);
    disp(['R squared = ', num2str(mdl1.Rsquared.Adjusted)] );
    disp(' ');
    disp(' ');
% figure(3);
plot(t/3600,(u(1)+u(2)*t+u(3)*tempInput'), t/3600, (u(1)+u(2)*t)); grid on; 
legend('Gauge', 'Linear', 'Linear + Tempco', 'Revised', 'Location', 'NorthWest');
% 
% % mdl1 =  LinearModel.fit(1:size(PressureMean1), volume * PressureMean1, 'linear'); % linear regression
% % disp(['previous hour: ', num2str(mdl1.Coefficients.Estimate(2)), ...
% %     ' +/- ' , num2str(mdl1.Coefficients.SE(2)), ...
% %     ' Tl/s']) % report slope and standard error
% % 
hold off;
