% simpcoag.m
% program to provide simple coagualation calculation
% simplified version of llayer3, coord with simpbeta.m
% work with constant shear, no disaggregation
% copyright George A. Jackson, Texas A&M Univ. 2001
% Input:  keryboard and a file of betas computed by simpbeta.m, default 

clear;%   start with a clean slate

cprog   ='simpcoag.m';                  % this file
cbeta   =  'betaf';                     % default input file with betas
cxrate  = @simpcoagr;

global beta25p beta1p    nsect linmat   % feed to ratechange for ode45  
global setcon	                          % use for set rate calc   

dum         =input('Accept the defaults parameters[cr] or customize [1]\n ');
linp        =~isempty(dum);

                % constants
                % the following are defaults
rt      = 3e-14;                        % relative and absolute tolerances 
at      = 1e-16;
cbeta   = 'betaf';                      % default name of beta file	
tfinal  = 20.;	                        % start and end times
alpha   = 1;                            % set the stickiness to 0.8
Lz      = 65;                           % thickness of the layers(m)
Lgam    = 0.1;                          % average shear, 1/s
Lmu     = 0.15;                         % specific growth rate in first section
num_1   = 40;                           % number of particles/cm^3 in first section
nsecgro = 4;                            % section at which start growth in aggreg.                    

if linp
	display('to accept the default value, simply press the cr');   
	
	dum =input(['Input name of file with betas  [',cbeta,']\n cbeta= '],'s');
	if ~isempty(dum),cbeta=dum;end;
		
	dum =input(['Input final time [',num2str(tfinal),']\n tfinal= ']);
	if ~isempty(dum),tfinal=dum;end;
	
	dum =input(['Input stickiness [',num2str(alpha),']\n alpha= ']);
	if ~isempty(dum),alpha=dum;end;
	
	dum =input(['Input layer thickness onelayer(m) [',num2str(Lz),']',...
         ' \n zthick= ']);
	if ~isempty(dum),Lz=dum;end;
	
	dum =input(['Input shear (1/s) [',num2str(Lgam),']\n',...
         ' Lgam= ']);
	if ~isempty(dum),Lgam=dum;end;
  
  dum = input(['Input specific growth rate (1/d) [',num2str(Lmu),']\n',...
      ' Lmu= ']);
  if ~isempty(dum),Lmu=dum;end;
  
  dum =input(['Input initial number of monomers (#/cm^3) [',num2str(num_1),...
      ']\n num_1= ']);
  if ~isempty(dum),num_1=dum;end;
  
  dum = input(['Input the section at which allow growth in aggregates.\n',...
      '  Set to 0 if no growth in aggregates.\n Default is ',num2str(nsecgro)',...
      '\n nsecgro = ']);
  if ~isempty(dum)
    if dum>(nsect-1) | dum<0 then
      disp(' nsecgro outside allow allowable range. using default');
    else, nsecgro=dum;
    end
  end
  
	dum =input(['Input relative tolerance [',num2str(rt),']\n rt= ']);
	if ~isempty(dum),rt=dum;end;
	
	dum =input(['Input absolut tolerance [',num2str(at),']\n at= ']);
	if ~isempty(dum),at=dum;end;
end

                % load the good stuff
load(cbeta)                             % data in structures: 
                                        %    betasd, physcon, speccon,+ betaheader

nsect   = speccon.nsect;			    
setcon  = physcon.setcon;


tinitial=0.;
daysec  =8.64e4;		                    % seconds/d  % convert betas to day units

                %  need to calc the initial distributions, given by Vcon
                %  assume that have equal mass in each section at first,
Vcon    =ones(nsect,1)*1.5*speccon.Vlow(1);
Vcon(2:end)=Vcon(2:end)/10;        % decrease the initial concentration of others
Vcon    = num_1 * Vcon;     % use the initialization data

calcomp =1:nsect;

                % setup shear kernels
beta25sn = betasd.beta25s*daysec*alpha;      % change units to days, bring in alpha
beta1sn  = betasd.beta1s *daysec*alpha;

			    % load differential sedimentation kernels
beta25o = betasd.beta25d*daysec*alpha;	
beta1o  = betasd.beta1d *daysec*alpha;

                % add the diffsed, shear terms 
beta25p = beta25o + beta25sn*Lgam;
beta1p  = beta1o  + beta1sn *Lgam;

setrat  = speccon.vset / 100*physcon.daysec;	  % convert to m/d

tempset = +diag(setrat)/Lz;		                  % make a matrix for settling, 1/d
                                                % now build the  matrices

                                                % do growth in first section
grodiag=zeros(nsect,1);                          % set the diagonal terms first
grosdiag=zeros(nsect-1,1);                       % set up the subdiagonal term
grodiag(1,1)=1;                                 % the growth term for solitary
if nsecgro>0
  grodiag(nsecgro:(nsect-1))=-1;                     % lose by growth
  grosdiag(nsecgro:end)=+2;
end
gromat=diag(grodiag) + diag(grosdiag,-1);
temp2   = Lmu*gromat;                           % growth


linmat  = temp2 -tempset;	                      % linear term: growth+set

			                                          % try the new ode15s
atmat   = (at * 1.5 .^ (-(calcomp-1)))';        % calc abs tol
atvec   = zeros(1,nsect);
atvec(:) = atmat; 
options = odeset('reltol',rt,'refine',0,'abstol',atvec);
tspan   = tinitial:(tfinal-1);
ntspan  = length(tspan);

%   try to set up a ss spectrum at first
Vtry=fsolve(@finitial,Vcon(2:end)*1e7,[],Vcon(1));
Vcon(2:end)=Vtry/1e7;

tic;

[tout5,y3d] = ode15s(cxrate,tspan,Vcon,options);

			     % display some parameters
[mout,nout] = size(y3d);	
sumv        = sum(y3d(end,:))
flux    = y3d(end,:)' .* (setrat/Lz)  ;
sumflux = sum(flux);

save simpout                            % temporary storage

toc                                     % get the time it has taken