clear
close all

%%Load----------------------------------------------------------------------
fpath1='C:\Users\miriam\Documents\Uni\Diplomarbeit\Fortran\Final_versions\UWT\KANAL\laminar\300EW\RE1000\';
data11 = importdata([fpath1 'KL_1000_0015\Nusselt .dat'],' ',3);
%% Extract DATA
x = data11.data(1:end,1);
nu11 = data11.data(1:end,2);
Pe11=1.5;

X_all=x./(4*Pe11); % According to Shah (125)
Y_all=nu11(:,1)./nu11(end,1);
X=X_all(2:end,1);
Y=Y_all(2:end,1);
%% Asymptote
As11=0.03*X_all.^(-0.466);
gerade=0.0375*X_all.^(-0.45);
As2(1:10001,1)=1;

%%  n-calculation
nu_infty=nu11(end,1);

% Intersection point of the two asymptotes
z=((nu_infty-0.4)/1.233)^(-3); % z=x/D*1/Pe_h;

% Nu according to tabulated data of Shah HMT,Bomby 1975
nush=8.91+(8.91-8.5166)/(0.004-0.005)*(z-0.004);

% n according to equation 18/ paper of Awad
g = @(n) n-log(1+((1.233*z^(-1/3)+0.4)/nu_infty)^n)/log(nush/nu_infty);
n = fzero(g,3);
 
%% Modelling 

%Startpunkt
st2_ = [0.46 0.03];

%Definitionsintervall der Koeffizienten C und B
c=0.38:0.01:0.43;
b=0.2:0.01:0.6;

%Initialisierung
errav=3;
h=1;
i=1;
k=1;

%% Optimierungsschleife der Fit-Funktion
while (errav>0.001) && (k<size(b,2)) 

while (h<size(c,2)) 
fo_= fitoptions('method','NonlinearLeastSquares','Lower',[c(h) b(k)],'Upper',[c(h+1) b(k+1)],'Startpoint',st2_);
ft_= fittype('(1+(C+B*X.^(-0.9)).^7.0874).^(1/7.0874)','dependent',{'y'},'independent',{'X'},'coefficients',{'C','B'});
cf_= fit(X(:,1),Y(:,1),ft_,fo_);
cff2(:,i)= cf_(X);
coeffvals2(i,:) = coeffvalues(cf_);
Err2(:,i)=abs((Y(:,1)-cff2(:,i))./Y(:,1))*100;
   
%Hier Optimierung
m=1;
for j=1:1000
if ((Err2(j,i)>=1.7)) && ((Err2(j,i)< 2.3))
       xerr(m,i)=j;   
      m=m+1;
else xerr(m,i)=0; 
   m=m+1;
end
end
xmin=min(xerr(xerr~=0));
xmax=max(xerr(:,i));
errav_alle(:,i)=mean(Err2(xmin:xmax,i));
errav=mean(Err2(xmin:xmax,i));
%Ende Optimierung

i=i+1;
h=h+1;
end
h=1;
k=k+1;
end

%Ermittlung bester Lösung
avmin=min(errav_alle);
for i=1:size(errav_alle,2)
if errav_alle(1,i)==avmin
    m=i
end
end

C_best=coeffvals2(m,1);
B_best=coeffvals2(m,2);

figure(1)
plot(errav_alle(1,:))

figure(2)
hold on
plot(X,Err2(:,m))
% plot(X,Err2(:,42),'r')

test=(1+(C_best+B_best*X.^(-0.91)).^7.0874).^(1/7.0874);

figure(3)
hold on
plot(X,Y)
plot(X,test,'r')







%save
out.Pe_D = Pe11;
out.n = n;
out.B = B_best;
out.C = C_best;
out.m = -0.91;
save('RE1000_0015.mat', '-struct', 'out');