% PROGRAMA: porgrama principal de cálculo
if calor==1
  % Tasa de calor automática 
   Ang_adelanto=[-25 -15 -5];
   Qdd=[10 30 35];
   Qpp=[20 8 13]; 
   gamdd=[88 80 86];
   gampp=[8 8.5 15];
   Mpp=[2 3 4];
   Qp=INTERP1(Ang_adelanto,Qpp,IIE,'spline');
   Qd=INTERP1(Ang_adelanto,Qdd,IIE,'spline');
   gamd=(pi/180)*INTERP1(Ang_adelanto,gamdd,IIE,'spline');
   gamp=(pi/180)*INTERP1(Ang_adelanto,gampp,IIE,'spline');
   Mp=INTERP1(Ang_adelanto,Mpp,IIE,'spline');
   if IIE>-18
     Md=0.5;
   else
     Ang_ade=[-25 -18];
     Mdd=[0.9 0.5];
     Md=INTERP1(Ang_ade,Mdd,IIE,'spline');
  end
end
% Comienza el programa
C=CC*10^-6;
D=DD*10^-2;
R=RR*10^-2;
L=LL*10^-2;
TE=TEE+273;
PE=PEE*10^5;
TW=TWW+273;
F=FM+((Carga-100)/100)*(FM-dosado_min);
II=IIE*pi/180;
Cf=0.15+(0.4*(Edad/10));
Vc=(C/Z)/(RC-1);
Vs=(RC/(RC-1))*(C/Z);
n=L/R;
ang_comb=II;
Rend_vol=0.9-(0.2/4000)*(RPM-1000);
Mol=(PE*Vs/(8.314*TE))*Rend_vol;%Moles de aire admitido
N=RPM*2*pi/60;%rad/s
Cdf=0.9;%Coeficiente de descarga del inyector de combustible
dn=1;%diámetro del inyector(mm)
dve=0.38*D;%diametro de la válvula de escape(m)
dvi=0.44*D;%diametro dela válvula de admisión(m)
CDes=0.35;
CDin=0.35;
beta=pi/4;%ángulo de asiento de la válvulas
pc(1)=PE;
Tc(1)=pc(1)*Vs/(8.314*Mol);
TM=Tc(1);
PM=pc(1);
Fc(1)=pc(1)*((pi*D^2)/4);
Dc=D;
Lc=5*D;
RPM_int=linspace(0,rev_maxima,7);
p_int=[1 0.999 .99 .98 .975 .97 .95];
p_esc=linspace(1,1.05,7);
T_int=linspace(1,0,7);
m_entrada=0;
m_salida=0;
%  Coeficientes de la ecuación de combustión
%  Combustión -->  a1 CO2 + a2 H2O + a3 02 + a4 N2 + a5 CnHm
nn=10.8;
mm=18.7;
PMA=28.64*(1-(HR/100))+18*(HR/100);
Masa=Mol*PMA;%(g)Masa de aire admitido

% Humedad relativa
pvs=10^((7.5*TEE/(TEE+273.309))+2.7858);
ws=0.622*(pvs/(PE-pvs));
w=ws*(HR/100);
mas=Masa/(1+w);
mw=Masa-mas;
a1(1)=0;
a2(1)=mw/18;
a3(1)=0.21*(mas/28.64);
a4(1)=0.79*(mas/28.64);
a5(1)=0;

mf1=(mas*F/14.5);%(g)Masa de combustible a inyectar
mol1=mf1/(nn*12+mm);
mc(1)=Masa;
nm(1)=mc(1)/PMA;
mf1f=0;%(g)masa de combustible en fugas
mol1f=0;


%Coeficientes para el módulo de fugas
if FU==1
H=D/10; %Ancho de los segmentos
holgura=0.003*D;
A=(pi/4)*((D+2*holgura)^2)-pi*(D^2)/4;
Ap=0.004*A;
Vh=A*H;
mfu(1)=0;%kg/s
mffu(1)=0;%g/rad
p2(1)=PE;%presión entre los segmentos primero y segundo
p3(1)=PE;%presión entre los segmentos segundo y tercero
m2(1)=p2(1)*Vh/(287*TW);%masa retenida entre los segmentos 1º y 2º
m3(1)=p3(1)*Vh/(287*TW);%masa retenida entre los segmentos 2º y 3º
T2(1)=TW;%Temperatura entre segmentos
T3(1)=TW;%Temperatura entre segmentos
else
p2(1)=0;
p3(1)=0;
m2(1)=0;
m3(1)=0;
T2(1)=0;
T3(1)=0;
end

