%=====================================================
%=====================================================
% Steady-state pressure profile in beamtube
% and residual gas index phase noise 
% 
% mez 25 Apr 2015
%=====================================================
%=====================================================
%
% geometry
L = 3994.5;                 % length (m)
Reff = .597 ;               % effective 'conductance' tube radius incl. effects of baffles (m)
R = .625 ;                  % physical tube radius (m)
V = pi*R^2*L;               % volume (m^3)
A = 2*pi*R*L;               % area (m^2)
%=====================================================
% gas parameters
% %==hydrogen=============
% M = 2;                    % 
% indexatm = 1.0+132e-6;    % H2 index @ atmospheric pressure 293K
% %
% %==nitrogen=============
% % 
% M = 28;                   % molecular mass (AMU)
% indexatm = 1.0+274e-6;    % N2 index @ atmospheric pressure 293K CHECK!!
% %
%==water=============
M = 18;                     % molecular mass (AMU)
indexatm = 1.0+261e-6;      % H20 index @ atmospheric pressure 293K ** CHECK!!**
%=====================================================
%
T = 293 ;                   % temperature [K]
AMU = 1.66053892e-27 ;      % [kg]
kB = 1.3806488e-23 ;        % Boltzmann [m^2 kg s^-2 K^-1] 
rhoNatm = 101325/(kB*T);    % number density at 1 atm = 101325 Pa 
alpha = (3/(4*pi*rhoNatm))*(indexatm^2-1)/(indexatm^2+2);   % molecular polarizability
v0 = sqrt(2*kB*T/(M*AMU));     % most probable molecular speed
vm = sqrt(8*kB*T/(pi*M*AMU));  % mean molecular speed
%% =====================================================
% optics
lambda = 1.064e-6 ;         % wavelength [m]
R1 = 1934;                  % ITM ROC [m]
R2 = 2245;                  % ETM ROC [m] 
g1 = 1-L/R1;                
g2 = 1-L/R2;
zR = L * sqrt( g1*g2*(1-g1*g2)/(g1+g2-2*g1*g2)^2 );     % Rayleigh range [m]
z1 = L * g2*(1-g1)/(g1+g2-2*g1*g2) ;    % waist location from ITM [m]
z2 = L * g1*(1-g2)/(g1+g2-2*g1*g2) ;    % waist location from ETM [m]
w0 = sqrt(lambda*zR/pi);                % waist radius [m]
%% =====================================================
% conductance of beamtube 
% CL = 1e-3*12.1*sqrt(T/293)*sqrt(28/M)*(2*Reff*100)^3/(100*L) ; % [m^3/s]
CL = (pi/12)*v0*(2*Reff)^3/L ;     % O'Hanlon 3rd ed., eq'n 3.22 (2003) [m^3/s]
tau = V/(pi*CL);                   % time constant [s]
% surface outgassing flux
% JM = 1.5e-15 * 1e-3 * 1e4;    % orig. spec [Tl/s/cm^2]*1e-3m^3/l*1e4cm^2/m^2
% JM = 1.5e-14 * 1e-3 * 1e4;    % 10x spec
JM = 1.3e-15 * 1e-3 * 1e4 ;     % T1400723, 11/14/2014 accumulation
% QM = pi*R^2*L * JM;           % YOU IDIOT!!! SERIOUSLY???
QM = JM*A;                      % OK that's better
%% =====================================================
% computation parameters
N = 81;                 % number of nodes
dV = V/(N-1) ;          % volume element
y = 0:L/(N-1):L ;       % position vector
MaxIter = 1e6 ;         % max no. iterations
tol = .005 ;            % acceptable convergence tolerance  
ncases = 5;             % no. pump arrangements to consider
%% =====================================================
% initial conditions
P = 1e-8 * ones(ncases,N);  % pressure init
Pdot = zeros(ncases,N);     % time derivative init
Pu = Pdot;              % analytic init
O = QM/(N-1)*ones(ncases,N); % outgasflux per node [T*m^3*s^-1] 
C = CL*(N-1);           % conductance of segment [m^3*s^-1]
%% =====================================================
% define pump array
FIP = - 1e-3 * 1/(1/600 + 1/(12.1*20^3/20));  % 1200LX speed, derated for tubulation
FSHORT = -85 ;          % short cryo, 8.5e4 l/s = 85 m^3/s
FLONG = -85 ;           % long cryo, 1e5 l/s = 100 m^3/s
F = zeros(ncases,N); 	% init, no pumping
F(:,1) = FLONG ;           
F(:,N) = FSHORT ;           
% 
F(2,round(N*2018.778/L)) = FIP;     % MY
%
F(3,round(N*2018.778/L)) = FIP;     % MY
F(3,round(N*1703.687/L)) = FIP;     % Y1-7 
% %
F(4,round(N*2330.813/L)) = FIP; 	% Y2-1 
F(4,round(N*2018.778/L)) = FIP;     % MY
F(4,round(N*1703.687/L)) = FIP;     % Y1-7
F(4,round(N*1465.943/L)) = FIP;     % Y1-6
% %
% F(5,round(N*2573.023/L)) = FIP;     % Y2-2 
% F(5,round(N*2330.813/L)) = FIP; 	% Y2-1 
% F(5,round(N*2018.778/L)) = FIP;     % MY
% F(5,round(N*1703.687/L)) = FIP;     % Y1-7 
% F(5,round(N*1465.943/L)) = FIP;     % Y1-6 
% % 
% % eqi-spaced pumps
pp = 15 ;              % no. intermediate pumps between vertex and end
FPORTS = 1+round((N-1)/(pp+1)):round((N-1)/(pp+1)):(pp+1)*round((N-1)/(pp+1));
F(5,FPORTS) = FIP ;    % 
%% =====================================================
tic
% calculate steady-state pressure profile for each pumping case
for j = 1:ncases
    % relax pressure to equilibrium
    for k = 1:MaxIter ;
        % net gas flux into each node
        Pdot(j,1) = (1/dV)*( C*(P(j,2)-P(j,1)) + P(j,1)*F(j,1) + O(1) );
        Pdot(j,N) = (1/dV)*( C*(P(j,N-1)-P(j,N)) + P(j,N)*F(j,N) + O(N) );
            for m = 2:N-1 ;
                Pdot(j,m) = (1/dV)*( C*((P(j,m-1)-P(j,m)) + (P(j,m+1)-P(j,m))) + P(j,m)*F(j,m) + O(m) );
            end
        % exit if pressure is stationary everywhere
            if max(abs(Pdot(j,:)./P(j,:))) < tol/tau;
                break ;
            end
        % if still changing, continue to iterate P
            for m = 1:N ;
                P(j,m) = P(j,m) + Pdot(j,m); 
            end
    end
    if k < MaxIter
        disp(['case', num2str(j),': ', num2str(k), ' iterations to ', num2str(100*tol), '% residual']);
        disp(['out/in = ', num2str(-sum((F(j,:).*P(j,:)))/QM)]); % check mass conservation
        disp(['peak P = ', num2str(max(P(j,:)))]);
        disp(['mean P = ', num2str(mean(P(j,:)))]);
        % disp(['(end-pumped peak & mean = ', num2str(max(Pu)),'   ', num2str(mean(Pu)),')']);
    else
        disp('failed to converge');
    end
