

close all;  clear all;


% change to your file name
fileID = fopen('p5.prof');

%fileID = fopen('p7RSm.prof');




% 
% read in the Exported ANSYS profile
%
Header = textscan(fileID,'%s point %d)',1,'Delimiter',' \n')
while (~feof(fileID))                               % For each block:                         
    str = textscan(fileID, '(%s'    );  VarName = str{1};
    p= textscan(fileID, '%f'    );  
    
    if( strcmp ( VarName , 'x') ) 
        x = p{1};
    elseif( strcmp ( VarName , 'y') )
       y = p{1};
    elseif( strcmp ( VarName , 'z') )
       z = p{1};
    elseif( strcmp ( VarName , 'turb-kinetic-energy') )
       k = p{1};
    elseif( strcmp ( VarName , 'uu-reynolds-stress') )
       uu = p{1};
    elseif( strcmp ( VarName , 'uv-reynolds-stress') )
          uv=p{1};
    elseif( strcmp ( VarName , 'vv-reynolds-stress') )
          vv=p{1};
    elseif( strcmp ( VarName , 'ww-reynolds-stress') )
          ww=p{1};
    elseif( strcmp ( VarName , 'dx-velocity-dx') )
          dudx=p{1};
    elseif( strcmp ( VarName , 'dy-velocity-dy') )
          dvdy=p{1};
    elseif( strcmp ( VarName , 'dy-velocity-dx') )
          dvdx=p{1};
    elseif( strcmp ( VarName , 'dx-velocity-dy') )
          dudy=p{1};
    elseif( strcmp ( VarName , 'viscosity-turb') )
          nu_T=p{1};
    elseif( strcmp ( VarName , 'density') )
          rho=p{1};
    end
    p= textscan(fileID, ')'    );  
end 
fclose (fileID)





figure(1)
set(gcf,'Units','inches','Position',[0 4 6 4]);
plot( x, k, 'linewidth', 1); hold on;
title('k vs. x');

figure (2);




%% Example of single Rij
%Rij =[ 1, 0, 0 ;   0, 0 0 ;  0, 0,0] ;   % example : 2 component
%Rij =[ 1, 0, 0 ;   0, 1 0 ;  0, 0,0] ;   % example : 2 component
%Rij =[ 5, 0, 0 ;   1, 2 0 ;  1, 1, 3] ; % three component
%Rij =[ 1, 0, 0 ;   0, 1 0 ;  0, 0,1] ;  % isotropic
%Rij =[ 5, 0, 0 ;   1, 2 0 ;  1, 1, 3] ; % three component
%Rij =[ 1, 5, 0 ;   1, 2 0 ;  1, 1, 3] ; % not realizable viloating cauchy-schawz inequality, 
%% The inputting Rij must be symmetrical, copy the low triangular part of Rij to the upper trianglular part .
%Rij = Rij + triu(Rij', 1);   
%% the defnition of bij = (<ui'uj'>-2/3 k\delta_ij )/2k 
%bij = ( Rij - trace(Rij)/3 * eye(3)  ) /  trace(Rij); 


set(gcf,'Units','inches','Position',[6 4 6 4]);

%draw_LumelyTriangle_Border;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             set(0,'DefaultLineMarkerSize',10);
            set(groot,'defaultLineLineWidth',1.5)
            set(0,'DefaultTextFontSize', 15)
            set(0,'DefaultAxesFontSize',15);
             %set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 4]);
            % set(gcf,'Units','inches','Position',[1 1 6 4]);
            set(gca,'YMinorTick','on');
           set(gca,'XMinorTick','on');
            fontsize_lable = 18;
            box on;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % The following part draw the 3 borders of lumley triangle 
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % left border
            xx = [ -1/6   0];                yy = [ 1/6   0];                       plot ( xx, yy , '-k' ); hold on;
            % right border 
            xx = [  0  1/3];                 yy = [ 0   1/3];                        plot ( xx, yy , '-k' ); hold on;
            % top 2C border 
            xx =   -1/6 :0.001: 1/3;  yy = sqrt(1/27+2*xx.^3 );     plot ( xx, yy , '-k' ); hold on;

            % lable 
            hylabel=ylabel('$\eta $' , 'FontSize',  fontsize_lable,'interpreter','latex') ;  
            hxlabel=xlabel('$\xi  $' , 'FontSize',  fontsize_lable,'interpreter','latex') ;  
            ax = gca;                                   ax.YAxisLocation = 'origin'; 
            xticks([-1/6 0 1/6 1/3]);        xticklabels( {'-1/6' '0' '1/6',  '1/3' } );
            yticks([ 0 1/6 1/3]);                yticklabels( { '0' '1/6',  '1/3' } );
            title ( 'Lumley Triangle') ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%plot_bij_Inside_LumelyTriangle(bij); return;