% BLOQUE DE CREACIÓN DE VARIABLES 

%variables volumen de la cámara y área de transferencia relativas al lazo de combustión

angc=linspace(-pi+RCA,pi-AAE,1000);
for i=1:1000
   Vcc(i)=Vc + ((Vs-Vc)/2)*(1+n-(n^2-sin(angc(i))^2)^0.5-cos(angc(i)));
   Ac(i)=pi*(D^2)/2+pi*D*(((Vcc(i)-Vc)/(Vs-Vc))*2*R+(Vc/(0.25*pi*D^2)));
end
%variables relativas al lazo de bombeo
TI=TE+(TW-TE)*INTERP1(RPM_int,T_int,RPM,'spline');
PI=PE*INTERP1(RPM_int,p_int,RPM,'spline');
PES=PE*INTERP1(RPM_int,p_esc,RPM,'spline');

angb=linspace(pi-AAE,3*pi+RCA,1000);
for i=1:1000
   Vcb(i)=Vc + ((Vs-Vc)/2)*(1+n-(n^2-sin(angb(i))^2)^0.5-cos(angb(i)));
   Ab(i)=pi*(D^2)/2+pi*D*(((Vcb(i)-Vc)/(Vs-Vc))*2*R+(Vc/(0.25*pi*D^2)));
end

% Bucle de realimentación 
for j=1:2

% BLOQUE DE CREACIÓN DE VARIABLES DEL PROBLEMA 

%cálculo del angulo de retraso   
k=1;
while angc(k)<ang_comb
   k=k+1;
end
Rc_comb=Vs/((Vcc(k)+Vcc(k+1))/2);
vp=2*(RPM/60)*(2*R);
Ea=618840/(CN+25);
TC=TM*(Rc_comb)^0.3;
PC=(PM*10^-5)*(Rc_comb)^1.3;
ang_retraso=(0.36+0.22*vp)*exp(Ea*(1/(8.314*TC)-1/17190)*(21.2/(PC-12.4))^0.63);
ang_comb=II+ang_retraso*pi/180;
tid=ang_retraso*10^3*(pi/180)/N;
%Creación de la variable calor liberado en la combustión
cont=1;  
for i=1:1000
   if angc(i)<ang_comb
      Qr(i)=0;
      cont=cont+1;
   else
      Qr(i)=(6.9*((Qp/gamp)*(Mp+1)*((angc(i)-ang_comb)/gamp)^Mp)*exp(-6.9*((angc(i)-ang_comb)/gamp)^(Mp+1))+6.9*(Qd/gamd)*(Md+1)*(((angc(i)-ang_comb)/gamd)^Md)*exp(-6.9*(((angc(i)-ang_comb)/gamd)^(Md+1))));
   end
end
% Creación de la variable (%) de masa quemada
Q(1)=0;
for i=1:999
   Q(i+1)=Q(i)+Qr(i)*(angc(i+1)-angc(i));
end
Qrr=((mf1-mf1f)*Hc*1000/Q(1000))*Qr;%(J/rad)
Mff=(mol1-mol1f)*Qr/Q(1000);%rate fuel-burning
Mf=Q/Q(1000);%tanto por ciento de masa quemada
Qid=mf1*Hc*Mf;%energía liberada ideal (KJ)
%BLOQUE DE RESOLUCIÓN DE LA COMBUSTIÓN
 g=gamf1(a1(1),a2(1),a3(1),a4(1),a5(1),TM);
 gammma_c(1)=g(1);
 A1=g(2);

%Resolución de la ecuación diferencial
    mf(1)=0;
    mff(1)=0;
  for i=1:999
     if aceptar==1
        contador=contador+1;
        set(completado,'string',round(contador/40));
     elseif aceptar1==1
        contador=contador+1;
        set(completado,'string',round(contador/360));
     elseif aceptar2==1
        contador=contador+1;
        set(completado,'string',round(contador/360));
     end
     %inyección de combustible
   if angc(i+1)<II
      mff(i+1)=0;
      mf(i+1)=0;
   else
      if mf(i) < mf1*(1+(Agua/100))   
      mff(i+1)=round(2*(CC/1000))*Cdf*(pi/4)*(dn^2)*(10^-3)*((2*0.920*(450*10^5-pc(i)))^0.5)/N;
      mf(i+1)=mf(i)+mff(i+1)*(angc(i+1)-angc(i));
      else
      mff(i+1)=0;
      mf(i+1)=mf(i);   
   end
   end

