subtest3.m

% subtest.m
% 1D QELS anomalous transport test simulator
% Modified from subtest2.m
% Calculate the walk, time series and power spectrum
% Calls
%   Powerspec
%   anomtransport
%   changesign
%   logNsamplin
%   reversevec

clear
format compact
tic

compu = computer

if ~strcmp(compu,'MACI64')
	addpath /home/nolte/UserFunctions/ -BEGIN
	rng shuffle
end

lam = 0.84;  % 0.84 wavelength
vel = 0.31;  % 0.31
vel_dist = 0;  % 0 = delta   1 = uniform   2 = gaussian

% 0 = randexp  1 = powerlaw
topt = 0;   alpha = 0.5;    % 1.0 waiting time: subdiffusion (0,1,2) for (exp, scaling,fBm)
tau_w = 0.0;                  %  0 0.5 5 tau_w is waiting time between runs
sopt = 0;   beta = 0.85;       % 1.0 segment length: superdiffusion  (0,1,2) for (exp, scaling,fBm)
tau_p = 0.5;                 %  0 0.5 5 tau_p is persistence time of run

Npart = 200;  % number of particles in each ensemble (100 200 400 800  4000)
Nsuper = 50;      % number of averages (10 20 40 80  800)
Nang = 1;    % Number of angles (Nang=1 gives 1-D max Doppler) (1 20)

Nphase = 400;
lp = vel*tau_p;

delT = 0.04;    % 0.04
T = 2000*delT;  % 2000*delT
NumT = floor(T/delT);
NDt = round(tau_p/delT/2);

figflag = 1;    % Turn off in-loop figures


N = 2*round(T/(tau_w + tau_p));  % number of segments (1000)

q = 4*pi*1.33/lam;

DopFreq = q*vel/2/pi
DopNum = q*vel*tau_p

nam0 = strcat('ST',num2str(vel),'-',num2str(topt),num2str(tau_w));
if topt ~=0
    nam0 = strcat(nam0,num2str(alpha));
end
nam0 = strcat(nam0,'-',num2str(sopt),num2str(tau_p));
if sopt ~=0
   nam0 = strcat(nam0,num2str(beta)); 
end
nam = strcat('anomtest/',nam0);

% find best time basis (sample 100 cases and find max)
for loop = 1:100
    [Dis,Tim] = anomtransport(N,vel,tau_w,tau_p,topt,sopt,alpha,beta,0);
    ttst(loop) = max(Tim);
    mxDis(loop) = max(abs(Dis));
end
% Ttst = min(ttst);
% D = max(mxDis);
% [B,I] = sort(mxDis);
% D = 1.05*B(100);   % Max Distance
D = vel*T;


tindex(1) = floor(0.01*NumT);
tindex(2) = floor(0.03*NumT);
tindex(3) = floor(0.09*NumT);
tindex(4) = floor(0.27*NumT);
tindex(5) = NumT;

%tindex = tindex/10;

