%  simpbeta.m
%  simplified beta calc program
%  modify betacalc2a.m
%  Input:  keyboard options
%  Output:  results of beta calculations saved to file, betaf.mat is default
%  Functions called:  simptcalc, setvel
%     results are used by simpcoag.m to calculate results
% copyright George A. Jackson, Texas A&M Univ. 2001
%	    errors in delrho calc made, 8 mar 01
global setcon
cfile='simpbeta.m';

                        % load spectrum: D3E c5E c6E dcomb histspectra ncomb tsamp
                        % set up constants, put in structures
rhofl       = 1.0275;			              % g cm^3  fluid density
visc        = 0.01;			                  % cm^2 /s   kinem visc
grav        = 980;			                  % grav acc, cm/s^2
daysec      = 8.64e4;			              % seconds/d
gamma       = 1.;                             % default shear 1/s,  adjust in coag calc
ratrirg=1.36;	        % ratio of  interaction radius to radius of gyration
			            % value comes from rogak and flagan "outer radius"
			            % note, they get hydraulic radius/rg=1

                        %   setup the default constants

fr          = 2.33;                          % fractal dimension 
dia0        = 20e-4;                          % unit particle diameter
delrho      = (4.5*2.48)*visc*rhofl/grav* (dia0/2)^(-0.83);     % correct 8 mar 2001
nsect       = 10;
nkernel     = 3;                             % which kernel? rect-1; curv-2; frac-3
coutfile    = 'betaf';

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

if linp                                      % if unwilling to accept the defaults
	display('to accept the default value, simply press the cr');   
	
	dum =input(['Fractal dimension [',num2str(fr),']\n',...
         ' fr= ']);
	if ~isempty(dum),fr=dum;end;
  
	dum =input(['Number of sections  [',num2str(nsect),']\n',...
         ' nsect= ']);
	if ~isempty(dum),nsect=dum;end;

 	dum =input(['Minimum diameter (cm) [',num2str(dia0),']\n',...
         ' dia0= ']);
	if ~isempty(dum),dia0=dum;end;

                         % calc the default delrho (g/cm^3)
                         % stokes flow: vs=2 g a^2 delrho/9 visc rhofl
                         % smayda data: vs=2.48 a^1.17 (cgs)
  delrho=(4.5*2.48)*visc*rhofl/grav* (dia0/2)^(-0.83);
 	dum =input(['Unit particle density (g/cc) [',num2str(delrho),']\n',...
         ' delrho= ']);
	if ~isempty(dum),delrho=dum;end;
  
 	dum =input(['Which coagulation kernel? (input number: 1,2,3):\n',...
            'Rectilinear (1)\n','   Curvilinear (2)\n   Fractal (3)\n',...
            'Default is  [',num2str(nkernel),']\n',' nkernel= ']);
	if ~isempty(dum),nkernel=dum;end;
  
 	dum =input(['Name of output file for betas: Default is ',coutfile],'s');
	if ~isempty(dum),coutfile=dum;end;
  
end

a0=dia0/2;                               % calc the radius of the unit particle

                                         % calculate physical properties of sections
V0      = pi/6* dia0^3;                  % conserved volume at lowest boundary
Vlow    = V0 * 2 .^( 0: (nsect-1) )';    % con volume of all lower bounds, col vector
Vup     = 2. * Vlow;                     % con volume of all upper bounds
volstd  = 1.5*Vlow;                      % mid volume
dcomb   = (Vlow*6/pi).^ (1/3);           % conserved diameter of lower bound
dwidth  = (2^.33333-1)*dcomb;            % width of a compartment

                                         % need various radii settling vel calc     
amfrac  = (4. / 3. *pi)^(-1/fr) * a0^(1-3/fr);   % rad-f from mass/vol
amfrac  = amfrac*(0.6)^(.5);             % this adjusts to rg from r0 at rv=ro
bmfrac  = 1. / fr;		                 % see notes, 27 nov

					                     % calculate the setting velocity terms
setcon  = 2/9*delrho/rhofl*grav/visc; % constant used in calculating settling velocity
					                     % calc cons for stokes sett with rv, rg
                                         
                                         % mid-section conserve, frac rad
rf      = amfrac*volstd .^ bmfrac;       % calculate the fractal radius
rv      = (.75/pi*volstd) .^ 0.3333;	 % calculate the conserved radius
vset    = setvel(rf,rv);				 % calculate settling velocity


if nkernel==1                            % rectilinear kernel
  cxmodea ='betrsh';
  cxmodeaf=@betrsh;
  cxmodeb ='betrds';
  cxmodebf=@betrds;
elseif nkernel==2;                       % curvilinear kernel
  cxmodea ='betcsh';
  cxmodeaf=@betcsh;
  cxmodeb ='betcds';
  cxmodebf=@betcds;
elseif nkernel==3                        % fractal kernel
  cxmodea ='betfsh';
  cxmodeaf=@betfsh;
  cxmodeb ='betfds';
  cxmodebf=@betfds;
else
  error(['incorrect choice for coagulation kernel in ',cfile]);
end

                % consolidate the constants into two structures

speccon=struct('dia0',dia0,'a0',a0,'fr',fr,'dcomb',dcomb,'nsect',nsect,...
  'rf',rf,'rv',rv, 'vset',vset,'dwidth',dwidth,'amfrac',amfrac,'bmfrac',bmfrac,...
  'Vlow',Vlow,'Vup',Vup,'ratrirg',ratrirg);

physcon=struct('delrho',delrho,'rhofl',rhofl,'visc',visc,'grav',grav,'setcon',setcon,...
  'gamma',gamma,'daysec',daysec);


tic             % start the clock
                % do the calcs of the betas
[beta25s,beta1s,beta2s,beta3s,beta4s,beta5s]=simptcalc(cxmodea,speccon,physcon);

[beta25d,beta1d,beta2d,beta3d,beta4d,beta5d]=simptcalc(cxmodeb,speccon,physcon);

                % prepare output 
betasd = struct('beta25s',beta25s,'beta1s',beta1s,'beta25d',beta25d,'beta1d',beta1d,...
    'nkernel',nkernel);

betaheader = [date,' program=',cfile, ' output file=',coutfile];

cxcomd2 = ['save ',coutfile,' betasd speccon physcon betaheader'];
eval(cxcomd2);

toc