% % % % Boussineq eddy-visocosity asssumptions
% % % % 
% % % %   <ui'uj'> - 2/3 k \delta_ij = 2  ( 0.09 *k^2/ e ) Sij       
% % % % 
% % % %       note 1 :    Vis_turb =  ( C_nu *k^2/ e ), C_nu = 0.09
% % % %       note 2:              S_ij =  ( d <u_i>/d x_j + d <u_j>/d x_i  )/2 
% % % %      
% % % %                             => b_ij =  0.09 * k/e*Sij
% % % % Define Sij_ke   as   k/e*Sij,  it must be symmetrical
%
% the following is demonstrate the plotting of 100 radomly generated "u_ij *k/epislon"
% for ii = 1: 100
%       uij_ke = rand (3,3)*3;  
%       uij_ke = uij_ke-  trace(uij_ke)/3 * eye(3); 
%       Sij_ke= (uij_ke + uij_ke.' )/2; 
%       bij = 0.09 *Sij_ke ; 
%       plot_bij_Inside_LumelyTriangle(bij);
%  end
%  return;

if exist('dudx','var') % k-equation , 
% Eddy-visocosity based on turbulence model: 
%                        k-epsilon, (realizble), k-omega, 3 eq. , ...
% Required input:  (1) turbulent viscosity 
%                                (2) rho 
%                                (3) turbulent kinetic energy 
%                               (4-7) dudx, dudy,dvdx,dvdy
            for m = 1: size( dudx, 1)
                 %uij_ke = rand (3,3)*3;  
                 uij = [ dudx(m), dudy(m), 0;  dvdx(m) , dvdy(m) , 0 ;    0  , 0  , 0 ] ;
                 % a converged fluent  simulation should guaranty zero trace
                 % uij = uij_ke-  trace(uij_ke)/3 * eye(3);  
                 Sij= (uij + uij.' )/2; 

                 bij = 2*nu_T(m)/rho(m)*Sij / (2 * k(m) ); 
                
                 figure(2) ; 
                 %plot_bij_onto_LumelyTriangle(bij,m);     
                            %  (Pope book, p.394, eq. 11.29) 6*xi^3 = 3III_b = bii^3 = b_ij*b_jk*b_ki ; 
                            xi     = nthroot(trace ( bij*bij*bij ) /6  , 3);
                            %  (Pope book, p.393, eq. 11.28) 6*eta^2 = -2II_b = bii^2 = b_ij*b_ji; 
                            eta  = sqrt(trace ( bij*bij ) /6 ) ;
                            plot ( xi, eta , 'o' ); hold on;                 
                 %%%%%%%%%%%%%%%%%%%%%%%%   
                 
                 figure(1);  plot( x(m) , k(m),'o' ) ; hold on;

                 pause(.02)
            end
            
else
% RSM models
% Required input:  k, uu,uv,vv,ww
            for m = 1: size( uu, 1)
                Rij = [ uu(m), uv(m), 0;  uv(m) , vv(m) , 0 ;    0  , 0  , ww(m) ] ;
                bij = ( Rij - trace(Rij)/3 * eye(3)  ) /  trace(Rij); 
                %note:  It should be verified that    " trace(Rij)/2  ==     k(m)"
                
                
                 figure(2) ;
                 %plot_bij_onto_LumelyTriangle(bij, m);
                            %  (Pope book, p.394, eq. 11.29) 6*xi^3 = 3III_b = bii^3 = b_ij*b_jk*b_ki ; 
                            xi     = nthroot(trace ( bij*bij*bij ) /6  , 3);
                            %  (Pope book, p.393, eq. 11.28) 6*eta^2 = -2II_b = bii^2 = b_ij*b_ji; 
                            eta  = sqrt(trace ( bij*bij ) /6 ) ;
                            plot ( xi, eta , 'o' ); hold on;                 
                %%%%%%%%%%%%%%%%%%%%    
                 
                 figure(1);    plot(  x(m) , k(m) , 'o' ) ; hold on;
                 
                pause(.02)
            end

end 


% 
% 
% 
% function [] = plot_bij_onto_LumelyTriangle ( bij , m )     
%             %  (Pope book, p.394, eq. 11.29) 6*xi^3 = 3III_b = bii^3 = b_ij*b_jk*b_ki ; 
%             xi     = nthroot(trace ( bij*bij*bij ) /6  , 3);
%             %  (Pope book, p.393, eq. 11.28) 6*eta^2 = -2II_b = bii^2 = b_ij*b_ji; 
%             eta  = sqrt(trace ( bij*bij ) /6 ) ;
%             plot ( xi, eta , 'o' ); hold on;
% end
% 
% 
% 
% function   draw_LumelyTriangle_Border()     
%         set(0,'DefaultLineMarkerSize',10);
%         set(groot,'defaultLineLineWidth',1.5)
%         set(0,'DefaultTextFontSize', 15)
%         set(0,'DefaultAxesFontSize',15);
%          %set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 4]);
%         % set(gcf,'Units','inches','Position',[1 1 6 4]);
%         set(gca,'YMinorTick','on');
%        set(gca,'XMinorTick','on');
%         fontsize_lable = 18;
%         box on;
% 
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         % The following part draw the 3 borders of lumley triangle 
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         % left border
%         x = [ -1/6   0];                y = [ 1/6   0];                       plot ( x, y , '-k' ); hold on;
%         % right border 
%         x = [  0  1/3];                 y = [ 0   1/3];                        plot ( x, y , '-k' ); hold on;
%         % top 2C border 
%         x =   -1/6 :0.001: 1/3;  y = sqrt(1/27+2*x.^3 );     plot ( x, y , '-k' ); hold on;
% 
%         % lable 
%         hylabel=ylabel('$\eta $' , 'FontSize',  fontsize_lable,'interpreter','latex') ;  
%         hxlabel=xlabel('$\xi  $' , 'FontSize',  fontsize_lable,'interpreter','latex') ;  
%         ax = gca;                                   ax.YAxisLocation = 'origin'; 
%         xticks([-1/6 0 1/6 1/3]);        xticklabels( {'-1/6' '0' '1/6',  '1/3' } );
%         yticks([ 0 1/6 1/3]);                yticklabels( { '0' '1/6',  '1/3' } );
%         title ( 'Lumley Triangle') ;
% end