% Launching main simulator
tic
ssum = zeros(1,101); Ssum = zeros(1,101); Ahetsum = zeros(1,2000); Ahomsum = Ahetsum;
phsum = zeros(1,Nphase);
Wsum = zeros(1,NumT);
hsum = zeros(5,100); hsum2 = zeros(5,100);
mnssum = 0; mntsum = 0;
for sloop = 1:Nsuper  % numer of averages
    displine('sloop =',sloop)
    
    ysum = zeros(1,NumT);   % The superposed field from Npart particles for each of Nangle angles
    if figflag == 1
        figure(1)
        clf
        axis off
        hold on
    end
    
    for eloop = 1:Npart
        for angloop = 1:Nang
            if angloop ==1 % Do the max Doppler shift
                ang = 0;
                angfac = 1;
            else  % Average the Doppler shifts over angle
                ang = angloop*(pi/2)/(Nang+1);
                angfac = sin(ang);
            end
            
            %[~,~,Dis,Tim] = randflightsub(N,alpha);
            if vel_dist == 0
                velocity_in = vel*cos(ang);
            elseif vel_dist == 1
                velocity_in = abs(rand*vel*cos(ang));
            else
                velocity_in = abs(randn*vel*cos(ang));
            end
            
            [Dis,Tim,val,tval,mean_s,mean_t] = anomtransport(N,velocity_in,tau_w,tau_p,topt,sopt,alpha,beta,0);  % Set contflag = 0
            maxT(eloop) = max(Tim);
            mnssum = mnssum + mean_s/Npart/Nsuper;
            mntsum = mntsum + mean_t/Npart/Nsuper;
            
            %Nsamp = N;
            %T = floor(5000000/(alpha^1.3));          % alpha = 0.5-25000   alpha = 0.95-10000
            %delT = T/Nsamp;
            %         flag = 0;
            for loop = 1:NumT
                tind = (loop-1)*delT;
                [I, ~] = changesign(tind - Tim);
                Ileft = I-1;
                Iright = I;
                if Ileft == 0
                    %                 flag = 1;
                    tlength = length(Tim);
                    Ileft = tlength-1;
                    Iright = tlength;
                end
                delt = tind - Tim(Ileft);
                Del = delt/(Tim(Iright)-Tim(Ileft) + 1e-3);
                tmp = Dis(Ileft) + Del*(Dis(Iright)-Dis(Ileft));
                if isnan(tmp)
                    keyboard
                end
                val(loop) = tmp;
                tval(loop) = tind;
            end
            
            if figflag == 1
                figure(1)
                plot(tval,val,'LineWidth',1.5)
                %keyboard
                pause(0.05)
            end
            
            A(eloop,:) = val;
            
            ysum = ysum + angfac*exp(i*2*pi*rand)*exp(i*q*val);     % This is the complex field
            
            % Particle densities
            for tloop = 1:5
                fac = tindex(tloop)/NumT;
                dat = A(:,tindex(tloop));
                [ytmp,xtmp] = histfix(dat,100,-1.05*fac*D,1.05*fac*D);  % scaled
                hs(tloop,:) = ytmp;
                [ytmp2,xtmp2] = histfix(dat,100,-1.05*D,1.05*D);  % non-scaled
                hs2(tloop,:) = ytmp2;
            end
            
        end     % end angloop
    end     % end eloop for Npart
    
    %keyboard
    
    if figflag == 1
        hold off
        
        figure(2)
        imagesc(A)
        ylabel('Particle Number')
        xlabel('Time (sec)')
        colormap(jet)
    end
    
    width = std(A);
    
    if figflag == 1
        figure(3)
        loglog(tval,width)
        xlabel('Time (sec)')
        ylabel('RMS displacement')
        title('RMS displacement')
    end
    
    Wsum = Wsum + width/Nsuper;
    
    if figflag == 1
        figure(4)
        plot(real(ysum))
        xlabel('Time')
        title('Real Field')
    end
    
    szysum = length(ysum);
    [y,freq,yss,freqss] = Powerspec(real(ysum));    % Heterodyne
    [Y,~,Yss,~] = Powerspec(abs(ysum).^2);          % Homodyne
    [~,szs] = size(freqss);
    [xp, yp] = logNsamplin(freqss(2:szs),yss(2:szs),100); % Het
    [~, Yp] = logNsamplin(freqss(2:szs),Yss(2:szs),100);  % Hom
    frq = (1/delT)*exp(xp);
    
    ssum = ssum + yp/Nsuper;
    Ssum = Ssum + Yp/Nsuper;
    
    phseraw = atan2(imag(ysum),real(ysum));
    delphse = mod(pi + diffDt(phseraw,NDt),2*pi)-pi;
    [hprob,hphse] = histfix(delphse(1:szysum-NDt),Nphase,-pi,pi);
    
    phsum = phsum + hprob/Nsuper;
    
    
    % Test autocorrelation analysis
    Ahet = autocorrel(real(ysum));
    Fhet = fft(fftshift(Ahet));
    Ahetsum = Ahetsum + abs(Fhet)/Nsuper;

    Ahom = autocorrel(abs(ysum).^2);
    Fhom = fft(fftshift(Ahom));
    Ahomsum = Ahomsum + abs(Fhom)/Nsuper;

    [~, Ayp] = logNsamplin(freqss(2:szs),Ahetsum(3:szs+1),100); % Het
    [~, AYp] = logNsamplin(freqss(2:szs),Ahomsum(3:szs+1),100);  % Hom
    
    if figflag == 1
       figure(14)
        loglog(frq,Ayp,frq,AYp)
        xlabel('Frequency (Hz)')
        title('Spectra')
        legend('Het','Hom')
    end

    
    if figflag == 1
        figure(5)
        loglog(frq,yp,frq,Yp)
        xlabel('Frequency (Hz)')
        title('Spectra')
        legend('Het','Hom')
    end
    
    for tloop = 1:5
        hsum(tloop,:) = hsum(tloop,:) + hs(tloop,:)/Nsuper;
        hsum2(tloop,:) = hsum2(tloop,:) + hs2(tloop,:)/Nsuper;
    end
    
    
end     % Nsuper
duration = toc

displine('Effective Doppler Number = ',q*mnssum^2/(mnssum + vel*mntsum))

figure(6)
loglog(tval,Wsum)
xlabel('Time (sec)')
ylabel('RMS Displacement')
title('Averaged RMS displacement')

%freq = DopFreq*exp(xp);
figure(7)
loglog(frq,ssum,frq,Ssum)
xlabel('Frequency (Hz)')
title('Average Spectrum')
legend('Het','Hom')

figure(8)
%plot(xtmp,hsum(1,:),xtmp,hsum(2,:),xtmp,hsum(3,:),xtmp,hsum(4,:),xtmp,hsum(5,:))
semilogy(xtmp,hsum(1,:),xtmp,hsum(2,:),xtmp,hsum(3,:),xtmp,hsum(4,:),xtmp,hsum(5,:))
title('Particle Density Function (PDF)')
xlabel('Position/vt')
ylabel('Density of Particles')

