using Roots

export ng_zfactor
"""
    ng_zfactor(p,T;CH4=0.965,N2=0.003,CO2=0.006,C2H6=0.018,C3H8=0.045,H2O=0.,H2S=0.,
    H2=0.,CO=0.,O2=0.,iC4H10=0.001,nC4H10=0.001,iC5H12=0.0005,nC5H12=0.0003,nC6H14=0.0007,
    nC7H16=0.,nC8H18=0.,nNonane=0.,nDecane=0.,He=0.,Ar=0.)
Calculate the compress factor of specific pressure and temperature nature gas
    
# Auguments
...
- "p::Float64" : thre pressure,unit MPa
- "T::Float64" : the temperature,unit K
...
# Example
```
julia>ng_z_factor(6,300;CO2=0.006,N2=0.003,CH4=0.965,C2H6=0.018,C3H8=0.0045,
iC4H10=0.001,nC4H10=0.001,iC5H12=0.0005,nC5H12=0.0003,nC6H14=0.0007))
0.8953514530758864
```
"""
function ng_zfactor(p,T;
    CH4=0.965,
    N2=0.003,
    CO2=0.006,
    C2H6=0.018,
    C3H8=0.0045,
    H2O=0.,
    H2S=0.,
    H2=0.,
    CO=0.,
    O2=0.,
    iC4H10=0.001,
    nC4H10=0.001,
    iC5H12=0.0005,
    nC5H12=0.0003,
    nC6H14=0.0007,
    nC7H16=0.,
    nC8H18=0.,
    nNonane=0.,
    nDecane=0.,
    HE=0.,
    AR=0.)
    #参数检查
    #=
    ul=[12,338,1,0.2,0.2,0.1,0.035,0.015,0.005,0.001,0.0005,0.0005,0.1,0.03,0.005,0.00015]
    ll=[0;263;0.7;zeros(13)]
    pars=["pressure"=>p,"Temperature"=>T,"CH4"=>CH4,"N2"=>N2,"CO2"=>CO2,"C2H6"=>C2H6,"C3H8"=>C3H8,
          "C4H10"=>iC4H10+nC4H10,"C5H12"=>iC5H12+nC5H12,"C6"=>nC6H14,"C7"=>nC7H16,"C8+"=>nC8H18+nNonane+nDecane,
          "H2"=>H2,"CO"=>CO,"HE"=>HE,"H2O"=>H2O]
    for i=1:16
        !(ll[i]<pars[i].second<ul[i]) && error("$(pars[i].first) should be in range of $(ll[i]) - $(ul[i])")
    end
    =#


    x=[CH4,N2,CO2,C2H6,C3H8,H2O,H2S,H2,CO,O2,iC4H10,nC4H10,iC5H12,nC5H12,nC6H14,nC7H16,nC8H18,nNonane,nDecane,HE,AR]
    x=x./sum(x)
    #sum(x)!=100. && error("all componnet sum should be 100")
    N=21
    R=8.31451e-3
    #状态参数值
    a=[  0.153832600, 1.341953000,-2.998583000, -0.048312280, 0.375796500, -1.589575000,-0.053588470, 0.886594630,
        -0.710237040,-1.471722000, 1.321850350, -0.786659250, 0.2291290e-08,0.157672400,-0.436386400,-0.044081590,
        -0.003433888, 0.032059050, 0.024873550,  0.073322790,-0.001600573,  0.642470600,-0.416260100,-0.066899570,
         0.279179500,-0.696605100,-0.002860589, -0.008098836, 3.150547000,  0.007224479,-0.705752900, 0.534979200,
        -0.079314910,-1.418465000,-0.5999050e-16,0.105840200, 0.034317290, -0.007022847, 0.024955870, 0.042968180,
         0.746545300,-0.291961300, 7.294616000, -9.936757000,-0.005399808, -0.243256700, 0.049870160, 0.003733797,
         1.874951000, 0.002168144,-0.658716400,  0.000205518, 0.009776195, -0.020487080, 0.015573220, 0.006862415,
        -0.001226752, 0.002850908]
    b=Float64[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,5,5,5,5,5,6,6,7,7,8,8,8,9,9]
    k=Float64[0,0,0,0,0,0,0,0,0,0,0,0,3,2,2,2,4,4,0,0,2,2,2,4,4,4,4,0,1,1,2,2,3,3,4,4,4,0,0,2,2,2,4,4,0,2,2,4,4,0,2,0,2,1,2,2,2,2]
    c=Float64[0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,0,1,1,1,1,0,1,0,1,1,1,1,1,1]
    g=Float64[0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,0,0]
    f=Float64[0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    q=Float64[0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0,1]
    s=Float64[0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    w=Float64[0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
    u=Float64[0,0.5 ,1,3.5 ,-0.5 ,4.5 ,0.5 ,7.5 ,9.5 ,6,12,12.5 ,-6,2,3,2,2,11,-0.5 ,0.5 ,0,4,6,21,23,22,-1,-0.5 ,7,-1,6,
              4,1,9,-13,21,8,-0.5 ,0,2,7,9,22,23,1,9,3,8,23,1.5 ,5,-0.5 ,4,7,3,0,1,0]
    #特征参数值
    M=[16.0430 , 28.0135 , 44.0100 , 30.0700 , 44.0970 , 18.0153 , 34.0820 ,  2.0159 , 28.0100 , 31.9988 ,
       58.1230 , 58.1230 , 72.1500 , 72.1500 , 86.1770 ,100.2040 ,114.2310 ,128.2580 ,142.2850 ,  4.0026 ,
       39.9480 ]
    E=[151.318300 ,  99.737780 , 241.960600 , 244.166700 , 298.118300 , 514.015600 , 296.355000 ,  26.957940 ,
       105.534800 , 122.766700 , 324.068900 , 337.638900 , 365.599900 , 370.682300 , 402.636293 , 427.722630 ,
       450.325022 , 470.840891 , 489.558373 ,   2.610111 , 119.629900 ]
    K=[ 0.4619255 ,  0.4479153 , 0.4557489 , 0.5279209 , 0.5837490 , 0.3825868 , 0.4618263 , 0.3514916 ,
        0.4533894 ,  0.4186954 , 0.6406937 , 0.6341423 , 0.6738577 , 0.6798307 , 0.7175118 , 0.7525189 ,
        0.7849550 ,  0.8152731 , 0.8437826 , 0.3589888 , 0.4216551 ]
    G=[0, 0.027815 , 0.189065 , 0.079300 , 0.141239 ,  0.332500 , 0.088500 , 0.034369 , 0.038953 , 0.021000 ,
          0.256692 , 0.281835 , 0.332267 , 0.366911 , 0.289731 , 0.337542 , 0.383381 , 0.427354 , 0.469659 , 0, 0]
    Q=[zeros(2);0.69; zeros(2); 1.06775; 0.633276; zeros(14)]
    F=[zeros(7);1.0;zeros(13,1)]
    S=[zeros(5); 1.5822; 0.390; zeros(14)]
    W=[zeros(5); 1; zeros(15)]
    
    #交互作用参数值
    Ex=[1.0 0.971640 0.960644 1.0 0.994635 0.708218 0.931484 1.170520 0.990126 1.0 1.019530 0.989844  1.002350   0.999268   1.107274   0.880880  0.880973   0.881067   0.881161   1.0  1.0;
    	ones(1,2) 1.022740   0.970120   0.945939   0.746954   0.902271  1.086320   1.005710   1.021000   0.946914   0.973384  0.959340   0.945520  ones(1,7);
        ones(1,3) 0.925053   0.960237   0.849408   0.955052   1.281790  1.5   1  0.906849   0.897362   0.726255  0.859764   0.855134   0.831229   0.808310   0.786323  0.765171  ones(1,2);
        ones(1,4) 1.022560   0.693168   0.946871  1.164460   ones(1,3)  1.013060   1  1.00532  ones(1,7);
        ones(1,7) 1.034787   ones(1,3)  1.0049  ones(1,9);
        ones(1,21);
        ones(1,14) 1.008692   1.010126   1.011501  1.012821   1.014089  ones(1,2);
        ones(1,8) 1.1   1  1.3   1.3  ones(1,9);ones(13,21)]
 	
    Gx=ones(21,21);Gx[1,3]=0.807653;Gx[1,8]=1.957310;Gx[2,3]=0.982746; Gx[3,4]=0.370296; Gx[3,6]=1.673090 
    
    Ux=[1.0 0.886106   0.963827   1  0.990877   1.0  0.736833  1.156390   ones(1,3)  0.992291   1  1.003670   1.302576  1.191904   1.205769   1.219634   1.233498  ones(1,2);
    ones(1,2) 0.835058   0.816431   0.915502   1.0 0.993476   0.408838   ones(1,3)  0.993556  ones(1,9);
    ones(1,3) 0.969870   ones(1,2)  1.045290   1.0  0.9  ones(1,5)  1.066638   1.077634   1.088178   1.098291  1.108021  ones(1,2);
    ones(1,4) 1.065173   1.0  0.971926   1.616660   ones(1,2)  ones(1,4)*1.25  ones(1,7);
    ones(2,21);
    ones(1,14)  1.028973   1.033754  1.038338   1.042735   1.046966  ones(1,2);
    ones(14,21)]

    Kx=[ones(1,1) 1.003630   0.995933   1  1.007619   1  1.000080  1.023260   ones(1,3)  0.997596   1  1.002529   0.982962  0.983565   0.982707   0.981849   0.980991  ones(1,2);
    ones(1,2) 0.982361   1.007960   1  1  0.9425960  1.032270  ones(1,13);
    ones(1,3) 1.008510   ones(1,2)  1.00779  ones(1,7) 0.910183   0.895362   0.881152   0.867520   0.854406  ones(1,2);
    ones(1,4) 0.986893   1  0.999969   1.020340  ones(1,13);
    ones(2,21);
    ones(1,16)  0.968130   0.962870   0.957828   0.952441   0.948338;
    ones(14,21)]

    Eij=[Ex[i,j]*sqrt(E[i]*E[j]) for i=1:N,j=1:N]
    Gij=[Gx[i,j]*(G[i]+G[j])/2.0 for i=1:N,j=1:N]
    Bnij=[(Gij[i,j]+1.0-g[n])^g[n]*
          (Q[i]*Q[j]+1.0-q[n])^q[n]*
          (sqrt(F[i]*F[j])+1.0-f[n])^f[n]*
          (S[i]*S[j]+1.0-s[n])^s[n]*
          (W[i]*W[j]+1.0-w[n])^w[n] for n=1:18,i=1:N,j=1:N]
    B=sum(a[n]*T^(-u[n])*sum(x[i]*x[j]*Bnij[n,i,j]*Eij[i,j]^u[n]*(K[i]*K[j])^1.5 for i=1:N,j=1:N) for n=1:18)

    U0=(sum(x.*(E.^2.5))^2.0+sum(x[i]*x[j]*(Ux[i,j]^5.0-1.0)*(E[i]*E[j])^2.5 for i=1:N-1 for j=i+1:N))^0.2
    G0=sum(x.*G)+sum(x[i]*x[j]*(Gx[i,j]-1.0)*(G[i]+G[j]) for i=1:N-1 for j=i+1:N)
    Q0=sum(x.*Q)
    F0=sum(x.^2.0.*F)
    K0=(sum(x.*K.^2.5)^2.0+2*sum(x[i]*x[j]*(Kx[i,j]^5.0-1.0)*(K[i]*K[j])^2.5  for i=1:N-1 for j=i+1:N))^0.2

    Cn=[a[n]*(G0+1-g[n])^g[n]*(Q0^2+1.0-q[n])^q[n]*(F0+1.0-f[n])^f[n]*U0^u[n]*T^(-u[n]) for n=13:58]
    Cn=[zeros(12);Cn]

    function fun(ρm)
        ρr=K0^3.0*ρm
        z=1.0+B*ρm-ρr*sum(Cn[n] for n=13:18)+sum(Cn[n]*(b[n]-c[n]*k[n]*ρr^k[n])*ρr^b[n]*exp(-c[n]*ρr^k[n]) for n=13:58)
        ρm*R*T*z-p
    end

    ρm=find_zero(fun,40,rtol=1e-6)
    p/(ρm*R*T)
end
