
%%This script loads the self-similar solutions, as well as their respective
%%high temperature transport properties

clc

files = dir('C:\Users\sedm5652\Desktop\Matlab_Calculations\Hypersonic_BL\Theory_Journal_Paper\Journal_Paper_const_cw_calcs\const_cw*');

C = {'k','b','r','g','y','m','c',[.5 .6 .7],'k','b','r','g','y','m','c',[.5 .6 .7],[.5 .6 .1]}; % Cell array of colros.

for i=1:length(files)
    i
    cd(['C:\Users\sedm5652\Desktop\Matlab_Calculations\Hypersonic_BL\Theory_Journal_Paper\Journal_Paper_const_cw_calcs\',files(i).name])
    addpath('C:\Users\sedm5652\Desktop\Matlab_Calculations\Hypersonic_BL\')
    if i==1
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(1500,11927,'Air','Air');
    elseif i==2
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(1500,11927,'Air','Air');
    elseif i==3
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(2000,11927*5,'Air','Air');
    elseif i==4
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(2000,11927*5,'He','Air');
    elseif i==5
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'Air','Air');
    elseif i==6
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'Air','He');
    elseif i==7
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'Ar','Air');
    elseif i==8
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'Ar','Ar');
    elseif i==9
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'He','Air');
    elseif i==10
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'He','He');
    elseif i==11
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'Ne','Air');
    elseif i==12
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'Ne','Ne');
    elseif i==13
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'Xe','Air');   
    elseif i==14
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'Xe','Xe');   
    elseif i==15
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927,'Air','Air');
    elseif i==16
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(300,11927*2,'Air','Air');
    elseif i==17
        [alpha,delta,d,l,tau,Gamma,mue,Pre,Sce,rhoe,D12,Mu]=High_Temperature_Transport_Properties(515,11927,'N2','Air');
    end
    
    for conc= [0:5:100] %load the individual files and save them.
        
    Files=dir(fullfile(cd,'*.mat'));
    pieces = regexp(Files(1).name, '\_', 'split');  %divide by _
    coolant = pieces{1};
    load([coolant,'_sols_cw',num2str(conc),'.mat']);
    
    cw_val(i,conc/5+1) = soly(6,1); 
    dcw(i,conc/5+1) = -soly(7,1);
    ddcw(i,conc/5+1) = -diff(soly(7,1:2))/diff(solx(1:2));
    cw_alpha(i,conc/5+1) = soly(6,1)./alpha(soly(6,1),soly(4,1));
    alpha_w(i,conc/5+1) = alpha(soly(6,1),soly(4,1));
    l_w(i,conc/5+1) = l(soly(6,1),soly(4,1));
    dalpha_w(i,conc/5+1) = diff(alpha(soly(6,1:2),soly(4,1:2)))./diff(solx(1:2));
    Dcoff(i,conc/5+1) = D12(1,1);
    end
    
    %plots
    if i < 9
    figure(300)
    plot(cw_alpha(i,:)',dcw(i,:)','x','color',C{i}), hold on
    figure(400)
    plot(1./alpha_w(i,:)',dcw(i,:)','x','color',C{i}), hold on
    figure(500)
    plot(dalpha_w(i,:)',dcw(i,:)','x','color',C{i}), hold on
    else 
    figure(300)
    plot(cw_alpha(i,:)',dcw(i,:)','o','color',C{i}), hold on
    figure(400)
    plot(1./alpha_w(i,:)',dcw(i,:)','o','color',C{i}), hold on
    figure(500)
    plot(dalpha_w(i,:)',dcw(i,:)','o','color',C{i}), hold on
    end
end

figure(300)
xlabel('$c_w\big(\frac{Sc}{l}\big)_w$','interpreter','latex')
ylabel('$\big(\frac{\partial c}{\partial \eta}\big)_w$','interpreter','latex')
legend(files.name)

figure(400)
xlabel('$\big(\frac{Sc}{l}\big)_w$','interpreter','latex')
ylabel('$\big(\frac{\partial c}{\partial \eta}\big)_w$','interpreter','latex')
legend(files.name)

figure(500)
xlabel('$\partial \big(\frac{Sc}{l}\big)_w$','interpreter','latex')
ylabel('$\big(\frac{\partial c}{\partial \eta}\big)_w$','interpreter','latex')
legend(files.name)

%% Calculate some additional quantities

%Molar mass normalised by molar mass of external gas.
Minj = [29,29,29,4,29,29,40,40,4,4,20,20,131,131,29,29,29];
Mext = [29,29,29,29,29,4,29,40,29,4,29,20,29,131,29,29,29];
M = Minj./Mext.*ones(21,1).*cw_val' + Mext./Mext.*ones(21,1).*(1-cw_val')

%Diffusion coefficient normalised by Dcoff if coolant = freestream gas.
Dcoff_ref=Dcoff;
Dcoff_ref(7,:) = Dcoff_ref(5,:);
Dcoff_ref(9,:) = Dcoff_ref(5,:);
Dcoff_ref(11,:) = Dcoff_ref(5,:);
Dcoff_ref(13,:) = Dcoff_ref(5,:);
Dcoff_ref(4,:) = Dcoff_ref(3,:);
Dcoff_ref(6,:) = Dcoff_ref(10,:);

