clear all
close all
clc

load('Figure8.mat');

%% goal for this figure is to give density and phase profiles for two spiral beams...  
%w0 = 8 um.
%N = 2e6
%tau for spiral is 4 s
%N for spiral is 3.

%% Constants

a = 100.44*5.29e-11;           % approx scattering length for Rb87 in meters.
m = 1.5e-25;                % m = mass of Rb87 in kg
hbar = 1.05457e-34;         % Planck's constant in J*s


%
%   Set up grids and time steps
%   All units are in mks.
%

%% space
%
xmax    = 120e-6;             % x grid size in meters
ymax    = xmax;
nxy      = 256;             % number of grid points
nmid  = floor(nxy/2);
v = 0:1:nxy-1;
dx = xmax/(nxy-1);
dy = ymax/(nxy-1);
[x,y] = meshgrid(v,v);
x = -xmax/2 + x*dx;
y = -ymax/2 + y*dy;
xar = -xmax/2 + v*dx;
yar = -ymax/2 + v*dy;

%% time
%
tspiral = ts;                %time for spiral
dt = t(2)-t(1)
nt = length(t);
ncycle = floor(nt/500);

%% make figure
%

[sy sx sz] = size(mov2D_psi);
 
% choose times to plot
ttt = [.5, 1, 1.55, 3]/tspiral;

% choosepts = 0;
%     x1_ar = 1.0e-04 * [0.4    0.15   -0.2   -0.25];
%     y1_ar = 1.0e-04 * [-.01    0.25   -0.25    0.1];

qctr = 0;   

%set up figure
Bshift = -20;
Bshift1 = -10;
Ashift = 5;  %shift for axis  old version was Ashift = 5
Lshift = 45;

%fontsizes
f1 = 21;  % fontsize for timing labels.
f2 = 20;  % fontsize for colormap numbers.   
f3 = 26;  % fontsize for colormap labels.    

h3 = figure(103); 
set(h3, 'position', [20        20    130*length(ttt)+5+Lshift 130*2]);
set(h3,'color','w');

for qq = round(ttt*ts/dt/ncycle)+1    
    qctr = qctr+1;

    % atomic density profile
    n =  abs(mov2D_psi(:,:,qq)).^2;   
    nmax = max(max(n));
    
    % create masks to clean up the edges where the atom density is very low.
    phmask = ones(size(n)).*(n/nmax > 5e-10).*(x.^2 + y.^2 < (1.05*RrTF)^2);
    nmask = ones(size(n)).*(n/nmax > 5e-4); 
    phmask2 = ones(size(n)).*(n/nmax > 5e-5);

    % now make black background for densities less than 5e-4*nmax 
    n = n.*nmask;
    n(n == 0) = NaN;
    imAlpha=ones(size(n));
    imAlpha(isnan(n))=0;

    % now make black background for phase ...
    ph = (angle(mov2D_psi(:,:,qq))+pi()).*phmask2;
    ph2 = ph.*phmask;
    ph2(ph2 == 0) = NaN;
    imAlpha1=ones(size(ph2));
    imAlpha1(isnan(ph2))=0;


    %make figure      
    figure(103)

    % density profile      
    axes('Units', 'pixels', 'position', [5+130*(qctr-1), 20+130+5+Bshift, 125, 125]) % 121 is for example
    imagesc(xar,yar, n/nmax,'AlphaData',imAlpha)    
    set(gca,'color',0*[1 1 1]);
    set(gca, 'CLim', [0 1])
    axis image
    axis xy
    set(gca,'XTickLabel','')
    set(gca,'YTickLabel','')

    hold on

    str = ['\bf{$t = ' num2str(round(ttt(qctr)*tspiral*100)/100) '\, \mathrm{s}$}'];

    text(xar(4), yar(sy-24),str, 'Interpreter', 'latex', 'fontsize', f1, 'color', 'white')
    hold off

    % add label
    if qctr == 1
        cax = colorbar('EastOutside');
        set(cax, 'Units', 'pixels', 'position', [Ashift+length(ttt)*130, 20+130+5+Bshift1, 10, 105], 'fontsize', f2, 'XTick', 0:1:1);
        cax.Label.Interpreter = 'latex';
        cax.TickLabelInterpreter = 'latex';
        cax.Label.String = '{$n/n_\mathrm{max}$}';
        cax.Label.Position = [4   0.5   0];   
        set(cax.Label,'Rotation',270);
        cax.Label.FontSize = f3;
    end

    % corresponding phase profile
    axes('Units', 'pixels', 'position', [5+130*(qctr-1), 20+5+Bshift, 125, 125]) 
    imagesc(ph2/2/pi,'AlphaData',imAlpha1)
    set(gca, 'CLim', [0 1])
    axis xy
    axis image
    set(gca,'color',0*[1 1 1]);
    set(gca,'XTickLabel','')
    set(gca,'YTickLabel','')

    % add labels
    if qctr == 1
        cax = colorbar('EastOutside');
        set(cax, 'Units', 'pixels', 'position', [Ashift+length(ttt)*130, 20+5+Bshift1, 10, 105], 'fontsize', f2, 'XTick', 0:1:1);
        cax.Label.Interpreter = 'latex';
        cax.TickLabelInterpreter = 'latex';
        cax.Label.String = '{$\phi/2\pi$}';
        cax.Label.Position = [3.9    0.5000         0];
        set(cax.Label,'Rotation',-90);
        cax.Label.FontSize = f3;
    end

end


%%
% print(h3,'-dpng', ['Figure_SimsImages_v0'])
% print(h3,'-dsvg', ['Figure_SimsImages_v0'])
% print(h3,'-dpdf', ['Figure_SimsImages_v0'])