anomtransport.m

% function [val,tval,Dis,Tim] = anomtransport(N,vel,tau_w,tau_p,topt,sopt,alpha,beta,contflag)
% anomalous transport simulator (1D)
% Inputs:
%  N is number of segments
%  vel is velocity in microns/sec
%  tau_w is waiting time between runs
%  tau_p is persistence time of run
%  topt = (0,1,2) is for waiting time (exp, scaling,fBm)
%  sopt = (0,1,2) is for segment length (exp, scaling,fBm)
%  alpha is anom exponent on waiting time
%  beta is anom exponent on segment length (beta = H for fBm)
%  contflag = (0,1) continuity flag for continuous val and tval
% Outputs:
%  Dis = discrete distances
%  Tim = discrete times
%  val = continuous distances
%  tval = continuous times
%
% Example
%     [Dis,Tim,val,tval] = anomtransport(N,tau_w,tau_p,topt,sopt,alpha,beta,1)

function [Dis,Tim,val,tval,mean_s,mean_t] = anomtransport(N,vel,tau_w,tau_p,topt,sopt,alpha,beta,contflag)

%vel = vel*(randexp(1,N)+0.01);
%vel = vel*randn(1,N);
vel = ones(1,N)*vel;

if topt == 0
    t = tau_w*randexp(1,N);
elseif topt == 1
    %t = tau_w*randscale(1,N,alpha);
    t = abs(tau_w*randLevy(N,alpha,1));
elseif topt == 2
    t = abs(tau_w*randHurst(N,alpha));    
else
    t = ones(N)*1e-3;
end
mean_t = mean(t);

if sopt == 0
    s = vel.*tau_p.*sign(randn(1,N)).*randexp(1,N);
elseif sopt == 1
    %s = vel*tau_p*sign(randn(1,N)).*randscale(1,N,beta);
    %s = vel*tau_p*sign(randn(1,N)).*randLevy(N,beta,1);
    s = sign(randn(1,N)).*vel.*tau_p.*randLevy(N,beta,1);
elseif sopt == 2
    s = vel.*tau_p.*randHurst(N,beta);
elseif sopt == 3
    dd = 1:N; tauage = N/4;
    %vel = vel.*exp(-dd/tauage);
    vel = vel./(dd.^0.25);
    s = vel.*tau_p.*sign(randn(1,N)).*randexp(1,N);
    %keyboard
else 
    s = ones(N)*1e-1;
end
mean_s = mean(abs(s));

Dis(1) = 0;
Tim(1) = 0;
Dis(2) = s(1);
Tim(2) = abs(s(1))/abs(vel(1));
Dis(3) = Dis(2);
Tim(3) = Tim(2) + t(1);
eloop = 3;
for loop = 2:N
    eloop = eloop + 1;
    Dis(eloop) = Dis(eloop-1) + s(loop);
    Tim(eloop) = Tim(eloop-1) + abs(s(loop))/abs(vel(loop));
    eloop = eloop + 1;
    Dis(eloop) = Dis(eloop-1);
    Tim(eloop) = Tim(eloop-1) + t(loop);
end

if contflag == 1
    Nsamp = N;
    T = max(Tim);
    delT = T/Nsamp;
    for loop = 1:Nsamp-1
        tind = loop*delT;
        [I, ~] = changesign(tind - Tim);
        Ileft = I-1;
        Iright = I;
        delt = tind - Tim(Ileft);
        Del = delt/(Tim(Iright)-Tim(Ileft));
        val(loop) = Dis(Ileft) + Del*(Dis(Iright)-Dis(Ileft));
        tval(loop) = tind;
    end
else
    val = 0;
    tval = 0;
end

%keyboard