%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;

fprintf(1, '\n\n== start == \n\n');

data_dir  = '../raw/';
%data_set  = 'spr_phi090_el7'

filestub =  'crd';
frames   =  999;
    
%%%%%%%%%%%%%%%%%%%%

%for run = ['B','D','J','L','S','T'];
for run = ['L', 'S', 'T' ]
    if run=='B'
        data_set     = 'B_frames'
        aspect_ratio = 9;
        filling      = 0.55;
        u0           = 250;
        margin       = 200;
    end
    
    if run=='D'
        data_set     = 'D_frames'
        aspect_ratio = 4;
        filling      = 0.063;
        u0           = 250;
        margin       = 300;
    end
    
    if run=='J'
        data_set     = 'J_frames'
        aspect_ratio = 3;
        filling      = 0.97;
        u0           = 250;
        margin       = 200;
    end
    
    if run=='L'
        data_set     = 'L_frames'
        aspect_ratio = 16;
        filling      = 0.53;
        u0           = 250;
        margin       = 300;
    end
    
    if run=='S'
        data_set     = 'S_frames'
        aspect_ratio = 16;
        filling      = 0.21;
        u0           = 250;
        margin       = 200;
    end
    
    if run=='T'
        data_set     = 'T_frames'
        aspect_ratio = 5;
        filling      = 0.84;
        u0           = 250;
        margin       = 200;
    end
    
    
    
    % contrast/figure settings
    
    total_scale = 0.045*aspect_ratio/(filling);
    margin = round(250);                
    
    %%%%%%%%%%%%%%%%%%%%
    
    
   
    input_dir  = [data_dir, data_set, '/'];
    
    % generate output directory
    output_dir = [data_set, '/'];
    %rmdir(output_dir);
    %mkdir(output_dir);
    
    
    
    
    for f= 1:frames
        
        fprintf(1, 'frame = %d \n', f);
        
        % load data
        
        if f<1000
            filename = [filestub,'_' num2str(f),'.dat'];
        end
        if f<100
            filename = [filestub, '_0', num2str(f),'.dat'];
        end
        if f<10
            filename = [filestub, '_00', num2str(f),'.dat'];
        end
        
        file =[input_dir , filename];
        A = load(file);
        
        [particles, cols] = size(A);
        
        
        
        X = A(:,1);
        Y = A(:,2);
        
        Nx = A(:,3);
        Ny = A(:,4);
        
        %   Vx = A(:,5);
        %   Vy = A(:,6);
        
        
        % compute intensity profile
        
        pix = round(sqrt(particles));
        pix = 1000;
        
        Vol          = pix * pix;       % total volume in pixels
        vol_per_part = filling * Vol / particles;   % total volume in pixels
        
        intensity = zeros( pix, pix );
        
        cut_off =0.2;
        scale  = 1 - 2 * cut_off;
        
        for p=1:particles
            
            if (X(p,1) > cut_off &&  X(p,1) < 1-cut_off )
                if (Y(p,1) > cut_off &&  Y(p,1) < 1-cut_off )
                    
                    % positions in pixels
                    x = (X(p,1)-cut_off)/scale * pix;
                    y = (Y(p,1)-cut_off)/scale * pix;
                    
                    nx = Nx(p,1) ;
                    ny = Ny(p,1) ;
                    
                    % nearest grid positions
                    xpos = round(x);
                    ypos = round(y);
                    
                    xpos = max(xpos, 1);
                    ypos = max(ypos, 1);
                    
                    xpos = min(xpos, pix);
                    ypos = min(ypos, pix);
                    
                    
                    % explore region around
                    
                    scan=2.0;
                    x1 = round(xpos - scan * aspect_ratio * nx * sqrt(vol_per_part));
                    x2 = round(xpos + scan * aspect_ratio * nx * sqrt(vol_per_part));
                    
                    y1 = round(ypos - scan * aspect_ratio * ny * sqrt(vol_per_part));
                    y2 = round(ypos + scan * aspect_ratio * ny * sqrt(vol_per_part));
                    
                    xmin = max( min(x1,x2), 3);
                    xmax = min( max(x1,x2), pix - 2);
                    
                    ymin = max( min(y1,y2), 3);
                    ymax = min( max(y1,y2), pix - 2);
                    
                    
                    for i= xmin-2 : xmax+2
                        for j= ymin-2 : ymax+2
                            
                            para2  = (nx*(x-i) + ny*(y-j))^2;
                            perp2  = (x-i)^2+ (y-j)^2 - para2;
                            scale2 = (aspect_ratio)^2;
                            
                            delta_intensity = u0*exp( - total_scale*( (para2/scale2)  + perp2) );
                            
                            intensity( i, j ) = intensity( i, j ) + delta_intensity;
                        end
                    end
                end
                
                
            end
        end
        
        
        subwindow = intensity(margin:round(pix-margin),margin: round(pix-margin));
        %subwindow = intensity;
        
        % generate figure
        figure(f);
        colormap('Gray');
        
        ip = pcolor(subwindow);
        
        set(ip, 'EdgeColor', 'none');
        set(gca,'YTick', [],'YTickLabel', [],'XTick', [],'XTickLabel', []);
        
        pbaspect([1 1 1]);
        
        set(gcf, 'PaperUnits', 'points');
        set(gcf, 'PaperSize', [pix pix]);
        set(gcf, 'PaperPositionMode', 'manual');
        set(gcf, 'PaperPosition', [0 0 pix pix]);
        
        set(gca, 'Position', get(gca, 'OuterPosition') - ...
            get(gca, 'TightInset') * [-1 0 1 0; 0 -1 0 1; 0 0 1 0; 0 0 0 1]);
        
        
        
        saveas(figure(f), [output_dir, 'f_', sprintf('%04.0f', f), '.tif']);
        close(figure(f));
        
        
    end
    
end
fprintf(1, '\n\n== end == \n\n');