%Cálculo de gamma para el siguiente paso
   a1(i+1)=a1(i)+nn*Mff(i)*(angc(i+1)-angc(i));a11=a1(i+1);
   a2(i+1)=a2(i)+(Agua/100)*(angc(i+1)-angc(i))*mff(i)/18+(mm/2)*Mff(i)*(angc(i+1)-angc(i));a22=a2(i+1);
   a3(i+1)=a3(i)-(nn+mm/4)*Mff(i)*(angc(i+1)-angc(i));a33=a3(i+1);
   a4(i+1)=a4(i);a44=a4(i+1);
   a5(i+1)=a5(i)+(1-(Agua/100))*mff(i)*(angc(i+1)-angc(i))/(nn*12+mm)-Mff(i)*(angc(i+1)-angc(i)); a55=a5(i+1);
   TT=Tc(i);
   g=gamf1(a11,a22,a33,a44,a55,TT);
   gamma=g(1);
   A1=g(2);
   gamma_c(i+1)=gamma;
   Mc=(a11*28+a22*18+a33*16+a44*32+a55*(nn*12+mm))/(a11+a22+a33+a44+a55);
   hf=(4180/(nn*12+mm))*(-9.1063*(Tc(i)/1000)+246.97*((Tc(i)/1000)^2)/2-143.74*((Tc(i)/1000)^3)/3+32.329*((Tc(i)/1000)^4)/4-0.0518/(Tc(i)/1000)-50.128+23.514);%entalpía del combustible
 %Cálculo de las variables del problema
   if FU==0
      mc(i+1)=mc(i)+mff(i+1)*(angc(i+1)-angc(i)); 
      mfu(i+1)=0;
      mffu(i+1)=0;
      p2(i+1)=0;
      p3(i+1)=0;
      m2(i+1)=0;
      m3(i+1)=0;
      T2(i+1)=0; 
      T3(i+1)=0;
   else
      res=fuga1(pc(i),Tc(i),Vcc(i),p2(i),m2(i),p3(i),m3(i),PE,TW,(angc(2)-angc(1)),N,gamma,Mc,A,Ap,Vh,D,vp,Za,Cf);
      mfu(i+1)=res(1);
      p2(i+1)=res(2);
      p3(i+1)=res(3);
      m2(i+1)=res(4);
      m3(i+1)=res(5);
      T2(i+1)=res(6);
      T3(i+1)=res(7);
      mffu(i+1)=1000*mfu(i+1)/(RPM/60);
      mc(i+1)=mc(i)+mff(i+1)*(angc(i+1)-angc(i))-mffu(i+1)*(angc(i+1)-angc(i));
   end
   if angc(i+1)<ang_comb;
      qc(i+1)=heat_loses(Tc(i),pc(i),TW,Vcc(i),Ac(i),gamma,vp,D,0);
   else
      qc(i+1)=heat_loses1(Tc(i),pc(i),Tc(cont-1),pc(cont-1),Vcc(cont-1),pc(1),Vcc(1),TW,Vcc(i),Ac(i),Vs,Vc,gamma,vp,D);
   end
   nm(i+1)=mc(i+1)/Mc;
   C1=1-(A1*Vcc(i)*pc(i))/(nm(i+1)*8.314*(gamma-1));
   C2=A1*pc(i)^2/(nm(i+1)*8.314*(gamma-1))-(gamma*pc(i)/Vcc(i));
   C3=(gamma-1)/Vcc(i);
   C4=A1*(pc(i)^2)*Vcc(i)/((nm(i+1)^2)*8.314*(gamma-1));
   pc(i+1)=pc(i)+(C2/C1)*(Vcc(i+1)-Vcc(i))+(C3/C1)*(angc(i+1)-angc(i))*Qrr(i)-(C3/C1)*(angc(i+1)-angc(i))*qc(i+1)/N+(C3/C1)*mff(i+1)*(angc(i+1)-angc(i))*hf-(C4/C1)*(nm(i+1)-nm(i));
   Tc(i+1)=pc(i+1)*Vcc(i+1)/(mc(i+1)*(8.314/Mc));
