#include "math.h"
#include "stdio.h"
#define PI 3.1415926

double D,S,L,ge,e,n,Ne,gnum,ysb,Hu,Lo,M,R,P_SC,T_SC;
/*缸径、行程、连杆长、比油耗、转速、功率、气缸数、压缩比、低热值、理论空燃比、Weibe函数常数m、气体常数、进气门关闭时压力和温度 */
    double m_T[720],m_P[720],m_a[720];
/*气缸温度、压力、瞬时过量空气系数 */
    int m_CA[720],A_SC,A_CB,A_CE,A_EO,A_SO,A_EC;
    double Pa,Ta,hs,P,V,T,CA,m,mf,nmd,Pe,mBo,mL,dXY;
    double Cv();
    double dV_dCA();
    double dQw_dCA();
    double dQB_dCA();
    double U();
    double V_CA();
    double dU_dnmd();
    double dT_Ys_dCA();
    double dT_Pz_dCA();
    double dT_Rs_dCA();
    double dT_Hq_dCA();
    double dT_dCA(double an);
    void Calculate();

/* A编写各子函数程序*/
    double Cv() /*求定容比热*/
{double Cv,x,y;
    x=-3*(0.0975+0.0485/pow(a,0.75))*(T-273)*(T-273)*0.000001;
    y=2*(7.768+3.36/pow(a,0.8))*(T-273)*0.0001;
    Cv=0.14455*(x+y+(489.6+46.4/pow(a,0.93)*0.01));
    return (Cv);
    }
