% Symbolic integration of transient 1-d diffusion impulse response 
% to give step and finite square pulse response
% Requires Symbolic Math toolbox
% mez 7/28/2012
%
% 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]
%
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 = 1e8; % 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]
W = 60;         % duration of injection in [s]
R = 0.6 ;       % radius of beam tube [m]
m = 4   ;       % molecular weight [AMU]
T = 300;        % temperature [K] 
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]
%
% 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
ezmeshc(subs(prect),[0 150 0 300], 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
%
% figure(3);
% offset = 10;
% sliced=subs(W*p);
% sliced=subs(sliced, z, offset);
% slicer=subs(prect, z, offset);
% slicer=subs(slicer);
% ezplot(sliced, [1,200]);
% hold on;
% ezplot(slicer, [1, 200]);
% % ezplot(log10(slicer),3);
% hold off;
% %