end

%BLOQUE DE CÁLCULO DEL LAZO DE BOMBEO

Tb(1)=Tc(1000);
pb(1)=pc(1000);
mb(1)=mc(1000);
qb(1)=qc(1000);
gamma_b(1)=gamma;
pa1=a11/(a11+a22+a33+a44+a55);a1(1001)=a11;
pa2=a22/(a11+a22+a33+a44+a55);a2(1001)=a22;
pa3=a33/(a11+a22+a33+a44+a55);a3(1001)=a33;
pa4=a44/(a11+a22+a33+a44+a55);a4(1001)=a44;
pa5=a55/(a11+a22+a33+a44+a55);a5(1001)=a55;
for i=1:999
       if aceptar==1
       contador=contador+1;
       set(completado,'string',round(contador/40));
       elseif aceptar1==1
       contador=contador+1;
       set(completado,'string',round(contador/360));
       elseif aceptar2==1
       contador=contador+1;
       set(completado,'string',round(contador/360));
       end
       % periodo de escape  
  if angb(i)<(2*pi-AAA)
     cont=cont+1;
   g=gamf1(a11,a22,a33,a44,a55,TT);
   gamma=g(1);
   A1=g(2);
   gamma_b(i+1)=gamma;
     lii=0.3*sin(pi*(angb(i)-(pi-AAE))/(AAE+pi+RCE));
     Ave=area(lii,beta,dve);
     % flujo bloqueado   
   if (PES/pb(i))<(2/(gamma+1))^(gamma/(gamma-1));   
   mb(i+1)=mb(i)-CDes*Ave*(1000/N)*(angb(i+1)-angb(i))*pb(i)*((287*Tb(i))^(-0.5))*(gamma^0.5)*((2/(gamma+1))^((gamma+1)/(2*(gamma-1))));
   else
     % flujo sin bloquear 
   mb(i+1)=mb(i)-CDes*Ave*(1000/N)*(angb(i+1)-angb(i))*pb(i)*((287*Tb(i))^(-0.5))*((PES/pb(i))^(1/gamma))*((2*gamma/(gamma-1))*(1-((PES/pb(i))^((gamma-1)/gamma))))^0.5;
   end 
   qb(i+1)=heat_loses(Tb(i),pb(i),TW,Vcb(i),Ab(i),gamma,vp,D,1);
   pb(i+1)=pb(i)-gamma*pb(i)*(Vcb(i+1)-Vcb(i))/Vcb(i)+pb(i)*(gamma_b(i+1)-gamma_b(i))/(gamma-1)+gamma*pb(i)*(mb(i+1)-mb(i))/mb(i)-((gamma-1)/Vcb(i))*(angb(i+1)-angb(i))*qb(i+1)/N;   
   Tb(i+1)=pb(i+1)*Vcb(i+1)/(mb(i+1)*(8.314/Mc));
   a1(1000+i+1)=a1(1000+i)+pa1*(mb(i+1)-mb(i))/Mc;
   a2(1000+i+1)=a2(1000+i)+pa2*(mb(i+1)-mb(i))/Mc;
   a3(1000+i+1)=a3(1000+i)+pa3*(mb(i+1)-mb(i))/Mc;
   a4(1000+i+1)=a4(1000+i)+pa4*(mb(i+1)-mb(i))/Mc;
   a5(1000+i+1)=a5(1000+i)+pa5*(mb(i+1)-mb(i))/Mc;