end
%
toc
%% =====================================================
% 
figure(1);
for j=1:ncases;
    plot(y, P(j,:)); hold all;
end
hold off; 
ylim(1e-9*[0 .5]);
title({'Pressure Distribution',['H_{2}O Flux = ', num2str(JM/10), ' Tl/s/cm^{2}']} );
ylabel('P (Torr)');
xlabel('y (m)');
legend('none','MY','MY, Y1-7','Y1-6, Y1-7, MY, Y2-1','All');
%
%% =====================================================
% calculate phase noise 
rhoN = zeros(ncases,N);
S_y = zeros(ncases,1);
ytilde = zeros(ncases,1);
htilde = zeros(ncases,1);
% 
rhoN = (P./760)*rhoNatm;                    % number density vs. y (each case)
w = w0 * sqrt( 1+(y-z1).*(y-z1)./zR^2 ) ;   % beam radius vs. y
sc = max(w)/max(max(rhoN));                 % plot scale factor
figure(2);
for j=1:ncases;
    plot(y,sc*rhoN(j,:)./w); hold all;
end
plot(y,w./w0, 'o');
title('Phase Noise Power Integrand');
ylabel('A.U.');
xlabel('y (m)');
legend('none','MY','MY, Y1-7','Y1-6, Y1-7, MY, Y2-1','all','(beam size)');
hold off;
% DC approximation f << v0/w
for j=1:ncases;
    S_y(j)=(4*(2*pi*alpha)^2/v0)*trapz(y,rhoN(j,:)./w); % LIGO-P940008 eq.(1); rho moved inside 
end;
ytilde = sqrt(S_y) ;
htilde = ytilde/L ;
%
figure(3);
for j=1:ncases;
    plot(j,htilde(j,1), 'O'); hold all;
end
xlim([.7 ncases+.3]);
set(gca,'XTick', 1:ncases);
ylim(1e-25*[0 3.5]);
title('Strain Noise (1 arm) vs. Pumping Arrangement');
ylabel('Strain (Hz^{-1/2})');
xlabel('Case');
legend('none','MY','MY,Y1-7','Y1-6,Y1-7,MY,Y2-1','all');
hold off;
%
%% =====================================================
figure(4);
for j=1:ncases;
    plot(j,htilde(j,1)/htilde(1,1), 'x'); hold all;
end
xlim([.7 ncases+.3]);
set(gca,'XTick', 1:ncases);
ylim([0 1.05]);
title('Normalized noise');
ylabel('Fraction');
xlabel('Case');
legend('none','MY','MY,Y1-7','Y1-6,Y1-7,MY,Y2-1','all');
hold off;

%% =====================================================
% frequency dependence
f = logspace(1, 4, 50) ; 
Sff = zeros(ncases,length(f));
figure(5);
for j = 1:ncases ; % select which case to do frequency dependence
    for k = 1:length(f);
         Sff(j,k) = (4*(2*pi*alpha)^2/v0) * trapz(y,(rhoN(j,:)./w).*exp((-2*pi*f(k)/v0)*w)); 
         loglog(f,sqrt(Sff(j,:))/L); hold all;
    end
end
title('Strain Noise (1 arm) vs. Frequency');
ylabel('Strain (Hz^{-1/2})');
xlabel('f (Hz)');
legend('none','MY','MY,Y1-7','Y1-6,Y1-7,MY,Y2-1','all');
hold off
%% =====================================================
% tabulate numbers
ctable = ones(ncases,5);
ctable(:,1)= 1:ncases ;                 % case number
ctable(:,2)= 1e10*mean(P,2);            % mean pressure
ctable(:,3)= 1e10*max(P');              % peak pressure
ctable(:,4)= 1e25*htilde(:);            % strain noise 
ctable(:,5)= htilde(:)/htilde(1);       % normalized strain noise
ctable
%% =====================================================
% % % analytic form, pressure distribution, end-pumped (MEZ slides, February 1995)
% figure(1);
% hold all;
% Pu = (QM/(2*CL))*(y/L-y.^2/L^2) + QM/sum(F(1,:));
% plot(y,Pu, 'mo');
% % 
% %





