% Symbolic integration of transient 1-d diffusion impulse response 
% to give finite square pulse response for test injection in Y tube
% performed 5/15/2013 
% Requires Symbolic Math toolbox
% mez 5/19/2013
% MEZ 6/14/2013 include error propagation for measured points
err=sqrt(.07^2 + .1^2 + .17^2 + .25^2 + .05^2);
% leak cal tag err, leak checker cal err,  
% parameter list
% Q =  injection rate [torr*m^3/s]
% W =  duration (width) of injection in [s]
% R =  radius of beam tube [m]
% m =  molecular weight [AMU]
% T =  temperature [K] 
% v =  mean molecular velocity [m/s]
% c =  diffusion constant for capillary of radius R [m^2/s]
load('testleak.mat', 'rate1');
%
syms p Q W R T v m c z t;
assume(t>0);
assume(c>0);
assume(z>0);
pi=sym('pi');
% impulse response of tube 
p = Q*(1/(2*pi^(3/2)*R^2))*(1/sqrt(c*t))*exp(-z^2/(4*c*t));
% indefinite integral on time, ignore singularities & unphys branches 
pint = int(p, t, 'IgnoreAnalyticConstraints', true);
% same as above delayed by W
pintw = heaviside(t-W) * subs(pint, t, abs(t-W));
% difference, i.e., rectangular pulse for duration W
prect = pint - pintw ;
%
% load numerical substitutions and plot scaling factor
%sc = 1; % plot scale factor (ez plot routine doesn't like small nos.)
% Q = sc * 1.15e-4 * 1e-3 * sqrt(28/4);  % injection flow [1/sc torr*m^3/s]
Q = 5.1e-6 * 0.759 * 1e-3 * (1-.032*.5);  % injection flow [torr*m^3/s]
% .759 Tl/s per atm-cc/s, 1e-3 m^3 per liter, 3.2% per year 
% age depletion, 6 months since factory calibration
W = 26;         % duration of injection in [s]
R = 0.6 ;       % radius of beam tube [m]
% R=.6;
m = 4   ;       % molecular weight [AMU]
% T = 294;      % temperature [K]
T = 310;        % bte was hot on May 15th 2013
v = 1197*sqrt(T/273)*sqrt(4/m); % mean molecular velocity [m/s]
c = (2/3)*R*v ; % diffusion constant for capillary radius R [m^2/s]
% cond = sqrt(28/4) * 12.1 * (1.8)^3 / 100 ; % conductance of leak check hose for He (l/s)
cond = (34/2) * sqrt(29/4); % conductance of angle valve pair into turbopump feeding LD
%
%
figure(1);
% offset = 253; % actual distance from test injection to sample point (m)
cases = 1:250:2001;
% cases = 250;
%%
for(k=1:length(cases));
    %%
    offset = cases(k);
    %%
sliced=subs(W*p);
sliced=subs(sliced, z, offset);
slicer=subs(prect, z, offset);
slicer=subs(slicer);
tt = 2:2:10000;
poft = subs(slicer, t, tt);
% cond = sqrt(28/4) * 12.1 * (1.8)^3 / 100 ; 
% conductance of leak check hose for He (l/s)
cond = (34/2) * sqrt(29/4); % conductance of valve pair into turbopump feeding LD
semilogy(tt,cond*poft);
% semilogy(tt,cond*poft, tt, (1+err)*cond*poft, '-.',  tt, (1-err)*cond*poft , '-.') ;
% plot(tt,cond*poft);
% loglog(tt,poft);
% legend(['Model, z = ', num2str(offset), '  m']);
hold all;
end;
% hold;
%
% % semilogy(rate1(:,1), (1/cond)* rate1(:,2), 'mo');
% plot(rate1(:,1)+10, rate1(:,2), 'mo');
% % semilogy(rate1(:,1), rate1(:,2), 'mo');
% % loglog(rate1(:,1), (1/cond)* rate1(:,2), 'mo');
ylim([1e-10  3e-8]);
% % ylim([0  1e-8]);
% % xlim([0 6000]);
% xlim([-100 1900]);
% title({'Diffusion model vs. measured points, beamtube cal leak He injection test' 'MEZ, 15 May 2013'} );
ylabel('Apparent Leak Rate (Tl/s)');
xlabel('t (sec)');
% hold off;
% legend(' Model', ' Measurement');
% legend(' Measurement');
legend(num2str(cases'));
% legend('Meas');
%
hold off;
%%  OLD STUFF
% ezplot(sliced, [1,6000]);
% hold on;
% ezplot(slicer, [1, 6000]);
% title('pressure at sampling port');
% % ezplot(log10(slicer),3);
% hold off;
% %
% %
% % plot impulse response scaled to Q*W total injection
% figure(1);
% % ezmesh(log10(subs(p)), [0 10 -100 100], 240);
% ezmeshc(W*subs(p), [0 150 0 300], 120);
% title(['P(t,z) for '...
%     num2str(W*Q*1e3/sc, '%0.1e') ' torr-liter injection at t=0 '], ...
%     'FontSize', 14);
% zlabel(['P(' num2str(1./sc, '%0.2g') ' torr)'] );
% %
% figure(2);
% %
% plot continuous injection of duration W at rate Q
% figure(1);
% ezmeshc(subs(prect),[0 6000 0 120 ], 120 );
% title(['P(t,z)  for ' ...
%     num2str(W,'%0.2g') ' sec injection at '...
%     num2str(Q*1e3/sc, '%0.1e') ' torr-liter/sec rate'], ...
%     'FontSize', 14 );
% zlabel(['P (' num2str(1./sc, '%0.2g') ' torr)'] );
% hold off
%
% figure(3);
% ezmesh(subs(pintw), [0 60 0 100], 60);
% hold on;
% ezmesh(subs(pint),[0 60 0 100], 60);
% hold off
%