%periodo de solape   
elseif angb(i)<(2*pi+RCE) 
   qb(i+1)=heat_loses(Tb(i),pb(i),TW,Vcb(i),Ab(i),gamma,vp,D,1);
   if pb(i)>PES %vaciado & puente  
   lie=0.3*sin(pi*(angb(i)-(pi-AAE))/(AAE+pi+RCE));
   Ave=area(lie,beta,dve);
   liin=0.3*sin(pi*(angb(i)-(2*pi-AAA))/(pi+RCA+AAA));  
   Avi=area(liin,beta,dvi);
   mbe=CDes*Ave*(1000/N)*(angb(i+1)-angb(i))*pb(i)*((287*Tb(i))^(-0.5))*((PES/pb(i))^(1/gamma))*((2*gamma/(gamma-1))*(1-((PES/pb(i))^((gamma-1)/gamma))))^0.5;
   mbin=CDin*Avi*(1000/N)*(angb(i+1)-angb(i))*pb(i)*((287*Tb(i))^(-0.5))*((PI/pb(i))^(1/gamma))*((2*gamma/(gamma-1))*(1-((PI/pb(i))^((gamma-1)/gamma))))^0.5;
   mb(i+1)=mb(i)-mbe-mbin;
   a1(1000+i+1)=a1(1000+i)+pa1*(mb(i+1)-mb(i))/Mc;a11=a1(1000+i+1);
   a2(1000+i+1)=a2(1000+i)+pa2*(mb(i+1)-mb(i))/Mc;a22=a2(1000+i+1);
   a3(1000+i+1)=a3(1000+i)+pa3*(mb(i+1)-mb(i))/Mc;a33=a3(1000+i+1);
   a4(1000+i+1)=a4(1000+i)+pa4*(mb(i+1)-mb(i))/Mc;a44=a4(1000+i+1);
   a5(1000+i+1)=a5(1000+i)+pa5*(mb(i+1)-mb(i))/Mc;a55=a5(1000+i+1);
   g=gamf1(a11,a22,a33,a44,a55,TT);
   gamma=g(1);
   A1=g(2);
   gamma_b(i+1)=gamma;
   Mb=Mc;
   pb(i+1)=pb(i)-gamma*pb(i)*(Vcb(i+1)-Vcb(i))/Vcb(i)+pb(i)*(gamma_b(i+1)-gamma_b(i))/(gamma-1)+gamma*pb(i)*(mb(i+1)-mb(i))/mb(i)-((gamma-1)/Vcb(i))*(angb(i+1)-angb(i))*qb(i+1)/N;   
   %backflow
   elseif pb(i)>PI 
   lie=0.3*sin(pi*(angb(i)-(pi-AAE))/(AAE+pi+RCE));
   Ave=area(lie,beta,dve);
   liin=0.3*sin(pi*(angb(i)-(2*pi-AAA))/(pi+RCA+AAA));  
   Avi=area(liin,beta,dvi);
   mbin=CDin*Avi*(1000/N)*(angb(i+1)-angb(i))*pb(i)*((287*Tb(i))^(-0.5))*((PI/pb(i))^(1/gamma))*((2*gamma/(gamma-1))*(1-((PI/pb(i))^((gamma-1)/gamma))))^0.5;
   mbe=CDes*Ave*(1000/N)*(angb(i+1)-angb(i))*PES*((287*TE)^(-0.5))*((pb(i)/PES)^(1/gamma))*((2*gamma/(gamma-1))*(1-((pb(i)/PES)^((gamma-1)/gamma))))^0.5;
   mb(i+1)=mb(i)-mbin+mbe;
   M=(1-(EGR/100))*PMA+(EGR/100)*Mc;
   a1(1000+i+1)=a1(1000+i)+pa1*(mb(i+1)-mb(i))/Mc;a11=a1(1000+i+1);
   a2(1000+i+1)=a2(1000+i)+pa2*(mb(i+1)-mb(i))/Mc;a22=a2(1000+i+1);
   a3(1000+i+1)=a3(1000+i)+pa3*(mb(i+1)-mb(i))/Mc;a33=a3(1000+i+1);
   a4(1000+i+1)=a4(1000+i)+pa4*(mb(i+1)-mb(i))/Mc;a44=a4(1000+i+1);
   a5(1000+i+1)=a5(1000+i)+pa5*(mb(i+1)-mb(i))/Mc;a55=a5(1000+i+1);
   g=gamf1(a11,a22,a33,a44,a55,TT);
   gamma=g(1);
   A1=g(2);
   gamma_b(i+1)=gamma;
   Mb=(a11*28+a22*18+a33*16+a44*32+a55*(nn*12+mm))/(a11+a22+a33+a44+a55);
   pb(i+1)=pb(i)-gamma*pb(i)*(Vcb(i+1)-Vcb(i))/Vcb(i)+pb(i)*(gamma_b(i+1)-gamma_b(i))/(gamma-1)+gamma*PE*(mb(i+1)-mb(i))/(1130*Vcb(i))-((gamma-1)/Vcb(i))*(angb(i+1)-angb(i))*qb(i+1)/N;   
   else%llenado
   lie=0.3*sin(pi*(angb(i)-(pi-AAE))/(AAE+pi+RCE));
   Ave=area(lie,beta,dve);
   liin=0.3*sin(pi*(angb(i)-(2*pi-AAA))/(pi+RCA+AAA));  
   Avi=area(liin,beta,dvi);
   mbe=CDin*Ave*(1000/N)*(angb(i+1)-angb(i))*PES*((287*TE)^(-0.5))*((pb(i)/PES)^(1/gamma))*((2*gamma/(gamma-1))*(1-((pb(i)/PES)^((gamma-1)/gamma))))^0.5;
   mb(i+1)=mb(i)+mbe+CDes*Avi*(1000/N)*(angb(i+1)-angb(i))*PI*((287*TI)^(-0.5))*((pb(i)/PI)^(1/gamma))*((2*gamma/(gamma-1))*(1-((pb(i)/PI)^((gamma-1)/gamma))))^0.5;
   M=(1-(EGR/100))*PMA+(EGR/100)*Mc;
   a1(1000+i+1)=a1(1000+i)+pa1*(EGR/100)*(mb(i+1)-mb(i))/M;a11=a1(1000+i+1);
   a2(1000+i+1)=a2(1000+i)+pa2*(EGR/100)*(mb(i+1)-mb(i))/M+(1/(1+(1/w)))*(1-(EGR/100))*(mb(i+1)-mb(i))*(PMA/(18*M));a22=a2(1000+i+1);
   a3(1000+i+1)=a3(1000+i)+pa3*(EGR/100)*(mb(i+1)-mb(i))/M+(0.21/(1+w))*(1-(EGR/100))*(mb(i+1)-mb(i))*(PMA/(28.64*M));a33=a3(1000+i+1);
   a4(1000+i+1)=a4(1000+i)+pa4*(EGR/100)*(mb(i+1)-mb(i))/M+(0.79/(1+w))*(1-(EGR/100))*(mb(i+1)-mb(i))*(PMA/(28.64*M));a44=a4(1000+i+1);
   a5(1000+i+1)=a5(1000+i)+pa5*(EGR/100)*(mb(i+1)-mb(i))/M;a55=a5(1000+i+1);
   g=gamf1(a11,a22,a33,a44,a55,TT);
   gamma=g(1);
   A1=g(2);
   gamma_b(i+1)=gamma;
   Mb=(a11*28+a22*18+a33*16+a44*32+a55*(nn*12+mm))/(a11+a22+a33+a44+a55);
   pb(i+1)=pb(i)-gamma*pb(i)*(Vcb(i+1)-Vcb(i))/Vcb(i)+pb(i)*(gamma_b(i+1)-gamma_b(i))/(gamma-1)+gamma*PI*(mb(i+1)-mb(i))/(1130*Vcb(i))-((gamma-1)/Vcb(i))*(angb(i+1)-angb(i))*qb(i+1)/N;   
   end
   Tb(i+1)=pb(i+1)*Vcb(i+1)/(mb(i+1)*(8.314/Mb));
 elseif angb(i)<(3*pi+RCA)
