function []=GlucInsModel(Gop)

%   Glucose-Dependent Insulin Secretion Model Comparison
%        with a Go variable
%   Carol Lucas
%   Calls glucfun4, glucfun2, and glucfunnorm
%%  function []=GlucinsModel(Gop);


global Vm Km Vs Ks Vmax2 K2 Vmax1 K1 L ko Io Go KI Ki tau n m

Go=Gop;

[t,x]=ode45('glucfun4',[0:80],[20;10;6;.3]);
figure(1);clf;plot(t,x);xlabel('time (secs)');title('Full 4 Equation Model');
pause;
[ts,xs]=ode45('glucfun2',[0:80],[6;.3]);
figure(2);clf;plot(ts,xs);xlabel('time (secs)');title('Two Equation Model');
pause;
figure(3);clf;plot(t,x(:,3),t,x(:,4),ts,xs(:,1),'ro',ts,xs(:,2),'m+');
xlabel('time (secs)');title('Comparison');pause;
[tna,xna]=ode45('glucfunnorm',[0 4],[6/9.8;.3]);
figure(4);clf;subplot(211);plot(tna,xna);xlabel('time (secs)');
title('Normalized Simulation');
subplot(212);plot(ts,xs(:,1)/9.8,ts,xs(:,2),20*tna,xna,'o');
xlabel('time (secs)');title('Comparison');pause;

Vm=tau*Vm/Km;
Vs=tau*Vs/KI;
Ks=tau*Ks/Km;
Vmax2=tau*Vmax2/Km;
K2=K2/Km;
Vmax1=tau*Vmax1/Km;
K1=K1/Km;
L=tau*L/Km;
ko=tau*ko;
Go=Go/Km;
Ki=Ki/KI;
G=Go;

Gis=[0:.01:4];
N4=Vm*Gis;
D4=1+Gis;
Rm=N4./D4;

Rs=(Vs*(Rm.^4+L^4))./(Ks^4+Rm.^4+L^4);
I=Rs/ko;
JisJN=1./(I+1);
A=(I.^n)./(Ki^n+I.^n);
N1=K1*Vmax1*(G-Gis);
D2=(K2+G)*(K2+Gis);
D1=(K1+G)*(K1+Gis);
N2=K2*Vmax2*Gis;
N3=K2*Vmax2*G;
R1=(N1./D1).*A;
R2=(N3*JisJN.^m-N2)./D2;

figure(5);clf;subplot(231);plot(Gis,Rm);xlabel('Gi');title('Rm');
subplot(232);plot(Gis,Rs);xlabel('Gi');title('Rs');
subplot(234);plot(Gis,I);xlabel('Gi');title('I');
subplot(235);plot(Gis,A);xlabel('Gi');title('A');
subplot(233);plot(Gis,R1,Gis,N1./D1);
xlabel('Gi');title('R1 with and without I factor');axis([0 4 -20 50]);
subplot(236);plot(Gis,R2);xlabel('Gi');title('R2 at Jinf');pause;


JisGN=real(((D2.*D1.*N4-A.*D4.*D2.*N1+N2.*D4.*D1)./(N3*D4.*D1)).^(1/m));


figure(6);clf;plot(xna(:,1),xna(:,2),'o',Gis,JisJN,Gis,JisGN);
xlabel('GiHat'),ylabel('JHat');
axis([0 .8 0 1]);pause;