function gas( t, p, w)
tin=t;pin=p;win=w;dis=0;%s1=0;s2=0;s3=0;
k=2;t0=273;nd=0.00955;dd=1000;g=9.8;theta=0;dx=100;%nd=labuda  dd=D dx=distance
s=[0 0 0 0;0 0 0 0;0 0 0 0];
q=0.7;i=1;
%input constant
a0 = 4.064815074200973e+000;b0 =  1.290337412215489e-002;c0 = 4.478469090956744e+005;d0 = 9.602660731459433e+006;
e0 = 5.690616464391430e+006;a =   1.128423965219249e+000;b =  6.437697890312933e-004;c =  6.426612444129175e+004;
d =  4.589139032818672e+001;ar =  1.644814633272486e-006;gm = 4.436194791067355e-004;
y1=0.9515  ; y2=0.485 ;
w1=0.0126 ;w2=0.0978;
Tc1=190.55; Tc2=305.43;
pc1=46.6232; pc2=57.7851;
r=8.3143;
%loop
while(dis<800)
for(cont=1:1:4)
%qiu q  
%zj1=exp(-gm*q^2);
zj2=c*q^2*(1+gm*q^2)*exp(-gm*q^2)/t^2;
mm=q*r+(b0*r+2*c0/(t^3)-3*d0/t^4+4*e0/t^5)*q^2+(b*r+d/t^2)*q^3+ar*(-d/t^2)*q^6-2*q*zj2/t;%*c*q^3*(1+gm*q^2)*exp(-gm*q^2)/t^3;
nn=2*q*(b0*r*t-a0-c0/t^2+d0/t^3-e0/t^4)+3*(b*r*t-a-d/t)*q^2+6*ar*(a+d/t)*q^5+(3-2*gm*q^2+2*gm*q/(1+gm*q^2))*zj2;  %+2*c*q^2*gm*q+c*q^3*(1+gm*q^2)*(-2*gm*q))*exp(-gm*q^2)/(t^2);
cons1=-mm/nn;
aa=r*t+nn +2*(q^2-q)*gm*zj2/(1+gm*q^2);
cons2=1.0/aa;
%  su=sum( y[i]*(b[i]+2*c[i]*t+3*d[i]*t^2+4*e[i]*t^3+5*f[i]*t^4);
%Ae1=135.84210  ;% Be1=2.39359 ;  Ce1=-22.18007*10^(-4) ;  De1=57.40220*10^(-7) ;  Ee1=-372.79050*10^(-11);   Fe1=85.49685*10^(-12);
%Ae2=379.27660  ;%Be2=1.10899  ;  Ce2=-1.88512*10^(-4)  ;  De2=39.65580*10^(-7)  ; Ee2=-314.02090*10^(-11)  ; Fe2=80.08189*10^(-12);
%su=0.9515*(2.39359-2*22.18007*10^(-4)*t+3*57.40220*10^(-7)*t^2-4*372.79050*10^(-11)*t^3+5*85.49685*10^(-12)*t^4)+0.485*(1.10899-2*1.88512*10^(-4)*t+3*39.65580*10^(-7)*t^2-4*314.02090*10^(-11)*t^3+5*80.08189*10^(-12)*t^4);
su=2.815361035000001e+000 -4.403723961000001e-003*t+ 2.215537689000000e-005*t^2-2.028041189000000e-008*t^3+6.009498471250000e-010*t^4;
ee=su+(b0*r+8*c0/t^3-15*d0/t^4+24*e0/t^5)*q+0.5*(2*b*r+4*d/t^2)*q^2+0.2*ar*q^5*(-7*d/t^2)-2*c/(gm*t^3)*(3-(3+gm*q^2/2-gm^2*q^4)*exp(-gm*q^2));
gg=(b0*r*t-2*a0-4*c0/t^2+5*d0/t^3-6*e0/t^4)+q*(2*b*r*t-3*a-4*d/t)+ar*q^4*(6*ar+7*d/t);
hh=(5*gm*q+gm^2*q^3*(5-2*gm*q^2))*exp(-gm*q^2);
ff=gg+c*hh/(gm*t^2);
cc=gg+hh;
a11=w/q*cons1;
a12=w/q*cons2;
a13=1;b1=0;a21=0;a22=1.0/q;a23=w;
b2=-nd/dd*w^2/2-g*sin(theta);
a31=ee+ff*cons1;
a32=cc*cons2;
a33=w;
b3=-4*k*(t-t0)/(q*w*dd)-g*sin(theta);
%eff=[a11,a12,a13,b1;a21,a22,a23,b2;a31,a32,a33,b3]
%s1=0;s2=0;s3=0;
while(i)
x1=(b3-a32*s(2,cont)-a31*s(3,cont))/a33;
x2=(b2-a23*x1-a21*s(3,cont))/a22;
x3=(b1-a13*x1-a12*x2)/a11;
if((abs(x3-s(3,cont))<10^(-4))&(abs(x1-s(1,cont))<10^(-4))&(abs(x2-s(2,cont))<10^(-4))) break;
end
s(1,cont)=x1;s(2,cont)=x2;s(3,cont)=x3;%s1=x1=w,s2=x2=p,s3=x3=t
end
s(1,cont)=x1;s(2,cont)=x2;s(3,cont)=x3;
%s(1,cont)=s1;s(2,cont)=s2;s(3,cont)=s3;
%format long ,s1,s2,s3
%second
t=t+dx/2*s(3,cont);p=p+dx/2*s(2,cont);w=w+dx/2*s(1,cont);
cont=cont+1;
end
tnext=tin+(dx/6)*(s(3,1)+2*s(3,2)+2*s(3,3)+s(3,4));
pnext=pin+(dx/6)*(s(2,1)+2*s(2,2)+2*s(2,3)+s(2,4));
wnext=win+(dx/6)*(s(1,1)+2*s(1,2)+2*s(1,3)+s(1,4));
tin=tnext;pin=pnext;win=wnext;
dis=dis+dx;
end
format long e,tout=tnext,pout=pnext,wout=wnext