double  dV_dCA() /*求dv/dφ*/
{double V_CA,x,y;
       x=sin(CA);
       y=S*sin(CA)*cos(CA)/2*sqrt(1-S*S/(4*L*L)*sin(CA)*sin(CA));
       V_CA=(D*D*S/8)*PI*(x+y)*PI/180;
       return(V_CA);
}
double V_CA() /*求气缸容积*/
{double V,x,y;
x=2/(e-1)+1-cos(CA);
y=(2*L/S)*(1-sqrt(1-S*S/(4*L*L)*sin(CA)*sin(CA)));
V=(D*D*S/8)*PI*(x+y);
return(V);
}
double dQw_dCA() /*求dQw/dφ*/
{double dQw_dCA,K,Tw1,Tw2,Tw3,Pe;
Tw1=273+120+30*Pe;
Tw2=273+100+7*Pe;
Tw3=273+100+4*Pe;
K=130*pow(D,-2)*pow(P,0.8)*pow(T,-0.53)*pow(2.28*S*n/30,0.8);
dQw_dCA=1/(60*n)*K*((PI*D*D/4)*(Tw1+Tw2-2*T)+(4/D)*V*(Tw3-T));
return(dQw_dCA);
}
double dQB_dCA() /*求dQB/dφ*/
{double dQB_dCA,y;
y=(CA-A_CB)/60;
dQB_dCA=Hu*ge*6.908*pow(y,M)*exp(-6.908*pow(y,(1.0+M)))*((1.0+M)/60.0);
return (dQB_dCA);
}
double U() /*求u*/
{double x,y,z,U;
x=-(0.0975+0.0485/pow(a,0.75))*(T-273)*(T-273)*(T-273)*0.000001;
y=(7.768+3.36/pow(a,0.8))*(T-273)*(T-273)*0.0001;
z=(489.6+46.4/pow(a,0.93))*(T-273)*0.01;
U=0.14455*(x+y+z+1356.8);
return (U);
}
double dU_dnmd()  /*求du/dλ*/
{double x,y,z,dU_dnmd;
x=0.0485*0.75/pow(a,0.75)*(T-273)*(T-273)*(T-273)*0.000001;
y=(3.36*0.8/pow(a,0.8))*(T-273)*(T-273)*0.0001;
z=(46.4*0.93/pow(a,1.93))*(T-273)*0.01;
dU_dnmd=0.14455*(x-y-z);
return (dU_dnmd);
}
double dT_Ys_dCA() /*压缩期*/
{double dT_Ys_dCA;
dT_Ys_dCA=(dQw_dCA-(mL*R*T/V)*V_CA)/(mL*Cv);
return (dT_Ys_dCA);
}
double dT_Pz_dCA() /*膨胀期*/
{doubledT_Pz_dCA();
dT_Pz_dCA()=(dQw_dCA-mL*R*T/V*V_CA)/(mL*Cv);
return (dT_Pz_dCA());
}
double dT_Rs_dCA() /*燃烧期*/
{double dT_Rs_dCA();
  dT_Rs_dCA()=(dQB_dCA+dQw_dCA-mL*R*T/V*V_CA-U/Hu*dQB_dCA-m*dU_dnmd*(-mL*dQB_dCA/(Lo*mBo*mBo*Hu)))/(mL*Cv);
  return ( dT_Rs_dCA());
}
double dT_dCA(double an) /*求dT/dφ*/
{
    CA=an*PI/180;
    V=V_CA();
    if(A_SC<an && an<A_CB)   return dT_Ys_dCA();
    if(A_CB<=an && an<=A_CE) return dT_Rs_dCA();
    if(A_CE<an && an<=A_EO)  return dT_Pz_dCA();
    if(A_EO<an && an<=A_SO)    /*换气过程不计算，线性赋值，避免零值 */
        P=110000+(m_P[(int)A_EO]-110000)*(A_SO-an)/(A_SO-A_EO);
    if(A_SO<an || an<A_EC)
        P=110000-20000*(an>A_SO?an-A_SO:an-A_SO+720)/(A_EC-A_SO+720);
    if(A_EC<=an && an<=A_SC) P=90000;
    return 0;
}
void Calculate(void)
{
    double Tpre,K1,K2,an;
    int index;
    Pe=30*4.0*Ne/(gnum*100000*n*S*PI*D*D/4);  /*平均有效压力bar*/
    mBo=Ne*ge/(30000000*n*gnum);   /*单缸循环供油量kg */
    CA=A_SC*PI/180; /*进气关闭时的角度转弧度*/
    V=V_CA();/*进气关闭时角度下的气缸容积*/
    P=P_SC;
    T=T_SC;
    index=(int)A_SC;
    m_T[index]=T;
    m_P[index]=P;
    m=P*V/(R*T);/*进气关闭时角度下的气缸内工质质量*/
    mL=m-0.01*mBo;/*进气关闭时角度下的气缸内空气质量*/
    mf=m-mL;
    Tpre=T;
    an=A_SC+1;
    do{if(k%2==0)
    {
        T[i+1]=T[i]+Tca[i];
        i++;
        j++;
        k++;
    }
    else
    {
        T[i]=T[i-1]+(Tca[i])/2;
        P=m*R*T[i]/V;
        printf("ca=%3lf P=%10.1lf T[%d]=%5.1lf QB=%4.1lf\n",ca,P,j,T[i],QB_CA);//
        fprintf(fp,"ca=%3.0lf, P=%10.1lf T[%d]=%5.1lf QB=%lf\n",ca,P,j,T[i],QB_CA);//
        ca--;
        k++;
    }
    }
    if(o==1)
    {j--;
    i--;
    k=0;
    }
      }while(an!=A_SC+1);
}
void main(void)
{
    int i;
    FILE *fp;
    D=0.11;   /*发动机数据按各自的任务书输入 */
    S=0.13;
    L=0.21;
    e=16;
    ge=236.0;
    n=2400;
    Ne=101568;
    gnum=6;  ysb=17;   M=0.7;
    Hu=41868000;
    R=287.08;
    Pa=101000; /*大气压力*/
    Ta=300;    /*大气温度*/
    hs=(735.0+R)*Ta;  /*大气焓值*/
    P_SC=90000;
    T_SC=313;
    A_SC=216; /*进气门关闭*/
    A_CB=350; /*燃烧开始*/
    A_CE=410; /*燃烧结束*/
    A_EO=484; /*排气门开启*/
    A_SO=705; /*进气门开启*/
    A_EC=15;  /*排气门关闭*/
    Lo=14.4;
  for(i=0;i<720;i++)
    {
        m_P[i]=0;
        m_CA[i]=i;
        m_T[i]=0;
    }
    Calculate();

    /*C编写数据输出程序，计算结束后按顺序输出0-720转角的“缸压，温度”数据到outdata.txt文件 */
   
}