      
  swain <- function(N,p,d,TML,TMLi,conf.level=0.90) {
  
  # N: sample size
  # p: number of observed variables
  # d: number of degrees of freedom for model M_j
  # TML: chi-square likelihood-ratio test statistic for model M_j, i.e., T_{ML_j}
  # TMLi: chi-square test statistic for the independence model M_i, T_{ML_i}
  # conf.level: confidence level for interval estimation (default 0.90)

  n = N - 1
  t = p*(p+1)/2 - d                       # number of parameters
  pTML <- 1 - pchisq(TML,d)               # p-value of T_ML_j
  q = (sqrt(1+8*t)-1)/2                   # (11) 
  num = p*(2*p^2+3*p-1)-q*(2*q^2+3*q-1)   # numerator in (10)
  den = 12*d*(N-1)                        # denominator in (10)
  s <- 1 - num/den                        # Swain's correction factor s (10)
  TMLs <- s*TML                           # Swain-corrected test statistic (12)
  pTMLs <- 1 - pchisq(TMLs,d)             # p-value of T_MLs
  di = p*(p-1)/2                          # degrees of freedom of model M_i
  RMSEA <- sqrt(max((TML-d)/(n*d),0))     # Root Mean Square Error of Approximation
  RMSEAs <- sqrt(max((TMLs-d)/(n*d),0))   # Swain's estimate RMSEAs
  uppers = uniroot(function(l) ifelse(l>0,pchisq(ncp=l,q=TMLs,df=d)-(1-conf.level)/2 ,1),c(0,5000))$root
  RMSEAsu <- sqrt(max((uppers)/(n*d),0))
  lowers = uniroot(function(l) ifelse(l>0,pchisq(ncp=l,q=TMLs,df=d)-(1-(1-conf.level)/2),1),c(0,5000))$root
  RMSEAsl <- sqrt(max((lowers)/(n*d),0))
  upper = uniroot(function(l) ifelse(l>0,pchisq(ncp=l,q=TML,df=d)-(1-conf.level)/2,1),c(0,5000))$root
  RMSEAu <- sqrt(max((upper)/(n*d),0))
  lower = uniroot(function(l) ifelse(l>0,pchisq(ncp=l,q=TML,df=d)-(1-(1-conf.level)/2),1),c(0,5000))$root
  RMSEAl <- sqrt(max((lower)/(n*d),0))
  L = 0.05^2*n*d
  pclose <- 1-pchisq(ncp=L,q=TML,df=d)
  pcloses <- 1-pchisq(ncp=L,q=TMLs,df=d)
  Gamma1 <- p/(p+2*(TML-d)/n)            # Gamma_1 (Steiger) 
  Gamma1l <- p/(p+2*(upper)/n)
  Gamma1u <- p/(p+2*(lower)/n)
  Gamma1s <- p/(p+2*(TMLs-d)/n)          # Swain's estimate of Gamma_1   
  Gamma1sl <- p/(p+2*(uppers)/n)
  Gamma1su <- p/(p+2*(lowers)/n)
  TLI <- (TMLi/di - TML/d)/(TMLi/di - 1)   # Tucker-Lewis Index
  TLIs <- (TMLi/di - TMLs/d)/(TMLi/di -1)  
  tau = max(TML-d, 0)                  
  taus = max(TMLs-d, 0)                    # Swain's estimate TLIs   
  taui = max(TMLi-di, TML-d, 0)          
  CFI <- 1 - tau/taui                      # Comparative Fit Index (Bentler)
  CFIs <- 1 - taus/taui                    # Swain's estimate CFIs
  CIV2=100*conf.level

 return(list(chi=TML,chipvalue=pTML,swain=s,swainchi=TMLs,swainchipvalue=pTMLs,
 RMSEA=RMSEA,RMSEAlower=RMSEAl,RMSEAupper=RMSEAu,pclose=pclose,swainRMSEA=RMSEAs,
 swainRMSEAlower=RMSEAsl,swainRMSEAupper=RMSEAsu,swainpclose=pcloses,tli=TLI,
 swaintli=TLIs,cfi=CFI,swaincfi=CFIs))

}



