function [FQS]=AnpStart
% Das ist die Rahmenfkt. fuer die ganze Parameteroptimierung

%Einlesen der Daten fuer die Anpassung
num=xlsread('E:\BPT zeuchs\MDCK_MC_data.xls');

% Mein Zeitvektor
tspan=num(:,1);

% Die Startwerte fuer die anzupassenden Parameter Y_... und mu_max
Y_Xglc=0.0369*10^3;
Y_Xgln=0.4159*10^3;
Y_Xlac=10^3/0.0141;
Y_Xamn=10^3/0.2847;
mu_max=0.0078;
para=[Y_Xglc Y_Xgln Y_Xlac Y_Xamn mu_max];

% Startwerte fuer das AWP des DGLS
X_sus=num(:,2);
X_at=num(:,3);
Glc=num(:,4);
Gln=num(:,6);
Lac=num(:,5);
Amn=num(:,8);
y=[X_sus(1);X_at(1);Glc(1);Gln(1);Lac(1);Amn(1)];

[t,y]=ode45(@(tspan,y) dffgl(tspan,y,para), tspan, y);
para=fminsearch(@(tspan,para) myode(tspan,y),tspan,para);


%x=fminsearch(@myode, (tspan,y,Y_Xglc,Y_Xgln,Y_Xlac,Y_Xamn,mu_max),[Y_Xglc,Y_Xgln,Y_Xlac,Y_Xamn,mu_max],[],t,X_sus,X_at,Glc,Gln,Lac,Amn);

myNorm=zeros(size(y));
xlsNorm=zeros(size(y));
FQS=zeros(size(y));
for i=1:6
    for j=1:length(t)
        myNorm(j,i)=y(j,i)/max(abs(y(:,i)));
        xlsNorm(j,i)=num(j,i)/max(abs(num(:,i)));
        FQS(j,i)=abs(myNorm(j,i)-xlsNorm(j,i));
    end
end

FQS=sum(sum(FQS));

function [y] = myode (tspan, y)
% Das ist die Fkt. in der 'ode' verwendet und das
% Differentialgleichungssystem aufgerufen wird.

% Schliessen evtl. offener Grafiken
for i=1:2
    close(gcf)
end

[t,y]=ode45(@(tspan,y) dffgl(tspan,y,para),tspan,y);

%{
Jetzt wird geplottet.
Dabei kommen die beiden Zellzahlen und die restlichen Groessen jeweils in 2
verschiedene Diagramme, weil sie unterschiedliche Groessenordnungen
besitzen.
%}

figure(1)
plot(t,y(:,1),'b','LineWidth',2)
hold on;
plot(t,y(:,2),'r')
legend('Zellen in Suspension','adhärente Zellen')
title('Zeitlicher Verlauf der Zellzahlen');

figure(2)
hold on
plot(t,y(:,3),'g','LineWidth',2)
plot(t,y(:,4),'m','LineWidth',2)
plot(t,y(:,5),'k','LineWidth',2)
plot(t,y(:,6),'b','LineWidth',2)
legend('Glucose','Glutamine','Lactate','Ammonia')
title('Zeitliche Aenderung der Substrate');


function [dy]=dffgl(~,y,para)
% Diese Fkt. enthaelt die DGLn

%function [dy]=myfunsys(t,y)
%{
Y_Xglc=0.3*10^3;
Y_Xgln=1.4*10^3;
Y_Xlac=10^3/2.1;
Y_Xamn=10^3/0.51;
mu_max=0.03;
%}

% Setzen der Parameter
K_glc=0.4;
K_gln=0.14;
K_lac=50;
K_amn=10;
k_dQ=0.0036;
k_d=0.007;
k_at=0.1;
k_det=0.002;
X_max=2.3;
%t_lag=24;

blubb(1)=y(3,1)/(K_glc+y(3,1));
blubb(2)=y(4,1)/(K_gln+y(4,1));
blubb(3)=K_lac/(K_lac+y(5,1));
blubb(4)=K_amn/(K_amn+y(6,1));
mu=para(1,5)*(X_max-y(2,1))/X_max*min(blubb);
%y(2,1)=0;


dy(1,1)=-k_at*(X_max-y(2,1))/X_max*y(1,1)+k_det*y(2,1)-k_d*y(1,1);

dy(2,1)=mu*y(2,1)+k_at*(X_max-y(2,1))/X_max*y(1,1)-k_det*y(2,1);

dy(3,1)=-1/para(1,1)*mu*y(2,1);

dy(4,1)=-1/para(1,2)*mu*y(2,1)-k_dQ*y(4,1);

dy(5,1)=mu*y(2,1)/para(1,3);

dy(6,1)=1/para(1,4)*mu*y(2,1)+k_dQ*y(4,1);