% Script

clear all,close all,cd bin,path(path,pwd),cd ..

% USER DEFINED PARAMETERS
% Specification of computational grid
Rmax = 100;
M = 300; % Number of radii
N = 300; % Number of angles
gamma_grid = 0.7; % Skewness of grid

% Specification af problem
Peclet = [ 0.0001 0.001 0.1 1 10 100 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 14000 18000 20000]/3;
gamma_turb = 1; 		% 1 Corresponds to pure molecular diffusion
Normalize = 'BC';	% 'BC' or 'Flow'
BC = 'Dirichlet';		% 'Dirichlet' or 'Neumann'

% Specification of flow files
FlowFileNames = { ...
      'FlowRe0R500M200N200' , ...
      'FlowRe1e-006R500M200N200' , ...
      'FlowRe0-0001R500M200N200' , ...
      'FlowRe0-001R500M200N200' , ...
      'FlowRe0-01R500M200N200' , ...
      'FlowRe0-1R500M200N200' , ...
      'FlowRe1R500M200N200' , ...
      'FlowRe2R500M200N200' , ...
      'FlowRe3R500M200N200' , ...
      'FlowRe4R500M200N200' , ... 
      'FlowRe5R500M200N200' , ...
      'FlowRe6R500M200N200' , ...
      'FlowRe7R500M200N200' , ...
      'FlowRe8R500M200N200' , ...
      'FlowRe9R500M200N200' , ...
      'FlowRe10R500M200N200' , ...
      'FlowRe14R500M200N200' , ...
		'FlowRe18R500M200N200' , ...
		'FlowRe20R500M200N200' , ...
	      
   };

% Output is stored in files
ConcFileNames = { ...
      'hCRe0Pe0-000033Rmax100M300N300DirBC', ...
      'hCRe1e-006Pe0-00033Rmax100M300N300DirBC', ...
      'hCRe0-0001Pe0-033Rmax100M300N300DirBC', ...
      'hCRe0-001Pe0.33Rmax100M300N300DirBC', ...
      'hCRe0-01Pe3.3Rmax100M300N300DirBC', ...
      'hCRe0-1Pe33Rmax100M300N300DirBC', ...
      'hCRe1Pe330Rmax100M300N300DirBC', ...
      'hCRe2Pe667Rmax100M300N300DirBC', ...
      'hCRe3Pe1000Rmax100M300N300DirBC', ...
      'hCRe4Pe1333000Rmax100M300N300DirBC', ...
      'hCRe5Pe1667Rmax100M300N300DirBC', ...
      'hCRe6Pe2000Rmax100M300N300DirBC', ...
      'hCRe7Pe2333Rmax100M300N300DirBC', ...
      'hCRe8Pe2667Rmax100M300N300DirBC', ...
      'hCRe9Pe3000Rmax100M300N300DirBC', ...
      'hCRe10Pe3333Rmax100M300N300DirBC', ...
      'hCRe14Pe4667Rmax100M300N300DirBC', ...
      'hCRe18Pe6000Rmax100M300N300DirBC', ...
      'hCRe20Pe6667Rmax100M300N300DirBC', ...
                  
      
   };

% No more user defined input below this
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

options.Method = '3-4 order upwind';

cd Userdata

for i = 1:length(FlowFileNames),
   if length(Peclet)>1,
      Pe = Peclet(i);
   else
      Pe = Peclet;
   end
   
   if isa(Normalize,'cell'),
      No = Normalize{i};
   else
      No = Normalize;
   end
   
   if isa(BC,'cell'),
      BCi = BC{i};
   else
      BCi = BC;
   end
   
   if length(M)>1,
      Mi = M(i);
   else
      Mi = M;
   end
   
   if length(Rmax)>1,
      Rm = Rmax(i);
   else
      Rm = Rmax;
   end
   
   if length(N)>1,
      Ni = N(i);
   else
      Ni = N;
   end
   
   if length(gamma_turb)>1,
      gt = gamma_turb(i);
   else
      gt = gamma_turb;
   end
   
   cd Flows;
   load(FlowFileNames{i});
   
   clear psi
   
   cd ../Concentrations;
      
   lg = loggrid(Rm,Mi,Ni,'Logarithmic grid',gamma_grid);
   
   options.BC = BCi;
   options.Peclet = Pe;
   options.Gamma = gt;
   
   ff = interp(ff,lg);
   TM = TensorMatrix(ff,options);
   
   nf = SolveEquation(TM,lg);
   
   if strcmpi(Normalize,'flow')
      nf = nf * (1/Flux(nf,ff,Pe,gt));
      Flow = 1;
   else
      Flow = Flux(nf,ff,Pe,gt);
   end
   
   save(ConcFileNames{i},'nf','Pe','Flow')
   cd ..
end

cd ..

return