% periodo de llenado  
   if pb(i)>=PI   
   liin=0.3*sin(pi*(angb(i)-(2*pi-AAA))/(pi+RCA+AAA));  
   Avi=area(liin,beta,dvi);
   mbin=CDin*Avi*(1000/N)*(angb(i+1)-angb(i))*pb(i)*((287*Tb(i))^(-0.5))*((PI/pb(i))^(1/gamma))*((2*gamma/(gamma-1))*(1-((PI/pb(i))^((gamma-1)/gamma))))^0.5;
   mb(i+1)=mb(i)-mbin;
   m_salida=m_salida+mbin;
   else
   liin=0.3*sin(pi*(angb(i)-(2*pi-AAA))/(pi+RCA+AAA));  
   Avi=area(liin,beta,dvi);
   if (pb(i)/PI)<(2/(gamma+1))^(gamma/(gamma-1));   
   mb(i+1)=mb(i)+CDin*Avi*(1000/N)*(angb(i+1)-angb(i))*PI*((287*TI)^(-0.5))*(gamma^0.5)*((2/(gamma+1))^((gamma+1)/(2*(gamma-1))));
   else
   % flujo sin bloquear 
   mb(i+1)=mb(i)+CDin*Avi*(1000/N)*(angb(i+1)-angb(i))*PI*((287*TI)^(-0.5))*((pb(i)/PI)^(1/gamma))*((2*gamma/(gamma-1))*(1-((pb(i)/PI)^((gamma-1)/gamma))))^0.5;
   end 
   end
   M=(1-(EGR/100))*PMA+(EGR/100)*Mc;
   a1(1000+i+1)=a1(1000+i)+pa1*(EGR/100)*(mb(i+1)-mb(i))/M;a11=a1(1000+i+1);
   a2(1000+i+1)=a2(1000+i)+pa2*(EGR/100)*(mb(i+1)-mb(i))/M+(1/(1+(1/w)))*(1-(EGR/100))*(mb(i+1)-mb(i))*(PMA/(18*M));a22=a2(1000+i+1);
   a3(1000+i+1)=a3(1000+i)+pa3*(EGR/100)*(mb(i+1)-mb(i))/M+(0.21/(1+w))*(1-(EGR/100))*(mb(i+1)-mb(i))*(PMA/(28.64*M));a33=a3(1000+i+1);
   a4(1000+i+1)=a4(1000+i)+pa4*(EGR/100)*(mb(i+1)-mb(i))/M+(0.79/(1+w))*(1-(EGR/100))*(mb(i+1)-mb(i))*(PMA/(28.64*M));a44=a4(1000+i+1);
   a5(1000+i+1)=a5(1000+i)+pa5*(EGR/100)*(mb(i+1)-mb(i))/M;a55=a5(1000+i+1);
   g=gamf1(a11,a22,a33,a44,a55,TT);
   gamma=g(1);
   A1=g(2);
   gamma_b(i+1)=gamma;
   Mb=(a11*28+a22*18+a33*16+a44*32+a55*(nn*12+mm))/(a11+a22+a33+a44+a55);
   qb(i+1)=heat_loses(Tb(i),pb(i),TW,Vcb(i),Ab(i),gamma,vp,D,1);
   pb(i+1)=pb(i)-gamma*pb(i)*(Vcb(i+1)-Vcb(i))/Vcb(i)+pb(i)*(gamma_b(i+1)-gamma_b(i))/(gamma-1)+gamma*PI*(mb(i+1)-mb(i))/(1130*Vcb(i))-((gamma-1)/Vcb(i))*(angb(i+1)-angb(i))*qb(i+1)/N;   
   Tb(i+1)=pb(i+1)*Vcb(i+1)/(mb(i+1)*(8.314/Mb));
 end 
