function [beta25,beta1,beta2,beta3,beta4,beta5]...
	=simptcalc(cxmode,speccon,physcon)
%simptcalc.m		 alter from tcalc-  include ratrig
% calculate the different interaction terms for the particle size spectrum
% do the numeric integrations for the 5 different terms, assemble into
% a matrix of interaction terms
%	INPUT:	ml, mu- vectors of lower, upper mass boundaries
%			mu should be 2*ml
%		xgamma: average shear
%		cxmode:  name of func to be used to calc bet  kern(betrect, betcurv)
%		xamfrac, xbmfrac:  coefs to calc interact radius from mass(vol)
%	OUTPUT: beta25:  matrix giving the effect of irow, jcol on jcol
%			beta1:	matrix giving the effect of irow, jcol-1 on jcol
%
% copyright George A. Jackson, Texas A&M Univ. 2001
global collint amfrac bmfrac cmode gamma
global ratrirg 
global collint2 mil miu mjl mju mjj rjj  rvjj

ldebug  =0;             % set to 1 to get intermediate results saved
ldebug2 =0;             % set to 1 to get timing on each integ

if ldebug2,tic;end      % start the clock ticking

ml    =speccon.Vlow;
mu    =speccon.Vup;
gamma =physcon.gamma;
amfrac=speccon.amfrac;
bmfrac=speccon.bmfrac;
cmode =cxmode;
ratrirg=speccon.ratrirg;
ncol  =speccon.nsect;
betamast=zeros(ncol,ncol);          %init the space

%			notice, have funny conventions
%			jcol is the compartment gaining, losing mass
%			icol is the other colliding particle
%			however, in notes, I used i for the affected compartment
%				j for the other interacting particle

             % calculate the loss from i to collision with larger j
beta5=betamast;
for jcol=1:(ncol-1)
	for irow=(jcol+1):ncol
		mil=ml(jcol);miu=mil*2;
		mjl=ml(irow);mju=mjl*2.;
		beta5(irow,jcol)=quadl(@integr5a,mjl,mju)/mil/mjl;
    %beta5(irow,jcol)=quad8('integr5a',mjl,mju)/mil/mjl;
	end
end
if ldebug2,toc5=toc,end
if ldebug, save tempbetcurv beta5 toc5,end;
             % calculate the loss in jcol from collision with self
beta4=betamast;
for jcol=1:ncol
	mil=ml(jcol);miu=mil*2;
	mjl=ml(jcol);mju=mjl*2;
	beta4(jcol,jcol)=quadl(@integr4a,mjl,mju)/mil/mjl;
	%beta4(jcol,jcol)=quad8('integr4a',mjl,mju)/mil/mjl;
end

if ldebug2,toc4=toc,end
if ldebug, save tempbetcurv beta5 toc5 beta4 toc4,end;
            % calculate the loss in jcol to collision with comp irow, irow<jcol
beta3=betamast;
for jcol=2:ncol
	for irow=1:(jcol-1)
		mil=ml(jcol);miu=mil*2;
		mjl=ml(irow);mju=mjl*2.;
		beta3(irow,jcol)=quadl(@integr3a,mjl,mju)/mil/mjl;
%		beta3(irow,jcol)=quad8('integr3a',mjl,mju)/mil/mjl;
	end
end

if ldebug2,toc3=toc,end
if ldebug, save tempbetcurv beta5 toc5 beta4 toc4 beta3 toc3,end;
                % calculate the gain in jcol from collision of irow with jcol, irow<jcol
warning off
beta2=betamast;
for jcol=2:ncol
	for irow=1:(jcol-1)
		mil=ml(jcol);miu=mil*2;
		mjl=ml(irow);mju=mjl*2.;
		beta2(irow,jcol)=quadl(@integr2a,mjl,mju)/mil/mjl;
%		beta2(irow,jcol)=quad8('integr2a',mjl,mju)/mil/mjl;
	end
end
warning on

if ldebug2,toc2=toc,end
if ldebug, save tempbetcurv beta5 toc5 beta4 toc4 beta3 toc3 beta2 toc2;end;

                % calculate the gain in jcol from collision of jcol-1 with irow, irow<jcol
beta1=betamast;
for jcol=2:ncol
	jtemp=jcol-1;
	for irow=1:jtemp
		mil=ml(jtemp);miu=mil*2;
		mjl=ml(irow);mju=mjl*2.;
		beta1(irow,jcol)=quadl(@integr1a,mjl,mju)/mil/mjl;
	%	beta1(irow,jcol)=quad8('integr1a',mjl,mju)/mil/mjl;
	end
end

if ldebug2,toc1=toc,end
if ldebug, save tempbetcurv beta5 toc5 beta4 toc4 beta3 toc3 beta2 toc2 beta1 toc1;end;

beta25=beta2-beta3-beta4-beta5;


%*******************end of simptcalc********************


function xint=integr1a(mj)
% this is the outer integral for case 1: gain mass to i+1 from comp j, i in
% collision with particle in comp i
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj

nj=length(mj);xint=mj*0;
rj=amfrac*mj.^bmfrac;
rvj=(.75/pi*mj).^ 0.3333333;
for iv=1:nj
	mjj=mj(iv);
	rjj=rj(iv);
	rvjj=rvj(iv);
	mlow=max([miu-mjj,mil]);	% handle same interaction  correctly
%	xint(iv)=quad8a('integr1b',mlow,miu);
	xint(iv)=quadl(@integr1b,mlow,miu);
end





function yint=integr1b(mi)
% this is the inner integral for case 1: gain mass to i+1 from comp j, i in
% collision with particle in comp i
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj

%
ri=amfrac*mi.^bmfrac;
ni=length(mi);rj=rjj*ones(1,ni);mj=mjj*ones(1,ni);
rvi=(0.75/pi*mi).^0.3333;rvj=rvjj*ones(1,ni);
beta=feval(cmode,ri,rj,rvi,rvj,gamma);
yint=beta.*(mi+mj)./mi./mj;





function xint=integr2a(mj)
% this is the outer integral for case 2: gain mass to i from comp j in
% collision with particle in comp j
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj
nj=length(mj);xint=mj*0;
rj=amfrac*mj.^bmfrac;
rvj=(.75/pi*mj).^ 0.3333333;
for iv=1:nj
	mjj=mj(iv);
	rjj=rj(iv);
	rvjj=rvj(iv);
%	xint(iv)=quad8a('integr2b',mil,miu-mjj);
	xint(iv)=quadl(@integr2b,mil,miu-mjj);
end





function yint=integr2b(mi)
% this is the inner integral for case 2: gain mass to j from comp i in
% collision with particle in comp j
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj
%
ri=amfrac*mi.^bmfrac;
ni=length(mi);rj=rjj*ones(1,ni);mj=mjj*ones(1,ni);
rvi=(0.75/pi*mi).^0.3333;rvj=rvjj*ones(1,ni);
beta=feval(cmode,ri,rj,rvi,rvj,gamma);
yint=beta./mi;





function xint=integr3a(mj)
% this is the outer integral for case 3: loss mass from comp i from
% collision with particle in comp j<i
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj
nj=length(mj);xint=mj*0;
rj=amfrac*mj.^bmfrac;
rvj=(.75/pi*mj).^ 0.3333333;
for iv=1:nj
	mjj=mj(iv);
	rjj=rj(iv);
	rvjj=rvj(iv);
%	xint(iv)=quad8a('integr3b',miu-mjj,miu);
	xint(iv)=quadl(@integr3b,miu-mjj,miu);
end
xint=xint ./ mj;






function yint=integr3b(mi)
% this is the inner integral for case 3: loss mass from comp i from
% collision with particle in comp j<i
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj
%
ri=amfrac*mi.^bmfrac;
ni=length(mi);rj=rjj*ones(1,ni);mj=mjj*ones(1,ni);
rvi=(0.75/pi*mi).^0.3333;rvj=rvjj*ones(1,ni);
beta=feval(cmode,ri,rj,rvi,rvj,gamma);
yint=beta;





function xint=integr4a(mj)
% this is the outer integral for case 4: loss mass from comp i from
% collision with particle in comp i
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj
nj=length(mj);xint=mj*0;
rj=amfrac*mj.^bmfrac;
rvj=(.75/pi*mj).^ 0.3333333;
for iv=1:nj
	mjj=mj(iv);
	rjj=rj(iv);
	rvjj=rvj(iv);
%	xint(iv)=quad8a('integr4b',mil,miu);
	xint(iv)=quadl(@integr4b,mil,miu);
end





function yint=integr4b(mi)
% this is the inner integral for case : loss mass from comp i from
% collision with particle in comp i
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj
%
ri=amfrac*mi.^bmfrac;
ni=length(mi);rj=rjj*ones(1,ni);mj=mjj*ones(1,ni);
rvi=(0.75/pi*mi).^0.3333;rvj=rvjj*ones(1,ni);
beta=feval(cmode,ri,rj,rvi,rvj,gamma);
yint=(mi+mjj) ./ mi ./ mjj .* beta;





function xint=integr5a(mj)
% this is the outer integral for case 5: loss mass from comp i from
% collision with particle in comp j, i<j
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj

nj=length(mj);xint=mj*0;
rj=amfrac*mj.^bmfrac;
rvj=(.75/pi*mj).^ 0.3333333;
for iv=1:nj
	mjj=mj(iv);
	rjj=rj(iv);
	rvjj=rvj(iv);
%	xint(iv)=quad8a('integr5b',mil,miu);
	xint(iv)=quadl(@integr5b,mil,miu);
end
xint=xint ./ mj;




function yint=integr5b(mi)
% this is the inner integral for case 5: loss mass from comp i from
% collision with particle in comp j, i<j
%
global collint amfrac bmfrac cmode gamma
global collint2 mil miu mjl mju mjj rjj rvjj
%
ri=amfrac*mi .^ bmfrac;
ni=length(mi);rj=rjj*ones(1,ni);mj=mjj*ones(1,ni);
rvi=(0.75/pi*mi).^0.3333;rvj=rvjj*ones(1,ni);
beta=feval(cmode,ri,rj,rvi,rvj,gamma);
yint=beta;