%Sc number if coolant = freestream gas.
Scw =l_w./alpha_w;
Scw(7,:) = Scw(5,:);
Scw(9,:) = Scw(5,:);
Scw(11,:) = Scw(5,:);
Scw(13,:) = Scw(5,:);
Scw(4,:) = Scw(3,:);

%Interesting relationship between the gradient of alpha and dcw.
i_sel = [4,5,6,7:2:13]; %cases which have different coolant than freestream gas
figure(510)
plot(dcw(i_sel,:)',alpha_w(i_sel,1)'./alpha_w(i_sel,:)','o'), hold on
K = [-0.66,0.159,0.621,-0.093,1e9,-0.063];
Mratio = [1.38,0.138,0.69,4.53,1,7.25];
figure(520)
plot((Mratio.^0.5-1),1./K,'x')

%See if this new relationship can be plugged into alpha of the old equation (it can't/doesn't collapse the data).
figure(530)
plot(0.63*cw_alpha',(0.63*cw_val./alpha_w(:,1))'./(1-(0.63*cw_val./alpha_w(:,1))'.*9.6.*(1-sqrt(Minj./Mext.*ones(21,1)))),'x'), hold on
plot(0.63*cw_alpha',dcw','o'), hold on

%This sum should be zero.
figure(540)
i_sel2 = [1:17]; %just the cases with different coolant and freestream gas
plot([0:0.05:1],ddcw(i_sel2,:)'.*alpha_w(i_sel2,:)' + dcw(i_sel2,:)'.*dalpha_w(i_sel2,:)','x')

%This figure plots shows the relationship between cw, alpha to dcw. Note
%that the MD ratio correction did not quite work for Helium
figure(550)
i_sel2 = [1,2,3,5,8,10,12,14,15,16,17]; %just the cases with same coolant and freestream gas
%i_sel2 = [1:17];
plot([0:0.05:1]',dcw(i_sel2,:)'.*alpha_w(i_sel2,:)'.*M(:,i_sel2).^-0.75.*(Dcoff'./Dcoff_ref').^-0.75.*Scw'.^0,'o') %this one is better!!

%same as previous figure, but applying the correction to x, instead of y.
figure(570)
i_sel2 = [1,2,3,5,8,10,12,14,15,16,17]; %just the cases with same coolant and freestream gas
i_sel2 = [1:17]; %all cases

p = polyfit([0:0.05:1]'./(alpha_w(5,:)'.*M(:,5).^-0.75.*(Dcoff(5,:)'./Dcoff_ref(5,:)').^-0.75),dcw(5,:)',1)
y1 = polyval(p,[0:0.05:1]'./(alpha_w(5,:)'.*M(:,5).^-0.75.*(Dcoff(5,:)'./Dcoff_ref(5,:)').^-0.75));
plot([0:0.05:1]'./(alpha_w(5,:)'.*M(:,5).^-0.75.*(Dcoff(5,:)'./Dcoff_ref(5,:)').^-0.75),y1(:,1),'--r'), hold on

plot([0:0.05:1]'./(alpha_w(i_sel2,:)'.*M(:,i_sel2).^-0.75.*(Dcoff(i_sel2,:)'./Dcoff_ref(i_sel2,:)').^-0.75),dcw(i_sel2,:)','kx') %this one is better!!
%xlabel('$c_w\big(\frac{Sc}{l}\big)_w$','interpreter','latex')
xlabel('$c_w\big(\frac{Sc}{l}\big)_w\bigg(\frac{M_\mathrm{ext}D_{\mathrm{ext;ext}}}{M D_{\mathrm{inj;ext}}}\bigg)_w^{-0.75}$','interpreter','latex')
ylabel('$\big(\frac{\partial c}{\partial \eta}\big)_w$','interpreter','latex')
% legend(files(i_sel2).name)
legend('Linear fit','Numerical solutions')
plot([0:0.05:1]'./(alpha_w(5,:)'.*M(:,5).^-0.75.*(Dcoff(5,:)'./Dcoff_ref(5,:)').^-0.75),y1(:,1),'--r'), hold on




%% Check what happens if B and numerator is different to B of denominator


m=0:0.001:1;


for Gamma=[0.1,1,10,100];

beta=(exp(Gamma*m)-1)/(exp(Gamma)-1)

plot(m,beta), hold on

end

legend('\Gamma=0.1','\Gamma=1','\Gamma=10','\Gamma=100')
xlabel('m')
ylabel('$\beta=\frac{e^{\Gamma m}-1}{e^{\Gamma}-1}$','interpreter', 'latex')

% B1=[0.01:0.1:10];
% B2=B1*1.1;
% 
% plot(B1,B1./(exp(B1)-1)), hold on
% plot(B1,B1./(exp(B2)-1)), hold on
% plot(B1,B2./(exp(B1)-1)), hold on