Print[Style["Evacuate Virgo vessel","Section"]]
L = 1000; (* 1000 m long vessel *)
L = 100 L; (* use centimeters *)
Diam = 1.2; (* 1.2 m diameter *)
Diam = 100 Diam; (* use centimeters *)
Ntotal= 10; (* elements in length *)
dt = 3600; (* length of time step (default 1 hour) *)
Ntime = 6; (* time steps *)
Ts = 300; (* surface temperature *)
M = 18; (* atomic number of gas molecules *)
Pinitial = 1 (10^(-9)); (* initial pressure in mbar *)
N0 = 1 (10^15); (* molecules per cm2 needed to form monolayer *)
f = 0.6; (* sticking coefficient *)
tFrenkel = 1 (10^(-13)); (* Frenkel time *)
Ed = 78500; (* desorption energy Joule/mole *)
R0 = 8.314; (* Joule/(K mole) *)
T = 300; (* gas temperature *)

ts = tFrenkel Exp[Ed /(R0 Ts)]; (* time spent on surface [ s ] *)
b =  3.51 (10^22) f tFrenkel Exp[Ed /(R0 Ts)]/(N0 Sqrt[M T]);
Theta = b Pinitial / (1 + b Pinitial); (* monolayer coverage *)
Print["Time spent on surface: ", ts, "  [ s ]"];
Print["Surface covered by: ", Theta, "  monolayer at start"];

(* calculate element values *)
Ci= 3.81 Sqrt[ T / M ] Diam^3 Ntotal/ L; (* conductance for water per element *)
Vi=  1 (10^(-3))Pi Diam^2 L /(4 Ntotal ); (* volume of each element *)
Qi = 1 (10^(-12)) Pi Diam L /(Ntotal ); (* desorption of stainless steel *)
Q0 = 1 (10^(-15)); (* gas load from mirror vessel *)

(* Declare memory for vectors *)
Pvector = Range[ Ntotal ];
PvectorOld = Range[ Ntotal ];
Cvector = Range[ Ntotal ];
Svector = Range[ Ntotal ];
Qvector = Range[ Ntotal ];
QvectorOld = Range[ Ntotal ];
Vvector = Range[ Ntotal ];
Zerovector = Range[ Ntotal ];
PressureList = Range[ Ntime + 1 ];

(* Fill vectors *)
Do[Cvector[[i]] = Ci,{i,1,Ntotal}];
Do[Vvector[[i]] = Vi,{i,1,Ntotal}];
Do[Pvector[[i]] =Pinitial,{i,1,Ntotal}];
Do[PvectorOld[[i]] = Pinitial,{i,1,Ntotal}];
Do[Svector[[i]] = 0,{i,1,Ntotal}];
Svector[[Ntotal]] = 1000;

(* Create desorption vectors *)
Do[Qvector[[i]] = Qi + Vi Pinitial / dt,{i,1,Ntotal}];
Qvector[[1]] = Qvector[[1]] + Q0;

Print["pressure distribution [ mbar ]: ", Pvector]
Print["desorption distribution [ mbar l/s ]: ", Qvector]

(* Create matrix and calculate inverse *)
Do[Zerovector[[i]] = 0,{i,1,Ntotal}];
Smatrix = DiagonalMatrix[ Zerovector];
Smatrix[[1,1]] = Svector[[1]] + Cvector[[1]] + Vi / dt;
Do[Smatrix[[i,i-1]] = -Cvector[[i-1]],{i,2,Ntotal}];
Do[Smatrix[[i,i]] = Svector[[i]] + Cvector[[i-1]] + Cvector[[i]] + Vi / dt ,{i,2,Ntotal}];
Do[Smatrix[[i,i+1]] = -Cvector[[i]],{i,1,Ntotal -1}];

PressureList[[0]] =Pvector
For[k=1, k < Ntime+1, k++,{Pvector = LinearSolve[ Smatrix, Qvector ];
Print["item: ",k,"    pressure distribution [ mbar ]: ", Pvector];
Do[QvectorOld[[i]] = Qvector[[i]],{i,1,Ntotal}];
Do[PvectorOld[[i]] = Pvector[[i]],{i,1,Ntotal}];
Do[Qvector[[i]] = Qi + Vi PvectorOld[[i]] / dt,{i,1,Ntotal}];
Qvector[[1]] = Qvector[[1]] + Q0;
Print["desorption distribution [ mbar l/s ]: ", Qvector];
PressureList[[k+1]] =Pvector;
}]
PressureList
g0 = ListPlot[PressureList[[0]], PlotRange -> All, PlotJoined -> True]
g1 = ListPlot[PressureList[[Ntime/2]], PlotRange -> All, PlotJoined -> True]
g2 = ListPlot[PressureList[[Ntime]], PlotRange -> All, PlotJoined -> True]
gl=Show[g0,g1,g2]