%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% To be used for processing data from GIS Simulator v.1.5
% Vinzenz Friedli, EMPA, Switzerland
% June 2008
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read simulation parameter from <..>Simulation_Summary.txt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[fname,pname] = uigetfile('*.txt', 'Select a <...>Simulation_Summary.txt file');
bname = fname(1:end-22);
if fname == 0
    disp([' ']);
    disp('Nothing selected!');
    disp([' ']);
elseif length(fname)<22
    disp([' ']);
    disp('Select a <...>Simulation_Summary.txt file!');
    disp([' ']);
elseif sum(fname(end-21:end)~='Simulation_Summary.txt')
    disp([' ']);
    disp('Select a <...>Simulation_Summary.txt file!');
    disp([' ']);
else
    fid = fopen([pname '/' fname]);

    fgetl(fid); fgetl(fid); fgetl(fid); fgetl(fid); fgetl(fid); fgetl(fid); fgetl(fid); fgetl(fid);

    fgetl(fid); config_set = fscanf(fid,'%f');
    for i = 1:config_set
        fgetl(fid); angle1(i) = fscanf(fid,'%f');
        fgetl(fid); angle2(i) = fscanf(fid,'%f');
        fgetl(fid); distance(i) = fscanf(fid,'%f');
    end
    clear i;

    fgetl(fid); tube_length = fscanf(fid,'%f');
    fgetl(fid); exit_tilt = fscanf(fid,'%f');
    fgetl(fid); ra = fscanf(fid,'%f');
    fgetl(fid); check_hole = fscanf(fid,'%f');
    fgetl(fid); hole_z = fscanf(fid,'%f');
    fgetl(fid); hole_tilt = fscanf(fid,'%f');
    fgetl(fid); hole_radius = fscanf(fid,'%f');
    fgetl(fid); hole_TC = fscanf(fid,'%f');
    fgetl(fid); check_free = fscanf(fid,'%f');
    fgetl(fid); MFP_in = fscanf(fid,'%f');
    fgetl(fid); check_decay = fscanf(fid,'%f');
    if(check_free == 1)
        tmp_check = 1;
    elseif(check_free == 0 || check_decay == 0)
        tmp_check = 2;
    elseif(check_free == 0 || check_decay == 1)
        tmp_check = 3;
    end
    % tmp_check: 1 molecular flow, 2 transient flow, 3 transient flow with pressure decay
    clear tmp_check;

    fgetl(fid); DeltaP = fscanf(fid,'%f');
    fgetl(fid); uptake = fscanf(fid,'%f');
    fgetl(fid); gridsize = fscanf(fid,'%f');
    fgetl(fid); imagesize = fscanf(fid,'%f');
    fgetl(fid); virtual_planes = fscanf(fid,'%f');
    fgetl(fid); AngleDistrCenterX = fscanf(fid,'%f');
    fgetl(fid); AngleDistrCenterY = fscanf(fid,'%f');
    fgetl(fid); filtersize = fscanf(fid,'%f');
    fgetl(fid); ang_res = fscanf(fid,'%f');

    fclose(fid);
    
    disp(' ');
    disp(['"' fname '" loaded']);
    disp([' ']);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read data from files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (exist([pname '/' bname 'Configuration.txt']) ~= 2)
    disp(' ');
    disp('No data available. Make sure the specified directory exists.');
    disp(' ');
