% simpgraph.m
%  handle the graphics stuff for simpcoag
    clear 
    load simpout;           % get the results from last coag calc                

  total= sum(y3d');
  axisvec=[0,nsect,min(total),max(total)];

  % calculate the number, mass spectrum
  
% allocate some memory
nspec=zeros(ntspan,nsect);
massspec=nspec;
fluxspec=nspec;
fluxsect=nspec;
diamat=ones(ntspan,1)*speccon.rf' *2;
timmat=tspan' * ones(1,nsect);

for jind=1:ntspan
  yout=y3d(jind,:);
  nspec(jind,:)=yout ./ (1.5*speccon.Vlow') ./ speccon.dwidth';
  massspec(jind,:)=yout ./ speccon.dwidth';
                              % units: setrat is m/d, yout vol/vol=cm^3/cm^3
                              % set fluxes as cm^3 m^-2 d^-1, cm^-1 by *1e6 cm3/m3
  fluxsect(jind,:)=yout .* setrat' * 1e6;       
  fluxspec(jind,:)=massspec(jind,:) .* setrat' * 1e6;
end

fgstyle=0;          % parameter for output style: 0 2-d, 1 3-d

dum = input(['How would you like output: 2-d figs (0-default); 3-d (1)?  ']);
if ~isempty(dum),fgstyle=dum;end;

totmass=(sum(y3d'))';

figure
clf

if fgstyle
  
	subplot(2,2,1);
	surf(diamat,timmat,fluxspec);
	xlabel('Apparent particle diameter (cm)');
	ylabel('Time (d)');
	zlabel('Volume flux spectrum {\it V n v_s}(cm^3 cm^{-1} m^{-2} d^{-1})');
	%set(gca,'xscale','log','zscale','log');
	view(50.5,46);
	title([cprog,' ',date]);
	shading interp
		
	subplot(2,2,3);
	surf(diamat,timmat,nspec);
	xlabel('Apparent particle diameter (cm)');
	ylabel('Time (d)');
	zlabel('Number spectrum {\it n}( # cm^{-4})');
	set(gca,'xscale','log','zscale','log');
	view(50.5,46);
	title([cprog,' ',date]);
	
	subplot(2,2,4);
	surf(diamat,timmat,massspec)
	xlabel('Apparent particle diameter (cm)');
	ylabel('Time (d)');
	zlabel('Conserved volume spectrum {\it V n}( # cm^{-1})');
	view(50.5,46);
	title([cprog,' ',date]);
	shading interp

else
	subplot(2,2,1);
%	semilogy(diamat',nspec');
  loglog(diamat',nspec');
  xlabel('Apparent particle diameter (cm)');
	ylabel('Number spectrum {\it n}( # cm^{-4})');
	 
  subplot(2,2,3);
  plot(diamat',fluxspec');
	xlabel('Apparent particle diameter (cm)');
	hy2=ylabel(char('Volume flux spectrum {\it V n v_s}', ...
    '(cm^3 cm^{-1} m^{-2} d^{-1})'));	
  set(hy2,'units','norm','pos',[-.266,.49]);
  
% 	subplot(2,2,4);
% 	plot(diamat',massspec')
% 	xlabel('Apparent particle diameter (cm)');
% 	ylabel('Conserved volume spectrum {\it V n}( # cm^{-1})');
%   
 
subplot(2,4,7)
totflux=(sum(fluxsect'))';
plot(tout5,fluxsect,tout5,totflux,'*--');
xlabel('time(d)');
ylabel('sectional flux (cm^3 m^{-2} d^{-1}/sect)');

subplot(2,4,8)
plot(tout5,totflux./totmass/1e6);
xlabel('time(d)');
ylabel('Average v (m d^{-1})');


end

subplot(2,4,3);
semilogy(tout5,y3d,tout5,totmass,'*--');
xlabel('time(d)');
ylabel('Sectional concentration(vol/vol/sect)');

% print out the parameter values:
subplot(2,4,4);
plot(0,0);
axis([0 1 0 1]);
cinfo=char(['fractal dimension fr= ',num2str(speccon.fr)],...
  ['unit part diam (cm) dia0=',num2str(speccon.dia0)],...
  ['\Delta \rho (g cm^{-3}) = ',num2str(physcon.delrho)],...
  ['number of sections, nsect= ',num2str(nsect)],...
  ['nkernel = ',num2str(betasd.nkernel)],...
  ['file with beta coeffs, cbeta = ',cbeta],...
  ['length of run (d), tfinal= ',num2str(tfinal)],...
  ['stickiness, \alpha (alpha)= ',num2str(alpha)],...
  ['mixed layer thickness (m), Lz= ',num2str(Lz)],...
  ['average shear (s^{-1}), \gamma (Lgam)= ',num2str(Lgam)],...
  ['Spec growth(d^{-1}), \mu (Lmu)=', num2str(Lmu)],...
  ['t_0 part #, sect 1 (# cm^{-3}), num\_1= ',num2str(num_1)],...
  ['First sect w agg gro, nsecgro= ',num2str(nsecgro)],...
  ['at= ',num2str(at),'    rt= ',num2str(rt)],...
  ' ',...
  'dashed, + lines represent total');
 ht=text(-.1,1.1,cinfo);
 set(gca,'vis','off');
 set(ht,'verticalal','top');

stamp('simpgraph.m');

orient landscape;