%%% To run: erode( flag, amplitude, frequency, pedestal )
%%%   where, flag defines erosion rate:
%%% 	1 = zero (steady state)
%%% 	2 = constant
%%% 	3 = full-wave rectification
%%% 	4 = square-wave
%%% 	5 = saw-tooth
%%% 	6 = exponential
%%% 	7 = step function
%%%   and amplitude/frequency/pedastal control the periodic erorison 
%%%   sequences (3-7) where erosion units are in m/MYrs. Amplitude 
%%%   defines max erosion and pedastal defines min erosion
%%%

function[concentration] = erode( flag, AMPLITUDE, FREQUENCY, PEDESTAL, ATLANTICBOREAL, ANTHROPOGEN, RECENT  );

%%% DEFAULT PARAMETERS OF PERIODIC EROSIONS 
%%%
if( ~exist( 'AMPLITUDE' ) ) AMPLITUDE = 150; end;
if( ~exist( 'FREQUENCY' ) ) FREQUENCY = 2; end;
if( ~exist( 'PEDESTAL' ) ) PEDESTAL = 22; end;
if( ~exist( 'ATLANTICBOREAL' ) ) ATLANTICBOREAL = 22; end;
if( ~exist( 'ANTHROPOGEN' ) ) ANTHROPOGEN = 150; end;
if( ~exist( 'RECENT' ) ) RECENT = 22; end; 

%%% VARIOUS CONSTANTS
%%%
flag = 4;
timestep = 10;		% time step - this times the below defines the timestep
						% for high erosion rates, lower the time step
maxtime  = 10e4;			% max time - sets the run duration of the model
time     = [timestep : timestep : maxtime];	% time interval

P0       = 17*timestep;			% surface production
C0       = 100;           % initial concentration
mu       = 2.65/165;	% constant
lambda   = log(2)/1.5e6 * timestep;% decay rate
e        = erosionrate( flag, time, AMPLITUDE, FREQUENCY, PEDESTAL, ATLANTICBOREAL, ANTHROPOGEN, RECENT );
e        = 1e2/1e6 * e; % m/Myrs to cm/yrs
e 		 = e * timestep;

%%% ERODE IT...
%%%
%%%100 sets the sample interval for the data - reduce for finer resolution
ttt = 1 : 1 : length(time);
concentration = zeros( 1, length(ttt) );
count = 1;
for t = ttt
	C = C0;
	CS = fliplr( cumsum( fliplr(e(1:t)) ) );
	for k = 1 : t
		c = P0 * exp(-mu*CS(k)) - lambda*C;
		C = C + c;
	end
	concentration(count) = C;
	count = count + 1;
end

%%% ESTIMATE EROSION FROM ASSUMED STEADY-STATE CONCENTRATION
%%%
eEst = P0 ./ (concentration*mu) -lambda;
eEst = eEst/timestep * 1e6/1e2;


%%% PLOT EVERYTHING
%%%
set( gca, 'FontSize', 12 );
subplot(211); plot( time(ttt), concentration, 'bo-' );
xlabel( 'Time (y)' );
ylabel( '^{10}Be (atoms/gram)' );
set( gca, 'xlim', [1e3 maxtime]);
set( gca, 'ylim', [10 2e5] );			%change this to reflect expected max concentration
grid on;
subplot(212); plot( time(ttt), e(ttt)/timestep*1e6/1e2, 'r' );
hold on;
plot( time(ttt), eEst, 'k--' );
hold off;
%set( gca, 'Xscale', 'log' );
xlabel( 'Time (y)' );
ylabel( 'Erosion (m/My)' );
set( gca, 'xlim', [1e3 maxtime] );		%change this to 1e5 when time is >10^6 yrs-1e2 otherwise
set( gca, 'ylim', [10 4e2] );            %change this to be just over the max erosion
format long;
disp( [time(ttt)' concentration' e(ttt)'/timestep*1e6/1e2 eEst'] );


return;

%%% =======================================================================
%%% DEFINE ERORSION RATE AS A FUNCTION OF TIME
%%%
function[e] = erosionrate( flag, time, amp, freq, ped, atl, ant, rec )
   
   if( flag == 1 ); % zero (steady state)
      e = zeros(size(time));
   elseif( flag == 2 ) % constant
      e = amp * ones(size(time));
   elseif( flag == 3 ) % full-wave rectification
      e = sin( freq*2*pi*time/max(time) ); 
      e(find(e<0)) = 0; 
	  e = (amp-ped) * e + ped;
   elseif( flag == 4 ) % square-wave
      e = sin( freq/1.79*3*pi*time/max(time) );
      e( find(e<0) ) = 0;
      e( find(e>0) ) = 1;
	  e = (amp-ped) * e + ped;
         elseif( flag == 5 ) % saw-tooth
      N = length(time);
      n = N / freq;
      saw = [1:n]/n;
      e = repmat( saw, 1, freq );
      e = (amp-ped) * e + ped;
      if( length(e)<N ) % freq is a non-integer multiple of time
	    e( length(e)+1:N ) = e(1:N-length(e));
      end
   elseif( flag == 6 ) % exponential
      N = length(time);
      n = N / freq;
      expo = exp(-([n:-1:1]/n)*6);
      e = repmat( expo, 1, freq );
      e = (amp-ped) * e + ped;
      if( length(e)<N ) % freq is a non-integer multiple of time
	    e( length(e)+1:N ) = e(1:N-length(e));
      end
   elseif( flag == 7 ) % step 
      N = length(time);
      n = round(N/1.5);
      e = zeros( 1,n );
      e( length(e)+1:N ) = 1;
	  e = (ped-amp) * e + amp;
   elseif( flag == 8 ) % constant
      e = amp * ones(size(time));
   else
      fprintf( 'erosionrate(): invalid flag (%d)-using steady state\n', flag );
      e = zeros(size(time));
   end