else
    
    % reads <..>Configuration.txt
    [Config_xs Config_ys Config_zs] = textread([pname '/' bname 'Configuration.txt'], '%s %s %s', 'whitespace',' ');

    i=1;
    l=1;
    m=1;
    for k = 1:length(Config_xs)
        if m <= length(Config_xs)
            if Config_xs{m} ~= '%'
                Config_x2(i,l) = str2num(Config_xs{m});
                Config_z2(i,l) = str2num(Config_zs{m});
                l = l+1;
            else
                i = i+1;
                l = 1;
            end
            m = m+1;
        end
    end
    no_disp = i-1;

    for i  = 1:config_set
        
        % reads <..>Image_Config<..>.txt        
        [n_px] = textread([pname '/' bname 'Image_Config' num2str(i) '.txt']);
        n_px = n_px(:,[1:(end-1)]); %skips the last row (inserted by textread!?)
        n_px_stack(:,:,i) = n_px;
        
        % reads <..>Profiles_Config<..>.txt
        [xs_pos(:,i) Counts(:,i) CountsCoeff(:,i) Sterr(:,i) SterrRel(:,i) FluxCoeff(:,i)] ...
            = textread([pname '/' bname 'Profiles_Config' num2str(i) '.txt'], '%f %f %f %f %f %f', 'whitespace',' ', 'headerlines', 1);
        
        % reads <..>Statistics_export_Config<..>.txt
        [TR TRH BS UT TC MC MolNum LIS(i) LOS(i) TR_LMC TR_LTC TR_DFT DB(i)] ...
            = textread([pname '/' bname 'Statistics_export_Config' num2str(i) '.txt'], '%f %f %f %f %f %f %f %f %f %f %f %f %f', 'whitespace',' ', 'headerlines', 1);

    end
    % reads <..>AngularDistributionHistogram.txt
    [HistoAngle zenithLMC azimuthLMC zenithLTC azimuthLTC zenithDFT azimuthDFT zenithAll azimuthAll] = textread([pname '/' bname 'AngularDistributionHistogram.txt'], '%f %u %u %u %u %u %u %u %u', 'whitespace',' ', 'headerlines', 1);
    if tube_length > 0
        % reads <..>PressureDistribution.txt
        [z_pos TRflux BSflux Netflux BSfluxNorm MFPsamplingNorm MFPsectionCount TubeCollCounter CollisionCounter] = textread([pname '/' bname 'PressureDistribution.txt'], '%f %f %f %f %f %f %f %f %f', 'whitespace',' ', 'headerlines', 1);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plot preparation
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %plot global settings

    color = 'black';
    width = 1;

    R_i = input('Enter the inner tube radius in mm (for absolute scale): ');
    disp([' ']);
    R_a = ra*R_i;

    for i  = 1:config_set

        xs = R_i*([-imagesize/2:imagesize/(gridsize-1):imagesize/2]);
        ys = xs;
        alpha1 = angle1(i)*pi/180;
        alpha2 = angle2(i)*pi/180;
        alpha = alpha1 + alpha2;
        beta = exit_tilt*pi/180;
        sigma = hole_tilt*pi/180;

        % Normal projection of tube inner diameter
        xs_di = [-1:0.001:1]*R_i/cos(beta)*cos(alpha);
        ys_di = R_i*sqrt(1-xs_di.^2/(R_i/cos(beta)*cos(alpha))^2);
        xs_di = xs_di - distance(i)*R_i*sin(alpha2); %ellipse shift to projection center

        % Normal projection of tube outer diameter
        xs_da = [-1:0.001:1]*R_a/cos(beta)*cos(alpha);
        ys_da = R_a*sqrt(1-xs_da.^2/(R_a/cos(beta)*cos(alpha))^2);
        xs_da = xs_da - distance(i)*R_i*sin(alpha2); %ellipse shift to projection center

        % Normal projection of hole
        xs_dh = [-1:0.001:1]*hole_radius*R_i/cos(sigma)*sin(pi/2-alpha-beta+sigma);
        ys_dh = hole_radius*R_i*sqrt(1-xs_dh.^2/(hole_radius*R_i/cos(sigma)*sin(pi/2-alpha-beta+sigma))^2);
        xs_dh = xs_dh - (tube_length-hole_z)*R_i*sin(alpha+beta) - distance(i)*R_i*sin(alpha) + R_i/cos(sigma)*sin(pi/2-alpha-beta+sigma); %ellipse shift to projection center

        xs_pos_center2 = - distance(i)*R_i*sin(alpha2)*[1 1];
        ys_pos = [-1/2 1/2]*R_i*imagesize;
        ys_pos_short = [-1.1 1.1]*R_a;
        x_axis = [-1/2 1/2]*R_i*imagesize;

        %44444444444444444444444444444444444444444444444444444444444
        HistRes = ang_res;
        max_zenith = max(zenithAll)/sum(zenithAll)/HistRes;
        %CosineSum = sum(sin(2*HistoAngle(find(HistoAngle<=90))*pi/180));
        CosineSum = 180/pi/HistRes;

        %555555555555555555555555555555555555555555555555555555555555
        zenithAll_projected = zenithAll/2./sin(HistoAngle*pi/180);
        

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % save data
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        fnamesave = [pname '/' bname 'Config' num2str(i) '_workspace.mat'];
        save(fnamesave);
        disp([bname 'Config' num2str(i) '_workspace.mat saved!'])
        disp([' ']);

        yn = input('Save projected tube data? (y/n): ', 's');
        disp([' ']);
        if yn == 'y'
            out = [[xs_pos_center2 xs_pos_center2(1) xs_di fliplr(xs_di) xs_di(1) xs_da(1) xs_da fliplr(xs_da) xs_da(1) xs_da(end)]'...
            [ys_pos_short 0 ys_di fliplr(-ys_di) 0 0 ys_da fliplr(-ys_da) 0 0]'];
            fnamesave = [pname '/' bname 'Config' num2str(i) '_ProjectedTubeData.txt'];
            save(fnamesave,'out','-ASCII')
            clear out;
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plot data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%% Data of several configurations is plotted within one figure %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %DISPLAY PLOT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %1 1 1
    figure(1);
    plot(Config_x2'*R_i, Config_z2'*R_i, '-', 'Linewidth', 2);
    axis image;
    set(gca, 'Box', 'on', 'FontSize', 10)
    ylabel('z_s [mm]', 'FontSize', 12);
    xlabel('x_s [mm]', 'FontSize', 12);

%     yn = input('Save Config plot? (y/n): ', 's');
%     disp([' ']);
%     if yn == 'y'
%         figure(1);
%         saveas(gcf,[pname '/' bname 'Config.jpeg'])
%     end

    % PROFILE PLOTS %%%%%%%%%%%%%%%%%%%%%%%%

    %11 11 11
    figure(11);
    hold on;
    for i  = 1:config_set
        errorbar(xs_pos(:,i)*R_i, FluxCoeff(:,i), Sterr(:,i)*pi/(imagesize/gridsize)^2, '-k');
    end
    title(['Profile along y_s = 0 (px size ' num2str(1e3*imagesize/gridsize*R_i,3) 'x' num2str(1e3*imagesize/gridsize*R_i,3) ' \mum^2)']);
    set(gca, 'Box', 'on', 'FontSize', 10)
    ylabel('Flux J/J_{tot}', 'FontSize', 12);
    xlabel('x_s [mm]', 'FontSize', 12);
    hold off;

    %12 12 12
    figure(12);
    hold on
    for i  = 1:config_set
        plot(xs_pos(:,i)*R_i, FluxCoeff(:,i), '.');
        Angle1 = angle1(i)*pi/180;
        Angle2 = angle2(i)*pi/180;
        Angle = Angle1+Angle2;
        if distance(i) ~= 0 && exit_tilt == 0
            for ind = 1:length(xs_pos(:,i))
                if xs_pos(ind,i) == 0
                    sigma(ind) = 0;
                else
                    Delta = (distance(i)*cos(Angle2))/cos(Angle);
                    sigma(ind) = atan(cos(Angle)./(Delta./xs_pos(ind,i)+sin(Angle)));
                end
            end
            I_x = (cos(sigma+Angle)/cos(Angle)).^2.*cos(sigma).*cos(sigma+Angle);
            I_x(find(I_x < 0)) = 0;
            plot(xs_pos(:,i)*R_i + distance(i)*(cos(Angle2)*tan(Angle)-sin(Angle2))*R_i, I_x/max(I_x)*max(FluxCoeff(:,i)));
        end
    end
    set(gca, 'Box', 'on', 'FontSize', 10)
    title(['Profile along y_s = 0 (px size ' num2str(1e3*imagesize/gridsize*R_i,3) 'x' num2str(1e3*imagesize/gridsize*R_i,3) ' \mum^2)']);
    ylabel('Flux J/J_{tot}', 'FontSize', 12);
    xlabel('x_s [mm]', 'FontSize', 12);
    legend('Simulation','Cosine point source \newline (normed to max value)');
    hold off

%     yn = input('Save Profile plot? (y/n): ', 's');
%     disp([' ']);
%     if yn == 'y'
%       figure(11);
%       saveas(gcf,[pname '/' bname 'Profile.jpeg'])
%     end
%     

    
    % 2D AND 3D PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%% Data of several configurations is plotted in separated figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for Config_ind = 1:config_set

        load([pname '/' bname 'Config' num2str(Config_ind) '_workspace.mat']);

        % 10X
        figure(100+Config_ind);
        hold on;
        contour_lines = 8;
        contour(xs, ys, n_px_stack(:,:,(Config_ind))/TR*pi/(imagesize/gridsize)^2, contour_lines); colormap('jet');
        plot([0 0], ys_pos, color, 'Linewidth', width);
        plot(xs_di, ys_di, [':' color], 'Linewidth', width); plot(xs_di, -ys_di, [':' color], 'Linewidth', width);
        plot(xs_da, ys_da, color, 'Linewidth', width); plot(xs_da, -ys_da, color, 'Linewidth', width);
        if check_hole(1) == 1
            plot(xs_dh, ys_dh, color, 'Linewidth', width); plot(xs_dh, -ys_dh, color, 'Linewidth', width);
        end
        plot(xs_pos_center2, ys_pos_short, color, 'Linewidth', width)
        plot(x_axis, [0 0], color, 'Linewidth', width)
        axis image
        set(gca, 'Box', 'on', 'FontSize', 10)
        title(['Isoflux contour plot of J/J_{tot} (px size ' num2str(1e3*imagesize/gridsize*R_i,3) 'x' num2str(1e3*imagesize/gridsize*R_i,3)...
            ' \mum^2) \newline \alpha_1 = ' num2str(angle1(Config_ind)) '°, \alpha_2 = ' num2str(angle2(Config_ind)) '°, d = ' num2str(distance(Config_ind)) 'r'])
        xlabel('x_s [mm]', 'FontSize', 12)
        ylabel('y_s [mm]', 'FontSize', 12)
        colorbar
        hold off;

%         yn = input('Save 2D plot? (y/n): ', 's');
%         disp([' ']);
%         if yn == 'y'
%             figure(100+Config_ind);
%             saveas(gcf,pname '/' bname '2D.jpeg'])
%         end

    end


    % ANGULAR PROBABILITY PLOTS %%%%%%%%%%%%%%%%%

    %1001 1001 1001
    figure(1001);
    title('Zenith angular probability distribution \newline (solid angle)');
    polar(HistoAngle*pi/180, zenithAll/sum(zenithAll)*CosineSum,'.k');
    set(gca, 'Box', 'on', 'FontSize', 10)
    axis image;

    %1002 1002 1002
    figure(1002);
    title('Zenith angular probability distribution \newline (2D projection)');
    polar(HistoAngle*pi/180, zenithAll_projected/max(zenithAll_projected),'.k');
    axis image;
    set(gca, 'Box', 'on', 'FontSize', 10)

%     yn = input('Save angular probability plot? (y/n): ', 's');
%     disp([' ']);
%     if yn == 'y'
%       figure(1001);
%       saveas(gcf,[pname '/' bname 'AngProb.jpeg'])
%     end

    
    % PRESSURE PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%

    if tube_length > 0

        %10001 10001 10001
        figure(10001);
        hold on
        plot(z_pos*R_i, BSflux, '.black');
        plot(z_pos*R_i, TRflux, '.red');
        plot(z_pos*R_i, Netflux, '.blue');
        hold off
        set(gca, 'Box', 'on', 'FontSize', 10)
        legend('Backscatter flux','Transmission flux', 'Net flux');
        xlabel('x [mm]', 'FontSize', 12);
        ylabel('Pressure (normed)', 'FontSize', 12);

%         yn = input('Save Pressure plot? (y/n): ', 's');
%         disp([' ']);
%         if yn == 'y'
%             figure(10001);
%             saveas(gcf,handles.output_dname '/' handles.output_fname 'StaticFlux.jpeg'])
%         end
    end


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    yn = input('Close figures? (y/n): ', 's');
    if yn == 'y'
        close all;
    end

end