end
% Valores de variables de realimentación del programa completo
if j==1   
M=(a11*28+a22*18+a33*16+a44*32+a55*(nn*12+mm))/(a11+a22+a33+a44+a55);
Tc(1)=Tb(1000);
pc(1)=pb(1000);
mc(1)=mb(1000);
TM=Tc(1);
PM=pc(1);
a1(1)=a11;
a2(1)=a22;
a3(1)=a33;
a4(1)=a44;
a5(1)=a55;
Mol=((mb(1000)-min(mb))/M)*(1-EGR/100);
Masa=Mol*PMA;%(g)
mas=Masa/(1+w);
mf1=(mas*F/14.5);%(g)
mol1=mf1/(12*nn+mm);
mf1f=(sum(mffu)*(angc(1000)-angc(1))/1000)*(F/14.5);
mol1f=mf1f/(12*nn+mm);
end
end
%Resultados obtenidos
ANG=[angc angb];
V=[Vcc Vcb];
P=[pc pb];
T=[Tc Tb];
m=[mc mb];
q=[qc qb];
gam=[gamma_c gamma_b];
Masa_ideal=(PE*Vs/(8.314*TE))*28.64;
Mol=((mb(1000)-min(mb))/M)*(1-EGR/100);
Masa=Mol*PMA;%(g)
mas=Masa/(1+w);
Rend_vol=mas/Masa_ideal;
pres_max_ciclo=max(P)/10^5;
 %potencia indicada
 trab=0;
 trab_bomb=0;
 cont1=1;
 cont2=1;
 cont3=1;
 cont4=1;
 while ANG(cont1)<=3*pi
    cont1=cont1+1;
 end
 while ANG(cont3)<=pi
    cont3=cont3+1;
 end
 for j=cont1:2000
    VV(cont2)=V(j);
    pp(cont2)=P(j);
    cont2=cont2+1;
 end
 for j=1:cont1-1
    VV(cont2)=V(j);
    pp(cont2)=P(j);
    cont2=cont2+1;
 end
 for i=1:(1000+(2000-cont1)+(cont3-1000))
    trab=trab+(VV(i+1)-VV(i))*(pp(i+1)+pp(i))/2;
 end
 for j=(1000+(2000-cont1)+(cont3-1000)):1999
    trab_bomb=trab_bomb-(VV(j+1)-VV(j))*(pp(j+1)+pp(j))/2;
 end
 
 trabajo_indicado=trab;
 trabajo_bombeo=trab_bomb;
 presion_indicada=(trab*10^-5)/(Vs-Vc);    
 pm1=(0.701*10^-3)*(30-4*(RPM/60)/1000)*(dvi^1.75)/((D^2)*2*R);
 pm2=(5.5*10^-2)*(D/(2*R))*(RPM/60000);
 pm3=(3.69*10^-3)*2*2*R*3/D^2+(7.56*10^-3)*(0.5*D)*vp/(2*R*D);
 presion_perdidas_mecanicas=pm1+pm2+pm3;
 presion_perdidas_bombeo=0.35*(trab_bomb*10^-5)/(Vs-Vc);
 presion_efectiva=presion_indicada-presion_perdidas_bombeo-presion_perdidas_mecanicas;
 potencia_indicada=presion_indicada*10^2*(Vs-Vc)*Z*(RPM/60)*0.5;
 potencia_efectiva=presion_efectiva*10^2*(Vs-Vc)*Z*(RPM/60)*0.5;
 par_indicado=potencia_indicada*1000/N;
 par_efectivo=potencia_efectiva*1000/N;
 cons_especifico=(Z/2)*(0.85*Masa*F/14.5)*RPM*60/(potencia_efectiva);
 rend_efectivo=(1/(cons_especifico*Hc))*3600;
 residuales_ciclo=100*(min(mb)+m_entrada+m_salida)/Masa;
 fugas_porcien=100*(sum(mffu)*(angc(1000)-angc(1))/1000)/Masa;
 fr(1)=0;
 parins=Fc*R;
 for i=1:999
    fr(i+1)=fr(i)+mffu(i+1)*(angc(i+1)-angc(i));
 end
 Qfuga=fr*(Hc*F/14.5);
 Qheat(1)=0;
 for i=1:999
    Qheat(i+1)=Qheat(i)+(angc(i+1)-angc(i))*qc(i+1)/N;
 end
 Qreal1=Qid-Qheat/1000;
 Qreal2=Qid-Qfuga-Qheat/1000;