figure(9)
for tloop = 1:5
    fac = tindex(tloop)/NumT;
    semilogy(xtmp*fac,hsum(tloop,:)/fac)
    HSumx(:,tloop) = xtmp*fac;
    HSumy(:,tloop) = hsum(tloop,:)/fac;
    hold on
end
hold off
title('Particle Density Function (PDF)')
xlabel('Position/vt')
ylabel('Density of Particles')

figure(10)
for tloop = 1:5
    fac = tindex(tloop)/NumT;
    mnhsum = 0.5*(hsum(tloop,51:100)/fac + reversevec(hsum(tloop,1:50)/fac));
    loglog(xtmp(51:100)*fac,mnhsum)
    hold on
end
hold off
xlabel('Position')
ylabel('Density of Particles')

figure(11)
plot(xtmp2,hsum2(1,:),xtmp2,hsum2(2,:),xtmp2,hsum2(3,:),xtmp2,hsum2(4,:),xtmp2,hsum2(5,:))
title('Non-ballistic PDF')

if topt == 1
    nam = strcat(nam,'-b',num2str(alpha));
end
if sopt == 1
    nam = strcat(nam,'-a',num2str(beta));
end

ind = 0;
for loop = 1:5
    eps = 1e-4;
    fac = tindex(loop)/NumT;
    
    ind = ind + 1;
    tmp = xtmp*fac;
    hsump(:,ind) = tmp';
    
    ind = ind + 1;
    tmp = hsum(loop,:)/fac + eps;
    hsump(:,ind) = tmp';
end

savefile = strcat(nam,'.mat');
save(savefile,'topt','sopt','tau_w','tau_p','alpha','beta','tval','Wsum','freq','Ssum','xtmp','hsump','vel','vel_dist')

pp1 = hphse(Nphase/2+1:Nphase);
pp2 = 0.5*(phsum(Nphase/2+1:Nphase) + reversevec(phsum(1:Nphase/2))) + 1e-6;

Printfile2(strcat(nam,'-MSD.txt'),tval,Wsum)  % MSD
Printfile3(strcat(nam,'-Spec.txt'),frq,ssum,Ssum)  % Spectrum
Printfile2(strcat(nam,'-Phase.txt'),pp1,pp2)  % Phase

hom = pwd;
cd('anomtest')
fid = fopen(strcat(nam0,'subtest3_Log.txt'),'w');
fprintf(fid,'%s\n','subtest3Log');
fprintf(fid,'%s\n','Logged from subtest3.m');
datein = date;
fprintf(fid,'%s\n',datein);
clockin = clock;
fprintf(fid,'%s\n',strcat(num2str(clockin(4)),':',num2str(clockin(5))));
fprintf(fid,'%s\n',strcat('duration = ',num2str(duration/3600)));
fprintf(fid,'%s\n',' ');
fprintf(fid,'%s\n',strcat('Npart = ',num2str(Npart)));
fprintf(fid,'%s\n',strcat('Nsuper = ',num2str(Nsuper)));
fprintf(fid,'%s\n',strcat('Nang = ',num2str(Nang)));
fprintf(fid,'%s\n',strcat('vel_dist = ',num2str(vel_dist)));
fprintf(fid,'%s\n',' ');
fprintf(fid,'%s\n',strcat('topt = ',num2str(topt)));
if topt ~= 0
    fprintf(fid,'%s\n',strcat('alpha = ',num2str(alpha)));
end
fprintf(fid,'%s\n',strcat('tau_w = ',num2str(tau_w)));
fprintf(fid,'%s\n',' ');
fprintf(fid,'%s\n',strcat('sopt = ',num2str(sopt)));
if sopt ~= 0
    fprintf(fid,'%s\n',strcat('beta = ',num2str(beta)));
end
fprintf(fid,'%s\n',strcat('tau_p = ',num2str(tau_p)));
fprintf(fid,'%s\n',' ');
fprintf(fid,'%s\n',strcat('DopFreq = ',num2str(DopFreq)));
fprintf(fid,'%s\n',strcat('DopNum = ',num2str(DopNum)));
fprintf(fid,'%s\n',strcat('mean s = ',num2str(mnssum)));
fprintf(fid,'%s\n',strcat('mean t = ',num2str(mntsum)));
fprintf(fid,'%s\n',strcat('Eff Dop/# = ',num2str(q*mnssum^2/(mnssum + vel*mntsum))));
fprintf(fid,'%s\n',strcat('lp (persistence length) = ',num2str(lp)));
fprintf(fid,'%s\n',strcat('lam = ',num2str(lam)));
fprintf(fid,'%s\n',strcat('vel = ',num2str(vel)));
fprintf(fid,'%s\n',strcat('delT = ',num2str(delT)));
fprintf(fid,'%s\n',strcat('NumT = ',num2str(NumT)));
fprintf(fid,'%s\n',strcat('N = ',num2str(N)));
fclose(fid);
cd(hom)

if strcmp(compu,'MACI64')
    chirpData = load('chirp.mat');
    chirpObj = audioplayer(chirpData.y,chirpData.Fs);
    playblocking(chirpObj);
end

xout = frq';
yout = ssum';

xout2 = pp1';
yout2 = pp2';