if aceptar==1 
%ciclo ideal 
gamid=median(gamma_c);
angid=linspace(-pi,pi,1000);
for i=1:1000
   Videal(i)=Vc + ((Vs-Vc)/2)*(1+n-(n^2-sin(angid(i))^2)^0.5-cos(angid(i)));
end

for i=1:500
   Vid(i)=Vc + ((Vs-Vc)/2)*(1+n-(n^2-sin(angid(i))^2)^0.5-cos(angid(i)));
end
pid(1)=pc(1);
for i=1:499
   pid(i+1)=pid(i)*(Vid(i)/Vid(i+1))^gamid;
end
pmax=max(pc);
T22=pid(500)*Vc/(8.314*mb(1000)/M);
T33=T22*(pmax/pid(500));
pid(501)=pmax;
Vid(501)=Vc;
T44=(1/1.005)*(mf1*Hc*1000-0.718*(T33-T22)+1.005*T33);
Vid(502)=(mb(1000)/M)*8.314*T44/pmax;
pid(502)=pmax;
k=1;
while Vid(502)>Videal(500+k)
   k=k+1;
end
for i=1:500-k
   Vid(502+i)=Videal(500+k+i);
end
for i=502:1002-k-1
   pid(i+1)=pid(i)*(Vid(i)/Vid(i+1))^gamid;
end
end   
if aceptar==1
 if presion_efectiva<0
 Hd_error=errordlg('El motor no puede vencer las perdidas mecánicas','Error');
 waitfor(Hd_error);
 aceptar=0;
 end
end

