text:
# unused

# I feel that dot fusion may be the culprit of slowness.
# Let me try loops instead.

#iMin,iMax,jMin,jMax,kMin,kMax = param.iMin, param.iMax, param.jMin, param.jMax, param.kMin, param.kMax

# face value
LState_XV::SubArray
RState_XV::SubArray
LState_YV::SubArray
RState_YV::SubArray
LState_ZV::SubArray
RState_ZV::SubArray

LState_XV = @view state_GV[iMin-1:iMax,jMin:jMax,kMin:kMax,:]
RState_XV = @view state_GV[iMin:iMax+1,jMin:jMax,kMin:kMax,:]
LState_YV = @view state_GV[iMin:iMax,jMin-1:jMax,kMin:kMax,:]
RState_YV = @view state_GV[iMin:iMax,jMin:jMax+1,kMin:kMax,:]
LState_ZV = @view state_GV[iMin:iMax,jMin:jMax,kMin-1:kMax,:]
RState_ZV = @view state_GV[iMin:iMax,jMin:jMax,kMin:kMax+1,:]

# Find dq_j = minmod{fwd diff, bwd diff, cntrl diff}
@views dqR_X =
   state_GV[iMin:iMax+2,  jMin:jMax,kMin:kMax,:] .-
   state_GV[iMin-1:iMax+1,jMin:jMax,kMin:kMax,:]
@views dqL_X =
   state_GV[iMin-1:iMax+1,jMin:jMax,kMin:kMax,:] .-
   state_GV[iMin-2:iMax  ,jMin:jMax,kMin:kMax,:]
@views dqC_X =
   state_GV[iMin:iMax+2,  jMin:jMax,kMin:kMax,:] .-
   state_GV[iMin-2:iMax  ,jMin:jMax,kMin:kMax,:]

dq_X = minmod.(dqR_X,dqL_X,dqC_X)

@views dqR_Y =
   state_GV[iMin:iMax,jMin:jMax+2,  kMin:kMax,:] .-
   state_GV[iMin:iMax,jMin-1:jMax+1,kMin:kMax,:]
@views dqL_Y =
   state_GV[iMin:iMax,jMin-1:jMax+1,kMin:kMax,:] .-
   state_GV[iMin:iMax,jMin-2:jMax  ,kMin:kMax,:]
@views dqC_Y =
   state_GV[iMin:iMax,jMin:jMax+2  ,kMin:kMax,:] .-
   state_GV[iMin:iMax,jMin-2:jMax  ,kMin:kMax,:]

dq_Y = minmod.(dqR_Y,dqL_Y,dqC_Y)

@views dqR_Z =
   state_GV[iMin:iMax,jMin:jMax,kMin:kMax+2,  :] .-
   state_GV[iMin:iMax,jMin:jMax,kMin-1:kMax+1,:]
@views dqL_Z =
   state_GV[iMin:iMax,jMin:jMax,kMin-1:kMax+1,:] .-
   state_GV[iMin:iMax,jMin:jMax,kMin-2:kMax  ,:]
@views dqC_Z =
   state_GV[iMin:iMax,jMin:jMax,kMin:kMax+2,  :] .-
   state_GV[iMin:iMax,jMin:jMax,kMin-2:kMax,  :]

dq_Z = minmod.(dqR_Z,dqL_Z,dqC_Z)

# Find dq_j = minmod{fwd diff, bwd diff}

@views dqR_X =
   state_GV[iMin:iMax+2  ,jMin:jMax,kMin:kMax,:] .-
   state_GV[iMin-1:iMax+1,jMin:jMax,kMin:kMax,:]
@views dqL_X =
   state_GV[iMin-1:iMax+1,jMin:jMax,kMin:kMax,:] .-
   state_GV[iMin-2:iMax  ,jMin:jMax,kMin:kMax,:]

dq_X = minmod.(dqR_X, dqL_X)

@views dqR_Y =
   state_GV[iMin:iMax,jMin:jMax+2,  kMin:kMax,:] .-
   state_GV[iMin:iMax,jMin-1:jMax+1,kMin:kMax,:]
@views dqL_Y =
   state_GV[iMin:iMax,jMin-1:jMax+1,kMin:kMax,:] .-
   state_GV[iMin:iMax,jMin-2:jMax  ,kMin:kMax,:]

dq_Y = minmod.(dqR_Y, dqL_Y)

@views dqR_Z =
   state_GV[iMin:iMax,jMin:jMax,kMin:kMax+2  ,:] .-
   state_GV[iMin:iMax,jMin:jMax,kMin-1:kMax+1,:]
@views dqL_Z =
   state_GV[iMin:iMax,jMin:jMax,kMin-1:kMax+1,:] .-
   state_GV[iMin:iMax,jMin:jMax,kMin-2:kMax  ,:]

dq_Z = minmod.(dqR_Z, dqL_Z)

# Linear interpolation onto edge centers
@views LState_XV = state_GV[iMin-1:iMax,jMin:jMax,kMin:kMax,:] .+
   0.5*dq_X[1:end-1,:,:,:]
@views RState_XV = state_GV[iMin:iMax+1,jMin:jMax,kMin:kMax,:] .-
   0.5*dq_X[2:end  ,:,:,:]
@views LState_YV = state_GV[iMin:iMax,jMin-1:jMax,kMin:kMax,:] .+
   0.5*dq_Y[:,1:end-1,:,:]
@views RState_YV = state_GV[iMin:iMax,jMin:jMax+1,kMin:kMax,:] .-
   0.5*dq_Y[:,2:end  ,:,:]
@views LState_ZV = state_GV[iMin:iMax,jMin:jMax,kMin-1:kMax,:] .+
   0.5*dq_Z[:,:,1:end-1,:]
@views RState_ZV = state_GV[iMin:iMax,jMin:jMax,kMin:kMax+1,:] .-
   0.5*dq_Z[:,:,2:end  ,:]

# face flux
bL = zeros(nI+1,nJ+1,nK+1)
bR = zeros(nI+1,nJ+1,nK+1)

bL = zeros(nI+1,nJ+1,nK+1)
bR = zeros(nI+1,nJ+1,nK+1)
uL = zeros(nI+1,nJ+1,nK+1)
uR = zeros(nI+1,nJ+1,nK+1)

bL = sum(LState_XV[:,:,:,B_].^2,dims=4)
bR = sum(RState_XV[:,:,:,B_].^2,dims=4)

bL = sum(LState_YV[:,:,:,B_].^2,dims=4) # takes time
bR = sum(RState_YV[:,:,:,B_].^2,dims=4)

bL = sum(LState_ZV[:,:,:,B_].^2,dims=4)
bR = sum(RState_ZV[:,:,:,B_].^2,dims=4) # takes time

@. LFlux_XV[:,:,:,Rho_] = LState_XV[:,:,:,Ux_]
@. RFlux_XV[:,:,:,Rho_] = RState_XV[:,:,:,Ux_]
@. LFlux_YV[:,:,:,Rho_] = LState_YV[:,:,:,Uy_]
@. RFlux_YV[:,:,:,Rho_] = RState_YV[:,:,:,Uy_]
@. LFlux_ZV[:,:,:,Rho_] = LState_ZV[:,:,:,Uz_]
@. RFlux_ZV[:,:,:,Rho_] = RState_ZV[:,:,:,Uz_]

@. LFlux_XV[:,:,:,Ux_] = LState_XV[:,:,:,Ux_]^2 / LState_XV[:,:,:,Rho_] + LState_XV[:,:,:,P_] + 0.5*$dropdims($sum(LState_XV[:,:,:,B_].^2,dims=4);dims=4) - LState_XV[:,:,:,Bx_]^2
@. RFlux_XV[:,:,:,Ux_] = RState_XV[:,:,:,Ux_]^2 / RState_XV[:,:,:,Rho_] + RState_XV[:,:,:,P_] + 0.5*$dropdims($sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4) - RState_XV[:,:,:,Bx_]^2
@. LFlux_XV[:,:,:,Uy_] = LState_XV[:,:,:,Ux_] * LState_XV[:,:,:,Uy_] / LState_XV[:,:,:,Rho_] - LState_XV[:,:,:,Bx_]*LState_XV[:,:,:,By_]
@. RFlux_XV[:,:,:,Uy_] = RState_XV[:,:,:,Ux_] * RState_XV[:,:,:,Uy_] / RState_XV[:,:,:,Rho_] - RState_XV[:,:,:,Bx_]*RState_XV[:,:,:,By_]
@. LFlux_XV[:,:,:,Uz_] = LState_XV[:,:,:,Ux_] * LState_XV[:,:,:,Uz_] / LState_XV[:,:,:,Rho_] - LState_XV[:,:,:,Bx_]*LState_XV[:,:,:,Bz_]
@. RFlux_XV[:,:,:,Uz_] = RState_XV[:,:,:,Ux_] * RState_XV[:,:,:,Uz_] / RState_XV[:,:,:,Rho_] - RState_XV[:,:,:,Bx_]*RState_XV[:,:,:,Bz_]

@. LFlux_YV[:,:,:,Ux_] = LState_YV[:,:,:,Uy_] * LState_YV[:,:,:,Ux_] / LState_YV[:,:,:,Rho_] - LState_YV[:,:,:,By_]*LState_YV[:,:,:,Bx_]
@. RFlux_YV[:,:,:,Ux_] = RState_YV[:,:,:,Uy_] * RState_YV[:,:,:,Ux_] / RState_YV[:,:,:,Rho_] - RState_YV[:,:,:,By_]*RState_YV[:,:,:,Bx_]
@. LFlux_YV[:,:,:,Uy_] = LState_YV[:,:,:,Uy_]^2 / LState_YV[:,:,:,Rho_] + LState_YV[:,:,:,P_] + 0.5*$dropdims($sum(LState_YV[:,:,:,B_].^2,dims=4);dims=4) - LState_YV[:,:,:,By_]^2
@. RFlux_YV[:,:,:,Uy_] = RState_YV[:,:,:,Uy_]^2 / RState_YV[:,:,:,Rho_] + RState_YV[:,:,:,P_] + 0.5*$dropdims($sum(RState_YV[:,:,:,B_].^2,dims=4);dims=4) - RState_YV[:,:,:,By_]^2
@. LFlux_YV[:,:,:,Uz_] = LState_YV[:,:,:,Uy_] * LState_YV[:,:,:,Uz_] / LState_YV[:,:,:,Rho_] - LState_YV[:,:,:,Bx_]*LState_YV[:,:,:,Bz_]
@. RFlux_YV[:,:,:,Uz_] = RState_YV[:,:,:,Uy_] * RState_YV[:,:,:,Uz_] / RState_YV[:,:,:,Rho_] - RState_YV[:,:,:,Bx_]*RState_YV[:,:,:,Bz_]

@. LFlux_ZV[:,:,:,Ux_] = LState_ZV[:,:,:,Uz_] * LState_ZV[:,:,:,Ux_] / LState_ZV[:,:,:,Rho_] - LState_ZV[:,:,:,Bz_]*LState_ZV[:,:,:,Bx_]
@. RFlux_ZV[:,:,:,Ux_] = RState_ZV[:,:,:,Uz_] * RState_ZV[:,:,:,Ux_] / RState_ZV[:,:,:,Rho_] - RState_ZV[:,:,:,Bz_]*RState_ZV[:,:,:,Bx_]
@. LFlux_ZV[:,:,:,Uy_] = LState_ZV[:,:,:,Uz_] * LState_ZV[:,:,:,Uy_] / LState_ZV[:,:,:,Rho_] - LState_ZV[:,:,:,Bz_]*LState_ZV[:,:,:,By_]
@. RFlux_ZV[:,:,:,Uy_] = RState_ZV[:,:,:,Uz_] * RState_ZV[:,:,:,Uy_] / RState_ZV[:,:,:,Rho_] - RState_ZV[:,:,:,Bz_]*RState_ZV[:,:,:,By_]
@. LFlux_ZV[:,:,:,Uz_] = LState_ZV[:,:,:,Uz_]^2 / LState_ZV[:,:,:,Rho_] + LState_ZV[:,:,:,P_] + 0.5*$dropdims($sum(LState_ZV[:,:,:,B_].^2,dims=4);dims=4) - LState_ZV[:,:,:,Bz_]^2
@. RFlux_ZV[:,:,:,Uz_] = RState_ZV[:,:,:,Uz_]^2 / RState_ZV[:,:,:,Rho_] + RState_ZV[:,:,:,P_] + 0.5*$dropdims($sum(RState_ZV[:,:,:,B_].^2,dims=4);dims=4) - RState_ZV[:,:,:,Bz_]^2

@. LFlux_XV[:,:,:,Bx_] = 0.0
@. RFlux_XV[:,:,:,Bx_] = 0.0
@. LFlux_XV[:,:,:,By_] = LState_XV[:,:,:,Ux_]*LState_XV[:,:,:,By_] - LState_XV[:,:,:,Bx_]*LState_XV[:,:,:,Uy_]
@. RFlux_XV[:,:,:,By_] = RState_XV[:,:,:,Ux_]*RState_XV[:,:,:,By_] - RState_XV[:,:,:,Bx_]*RState_XV[:,:,:,Uy_]
@. LFlux_XV[:,:,:,Bz_] = LState_XV[:,:,:,Ux_]*LState_XV[:,:,:,Bz_] - LState_XV[:,:,:,Bx_]*LState_XV[:,:,:,Uz_]
@. RFlux_XV[:,:,:,Bz_] = RState_XV[:,:,:,Ux_]*RState_XV[:,:,:,Bz_] - RState_XV[:,:,:,Bx_]*RState_XV[:,:,:,Uz_]

@. LFlux_YV[:,:,:,Bx_] = LState_YV[:,:,:,Uy_]*LState_YV[:,:,:,Bx_] - LState_YV[:,:,:,By_]*LState_YV[:,:,:,Ux_]
@. RFlux_YV[:,:,:,Bx_] = RState_YV[:,:,:,Uy_]*RState_YV[:,:,:,Bx_] - RState_YV[:,:,:,By_]*RState_YV[:,:,:,Ux_]
@. LFlux_YV[:,:,:,By_] = 0.0
@. RFlux_YV[:,:,:,By_] = 0.0
@. LFlux_YV[:,:,:,Bz_] = LState_YV[:,:,:,Uy_]*LState_YV[:,:,:,Bz_] - LState_YV[:,:,:,By_]*LState_YV[:,:,:,Uz_]
@. RFlux_YV[:,:,:,Bz_] = RState_YV[:,:,:,Uy_]*RState_YV[:,:,:,Bz_] - RState_YV[:,:,:,By_]*RState_YV[:,:,:,Uz_]

@. LFlux_ZV[:,:,:,Bx_] = LState_ZV[:,:,:,Uz_]*LState_ZV[:,:,:,Bx_] - LState_ZV[:,:,:,Bz_]*LState_ZV[:,:,:,Ux_]
@. RFlux_ZV[:,:,:,Bx_] = RState_ZV[:,:,:,Uz_]*RState_ZV[:,:,:,Bx_] - RState_ZV[:,:,:,Bz_]*RState_ZV[:,:,:,Ux_]
@. LFlux_ZV[:,:,:,By_] = LState_ZV[:,:,:,Uz_]*LState_ZV[:,:,:,By_] - LState_ZV[:,:,:,Bz_]*LState_ZV[:,:,:,Uy_]
@. RFlux_ZV[:,:,:,By_] = RState_ZV[:,:,:,Uz_]*RState_ZV[:,:,:,By_] - RState_ZV[:,:,:,Bz_]*RState_ZV[:,:,:,Uy_]
@. LFlux_ZV[:,:,:,Bz_] = 0.0
@. RFlux_ZV[:,:,:,Bz_] = 0.0

@. LFlux_XV[:,:,:,P_] = LState_XV[:,:,:,P_]*LState_XV[:,:,:,Ux_]/LState_XV[:,:,:,Rho_]
@. RFlux_XV[:,:,:,P_] = RState_XV[:,:,:,P_]*RState_XV[:,:,:,Ux_]/RState_XV[:,:,:,Rho_]
@. LFlux_YV[:,:,:,P_] = LState_YV[:,:,:,P_]*LState_YV[:,:,:,Uy_]/LState_YV[:,:,:,Rho_]
@. RFlux_YV[:,:,:,P_] = RState_YV[:,:,:,P_]*RState_YV[:,:,:,Uy_]/RState_YV[:,:,:,Rho_]
@. LFlux_ZV[:,:,:,P_] = LState_ZV[:,:,:,P_]*LState_ZV[:,:,:,Uz_]/LState_ZV[:,:,:,Rho_]
@. RFlux_ZV[:,:,:,P_] = RState_ZV[:,:,:,P_]*RState_ZV[:,:,:,Uz_]/RState_ZV[:,:,:,Rho_]

@. Flux_XV = 0.5 * (LFlux_XV + RFlux_XV)
@. Flux_YV = 0.5 * (LFlux_YV + RFlux_YV)
@. Flux_ZV = 0.5 * (LFlux_ZV + RFlux_ZV)

u = zeros(nI+1,nJ+1,nK+1)
b = zeros(nI+1,nJ+1,nK+1)
ub= zeros(nI+1,nJ+1,nK+1)

u = sum(LState_XV[:,:,:,U_].^2,dims=4)
b = sum(LState_XV[:,:,:,B_].^2,dims=4)
ub = sum(LState_XV[:,:,:,U_].*LState_XV[:,:,:,B_],dims=4)

@. LFlux_XV[:,:,:,E_] = LState_XV[:,:,:,Ux_]/LState_XV[:,:,:,Rho_]*
      ((LState_XV[:,:,:,P_]/(γ-1.0) + 0.5/LState_XV[:,:,:,Rho_]*
      u + 0.5*b + LState_XV[:,:,:,P_] + 0.5*b)) - ub*LState_XV[:,:,:,Bx_]

u = sum(RState_XV[:,:,:,U_].^2,dims=4) # takes time!
b = sum(RState_XV[:,:,:,B_].^2,dims=4) # takes time!
ub = sum(RState_XV[:,:,:,U_].*RState_XV[:,:,:,B_],dims=4)

@. RFlux_XV[:,:,:,E_] = RState_XV[:,:,:,Ux_]/RState_XV[:,:,:,Rho_]*
      ((RState_XV[:,:,:,P_]/(γ-1.0) + 0.5/RState_XV[:,:,:,Rho_]*
      $dropdims($sum(RState_XV[:,:,:,U_].^2,dims=4);dims=4) +
      0.5*$dropdims($sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4) + RState_XV[:,:,:,P_] +
      0.5*$dropdims($sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4))) -
      $dropdims($sum(RState_XV[:,:,:,U_].*RState_XV[:,:,:,B_],dims=4);dims=4)*
      RState_XV[:,:,:,Bx_]

u = sum(LState_YV[:,:,:,U_].^2,dims=4)
b = sum(LState_YV[:,:,:,B_].^2,dims=4)
ub = sum(LState_YV[:,:,:,U_].*LState_YV[:,:,:,B_],dims=4)

@. LFlux_YV[:,:,:,E_] = LState_YV[:,:,:,Uy_]/LState_YV[:,:,:,Rho_]*
      ((LState_YV[:,:,:,P_]/(γ-1.0) + 0.5/LState_YV[:,:,:,Rho_]*
      $dropdims($sum(LState_YV[:,:,:,U_].^2,dims=4);dims=4) +
      0.5*$dropdims($sum(LState_YV[:,:,:,B_].^2,dims=4);dims=4) + LState_YV[:,:,:,P_] +
      0.5*$dropdims($sum(LState_YV[:,:,:,B_].^2,dims=4);dims=4))) -
      $dropdims($sum(LState_YV[:,:,:,U_].*LState_YV[:,:,:,B_],dims=4);dims=4)*
      LState_YV[:,:,:,By_]

u = sum(RState_YV[:,:,:,U_].^2,dims=4) # dropdims takes time
b = sum(RState_YV[:,:,:,B_].^2,dims=4)
ub = sum(RState_YV[:,:,:,U_].*RState_YV[:,:,:,B_],dims=4) # takes time

@. RFlux_YV[:,:,:,E_] = RState_YV[:,:,:,Uy_]/RState_YV[:,:,:,Rho_]*
      ((RState_YV[:,:,:,P_]/(γ-1.0) +
      0.5/RState_YV[:,:,:,Rho_]*$dropdims($sum(RState_YV[:,:,:,U_].^2,dims=4);dims=4) +
      0.5*$dropdims($sum(RState_YV[:,:,:,B_].^2,dims=4);dims=4) + RState_YV[:,:,:,P_] +
      0.5*$dropdims($sum(RState_YV[:,:,:,B_].^2,dims=4);dims=4))) -
      $dropdims($sum(RState_YV[:,:,:,U_].*RState_YV[:,:,:,B_],dims=4);dims=4)*
      RState_YV[:,:,:,By_]

u = sum(LState_ZV[:,:,:,U_].^2,dims=4)
b = sum(LState_ZV[:,:,:,B_].^2,dims=4)
ub = sum(LState_ZV[:,:,:,U_].*LState_ZV[:,:,:,B_],dims=4)

@. LFlux_ZV[:,:,:,E_] = LState_ZV[:,:,:,Uz_]/LState_ZV[:,:,:,Rho_]*
      ((LState_ZV[:,:,:,P_]/(γ-1.0) + 0.5/LState_ZV[:,:,:,Rho_]*
      $dropdims($sum(LState_ZV[:,:,:,U_].^2,dims=4);dims=4) +
      0.5*$dropdims($sum(LState_ZV[:,:,:,B_].^2,dims=4);dims=4) + LState_ZV[:,:,:,P_] +
      0.5*$dropdims($sum(LState_ZV[:,:,:,B_].^2,dims=4);dims=4))) -
      $dropdims($sum(LState_ZV[:,:,:,U_].*LState_ZV[:,:,:,B_],dims=4);dims=4)*
      LState_ZV[:,:,:,Bz_]

u = sum(RState_ZV[:,:,:,U_].^2,dims=4)
b = sum(RState_ZV[:,:,:,B_].^2,dims=4)
ub = sum(RState_ZV[:,:,:,U_].*RState_ZV[:,:,:,B_],dims=4)

@. RFlux_ZV[:,:,:,E_] = RState_ZV[:,:,:,Uz_]/RState_ZV[:,:,:,Rho_]*
   ((RState_ZV[:,:,:,P_]/(γ-1.0) + 0.5/RState_ZV[:,:,:,Rho_]*
   $dropdims($sum(RState_ZV[:,:,:,U_].^2,dims=4);dims=4) +
   0.5*$dropdims($sum(RState_ZV[:,:,:,B_].^2,dims=4);dims=4) + RState_ZV[:,:,:,P_] +
   0.5*$dropdims($sum(RState_ZV[:,:,:,B_].^2,dims=4);dims=4))) -
   $dropdims($sum(RState_ZV[:,:,:,U_].*RState_ZV[:,:,:,B_],dims=4);dims=4)*
   RState_ZV[:,:,:,Bz_]

@. Flux_XV -= 0.5*Cmax_XF*(RState_XV - LState_XV)
@. Flux_YV -= 0.5*Cmax_YF*(RState_YV - LState_YV)
@. Flux_ZV -= 0.5*Cmax_ZF*(RState_ZV - LState_ZV)

uL = sum(LState_XV[:,:,:,U_].^2,dims=4)
bL = sum(LState_XV[:,:,:,B_].^2,dims=4)
uR = sum(RState_XV[:,:,:,U_].^2,dims=4)
bR = sum(RState_XV[:,:,:,B_].^2,dims=4)

@. Flux_XV[:,:,:,Rho_:Bz_] -=
0.5*Cmax_XF*(RState_XV[:,:,:,Rho_:Bz_] - LState_XV[:,:,:,Rho_:Bz_])
@. Flux_XV[:,:,:,E_] -= 0.5*Cmax_XF*(
(RState_XV[:,:,:,P_]/(γ-1.0) + 0.5/RState_XV[:,:,:,Rho_]*
$dropdims($sum(RState_XV[:,:,:,U_].^2,dims=4);dims=4) +
0.5*$dropdims($sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4)) -
(LState_XV[:,:,:,P_]/(γ-1.0) +
0.5/LState_XV[:,:,:,Rho_]*$dropdims($sum(LState_XV[:,:,:,U_].^2,dims=4);dims=4) +
0.5*$dropdims($sum(LState_XV[:,:,:,B_].^2,dims=4);dims=4)))

uL = sum(LState_YV[:,:,:,U_].^2,dims=4)
bL = sum(LState_YV[:,:,:,B_].^2,dims=4)
uR = sum(RState_YV[:,:,:,U_].^2,dims=4)
bR = sum(RState_YV[:,:,:,B_].^2,dims=4)

@. Flux_YV[:,:,:,Rho_:Bz_] -=
0.5*Cmax_YF*(RState_YV[:,:,:,Rho_:Bz_] - LState_YV[:,:,:,Rho_:Bz_])
@. Flux_YV[:,:,:,E_] -= 0.5*Cmax_YF*(
(RState_YV[:,:,:,P_]/(γ-1.0) + 0.5/RState_YV[:,:,:,Rho_]*
$dropdims($sum(RState_YV[:,:,:,U_].^2,dims=4);dims=4) +
0.5*$dropdims($sum(RState_YV[:,:,:,B_].^2,dims=4);dims=4)) -
(LState_YV[:,:,:,P_]/(γ-1.0) +
0.5/LState_YV[:,:,:,Rho_]*$dropdims($sum(LState_YV[:,:,:,U_].^2,dims=4);dims=4) +
0.5*$dropdims($sum(LState_YV[:,:,:,B_].^2,dims=4);dims=4)))

uL = sum(LState_ZV[:,:,:,U_].^2,dims=4)
bL = sum(LState_ZV[:,:,:,B_].^2,dims=4)
uR = sum(RState_ZV[:,:,:,U_].^2,dims=4)
bR = sum(RState_ZV[:,:,:,B_].^2,dims=4)

@. Flux_ZV[:,:,:,Rho_:Bz_] -=
0.5*Cmax_ZF*(RState_ZV[:,:,:,Rho_:Bz_] - LState_ZV[:,:,:,Rho_:Bz_])
@. Flux_ZV[:,:,:,E_] -= 0.5*Cmax_ZF*(
(RState_ZV[:,:,:,P_]/(γ-1.0) + 0.5/RState_ZV[:,:,:,Rho_]*
$dropdims($sum(RState_ZV[:,:,:,U_].^2,dims=4);dims=4) +
0.5*$dropdims($sum(RState_ZV[:,:,:,B_].^2,dims=4);dims=4)) -
(LState_ZV[:,:,:,P_]/(γ-1.0) + 0.5/LState_ZV[:,:,:,Rho_]*
$dropdims($sum(LState_ZV[:,:,:,U_].^2,dims=4);dims=4) +
0.5*$dropdims($sum(LState_ZV[:,:,:,B_].^2,dims=4);dims=4)))

Cs2_XF = γ.*(LS_XV[:,:,:,P_] .+ RS_XV[:,:,:,P_]) ./
   (LS_XV[:,:,:,Rho_] .+ RS_XV[:,:,:,Rho_])
Cs2_YF = γ.*(LS_YV[:,:,:,P_] .+ RS_YV[:,:,:,P_]) ./
   (LS_YV[:,:,:,Rho_] .+ RS_YV[:,:,:,Rho_])
Cs2_ZF = γ.*(LS_ZV[:,:,:,P_] .+ RS_ZV[:,:,:,P_]) ./
   (LS_ZV[:,:,:,Rho_] .+ RS_ZV[:,:,:,Rho_])

Ca2_XF = ( (LS_XV[:,:,:,Bx_] .+ RS_XV[:,:,:,Bx_]).^2 .+
   (LS_XV[:,:,:,By_] .+ RS_XV[:,:,:,By_]).^2 .+
   (LS_XV[:,:,:,Bz_] .+ RS_XV[:,:,:,Bz_]).^2 ) ./
   (2.0 .*(LS_XV[:,:,:,Rho_] .+ RS_XV[:,:,:,Rho_]))
Ca2_YF = ( (LS_YV[:,:,:,Bx_] .+ RS_YV[:,:,:,Bx_]).^2 .+
   (LS_YV[:,:,:,By_] .+ RS_YV[:,:,:,By_]).^2 .+
   (LS_YV[:,:,:,Bz_] .+ RS_YV[:,:,:,Bz_]).^2 ) ./
   (2.0 .*(LS_YV[:,:,:,Rho_] .+ RS_YV[:,:,:,Rho_]))
Ca2_ZF = ( (LS_ZV[:,:,:,Bx_] .+ RS_ZV[:,:,:,Bx_]).^2 .+
   (LS_ZV[:,:,:,By_] .+ RS_ZV[:,:,:,By_]).^2 .+
   (LS_ZV[:,:,:,Bz_] .+ RS_ZV[:,:,:,Bz_]).^2 ) ./
   (2.0 .*(LS_ZV[:,:,:,Rho_] .+ RS_ZV[:,:,:,Rho_]))

Can2_XF = ( (LS_XV[:,:,:,Bx_] .+ RS_XV[:,:,:,Bx_]).^2 ) ./
   (2.0 .*(LS_XV[:,:,:,Rho_] .+ RS_XV[:,:,:,Rho_]))
Can2_YF = ( (LS_YV[:,:,:,By_] .+ RS_YV[:,:,:,By_]).^2 ) ./
   (2.0 .*(LS_YV[:,:,:,Rho_] + RS_YV[:,:,:,Rho_]))
Can2_ZF = ( (LS_ZV[:,:,:,Bz_] .+ RS_ZV[:,:,:,Bz_]).^2 ) ./
   (2.0 .*(LS_ZV[:,:,:,Rho_] .+ RS_ZV[:,:,:,Rho_]))

@. Cmax_XF = 0.5 * abs(LS_XV[:,:,:,Ux_]/LS_XV[:,:,:,Rho_] +
   RS_XV[:,:,:,Ux_]/RS_XV[:,:,:,Rho_]) + sqrt( 0.5*(Cs2_XF + Ca2_XF +
   sqrt((Cs2_XF + Ca2_XF)^2 - 4.0 *Cs2_XF*Can2_XF)) )

@. Cmax_YF = 0.5 * abs(LS_YV[:,:,:,Uy_]/LS_YV[:,:,:,Rho_] +
   RS_YV[:,:,:,Uy_]/RS_YV[:,:,:,Rho_]) + sqrt( 0.5*(Cs2_YF + Ca2_YF +
   sqrt((Cs2_YF + Ca2_YF)^2 - 4.0 *Cs2_YF*Can2_YF)) )

@. Cmax_ZF = 0.5 * abs(LS_ZV[:,:,:,Uz_]/LS_ZV[:,:,:,Rho_] +
   RS_ZV[:,:,:,Uz_]/RS_ZV[:,:,:,Rho_]) + sqrt( 0.5*(Cs2_ZF + Ca2_ZF +
   sqrt((Cs2_ZF + Ca2_ZF)^2 - 4.0 *Cs2_ZF*Can2_ZF)) )

# get_speed_maxmin
Cs2_LXF = @. γ*LS_XV[:,:,:,P_]/LS_XV[:,:,:,Rho_]
Cs2_LYF = @. γ*LS_YV[:,:,:,P_]/LS_YV[:,:,:,Rho_]
Cs2_LZF = @. γ*LS_ZV[:,:,:,P_]/LS_ZV[:,:,:,Rho_]
Cs2_RXF = @. γ*RS_XV[:,:,:,P_]/RS_XV[:,:,:,Rho_]
Cs2_RYF = @. γ*RS_YV[:,:,:,P_]/RS_YV[:,:,:,Rho_]
Cs2_RZF = @. γ*RS_ZV[:,:,:,P_]/RS_ZV[:,:,:,Rho_]

Ca2_LXF = @. (LS_XV[:,:,:,Bx_]^2 + LS_XV[:,:,:,By_]^2 +
LS_XV[:,:,:,Bz_]^2 ) / LS_XV[:,:,:,Rho_]
Ca2_LYF = @. (LS_YV[:,:,:,Bx_]^2 + LS_YV[:,:,:,By_]^2 +
LS_YV[:,:,:,Bz_]^2 ) / LS_YV[:,:,:,Rho_]
Ca2_LZF = @. (LS_ZV[:,:,:,Bx_]^2 + LS_ZV[:,:,:,By_]^2 +
LS_ZV[:,:,:,Bz_]^2 ) / LS_ZV[:,:,:,Rho_]
Ca2_RXF = @. (RS_XV[:,:,:,Bx_]^2 + RS_XV[:,:,:,By_]^2 +
RS_XV[:,:,:,Bz_]^2 ) / RS_XV[:,:,:,Rho_]
Ca2_RYF = @. (RS_YV[:,:,:,Bx_]^2 + RS_YV[:,:,:,By_]^2 +
RS_YV[:,:,:,Bz_]^2 ) / RS_YV[:,:,:,Rho_]
Ca2_RZF = @. (RS_ZV[:,:,:,Bx_]^2 + RS_ZV[:,:,:,By_]^2 +
RS_ZV[:,:,:,Bz_]^2 ) / RS_ZV[:,:,:,Rho_]

Can2_LXF = @. LS_XV[:,:,:,Bx_]^2 / LS_XV[:,:,:,Rho_]
Can2_LYF = @. LS_YV[:,:,:,By_]^2 / LS_YV[:,:,:,Rho_]
Can2_LZF = @. LS_ZV[:,:,:,Bz_]^2 / LS_ZV[:,:,:,Rho_]
Can2_RXF = @. RS_XV[:,:,:,Bx_]^2 / RS_XV[:,:,:,Rho_]
Can2_RYF = @. RS_YV[:,:,:,By_]^2 / RS_YV[:,:,:,Rho_]
Can2_RZF = @. RS_ZV[:,:,:,Bz_]^2 / RS_ZV[:,:,:,Rho_]

u_LXF = @. LS_XV[:,:,:,Ux_]/LS_XV[:,:,:,Rho_]
u_LYF = @. LS_YV[:,:,:,Uy_]/LS_YV[:,:,:,Rho_]
u_LZF = @. LS_ZV[:,:,:,Uz_]/LS_ZV[:,:,:,Rho_]
u_RXF = @. RS_XV[:,:,:,Ux_]/RS_XV[:,:,:,Rho_]
u_RYF = @. RS_YV[:,:,:,Uy_]/RS_YV[:,:,:,Rho_]
u_RZF = @. RS_ZV[:,:,:,Uz_]/RS_ZV[:,:,:,Rho_]

c_LXF = @. sqrt( 0.5*(Cs2_LXF + Ca2_LXF +
sqrt((Cs2_LXF + Ca2_LXF)^2 - 4.0*Cs2_LXF*Can2_LXF)) )
c_LYF = @. sqrt( 0.5*(Cs2_LYF + Ca2_LYF +
sqrt((Cs2_LYF + Ca2_LYF)^2 - 4.0*Cs2_LYF*Can2_LYF)) )
c_LZF = @. sqrt( 0.5*(Cs2_LZF + Ca2_LZF +
sqrt((Cs2_LZF + Ca2_LZF)^2 - 4.0*Cs2_LZF*Can2_LZF)) )
c_RXF = @. sqrt( 0.5*(Cs2_RXF + Ca2_RXF +
sqrt((Cs2_RXF + Ca2_RXF)^2 - 4.0*Cs2_RXF*Can2_RXF)) )
c_RYF = @. sqrt( 0.5*(Cs2_RYF + Ca2_RYF +
sqrt((Cs2_RYF + Ca2_RYF)^2 - 4.0*Cs2_RYF*Can2_RYF)) )
c_RZF = @. sqrt( 0.5*(Cs2_RZF + Ca2_RZF +
sqrt((Cs2_RZF + Ca2_RZF)^2 - 4.0*Cs2_RZF*Can2_RZF)) )

sLmax_XF = max(0,u_LXF+c_LXF)
sLmin_XF = min(0,u_LXF-c_LXF)
sRmax_XF = max(0,u_RXF+c_RXF)
sRmin_XF = min(0,u_RXF-c_RXF)
sLmax_YF = max(0,u_LYF+c_LYF)
sLmin_YF = min(0,u_LYF-c_LYF)
sRmax_YF = max(0,u_RYF+c_RYF)
sRmin_YF = min(0,u_RYF-c_RYF)
sLmax_ZF = max(0,u_LZF+c_LZF)
sLmin_ZF = min(0,u_LZF-c_LZF)
sRmax_ZF = max(0,u_RZF+c_RZF)
sRmin_ZF = min(0,u_RZF-c_RZF)

smax_XF = max(sLmax_XF,sRmax_XF)
smin_XF = min(sLmin_XF,sRmin_XF)
smax_YF = max(sLmax_YF,sRmax_YF)
smin_YF = min(sLmin_YF,sRmin_YF)
smax_ZF = max(sLmax_ZF,sRmax_ZF)
smin_ZF = min(sLmin_ZF,sRmin_ZF)

# remove dropdims?
bL = dropdims(sum(LState_XV[:,:,:,B_].^2,dims=4);dims=4)
bR = dropdims(sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4)

# HLLE
@. Flux_XV -=
   0.5*(Cmax_XF + Cmin_XF)/(Cmax_XF - Cmin_XF)*(RFlux_XV - LFlux_XV) -
   Cmax_XF*Cmin_XF/(Cmax_XF - Cmin_XF)*(RState_XV - LState_XV)
@. Flux_YV -=
   0.5*(Cmax_YF + Cmin_YF)/(Cmax_YF - Cmin_YF)*(RFlux_YV - LFlux_YV) -
   Cmax_YF*Cmin_YF/(Cmax_YF - Cmin_YF)*(RState_YV - LState_YV)
@. Flux_ZV -=
   0.5*(Cmax_ZF + Cmin_ZF)/(Cmax_ZF - Cmin_ZF)*(RFlux_ZV - LFlux_ZV) -
   Cmax_ZF*Cmin_ZF/(Cmax_ZF - Cmin_ZF)*(RState_ZV - LState_ZV)

 @. Flux_XV[:,:,:,Rho_:Bz_] -= 0.5*(Cmax_XF + Cmin_XF)/(Cmax_XF - Cmin_XF)*
 (RFlux_XV[:,:,:,Rho_:Bz_] - LFlux_XV[:,:,:,Rho_:Bz_]) -
 Cmax_XF*Cmin_XF/(Cmax_XF - Cmin_XF)*(RState_XV[:,:,:,Rho_:Bz_] - LState_XV[:,:,:,Rho_:Bz_])

 @. Flux_XV[:,:,:,E_] -= 0.5*(Cmax_XF + Cmin_XF)/(Cmax_XF - Cmin_XF)*
 (RFlux_XV[:,:,:,E_] - LFlux_XV[:,:,:,E_]) -
 Cmax_XF*Cmin_XF/(Cmax_XF - Cmin_XF)* (
 (RState_XV[:,:,:,P_] / (γ-1) +
 0.5/RState_XV[:,:,:,Rho_]*$dropdims($sum(RState_XV[:,:,:,U_].^2,dims=4);dims=4) +
 0.5*$dropdims($sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4)) -
 (LState_XV[:,:,:,P_] / (γ-1) +
 0.5/LState_XV[:,:,:,Rho_]*$dropdims($sum(LState_XV[:,:,:,U_].^2,dims=4);dims=4) +
 0.5*$dropdims($sum(LState_XV[:,:,:,B_].^2,dims=4);dims=4)))

 @. Flux_YV[:,:,:,Rho_:Bz_] -= 0.5*(Cmax_YF + Cmin_YF)/(Cmax_YF - Cmin_YF)*
 (RFlux_YV[:,:,:,Rho_:Bz_] - LFlux_YV[:,:,:,Rho_:Bz_]) -
 Cmax_YF*Cmin_YF/(Cmax_YF - Cmin_YF)*
 (RState_YV[:,:,:,Rho_:Bz_] - LState_YV[:,:,:,Rho_:Bz_])
 @. Flux_YV[:,:,:,E_] -=
 0.5*(Cmax_YF + Cmin_YF)/(Cmax_YF - Cmin_YF)*
 (RFlux_YV[:,:,:,E_] - LFlux_YV[:,:,:,E_]) -
 Cmax_YF*Cmin_YF/(Cmax_YF - Cmin_YF)* (
 (RState_YV[:,:,:,P_] / (γ-1) +
 0.5/RState_YV[:,:,:,Rho_]*$dropdims($sum(RState_YV[:,:,:,U_].^2,dims=4);dims=4) +
 0.5*$dropdims($sum(RState_YV[:,:,:,B_].^2,dims=4);dims=4)) -
 (LState_YV[:,:,:,P_] / (γ-1) +
 0.5/LState_YV[:,:,:,Rho_]*$dropdims($sum(LState_YV[:,:,:,U_].^2,dims=4);dims=4) +
 0.5*$dropdims($sum(LState_YV[:,:,:,B_].^2,dims=4);dims=4)))

 @. Flux_ZV[:,:,:,Rho_:Bz_] -= 0.5*(Cmax_ZF + Cmin_ZF)/(Cmax_ZF - Cmin_ZF)*
 (RFlux_ZV[:,:,:,Rho_:Bz_] - LFlux_ZV[:,:,:,Rho_:Bz_]) -
 Cmax_ZF*Cmin_ZF/(Cmax_ZF - Cmin_ZF)*
 (RState_ZV[:,:,:,Rho_:Bz_] - LState_ZV[:,:,:,Rho_:Bz_])
 @. Flux_ZV[:,:,:,E_] -=
 0.5*(Cmax_ZF + Cmin_ZF)/(Cmax_ZF - Cmin_ZF)*
 (RFlux_ZV[:,:,:,E_] - LFlux_ZV[:,:,:,E_]) -
 Cmax_ZF*Cmin_ZF/(Cmax_ZF - Cmin_ZF)* (
 (RState_ZV[:,:,:,P_] / (γ-1) +
 0.5/RState_ZV[:,:,:,Rho_]*$dropdims($sum(RState_ZV[:,:,:,U_].^2,dims=4);dims=4) +
 0.5*$dropdims($sum(RState_ZV[:,:,:,B_].^2,dims=4);dims=4)) -
 (LState_ZV[:,:,:,P_] / (γ-1) +
 0.5/LState_ZV[:,:,:,Rho_]*$dropdims($sum(LState_ZV[:,:,:,U_].^2,dims=4);dims=4) +
 0.5*$dropdims($sum(LState_ZV[:,:,:,B_].^2,dims=4);dims=4)))

# source
for k = 1:nK, j = 1:nJ, i = 1:nI
   source_GV[i,j,k,Rho_] = 0.0
end

# This is slower than dot fusion, because it is wrong in indexing!
for k = 1:nK, j = 1:nJ, i = 1:nI
   source_GV[i,j,k,Ux_] = -state_GV[i,j,k,Bx_]*DivB[i+1,j+1,k+1]
   source_GV[i,j,k,Uy_] = -state_GV[i,j,k,By_]*DivB[i+1,j+1,k+1]
   source_GV[i,j,k,Uz_] = -state_GV[i,j,k,Bz_]*DivB[i+1,j+1,k+1]
   source_GV[i,j,k,Bx_] =  state_GV[i,j,k,Ux_]*DivB[i+1,j+1,k+1]
   source_GV[i,j,k,By_] =  state_GV[i,j,k,Uy_]*DivB[i+1,j+1,k+1]
   source_GV[i,j,k,Bz_] =  state_GV[i,j,k,Uz_]*DivB[i+1,j+1,k+1]
end

@. source_GV[:,:,:,U_] = -state_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_] *
   DivB[iMin:iMax,jMin:jMax,kMin:kMax]

@. source_GV[:,:,:,B_] =  state_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_] *
   DivB[iMin:iMax,jMin:jMax,kMin:kMax]

@. source_GV[:,:,:,P_] =-(γ-1.0)*state_GV[iMin:iMax,jMin:jMax,kMin:kMax,P_]*
   DivU[iMin:iMax,jMin:jMax,kMin:kMax]

@. source_GV[:,:,:,E_] =
   -$dropdims($sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_] *
   state_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4);dims=4) *
   DivB[iMin:iMax,jMin:jMax,kMin:kMax]

ub = dropdims(sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_].*
   state_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4);dims=4)
for k = 1:nK, j = 1:nJ, i = 1:nI
   source_GV[i,j,k,E_] = -ub[i,j,k]*DivB[i,j,k]
end

# test for simd usage, this may be faster.
@inbounds for k = kMin:kMax
   for j = jMin:jMax
      @simd for i = iMin:iMax
         source_GV[i-nG,j-nG,k-nG,Ux_] = -state_GV[i,j,k,Bx_]*DivB[i,j,k]
         source_GV[i-nG,j-nG,k-nG,Uy_] = -state_GV[i,j,k,By_]*DivB[i,j,k]
         source_GV[i-nG,j-nG,k-nG,Uz_] = -state_GV[i,j,k,Bz_]*DivB[i,j,k]
         source_GV[i-nG,j-nG,k-nG,Bx_] =  state_GV[i,j,k,Ux_]*DivB[i,j,k]
         source_GV[i-nG,j-nG,k-nG,By_] =  state_GV[i,j,k,Uy_]*DivB[i,j,k]
         source_GV[i-nG,j-nG,k-nG,Bz_] =  state_GV[i,j,k,Uz_]*DivB[i,j,k]
      end
   end
end

function calc_source!(param::Param, state_GV::Array{Float64,4},
   source_GV::Array{Float64,4})

   nVar, nG = param.nVar, param.nG
   x, y, z = param.x, param.y, param.z

   iMin, iMax, jMin, jMax, kMin, kMax =
   param.iMin, param.iMax, param.jMin, param.jMax, param.kMin, param.kMax

   # I can preallocate DivB and DivU!
   nI,nJ,nK = param.nI, param.nJ, param.nK

   # Calculate divergence of B using central difference
   DivB = divergence_ndgrid(x,y,z,state_GV[:,:,:,B_])

   source_GV[:,:,:,Rho_] .= 0.0

   @inbounds for k = kMin:kMax, j = jMin:jMax, i = iMin:iMax
      source_GV[i-nG,j-nG,k-nG,Ux_] = -state_GV[i,j,k,Bx_]*DivB[i,j,k]
      source_GV[i-nG,j-nG,k-nG,Uy_] = -state_GV[i,j,k,By_]*DivB[i,j,k]
      source_GV[i-nG,j-nG,k-nG,Uz_] = -state_GV[i,j,k,Bz_]*DivB[i,j,k]
      source_GV[i-nG,j-nG,k-nG,Bx_] =  state_GV[i,j,k,Ux_]*DivB[i,j,k]
      source_GV[i-nG,j-nG,k-nG,By_] =  state_GV[i,j,k,Uy_]*DivB[i,j,k]
      source_GV[i-nG,j-nG,k-nG,Bz_] =  state_GV[i,j,k,Uz_]*DivB[i,j,k]
   end

   # Calculate divergence of U
   DivU = divergence_ndgrid(x,y,z,state_GV[:,:,:,U_]./state_GV[:,:,:,Rho_])

   if !param.UseConservative
      @inbounds for k = kMin:kMax, j = jMin:jMax, i = iMin:iMax
         source_GV[i-nG,j-nG,k-nG,P_] = -(γ-1.0)*state_GV[i,j,k,P_]*DivU[i,j,k]
      end
   else
      @inbounds for k = kMin:kMax, j = jMin:jMax, i = iMin:iMax
         ub = state_GV[i,j,k,Ux_]*state_GV[i,j,k,Bx_] +
            state_GV[i,j,k,Uy_]*state_GV[i,j,k,By_] +
            state_GV[i,j,k,Uz_]*state_GV[i,j,k,Bz_]
         source_GV[i-nG,j-nG,k-nG,E_] = -ub*DivB[i,j,k]
      end
   end

   return source_GV
end

# Update state
@. stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,:] =
   state_GV[iMin:iMax,jMin:jMax,kMin:kMax,:] - dt*(
   (Flux_XV[2:end,:,:,:] - Flux_XV[1:end-1,:,:,:])/CellSize_D[1] +
   (Flux_YV[:,2:end,:,:] - Flux_YV[:,1:end-1,:,:])/CellSize_D[2] +
   (Flux_ZV[:,:,2:end,:] - Flux_ZV[:,:,1:end-1,:])/CellSize_D[3] -
   source_GV)

@. stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,Rho_:Bz_] =
   state_GV[iMin:iMax,jMin:jMax,kMin:kMax,Rho_:Bz_] - dt*(
   (Flux_XV[2:end,:,:,Rho_:Bz_] - Flux_XV[1:end-1,:,:,Rho_:Bz_])
   /CellSize_D[1] +
   (Flux_YV[:,2:end,:,Rho_:Bz_] - Flux_YV[:,1:end-1,:,Rho_:Bz_])
   /CellSize_D[2] +
   (Flux_ZV[:,:,2:end,Rho_:Bz_] - Flux_ZV[:,:,1:end-1,Rho_:Bz_])
   /CellSize_D[3] +
   source_GV[:,:,:,Rho_:Bz_])

@. state_GV[iMin:iMax,jMin:jMax,kMin:kMax,E_] =
   state_GV[iMin:iMax,jMin:jMax,kMin:kMax,P_]/(γ-1) +
   0.5/state_GV[iMin:iMax,jMin:jMax,kMin:kMax,Rho_]*
   $dropdims($sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_].^2,dims=4)
   ;dims=4) + 0.5*
   $dropdims($sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4);
   dims=4)

@. stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,E_] =
   state_GV[iMin:iMax,jMin:jMax,kMin:kMax,E_] - dt*(
   (Flux_XV[2:end,:,:,E_] - Flux_XV[1:end-1,:,:,E_])/CellSize_D[1] +
   (Flux_YV[:,2:end,:,E_] - Flux_YV[:,1:end-1,:,E_])/CellSize_D[2] +
   (Flux_ZV[:,:,2:end,E_] - Flux_ZV[:,:,1:end-1,E_])/CellSize_D[3] +
   source_GV[:,:,:,E_])

@. stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,P_] = (γ-1.0)*
   (stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,E_] -
   0.5/stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,Rho_]*
   $dropdims($sum(stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_].^2,
   dims=4);dims=4) - 0.5*
   $dropdims($sum(stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4)
   ;dims=4) )

u = sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_].^2,dims=4)
b = sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4) # takes time

u = sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_].^2,dims=4)
b = sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4)

   function update_state!(param::Param, state_GV::Array{Float64,4}, dt::Float64,
      faceFlux::FaceFlux, source_GV::Array{Float64,4})

      # In local timestepping, dt can be an array!

      Flux_XV = faceFlux.Flux_XV
      Flux_YV = faceFlux.Flux_YV
      Flux_ZV = faceFlux.Flux_ZV

      #stateNew_GV = similar(state_GV)
      stateNew_GV = copy(state_GV)

      CellSize_D = param.CellSize_D

      iMin, iMax, jMin, jMax, kMin, kMax =
      param.iMin, param.iMax, param.jMin, param.jMax, param.kMin, param.kMax
      nVar = param.nVar
      nI,nJ,nK = param.nI, param.nJ, param.nK

      if param.TypeGrid == "Cartesian"
         # No need for volume and face if the grid is uniform Cartesian
         if !param.UseConservative
            for iVar=1:nVar, k=1:nK, j=1:nJ, i=1:nI
               stateNew_GV[i+1,j+1,k+1,iVar] = state_GV[i+1,j+1,k+1,iVar] - dt*(
               (Flux_XV[i+1,j,k,iVar] - Flux_XV[i,j,k,iVar])/CellSize_D[1] +
               (Flux_YV[i,j+1,k,iVar] - Flux_YV[i,j,k,iVar])/CellSize_D[2] +
               (Flux_ZV[i,j,k+1,iVar] - Flux_ZV[i,j,k,iVar])/CellSize_D[3] +
               source_GV[i,j,k,iVar])
            end
         else
            for iVar=Rho_:Bz_, k=1:nK, j=1:nJ, i=1:nI
               stateNew_GV[i+1,j+1,k+1,iVar] = state_GV[i+1,j+1,k+1,iVar] - dt*(
               (Flux_XV[i+1,j,k,iVar] - Flux_XV[i,j,k,iVar])/CellSize_D[1] +
               (Flux_YV[i,j+1,k,iVar] - Flux_YV[i,j,k,iVar])/CellSize_D[2] +
               (Flux_ZV[i,j,k+1,iVar] - Flux_ZV[i,j,k+1,iVar])/CellSize_D[3] +
               source_GV[i,j,k,iVar])
            end

            u = dropdims(sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_].^2,dims=4);dims=4)
            b = dropdims(sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4);dims=4)
            for k=1:nK, j=1:nJ, i=1:nI
               state_GV[i+1,j+1,k+1,E_] = state_GV[i+1,j+1,k+1,P_]/(γ-1.0) +
                  0.5/state_GV[i+1,j+1,k+1,Rho_]*u[i,j,k] + 0.5*b[i,j,k]
            end

            for k=1:nK, j=1:nJ, i=1:nI
               stateNew_GV[i+1,j+1,k+1,E_] = state_GV[i+1,j+1,k+1,E_] - dt*(
               (Flux_XV[i+1,j,k,E_] - Flux_XV[i,j,k,E_])/CellSize_D[1] +
               (Flux_YV[i,j+1,k,E_] - Flux_YV[i,j,k,E_])/CellSize_D[2] +
               (Flux_ZV[i,j,k+1,E_] - Flux_ZV[i,j,k,E_])/CellSize_D[3] +
               source_GV[i,j,k,E_])
            end

            u = dropdims(sum(stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_].^2,dims=4);dims=4)
            b = dropdims(sum(stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4);dims=4)

            for k=1:nK, j=1:nJ, i=1:nI
               stateNew_GV[i+1,j+1,k+1,P_] = (γ-1.0)*(stateNew_GV[i+1,j+1,k+1,E_] -
                  0.5/stateNew_GV[i+1,j+1,k+1,Rho_]*u[i,j,k] - 0.5*b[i,j,k])
            end
         end
      else
         # Need volume and face
         stateNew_GV .= 0.0
      end

      state_GV .= stateNew_GV

   end

   function update_state!(param::Param, state_GV::Array{Float64,4}, dt::Float64,
      faceFlux::FaceFlux, source_GV::Array{Float64,4})

      # In local timestepping, dt can be an array!

      Flux_XV = faceFlux.Flux_XV
      Flux_YV = faceFlux.Flux_YV
      Flux_ZV = faceFlux.Flux_ZV

      #stateNew_GV = similar(state_GV)
      stateNew_GV = copy(state_GV)

      CellSize_D = param.CellSize_D

      iMin, iMax, jMin, jMax, kMin, kMax =
      param.iMin, param.iMax, param.jMin, param.jMax, param.kMin, param.kMax
      nVar = param.nVar
      nI,nJ,nK = param.nI, param.nJ, param.nK

      if param.TypeGrid == "Cartesian"
         # No need for volume and face if the grid is uniform Cartesian
         if !param.UseConservative
            @. stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,:] =
               state_GV[iMin:iMax,jMin:jMax,kMin:kMax,:] - dt*(
               (Flux_XV[2:end,:,:,:] - Flux_XV[1:end-1,:,:,:])/CellSize_D[1] +
               (Flux_YV[:,2:end,:,:] - Flux_YV[:,1:end-1,:,:])/CellSize_D[2] +
               (Flux_ZV[:,:,2:end,:] - Flux_ZV[:,:,1:end-1,:])/CellSize_D[3] +
               source_GV)
         else
            @. stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,Rho_:Bz_] =
               state_GV[iMin:iMax,jMin:jMax,kMin:kMax,Rho_:Bz_] - dt*(
               (Flux_XV[2:end,:,:,Rho_:Bz_] - Flux_XV[1:end-1,:,:,Rho_:Bz_])
               /CellSize_D[1] +
               (Flux_YV[:,2:end,:,Rho_:Bz_] - Flux_YV[:,1:end-1,:,Rho_:Bz_])
               /CellSize_D[2] +
               (Flux_ZV[:,:,2:end,Rho_:Bz_] - Flux_ZV[:,:,1:end-1,Rho_:Bz_])
               /CellSize_D[3] +
               source_GV[:,:,:,Rho_:Bz_])

            @. state_GV[iMin:iMax,jMin:jMax,kMin:kMax,E_] =
               state_GV[iMin:iMax,jMin:jMax,kMin:kMax,P_]/(γ-1) +
               0.5/state_GV[iMin:iMax,jMin:jMax,kMin:kMax,Rho_]*
               $dropdims($sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_].^2,dims=4)
               ;dims=4) + 0.5*
               $dropdims($sum(state_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4);
               dims=4)

            @. stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,E_] =
               state_GV[iMin:iMax,jMin:jMax,kMin:kMax,E_] - dt*(
               (Flux_XV[2:end,:,:,E_] - Flux_XV[1:end-1,:,:,E_])/CellSize_D[1] +
               (Flux_YV[:,2:end,:,E_] - Flux_YV[:,1:end-1,:,E_])/CellSize_D[2] +
               (Flux_ZV[:,:,2:end,E_] - Flux_ZV[:,:,1:end-1,E_])/CellSize_D[3] +
               source_GV[:,:,:,E_])

            @. stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,P_] = (γ-1.0)*
               (stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,E_] -
               0.5/stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,Rho_]*
               $dropdims($sum(stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,U_].^2,
               dims=4);dims=4) - 0.5*
               $dropdims($sum(stateNew_GV[iMin:iMax,jMin:jMax,kMin:kMax,B_],dims=4)
               ;dims=4) )
         end
      else
         # Need volume and face
         stateNew_GV .= 0.0
      end

      state_GV .= stateNew_GV

   end

   time_G .= CFL ./ (
   (nI > 1).*Cmax_XF./CellSize_D[1] .+
   (nJ > 1).*Cmax_YF./CellSize_D[2] .+
   (nK > 1).*Cmax_ZF./CellSize_D[3])

# Divergence
@. px[2:n-1,:,:] = (var[3:n,:,:,1] - var[1:n-2,:,:,1])/(hx[3:n] - hx[1:n-2])
@. qy[:,2:n-1,:] = (var[:,3:n,:,2] - var[:,1:n-2,:,2])/(hy[3:n] - hy[1:n-2])
@. rz[:,:,2:n-1] = (var[:,:,3:n,3] - var[:,:,1:n-2,3])/(hz[3:n] - hz[1:n-2])

# Boundary

# periodic
state_GV[1:nG,jMin:jMax,kMin:kMax,:] .= state_GV[iMax-nG+1:iMax,jMin:jMax,kMin:kMax,:]
state_GV[iMax+1:iMaxAll,jMin:jMax,kMin:kMax,:] .= state_GV[iMin:iMin+nG-1,jMin:jMax,kMin:kMax,:]
state_GV[iMin:iMax,1:nG,kMin:kMax,:] .= state_GV[iMin:iMax,jMax-nG+1:jMax,kMin:kMax,:]
state_GV[iMin:iMax,jMax+1:jMaxAll,kMin:kMax,:] .= state_GV[iMin:iMax,jMin:jMin+nG-1,kMin:kMax,:]
state_GV[iMin:iMax,jMin:jMax,1:nG,:] .= state_GV[iMin:iMax,jMin:jMax,kMax-nG+1:kMax,:]
state_GV[iMin:iMax,jMin:jMax,kMax+1:kMaxAll,:] .= state_GV[iMin:iMax,jMin:jMax,kMin:kMin+nG-1,:]
# float
state_GV[1:nG,jMin:jMax,kMin:kMax,:] .=
   state_GV[iMin:iMin,jMin:jMax,kMin:kMax,:]
state_GV[iMax+1:iMaxAll,jMin:jMax,kMin:kMax,:] .=
   state_GV[iMax:iMax,jMin:jMax,kMin:kMax,:]

state_GV[iMin:iMax,1:nG,kMin:kMax,:] .=
   state_GV[iMin:iMax,jMin:jMin,kMin:kMax,:]
state_GV[iMin:iMax,jMax+1:jMaxAll,kMin:kMax,:] .=
   state_GV[iMin:iMax,jMax:jMax,kMin:kMax,:]

state_GV[iMin:iMax,jMin:jMax,1:nG,:] .=
   state_GV[iMin:iMax,jMin:jMax,kMin:kMin,:]
state_GV[iMin:iMax,jMin:jMax,kMax+1:kMaxAll,:] .=
   state_GV[iMin:iMax,jMin:jMax,kMax:kMax,:]

# calc_timestep
Cmax_XG = max.(Cmax_XF[2:end,:,:],Cmax_XF[1:end-1,:,:]) # these allocations are not needed!
Cmax_YG = max.(Cmax_YF[:,2:end,:],Cmax_YF[:,1:end-1,:])
Cmax_ZG = max.(Cmax_ZF[:,:,2:end],Cmax_ZF[:,:,1:end-1])

time_G = CFL ./ (
   (nI > 1).*Cmax_XG./CellSize_D[1] .+
   (nJ > 1).*Cmax_YG./CellSize_D[2] .+
   (nK > 1).*Cmax_ZG./CellSize_D[3])

   # flux module using dot fusion.
   module Flux

   export calc_face_flux!, init_flux, SpeedFlux, FaceFlux, FaceFluxLR

   using ..Parameters: Param, Rho_, Ux_, Uy_, Uz_, Bx_, By_, Bz_, P_, E_, U_, B_
   using ..Parameters: γ
   using ..FaceValue: FaceState

   struct SpeedFlux
      Cmax_XF::Array{Float64,3}
      Cmax_YF::Array{Float64,3}
      Cmax_ZF::Array{Float64,3}
   end

   struct FaceFlux
      Flux_XV::Array{Float64,4}
      Flux_YV::Array{Float64,4}
      Flux_ZV::Array{Float64,4}
   end

   struct FaceFluxLR
      LFlux_XV::Array{Float64,4}
      RFlux_XV::Array{Float64,4}
      LFlux_YV::Array{Float64,4}
      RFlux_YV::Array{Float64,4}
      LFlux_ZV::Array{Float64,4}
      RFlux_ZV::Array{Float64,4}
   end

   function init_flux(param::Param)

      GridSize = param.GridSize
      nVar = param.nVar

      Flux_XV = Array{Float64,4}(undef,GridSize+[1,0,0]...,nVar)
      Flux_YV = Array{Float64,4}(undef,GridSize+[0,1,0]...,nVar)
      Flux_ZV = Array{Float64,4}(undef,GridSize+[0,0,1]...,nVar)

      LFlux_XV = Array{Float64,4}(undef,GridSize+[1,0,0]...,nVar)
      RFlux_XV = Array{Float64,4}(undef,GridSize+[1,0,0]...,nVar)
      LFlux_YV = Array{Float64,4}(undef,GridSize+[0,1,0]...,nVar)
      RFlux_YV = Array{Float64,4}(undef,GridSize+[0,1,0]...,nVar)
      LFlux_ZV = Array{Float64,4}(undef,GridSize+[0,0,1]...,nVar)
      RFlux_ZV = Array{Float64,4}(undef,GridSize+[0,0,1]...,nVar)

      Cmax_XF = Array{Float64,3}(undef,GridSize+[1,0,0]...)
      Cmax_YF = Array{Float64,3}(undef,GridSize+[0,1,0]...)
      Cmax_ZF = Array{Float64,3}(undef,GridSize+[0,0,1]...)

      faceFluxLR = FaceFluxLR(LFlux_XV,RFlux_XV,LFlux_YV,RFlux_YV,LFlux_ZV, RFlux_ZV)
      faceFlux = FaceFlux(Flux_XV, Flux_YV, Flux_ZV)
      speedFlux = SpeedFlux(Cmax_XF, Cmax_YF, Cmax_ZF)

      return faceFluxLR, faceFlux, speedFlux
   end

   function calc_face_flux!(param::Param, faceValue::FaceState,
      faceFlux::FaceFlux, speedFlux::SpeedFlux, faceFluxLR::FaceFluxLR)

      nVar = param.nVar
      GridSize = param.GridSize

      if param.Scheme == "Rusanov"
         get_physical_flux!(param, faceValue, faceFlux, speedFlux, faceFluxLR)
      elseif param.Scheme == "HLLE"
         get_HLLE_flux(param, faceValue, faceFlux, speedFlux)
      end

      return
   end

   function get_physical_flux!(param::Param, faceValue::FaceState,
      faceFlux::FaceFlux, speedFlux::SpeedFlux, faceFluxLR::FaceFluxLR)

      LState_XV = faceValue.LState_XV
      RState_XV = faceValue.RState_XV
      LState_YV = faceValue.LState_YV
      RState_YV = faceValue.RState_YV
      LState_ZV = faceValue.LState_ZV
      RState_ZV = faceValue.RState_ZV

      LFlux_XV = faceFluxLR.LFlux_XV
      RFlux_XV = faceFluxLR.RFlux_XV
      LFlux_YV = faceFluxLR.LFlux_YV
      RFlux_YV = faceFluxLR.RFlux_YV
      LFlux_ZV = faceFluxLR.LFlux_ZV
      RFlux_ZV = faceFluxLR.RFlux_ZV

      Flux_XV = faceFlux.Flux_XV
      Flux_YV = faceFlux.Flux_YV
      Flux_ZV = faceFlux.Flux_ZV

      GridSize = param.GridSize
      nVar = param.nVar

      nI,nJ,nK = param.nI, param.nJ, param.nK

      # Density flux
      @. LFlux_XV[:,:,:,Rho_] = LState_XV[:,:,:,Ux_]
      @. RFlux_XV[:,:,:,Rho_] = RState_XV[:,:,:,Ux_]
      @. LFlux_YV[:,:,:,Rho_] = LState_YV[:,:,:,Uy_]
      @. RFlux_YV[:,:,:,Rho_] = RState_YV[:,:,:,Uy_]
      @. LFlux_ZV[:,:,:,Rho_] = LState_ZV[:,:,:,Uz_]
      @. RFlux_ZV[:,:,:,Rho_] = RState_ZV[:,:,:,Uz_]

      # Momentum flux
      @. LFlux_XV[:,:,:,Ux_] = LState_XV[:,:,:,Ux_]^2 / LState_XV[:,:,:,Rho_] + LState_XV[:,:,:,P_] + 0.5*$dropdims($sum(LState_XV[:,:,:,B_].^2,dims=4);dims=4) - LState_XV[:,:,:,Bx_]^2
      @. RFlux_XV[:,:,:,Ux_] = RState_XV[:,:,:,Ux_]^2 / RState_XV[:,:,:,Rho_] + RState_XV[:,:,:,P_] + 0.5*$dropdims($sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4) - RState_XV[:,:,:,Bx_]^2
      @. LFlux_XV[:,:,:,Uy_] = LState_XV[:,:,:,Ux_] * LState_XV[:,:,:,Uy_] / LState_XV[:,:,:,Rho_] - LState_XV[:,:,:,Bx_]*LState_XV[:,:,:,By_]
      @. RFlux_XV[:,:,:,Uy_] = RState_XV[:,:,:,Ux_] * RState_XV[:,:,:,Uy_] / RState_XV[:,:,:,Rho_] - RState_XV[:,:,:,Bx_]*RState_XV[:,:,:,By_]
      @. LFlux_XV[:,:,:,Uz_] = LState_XV[:,:,:,Ux_] * LState_XV[:,:,:,Uz_] / LState_XV[:,:,:,Rho_] - LState_XV[:,:,:,Bx_]*LState_XV[:,:,:,Bz_]
      @. RFlux_XV[:,:,:,Uz_] = RState_XV[:,:,:,Ux_] * RState_XV[:,:,:,Uz_] / RState_XV[:,:,:,Rho_] - RState_XV[:,:,:,Bx_]*RState_XV[:,:,:,Bz_]

      @. LFlux_YV[:,:,:,Ux_] = LState_YV[:,:,:,Uy_] * LState_YV[:,:,:,Ux_] / LState_YV[:,:,:,Rho_] - LState_YV[:,:,:,By_]*LState_YV[:,:,:,Bx_]
      @. RFlux_YV[:,:,:,Ux_] = RState_YV[:,:,:,Uy_] * RState_YV[:,:,:,Ux_] / RState_YV[:,:,:,Rho_] - RState_YV[:,:,:,By_]*RState_YV[:,:,:,Bx_]
      @. LFlux_YV[:,:,:,Uy_] = LState_YV[:,:,:,Uy_]^2 / LState_YV[:,:,:,Rho_] + LState_YV[:,:,:,P_] + 0.5*$dropdims($sum(LState_YV[:,:,:,B_].^2,dims=4);dims=4) - LState_YV[:,:,:,By_]^2
      @. RFlux_YV[:,:,:,Uy_] = RState_YV[:,:,:,Uy_]^2 / RState_YV[:,:,:,Rho_] + RState_YV[:,:,:,P_] + 0.5*$dropdims($sum(RState_YV[:,:,:,B_].^2,dims=4);dims=4) - RState_YV[:,:,:,By_]^2
      @. LFlux_YV[:,:,:,Uz_] = LState_YV[:,:,:,Uy_] * LState_YV[:,:,:,Uz_] / LState_YV[:,:,:,Rho_] - LState_YV[:,:,:,Bx_]*LState_YV[:,:,:,Bz_]
      @. RFlux_YV[:,:,:,Uz_] = RState_YV[:,:,:,Uy_] * RState_YV[:,:,:,Uz_] / RState_YV[:,:,:,Rho_] - RState_YV[:,:,:,Bx_]*RState_YV[:,:,:,Bz_]

      @. LFlux_ZV[:,:,:,Ux_] = LState_ZV[:,:,:,Uz_] * LState_ZV[:,:,:,Ux_] / LState_ZV[:,:,:,Rho_] - LState_ZV[:,:,:,Bz_]*LState_ZV[:,:,:,Bx_]
      @. RFlux_ZV[:,:,:,Ux_] = RState_ZV[:,:,:,Uz_] * RState_ZV[:,:,:,Ux_] / RState_ZV[:,:,:,Rho_] - RState_ZV[:,:,:,Bz_]*RState_ZV[:,:,:,Bx_]
      @. LFlux_ZV[:,:,:,Uy_] = LState_ZV[:,:,:,Uz_] * LState_ZV[:,:,:,Uy_] / LState_ZV[:,:,:,Rho_] - LState_ZV[:,:,:,Bz_]*LState_ZV[:,:,:,By_]
      @. RFlux_ZV[:,:,:,Uy_] = RState_ZV[:,:,:,Uz_] * RState_ZV[:,:,:,Uy_] / RState_ZV[:,:,:,Rho_] - RState_ZV[:,:,:,Bz_]*RState_ZV[:,:,:,By_]
      @. LFlux_ZV[:,:,:,Uz_] = LState_ZV[:,:,:,Uz_]^2 / LState_ZV[:,:,:,Rho_] + LState_ZV[:,:,:,P_] + 0.5*$dropdims($sum(LState_ZV[:,:,:,B_].^2,dims=4);dims=4) - LState_ZV[:,:,:,Bz_]^2
      @. RFlux_ZV[:,:,:,Uz_] = RState_ZV[:,:,:,Uz_]^2 / RState_ZV[:,:,:,Rho_] + RState_ZV[:,:,:,P_] + 0.5*$dropdims($sum(RState_ZV[:,:,:,B_].^2,dims=4);dims=4) - RState_ZV[:,:,:,Bz_]^2

      # Magnetic flux
      @. LFlux_XV[:,:,:,Bx_] = 0.0
      @. RFlux_XV[:,:,:,Bx_] = 0.0
      @. LFlux_XV[:,:,:,By_] = LState_XV[:,:,:,Ux_]*LState_XV[:,:,:,By_] - LState_XV[:,:,:,Bx_]*LState_XV[:,:,:,Uy_]
      @. RFlux_XV[:,:,:,By_] = RState_XV[:,:,:,Ux_]*RState_XV[:,:,:,By_] - RState_XV[:,:,:,Bx_]*RState_XV[:,:,:,Uy_]
      @. LFlux_XV[:,:,:,Bz_] = LState_XV[:,:,:,Ux_]*LState_XV[:,:,:,Bz_] - LState_XV[:,:,:,Bx_]*LState_XV[:,:,:,Uz_]
      @. RFlux_XV[:,:,:,Bz_] = RState_XV[:,:,:,Ux_]*RState_XV[:,:,:,Bz_] - RState_XV[:,:,:,Bx_]*RState_XV[:,:,:,Uz_]

      @. LFlux_YV[:,:,:,Bx_] = LState_YV[:,:,:,Uy_]*LState_YV[:,:,:,Bx_] - LState_YV[:,:,:,By_]*LState_YV[:,:,:,Ux_]
      @. RFlux_YV[:,:,:,Bx_] = RState_YV[:,:,:,Uy_]*RState_YV[:,:,:,Bx_] - RState_YV[:,:,:,By_]*RState_YV[:,:,:,Ux_]
      @. LFlux_YV[:,:,:,By_] = 0.0
      @. RFlux_YV[:,:,:,By_] = 0.0
      @. LFlux_YV[:,:,:,Bz_] = LState_YV[:,:,:,Uy_]*LState_YV[:,:,:,Bz_] - LState_YV[:,:,:,By_]*LState_YV[:,:,:,Uz_]
      @. RFlux_YV[:,:,:,Bz_] = RState_YV[:,:,:,Uy_]*RState_YV[:,:,:,Bz_] - RState_YV[:,:,:,By_]*RState_YV[:,:,:,Uz_]

      @. LFlux_ZV[:,:,:,Bx_] = LState_ZV[:,:,:,Uz_]*LState_ZV[:,:,:,Bx_] - LState_ZV[:,:,:,Bz_]*LState_ZV[:,:,:,Ux_]
      @. RFlux_ZV[:,:,:,Bx_] = RState_ZV[:,:,:,Uz_]*RState_ZV[:,:,:,Bx_] - RState_ZV[:,:,:,Bz_]*RState_ZV[:,:,:,Ux_]
      @. LFlux_ZV[:,:,:,By_] = LState_ZV[:,:,:,Uz_]*LState_ZV[:,:,:,By_] - LState_ZV[:,:,:,Bz_]*LState_ZV[:,:,:,Uy_]
      @. RFlux_ZV[:,:,:,By_] = RState_ZV[:,:,:,Uz_]*RState_ZV[:,:,:,By_] - RState_ZV[:,:,:,Bz_]*RState_ZV[:,:,:,Uy_]
      @. LFlux_ZV[:,:,:,Bz_] = 0.0
      @. RFlux_ZV[:,:,:,Bz_] = 0.0

      # Pressure flux / energy flux
      if !param.UseConservative
         for k = 1:nK+1, j = 1:nJ, i = 1:nI
            @. LFlux_XV[:,:,:,P_] = LState_XV[:,:,:,P_]*LState_XV[:,:,:,Ux_]/LState_XV[:,:,:,Rho_]
            @. RFlux_XV[:,:,:,P_] = RState_XV[:,:,:,P_]*RState_XV[:,:,:,Ux_]/RState_XV[:,:,:,Rho_]
            @. LFlux_YV[:,:,:,P_] = LState_YV[:,:,:,P_]*LState_YV[:,:,:,Uy_]/LState_YV[:,:,:,Rho_]
            @. RFlux_YV[:,:,:,P_] = RState_YV[:,:,:,P_]*RState_YV[:,:,:,Uy_]/RState_YV[:,:,:,Rho_]
            @. LFlux_ZV[:,:,:,P_] = LState_ZV[:,:,:,P_]*LState_ZV[:,:,:,Uz_]/LState_ZV[:,:,:,Rho_]
            @. RFlux_ZV[:,:,:,P_] = RState_ZV[:,:,:,P_]*RState_ZV[:,:,:,Uz_]/RState_ZV[:,:,:,Rho_]
         end
      else
         # Currently use the same index for pressure/energy
         @. LFlux_XV[:,:,:,E_] = LState_XV[:,:,:,Ux_]/LState_XV[:,:,:,Rho_]*
               ((LState_XV[:,:,:,P_]/(γ-1) + 0.5/LState_XV[:,:,:,Rho_]*
               $dropdims($sum(LState_XV[:,:,:,U_].^2,dims=4);dims=4) +
               0.5*$dropdims($sum(LState_XV[:,:,:,B_].^2,dims=4);dims=4) + LState_XV[:,:,:,P_] +
               0.5*$dropdims($sum(LState_XV[:,:,:,B_].^2,dims=4);dims=4))) -
               $dropdims($sum(LState_XV[:,:,:,U_].*LState_XV[:,:,:,B_],dims=4);dims=4)*
               LState_XV[:,:,:,Bx_]

         @. RFlux_XV[:,:,:,E_] = RState_XV[:,:,:,Ux_]/RState_XV[:,:,:,Rho_]*
               ((RState_XV[:,:,:,P_]/(γ-1) + 0.5/RState_XV[:,:,:,Rho_]*
               $dropdims($sum(RState_XV[:,:,:,U_].^2,dims=4);dims=4) +
               0.5*$dropdims($sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4) + RState_XV[:,:,:,P_] +
               0.5*$dropdims($sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4))) -
               $dropdims($sum(RState_XV[:,:,:,U_].*RState_XV[:,:,:,B_],dims=4);dims=4)*
               RState_XV[:,:,:,Bx_]

         @. LFlux_YV[:,:,:,E_] = LState_YV[:,:,:,Uy_]/LState_YV[:,:,:,Rho_]*
               ((LState_YV[:,:,:,P_]/(γ-1) + 0.5/LState_YV[:,:,:,Rho_]*
               $dropdims($sum(LState_YV[:,:,:,U_].^2,dims=4);dims=4) +
               0.5*$dropdims($sum(LState_YV[:,:,:,B_].^2,dims=4);dims=4) + LState_YV[:,:,:,P_] +
               0.5*$dropdims($sum(LState_YV[:,:,:,B_].^2,dims=4);dims=4))) -
               $dropdims($sum(LState_YV[:,:,:,U_].*LState_YV[:,:,:,B_],dims=4);dims=4)*
               LState_YV[:,:,:,By_]

         @. RFlux_YV[:,:,:,E_] = RState_YV[:,:,:,Uy_]/RState_YV[:,:,:,Rho_]*
               ((RState_YV[:,:,:,P_]/(γ-1) +
               0.5/RState_YV[:,:,:,Rho_]*$dropdims($sum(RState_YV[:,:,:,U_].^2,dims=4);dims=4) +
               0.5*$dropdims($sum(RState_YV[:,:,:,B_].^2,dims=4);dims=4) + RState_YV[:,:,:,P_] +
               0.5*$dropdims($sum(RState_YV[:,:,:,B_].^2,dims=4);dims=4))) -
               $dropdims($sum(RState_YV[:,:,:,U_].*RState_YV[:,:,:,B_],dims=4);dims=4)*
               RState_YV[:,:,:,By_]

         @. LFlux_ZV[:,:,:,E_] = LState_ZV[:,:,:,Uz_]/LState_ZV[:,:,:,Rho_]*
               ((LState_ZV[:,:,:,P_]/(γ-1) + 0.5/LState_ZV[:,:,:,Rho_]*
               $dropdims($sum(LState_ZV[:,:,:,U_].^2,dims=4);dims=4) +
               0.5*$dropdims($sum(LState_ZV[:,:,:,B_].^2,dims=4);dims=4) + LState_ZV[:,:,:,P_] +
               0.5*$dropdims($sum(LState_ZV[:,:,:,B_].^2,dims=4);dims=4))) -
               $dropdims($sum(LState_ZV[:,:,:,U_].*LState_ZV[:,:,:,B_],dims=4);dims=4)*
               LState_ZV[:,:,:,Bz_]

         @. RFlux_ZV[:,:,:,E_] = RState_ZV[:,:,:,Uz_]/RState_ZV[:,:,:,Rho_]*
            ((RState_ZV[:,:,:,P_]/(γ-1) + 0.5/RState_ZV[:,:,:,Rho_]*
            $dropdims($sum(RState_ZV[:,:,:,U_].^2,dims=4);dims=4) +
            0.5*$dropdims($sum(RState_ZV[:,:,:,B_].^2,dims=4);dims=4) + RState_ZV[:,:,:,P_] +
            0.5*$dropdims($sum(RState_ZV[:,:,:,B_].^2,dims=4);dims=4))) -
            $dropdims($sum(RState_ZV[:,:,:,U_].*RState_ZV[:,:,:,B_],dims=4);dims=4)*
            RState_ZV[:,:,:,Bz_]
      end

      # Collect all the physical fluxes
      @. Flux_XV = 0.5 * (LFlux_XV + RFlux_XV)
      @. Flux_YV = 0.5 * (LFlux_YV + RFlux_YV)
      @. Flux_ZV = 0.5 * (LFlux_ZV + RFlux_ZV)

      # Cmax_XF should be passed in?
      add_numerical_flux!(param, faceValue, faceFlux, speedFlux)
   end

   function add_numerical_flux!(param::Param, faceValue::FaceState,
      faceFlux::FaceFlux, speedFlux::SpeedFlux)

      nI,nJ,nK = param.nI, param.nJ, param.nK
      nVar = param.nVar

      LState_XV = faceValue.LState_XV
      RState_XV = faceValue.RState_XV
      LState_YV = faceValue.LState_YV
      RState_YV = faceValue.RState_YV
      LState_ZV = faceValue.LState_ZV
      RState_ZV = faceValue.RState_ZV

      Flux_XV = faceFlux.Flux_XV
      Flux_YV = faceFlux.Flux_YV
      Flux_ZV = faceFlux.Flux_ZV

      Cmax_XF = speedFlux.Cmax_XF
      Cmax_YF = speedFlux.Cmax_YF
      Cmax_ZF = speedFlux.Cmax_ZF

      if param.Scheme == "Rusanov"
         get_speed_max!(param, faceValue, speedFlux)

         if !param.UseConservative
            for iVar = 1:nVar, k = 1:nK, j = 1:nJ, i = 1:nI+1
               Flux_XV[i,j,k,iVar] -= 0.5*Cmax_XF[i,j,k]*(RState_XV[i,j,k,iVar] -
                  LState_XV[i,j,k,iVar])
               Flux_YV[i,j,k,iVar] -= 0.5*Cmax_YF[i,j,k]*(RState_YV[i,j,k,iVar] -
                  LState_YV[i,j,k,iVar])
               Flux_ZV[i,j,k,iVar] -= 0.5*Cmax_ZF[i,j,k]*(RState_ZV[i,j,k,iVar] -
                  LState_ZV[i,j,k,iVar])
            end
         else
            # If I solve energy equation instead of pressure, there's
            # duplicate calculation above, even though the expression
            # looks compact. That's why I use an if-else statement
            for iVar = Rho_:Bz_, k = 1:nK, j = 1:nJ, i = 1:nI+1
               Flux_XV[i,j,k,iVar] -= 0.5*Cmax_XF[i,j,k]*(RState_XV[i,j,k,iVar] -
                  LState_XV[i,j,k,iVar])
            end

            uL = dropdims(sum(LState_XV[:,:,:,U_].^2,dims=4);dims=4)
            bL = dropdims(sum(LState_XV[:,:,:,B_].^2,dims=4);dims=4)
            uR = dropdims(sum(RState_XV[:,:,:,U_].^2,dims=4);dims=4)
            bR = dropdims(sum(RState_XV[:,:,:,B_].^2,dims=4);dims=4)

            for k = 1:nK, j = 1:nJ, i = 1:nI+1
               Flux_XV[i,j,k,E_] -= 0.5*Cmax_XF[i,j,k]*(
                  (RState_XV[i,j,k,P_]/(γ-1) + 0.5/RState_XV[i,j,k,Rho_]*uR[i,j,k]
                  + 0.5*bR[i,j,k]) -
                  (LState_XV[i,j,k,P_]/(γ-1) + 0.5/LState_XV[i,j,k,Rho_]*uL[i,j,k]
                  + 0.5*bL[i,j,k]))
            end

            for iVar = Rho_:Bz_, k = 1:nK, j = 1:nJ+1, i = 1:nI
               Flux_YV[i,j,k,iVar] -= 0.5*Cmax_YF[i,j,k]*(RState_YV[i,j,k,iVar] -
                  LState_YV[i,j,k,iVar])
            end

            uL = dropdims(sum(LState_YV[:,:,:,U_].^2,dims=4);dims=4)
            bL = dropdims(sum(LState_YV[:,:,:,B_].^2,dims=4);dims=4)
            uR = dropdims(sum(RState_YV[:,:,:,U_].^2,dims=4);dims=4)
            bR = dropdims(sum(RState_YV[:,:,:,B_].^2,dims=4);dims=4)

            for k = 1:nK, j = 1:nJ+1, i = 1:nI
               Flux_YV[i,j,k,E_] -= 0.5*Cmax_YF[i,j,k]*(
                  (RState_YV[i,j,k,P_]/(γ-1) + 0.5/RState_YV[i,j,k,Rho_]*uR[i,j,k]
                  + 0.5*bR[i,j,k]) -
                  (LState_YV[i,j,k,P_]/(γ-1) + 0.5/LState_YV[i,j,k,Rho_]*uL[i,j,k]
                  + 0.5*bL[i,j,k]))
            end

            for iVar = Rho_:Bz_, k = 1:nK+1, j = 1:nJ, i = 1:nI
               Flux_ZV[i,j,k,iVar] -= 0.5*Cmax_ZF[i,j,k]*(RState_ZV[i,j,k,iVar] -
                  LState_ZV[i,j,k,iVar])
            end

            uL = dropdims(sum(LState_ZV[:,:,:,U_].^2,dims=4);dims=4)
            bL = dropdims(sum(LState_ZV[:,:,:,B_].^2,dims=4);dims=4)
            uR = dropdims(sum(RState_ZV[:,:,:,U_].^2,dims=4);dims=4)
            bR = dropdims(sum(RState_ZV[:,:,:,B_].^2,dims=4);dims=4)

            for k = 1:nK+1, j = 1:nJ, i = 1:nI
               Flux_ZV[i,j,k,E_] -= 0.5*Cmax_ZF[i,j,k]*(
                  (RState_ZV[i,j,k,P_]/(γ-1) + 0.5/RState_ZV[i,j,k,Rho_]*uR[i,j,k]
                  + 0.5*bR[i,j,k]) -
                  (LState_ZV[i,j,k,P_]/(γ-1) + 0.5/LState_ZV[i,j,k,Rho_]*uL[i,j,k]
                  + 0.5*bL[i,j,k]))
            end

         end
      end
   end

   function get_HLLE_flux()

   end

   ```Calculate the maximum speed in each direction.```
   function get_speed_max!(param::Param,faceValue::FaceState,speedFlux::SpeedFlux)
      nI, nJ, nK = param.nI, param.nJ, param.nK
      # Aliases
      LS_XV, RS_XV = faceValue.LState_XV, faceValue.RState_XV
      LS_YV, RS_YV = faceValue.LState_YV, faceValue.RState_YV
      LS_ZV, RS_ZV = faceValue.LState_ZV, faceValue.RState_ZV

      Cmax_XF = speedFlux.Cmax_XF
      Cmax_YF = speedFlux.Cmax_YF
      Cmax_ZF = speedFlux.Cmax_ZF

      for k = 1:nK, j = 1:nJ, i = 1:nI+1
         Cs2_XF = γ*(LS_XV[i,j,k,P_] + RS_XV[i,j,k,P_]) /
            (LS_XV[i,j,k,Rho_] + RS_XV[i,j,k,Rho_])
         Ca2_XF = ( (LS_XV[i,j,k,Bx_] + RS_XV[i,j,k,Bx_])^2 +
            (LS_XV[i,j,k,By_] + RS_XV[i,j,k,By_])^2 +
            (LS_XV[i,j,k,Bz_] + RS_XV[i,j,k,Bz_])^2 ) /
            (2.0*(LS_XV[i,j,k,Rho_] + RS_XV[i,j,k,Rho_]))
         Can2_XF = ( (LS_XV[i,j,k,Bx_] .+ RS_XV[i,j,k,Bx_])^2 ) /
            (2.0 *(LS_XV[i,j,k,Rho_] + RS_XV[i,j,k,Rho_]))
         Cmax_XF[i,j,k] = 0.5 * abs(LS_XV[i,j,k,Ux_]/LS_XV[i,j,k,Rho_] +
            RS_XV[i,j,k,Ux_]/RS_XV[i,j,k,Rho_]) + sqrt( 0.5*(Cs2_XF + Ca2_XF +
            sqrt((Cs2_XF + Ca2_XF)^2 - 4*Cs2_XF*Can2_XF)) )
      end

      for k = 1:nK, j = 1:nJ+1, i = 1:nI
         Cs2_YF = γ*(LS_YV[i,j,k,P_] + RS_YV[i,j,k,P_]) /
            (LS_YV[i,j,k,Rho_] + RS_YV[i,j,k,Rho_])
         Ca2_YF = ( (LS_YV[i,j,k,Bx_] + RS_YV[i,j,k,Bx_])^2 +
            (LS_YV[i,j,k,By_] + RS_YV[i,j,k,By_])^2 +
            (LS_YV[i,j,k,Bz_] + RS_YV[i,j,k,Bz_])^2 ) /
            (2.0 *(LS_YV[i,j,k,Rho_] + RS_YV[i,j,k,Rho_]))
         Can2_YF = ( (LS_YV[i,j,k,By_] + RS_YV[i,j,k,By_])^2 ) /
            (2.0 *(LS_YV[i,j,k,Rho_] + RS_YV[i,j,k,Rho_]))
         Cmax_YF[i,j,k] = 0.5 * abs(LS_YV[i,j,k,Uy_]/LS_YV[i,j,k,Rho_] +
            RS_YV[i,j,k,Uy_]/RS_YV[i,j,k,Rho_]) + sqrt( 0.5*(Cs2_YF + Ca2_YF +
            sqrt((Cs2_YF + Ca2_YF)^2 - 4*Cs2_YF*Can2_YF)) )
      end

      for k = 1:nK+1, j = 1:nJ, i = 1:nI
         Cs2_ZF = γ*(LS_ZV[i,j,k,P_] + RS_ZV[i,j,k,P_]) /
            (LS_ZV[i,j,k,Rho_] + RS_ZV[i,j,k,Rho_])
         Ca2_ZF = ( (LS_ZV[i,j,k,Bx_] + RS_ZV[i,j,k,Bx_])^2 +
            (LS_ZV[i,j,k,By_] + RS_ZV[i,j,k,By_])^2 +
            (LS_ZV[i,j,k,Bz_] + RS_ZV[i,j,k,Bz_])^2 ) /
            (2.0 *(LS_ZV[i,j,k,Rho_] + RS_ZV[i,j,k,Rho_]))
         Can2_ZF = ( (LS_ZV[i,j,k,Bz_] + RS_ZV[i,j,k,Bz_])^2 ) /
            (2.0 *(LS_ZV[i,j,k,Rho_] + RS_ZV[i,j,k,Rho_]))
         Cmax_ZF[i,j,k] = 0.5 * abs(LS_ZV[i,j,k,Uz_]/LS_ZV[i,j,k,Rho_] +
            RS_ZV[i,j,k,Uz_]/RS_ZV[i,j,k,Rho_]) + sqrt( 0.5*(Cs2_ZF + Ca2_ZF +
            sqrt((Cs2_ZF + Ca2_ZF)^2 - 4*Cs2_ZF*Can2_ZF)) )
      end

   end

   function get_speed_maxmin(faceValue::FaceState)

   end

   end


text:
# ---------------------------------
#        Humid Air Properties
# ---------------------------------

"""
    HAPropsSI(output::AbstractString, name1::AbstractString, value1::Real, name2::AbstractString, value2::Real, name3::AbstractString, value3::Real)

DLL wrapper of the HAPropsSI function.

# Arguments
* `output`: Output name for desired property, accepetd values are listed in the following table
* `name1`, `name2`, `name3`: Input name of  given state values, one must be "P"

# Note
Here, all outputs calculated as functions of these three: "Y"(Water mole fraction), "T" and "P", as "P" is mandatory so more perfomance achived when "Y" and "T" is given (or at least one of them).

Parameter(s) name     |Description                                 |Unit      |Formula
:---------------------|:-------------------------------------------|:---------|:-------------------------------
Omega HumRat W        |Humidity ratio                              |          |
psi_w Y               |Water mole fraction                         |mol_w/mol |
Tdp T_dp DewPoint D   |Dew point temperature                       |K         |
Twb T_wb WetBulb B    |Wet bulb temperature                        |K         |
Enthalpy H Hda        |Enthalpy                                    |J/kg_da   |
Hha                   |Enthalpy per kg of humid air                |J/kg_ha   |
InternalEnergy U Uda  |Internal energy                             |J/kg_da   |
Uha                   |Internal energy per kg of humid air         |J/kg_ha   |
Entropy S Sda         |Entropy                                     |J/kg_da/K |
Sha                   |Entropy per kg of humid air                 |J/kg_ha/K |
RH RelHum R           |Relative humidity                           |          |
Tdb T_db T            |Temperature                                 |K         |
P                     |Pressure                                    |Pa        |
V Vda                 |Specific volume                             |m^3/kg_da |``MolarVolume*(1+HumidityRatio)/M_ha``
Vha                   |Specific volume per kg of humid air         |m^3/kg_ha |``MolarVolume/M_ha``
mu Visc M             |Viscosity                                   |          |
k Conductivity K      |Conductivity                                |          |
C cp                  |Heat cap. const. press.                     |J/kg_da/K |
Cha cp_ha             |Heat cap. const. press. per kg of humid air |J/kg_ha/K |
CV                    |Heat cap. const. vol.                       |J/kg_da/K |
CVha cv_ha            |Heat cap. const. vol. per kg of humid air   |J/kg_ha/K |
P_w                   |Partial pressure of water                   |          |
isentropic_exponent   |Isentropic exponent                         |          |
speed_of_sound        |Speed of sound                              |          |``sqrt(1/M_ha*cp/cv*dpdrho__constT)``
Z                     |Compressibility factor                      |          |``P*MolarVolume/(R*T)``

# Example
```julia
# Enthalpy (J per kg dry air) as a function of temperature, pressure,
# and relative humidity at STP
julia> h = HAPropsSI("H", "T", 298.15, "P", 101325, "R", 0.5)
50423.45039247604
# Temperature of saturated air at the previous enthalpy
julia> T = HAPropsSI("T", "P", 101325, "H", h, "R", 1.0)
290.9620891952412
# Temperature of saturated air - order of inputs doesn't matter
julia> T = HAPropsSI("T", "H", h, "R", 1.0, "P", 101325)
290.9620891952412
```

# Ref
HumidAir::HAPropsSI(const char* OutputName, const char* Input1Name, double Input1, const char* Input2Name, double Input2, const char* Input3Name, double Input3);
"""
function HAPropsSI(output::AbstractString, name1::AbstractString, value1::Real, name2::AbstractString, value2::Real, name3::AbstractString, value3::Real)
    val = ccall( (:HAPropsSI, "CoolProp"), Cdouble, (Cstring, Cstring, Cdouble, Cstring, Cdouble, Cstring, Cdouble), output, name1, value1, name2, value2, name3, value3)
    if val == Inf
        error("CoolProp: ", get_global_param_string("errstring"))
    end
    return val
end

"""
    cair_sat(t::Real)

Humid air saturation specific in [kJ/kg-K] heat at 1 atmosphere, based on a correlation from EES.

# Arguments
* `t`: T [K] good from 250K to 300K, no error bound checking is carried out.

# Ref
HumidAir::cair_sat(double);

# Note
Equals partial derivative of enthalpy with respect to temperature at constant relative humidity of 100 percent and pressure of 1 atmosphere.
"""
function cair_sat(t::Real)
    val = ccall( (:cair_sat, "CoolProp"), Cdouble, (Cdouble, ), t)
    return val;
end

ext: jl


text:
module SatelliteToolboxTLE

using Crayons, Dates, Parameters, Printf

import Base: show

################################################################################
#                                  Constants
################################################################################

# Escape sequences related to the crayons.
const _d = Crayon(reset = true)
const _b = crayon"bold"
const _g = crayon"bold green"
const _u = crayon"bold blue"
const _y = crayon"bold yellow"

################################################################################
#                                    Types
################################################################################

include("types.jl")

################################################################################
#                                   Includes
################################################################################

include("main.jl")
include("private.jl")

end


text:
@testset "qm7b" begin
    raw = rawdata(QM7b())
    @test typeof(raw[:names]) == Vector{String}
    @test size(raw[:names]) == (14,)
    @test typeof(raw[:X]) == Array{Float32,3}
    @test size(raw[:X]) == (7211, 23, 23)
    @test typeof(raw[:T]) == Matrix{Float32}
    @test size(raw[:T]) == (7211, 14)
end


text:
module PolyJuMP

# MultivariatePolynomials extension

using MultivariatePolynomials
using MultivariateMoments
using SemialgebraicSets

include("basis.jl")

# MOI extension

using MathOptInterface
const MOI = MathOptInterface

include("attributes.jl")
include("zero_polynomial.jl")
include("functions.jl")

# Bridges
const MOIB = MOI.Bridges
include("zero_polynomial_bridge.jl")
include("zero_polynomial_in_algebraic_set_bridge.jl")
include("plus_minus_bridge.jl")

# JuMP extension

using JuMP

include("variable.jl")
include("constraint.jl")

include("data.jl")
include("default.jl")

end # module


text:
module MKS
    SpeedOfLight = 2.99792458e8;                            # m / s
    GravitationalConstant = 6.673e-11;                      # m^3 / kg s^2
    PlancksConstantH = 6.62606896e-34;                      # kg m^2 / s
    PlancksConstantHbar = 1.05457162825e-34;                # kg m^2 / s
    AstronomicalUnit = 1.49597870691e11;                    # m
    LightYear = 9.46053620707e15;                           # m
    Parsec = 3.08567758135e16;                              # m
    GravAccel = 9.80665e0;                                  # m / s^2
    ElectronVolt = 1.602176487e-19;                         # kg m^2 / s^2
    Rydberg = 2.17987196968e-18;                            # kg m^2 / s^2
    Boltzmann = 1.3806504e-23;                              # kg m^2 / K s^2
    MolarGas = 8.314472e0;                                  # kg m^2 / K ol s^2
    StandardGasVolume = 2.2710981e-2;                       # m^3 / mol
    Minute = 6e1;                                           # s
    Hour = 3.6e3;                                           # s
    Day = 8.64e4;                                           # s
    Week = 6.048e5;                                         # s
    Inch = 2.54e-2;                                         # m
    Foot = 3.048e-1;                                        # m
    Yard = 9.144e-1;                                        # m
    Mile = 1.609344e3;                                      # m
    NauticalMile = 1.852e3;                                 # m
    Fathom = 1.8288e0;                                      # m
    Mil = 2.54e-5;                                          # m
    Point = 3.52777777778e-4;                               # m
    Texpoint = 3.51459803515e-4;                            # m
    Micron = 1e-6;                                          # m
    Angstrom = 1e-10;                                       # m
    Hectare = 1e4;                                          # m^2
    Acre = 4.04685642241e3;                                 # m^2
    Barn = 1e-28;                                           # m^2
    Liter = 1e-3;                                           # m^3
    UsGallon = 3.78541178402e-3;                            # m^3
    Quart = 9.46352946004e-4;                               # m^3
    Pint = 4.73176473002e-4;                                # m^3
    Cup = 2.36588236501e-4;                                 # m^3
    FluidOunce = 2.95735295626e-5;                          # m^3
    Tablespoon = 1.47867647813e-5;                          # m^3
    Teaspoon = 4.92892159375e-6;                            # m^3
    CanadianGallon = 4.54609e-3;                            # m^3
    UkGallon = 4.546092e-3;                                 # m^3
    MilesPerHour = 4.4704e-1;                               # m / s
    KilometersPerHour = 2.77777777778e-1;                   # m / s
    Knot = 5.14444444444e-1;                                # m / s
    PoundMass = 4.5359237e-1;                               # kg
    OunceMass = 2.8349523125e-2;                            # kg
    Ton = 9.0718474e2;                                      # kg
    MetricTon = 1e3;                                        # kg
    UkTon = 1.0160469088e3;                                 # kg
    TroyOunce = 3.1103475e-2;                               # kg
    Carat = 2e-4;                                           # kg
    UnifiedAtomicMass = 1.660538782e-27;                    # kg
    GramForce = 9.80665e-3;                                 # kg m / s^2
    PoundForce = 4.44822161526e0;                           # kg m / s^2
    KilopoundForce = 4.44822161526e3;                       # kg m / s^2
    Poundal = 1.38255e-1;                                   # kg m / s^2
    Calorie = 4.1868e0;                                     # kg m^2 / s^2
    Btu = 1.05505585262e3;                                  # kg m^2 / s^2
    Therm = 1.05506e8;                                      # kg m^2 / s^2
    Horsepower = 7.457e2;                                   # kg m^2 / s^3
    Bar = 1e5;                                              # kg / m s^2
    StdAtmosphere = 1.01325e5;                              # kg / m s^2
    Torr = 1.33322368421e2;                                 # kg / m s^2
    MeterOfMercury = 1.33322368421e5;                       # kg / m s^2
    InchOfMercury = 3.38638815789e3;                        # kg / m s^2
    InchOfWater = 2.490889e2;                               # kg / m s^2
    Psi = 6.89475729317e3;                                  # kg / m s^2
    Poise = 1e-1;                                           # kg m^-1 s^-1
    Stokes = 1e-4;                                          # m^2 / s
    Stilb = 1e4;                                            # cd / m^2
    Lumen = 1e0;                                            # cd sr
    Lux = 1e0;                                              # cd sr / m^2
    Phot = 1e4;                                             # cd sr / m^2
    Footcandle = 1.076e1;                                   # cd sr / m^2
    Lambert = 1e4;                                          # cd sr / m^2
    Footlambert = 1.07639104e1;                             # cd sr / m^2
    Curie = 3.7e10;                                         # 1 / s
    Roentgen = 2.58e-4;                                     # A s / kg
    Rad = 1e-2;                                             # m^2 / s^2
    SolarMass = 1.98892e30;                                 # kg
    BohrRadius = 5.291772083e-11;                           # m
    Newton = 1e0;                                           # kg m / s^2
    Dyne = 1e-5;                                            # kg m / s^2
    Joule = 1e0;                                            # kg m^2 / s^2
    Erg = 1e-7;                                             # kg m^2 / s^2
    StefanBoltzmannConstant = 5.67040047374e-8;             # kg / K^4 s^3
    ThomsonCrossSection = 6.65245893699e-29;                # m^2
    BohrMagneton = 9.27400899e-24;                          # A m^2
    NuclearMagneton = 5.05078317e-27;                       # A m^2
    ElectronMagneticMoment = 9.28476362e-24;                # A m^2
    ProtonMagneticMoment = 1.410606633e-26;                 # A m^2
    Faraday = 9.64853429775e4;                              # A s / mol
    ElectronCharge = 1.602176487e-19;                       # A s
    VacuumPermittivity = 8.854187817e-12;                   # A^2 s^4 / kg m^3
    VacuumPermeability = 1.25663706144e-6;                  # kg m / a^2 s^2
    Debye = 3.33564095198e-30;                              # A s^2 / m^2
    Gauss = 1e-4;                                           # kg / A s^2
end


text:
abstract type AbstractMonitor end

include("./Monitors/FrequencyDomainMonitor.jl")
include("./Monitors/TimeDomainMonitor.jl")
include("./Monitors/ModeExpansionMonitor.jl")


text:
function BFsub(scen_data, gap, delayMode, filename_Node, filename_Generator, filename_Line, filename_Storage, filename_Station, filename_SMP)

    m = Model(solver = GurobiSolver(MIPGap=gap))
    buses, generators, lines, storages, stations, SMP = DataImport(filename_Node, filename_Generator, filename_Line, filename_Storage, filename_Station, filename_SMP)

    ### Code begins here ###
    N_bus = length(buses)
    N_line = length(lines)
    # N_gen = length(generators)
    N_storage = length(storages)

	B_G = Int64[]
	N_gen = 0
	for b in 1:N_bus
		if generators[b].Pgmax > 0
			N_gen += 1
			push!(B_G,b)
		end
	end
	N_microgrid = N_gen + N_storage

    N_station = Array{Int64}(N_bus,1)
    N_station[1:end] = 0
    for b = 1:N_bus
        for h = 1:length(stations)
            if stations[h].busidx == b
                N_station[b] += 1
            end
        end
    end

    N_time = 24
    N_scen = size(scen_data.contin,3)
    S_N = Int64[]
    S_E = Int64[]
    for s = 1:N_scen
        if sum(scen_data.contin,1:2)[s] == 0
            push!(S_N,s)
        elseif sum(scen_data.contin,1:2)[s] > 0
            push!(S_E,s)
        end
    end

    ### 24 hour load profile ###
    # load24 = readcsv("load24.csv")[2,:]
    load24 = Array{Float64}(N_scen,N_time)
    for s in S_N
        load24[s,:] = scen_data.load24[s,:]
    end
    for s in S_E
        load24[s,:] = mean(scen_data.load24,1)
    end

    delay = DelayMatrix(filename_Line)
    t_delay = Array{Int64}(N_bus,N_bus,N_time)
    t_delay_ini = Array{Int64}(N_bus,N_bus)
    for t = 1:N_time
        t_delay[:,:,t] = delay
    end
    t_delay_ini = delay
    #t_delay[1:end] = 2
    #t_delay_ini[1:end] = 2

    b = 1:N_bus
    l = 1:N_line
	i = 1:N_gen
	mg = 1:N_microgrid
    k = 1:N_storage
    h = 1:length(stations)
    t = 1:N_time
    s = 1:N_scen

    PF_ES = tan(acos(0.9))
    PF_LS = tan(acos(0.9))

    Prob = Array{Float64}(N_scen,1) # Scenario Probability
    Prob[S_N] = 0.9 / length(S_N)
    Prob[S_E] = 0.1 / length(S_E)

    LScost = Array{Float64}(1,N_bus)
    LScost[1:end] = 5000
    LScost[13:15] = 10000

    @variable(m, v[b,t,s] >= 0) # variable for voltage square, voltage^2
    @variable(m, a[b,t,s] >= 0) # variable for current

    @variable(m, fp[b,t,s])
    @variable(m, fq[b,t,s])

    @variable(m, pg[b,t,s] >= 0)
    @variable(m, qg[b,t,s])

    @variable(m, e[k,t,s] >= 0)
    @variable(m, SoC[k,t,s] >= 0)
    @variable(m, pch[k,b,t,s] >= 0) # P^ch_(k,b,t): Charging Power of Energy Storage
    @variable(m, pdisch[k,b,t,s] >= 0) # P^disch_(k,b,t): Discharging Power of Energy Storage
    @variable(m, qch[k,b,t,s] >= 0) # P^ch_(k,b,t): Charging Power of Energy Storage
    @variable(m, qdisch[k,b,t,s] >= 0) # P^disch_(k,b,t): Discharging Power of Energy Storage

    @variable(m, u[k,b,t,s], Bin)
    @variable(m, z5[k,t,s], Bin) # ||Journal|| z -> z5 # prevent simantaneous charging
    @variable(m, x[k], Bin)

    @variable(m, z[k,b], Bin) # ||Journal|| Incidence matrix variable of ES unit in stationary mode: 1 if ES unit k is located at bus b, 0 otherwise.

    @variable(m, delta[mg,b,t,s], Bin) # ||Journal||
    @variable(m, mu[mg,l,t,s], Bin) # ||Journal||
    @variable(m, sigma_l[l,t,s], Bin) # ||Journal||
    @variable(m, sigma_d[b,t,s], Bin) # ||Journal||

    @variable(m, IC >= 0)
    @variable(m, OC[s] >= 0)
    @variable(m, EC[s] >= 0) # ||Journal||
    @variable(m, OCgen[s] >= 0)
    @variable(m, OCstorage[s] >= 0)
    @variable(m, OCloadshed[s] >= 0)
	@variable(m, OCinter[s] >= 0)
    @variable(m, PreventMultiSol >= 0)

    #@objective(m, Min, 0)
    @objective(m, Min, IC + sum( Prob[s] * OC[s] for s in S_N) + sum( Prob[s] * EC[s] for s in S_E) + PreventMultiSol)
    #@objective(m, Min, IC + sum( Prob[s] * OC[s] for s in 1:N_scen))
    @constraint(m, PreventMultiSol == 1e-12*sum(sum(sum(sum(b*u[k,b,t,s] for t in 1:N_time) for k in 1:N_storage) for b in 1:N_bus )for s in 1:N_scen))

    DCRF=3.5481*1e-4 # h=10years
    #@constraint(m, IC == sum(DCRF*(storages[k].Cp * storages[k].Pkmax + storages[k].Ce * storages[k].Ekmax) for k=1:N_storage) )
    @constraint(m, IC == sum(DCRF*(storages[k].Cp * storages[k].Pkmax + storages[k].Ce * storages[k].Ekmax) * x[k] for k=1:N_storage) )
    #@constraint(m, IC == sum(storages[k].Cp * storages[k].Pkmax + storages[k].Ce * storages[k].Ekmax for k=1:N_storage) + sum(stations[h].Cf for h=1:length(stations)))
    @constraint(m, OperatingCost[s in S_N], OC[s] == OCgen[s] + OCstorage[s] + OCinter[s])
    @constraint(m, EmergencyCost[s in S_E], EC[s] == OCgen[s] + OCstorage[s] + OCinter[s] + OCloadshed[s] )
    @constraint(m, InterfaceLineCost[s in 1:N_scen], OCinter[s] == sum( SMP[t] * pg[1,t,s] for t = 1:N_time ))
    @constraint(m, GenerationOC[s in 1:N_scen], OCgen[s] == sum(sum(generators[b].cost * pg[b,t,s] for t =1:N_time) for b=2:N_bus))
    @constraint(m, StorageOC[s in 1:N_scen], OCstorage[s] == sum(sum(sum(storages[k].Co * (pch[k,b,t,s] + pdisch[k,b,t,s]) for k=1:N_storage) for b=1:N_bus) for t=1:N_time))
    @constraint(m, LoadShedOC[s in S_E], OCloadshed[s] == sum(sum(LScost[b] * (1-sigma_d[b,t,s]) * buses[b].Pd * load24[s,t] for t =1:N_time) for b=1:N_bus))

    # Energy Storage Initial & Final SoC constraints
    for k = 1:N_storage
        for s in S_N
            @constraint(m, e[k,1,s] == 0.5*storages[k].Ekmax + sum(pch[k,b,1,s] * u[k,b,1,s] * storages[k].Aleph for b=1:N_bus) - sum(pdisch[k,b,1,s] * u[k,b,1,s]/ storages[k].Aleph for b=1:N_bus))
            @constraint(m, e[k,24,s] >= 0.5*storages[k].Ekmax)
        end
        for s in S_E
            @constraint(m, e[k,1,s] == 0.5*storages[k].Ekmax + sum(pch[k,b,1,s] * u[k,b,1,s] * storages[k].Aleph for b=1:N_bus) - sum(pdisch[k,b,1,s] * u[k,b,1,s]/ storages[k].Aleph for b=1:N_bus))
            @constraint(m, e[k,24,s] >= 0*storages[k].Ekmax)
        end
    end

    # Energy Storage Dynamic Equation
    for s in 1:N_scen
        for k = 1:N_storage
            for t = 2:N_time
                @constraint(m, e[k,t,s] == e[k,t-1,s]+sum(pch[k,b,t,s] * u[k,b,t,s] * storages[k].Aleph for b=1:N_bus) - sum(pdisch[k,b,t,s] * u[k,b,t,s]/ storages[k].Aleph for b=1:N_bus))
            end
            for t = 1:N_time
                @constraint(m, SoC[k,t,s] == e[k,t,s] / storages[k].Ekmax)
                @constraint(m, e[k,t,s] <= storages[k].Ekmax)
                for b = 1:N_bus
                    @constraint(m, pch[k,b,t,s]*storages[k].Aleph <= storages[k].Pkmax*z5[k,t,s])
                    @constraint(m, pdisch[k,b,t,s]/storages[k].Aleph <= storages[k].Pkmax*(1-z5[k,t,s]))
                end
            end
        end

        for k = 1:N_storage
            for t = 1:N_time
                for b = 1:N_bus
                    @constraint(m, qch[k,b,t,s] >= - pch[k,b,t,s] * tan(acos(PF_ES)))
                    @constraint(m, qch[k,b,t,s] <= + pch[k,b,t,s] * tan(acos(PF_ES)))
                    @constraint(m, qdisch[k,b,t,s] >= - pdisch[k,b,t,s] * tan(acos(PF_ES)))
                    @constraint(m, qdisch[k,b,t,s] <= + pdisch[k,b,t,s] * tan(acos(PF_ES)))
                end
            end
        end

        for l = 1:N_line
            for t = 1:N_time
                @constraint(m, (fp[lines[l].head,t,s])^2 + (fq[lines[l].head,t,s])^2 <= ( lines[lines[l].arcID].u)^2 )
                @constraint(m, (fp[lines[l].head,t,s] - a[lines[l].head,t,s]*lines[l].r)^2 + (fq[lines[l].head,t,s] - a[lines[l].head,t,s]*lines[l].x)^2 <= ( lines[lines[l].arcID].u)^2 )
                if scen_data.contin[l,t,s] == 0
                    @constraint(m, v[lines[l].head,t,s] - 2*(lines[l].r*fp[lines[l].head,t,s] + lines[l].x*fq[lines[l].head,t,s]) + a[lines[l].head,t,s]*(lines[l].r^2 + lines[l].x^2) == v[lines[l].tail,t,s])
                end
                @constraint(m, ((fp[lines[l].head,t,s])^2 + (fq[lines[l].head,t,s])^2) <= v[lines[l].head,t,s] * a[lines[l].head,t,s] ) # SOCP CONSTRAINT
            end
        end

        for t = 1:N_time
            @constraint(m, v[buses[1].root,t,s] == 1) #Voltage constraint for root node
            @constraint(m, a[buses[1].root,t,s] == 0)
            @constraint(m, fp[buses[1].root,t,s] == 0) # one-way flow: only from Transmission to Distribution
            @constraint(m, fq[buses[1].root,t,s] == 0)
        end
    end

    for s in S_N
        for b = 1:N_bus
            for t = 1:N_time
                @constraint(m, fp[b,t,s]  - sum(fp[l,t,s] - buses[l].R*a[l,t,s] for l in buses[b].children) - pg[b,t,s] - sum(pdisch[k,b,t,s] * u[k,b,t,s] for k in 1:N_storage) + sum(pch[k,b,t,s] * u[k,b,t,s] for k in 1:N_storage) + buses[b].Pd * load24[s,t] == 0)
                @constraint(m, fq[b,t,s]  - sum(fq[l,t,s] - buses[l].X*a[l,t,s] for l in buses[b].children) - qg[b,t,s] - sum(qdisch[k,b,t,s] * u[k,b,t,s] for k in 1:N_storage) + sum(qch[k,b,t,s] * u[k,b,t,s] for k in 1:N_storage) + buses[b].Qd * load24[s,t] - v[b,t,s]*buses[b].B == 0)
            end
        end
    end
    for s in 1:N_scen
        for b = 1:N_bus
            for t = 1:N_time
                if s in S_N
                    @constraint(m, v[b,t,s] <= buses[b].Vmax) #upper limit constraint for voltage square
                    @constraint(m, v[b,t,s] >= buses[b].Vmin) #lower limit constraint for voltage square
                end
            end
        end

        for b = 1:N_bus
            for t = 1:N_time
                @constraint(m, pg[b,t,s] <= generators[b].Pgmax)
                @constraint(m, qg[b,t,s] >= -generators[b].Qgmax)
                @constraint(m, qg[b,t,s] <= generators[b].Qgmax)
            end
        end

        for k = 1:N_storage
            #@constraint(m, sum(u[k,b,1,s] for b in 1:N_bus) == 1)
            #@constraint(m, sum(u[k,b,1,s] for b in 1:N_bus) == x[k])
            for t = 1:N_time
                #@constraint(m, sum(u[k,b,t,s] for b in 1:N_bus) <= 1)
                @constraint(m, sum(u[k,b,t,s] for b in 1:N_bus) <= x[k])
            end
        end

        #**** index generation
        temp = collect(1:N_bus)
        for h = 1:length(stations)
            temp_idx = find(temp .== stations[h].busidx)
            if (!isempty(temp_idx))
                splice!(temp,temp_idx[1])
            end
        end
        B_not_ch = temp
        temp = collect(1:N_bus)
        B_ch = setdiff(temp,B_not_ch)

        # equation (19)
        for t in 1:N_time
            for b in B_ch
                @constraint(m, sum( u[k,b,t,s] for k in 1:N_storage) <= N_station[b] )
            end
        end
        # charging station installed location only
        for k in 1:N_storage
            for b in B_not_ch
                for t in 1:N_time
                    @constraint(m, u[k,b,t,s] == 0)
                end
            end
        end
    end

    for k in 1:N_storage
        for b in 1:N_bus
            for t in 1:N_time
                for s in S_N
                    @constraint(m, u[k,b,t,s] == z[k,b] )
                end
            end
        end
    end
    for k in 1:N_storage
        @constraint(m, sum(z[k,b] for b in 1:N_bus) == x[k])
    end

    ##  Linkage z[k,b] between Normal Operation and Emergency Operations
    for s in S_E
        for k in 1:N_storage
            for i in 1:N_bus
                for j in 1:N_bus
                    if (t_delay_ini[i,j]!=0)
                        for tt in 1:t_delay_ini[i,j]
                            if (i!=j) & (tt <= N_time)
                                @constraint(m, z[k,i] <= 1 - u[k,j,tt,s])
                            end
                            if (i!=j) & (1+tt <= N_time)
                                @constraint(m, z[k,i] - u[k,i,1,s] <= 1 - u[k,j,1+tt,s])
                            end
                        end
                    end
                end
            end
        end
    end

    # Delay Model
    for s in S_E
        for k in 1:N_storage
            for i in 1:N_bus
                for j in 1:N_bus
                    for t in 1:N_time
                        if (t_delay[i,j,t]!=0)
                            for tt in 1:t_delay[i,j,t]
                                if (i!=j) & (t+tt <= N_time) & (delayMode==1)
                                    @constraint(m, u[k,i,t,s] - u[k,i,t+1,s] <= 1 - u[k,j,t+tt,s])
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    # Microgrid operation mode
    for s in S_E
        for l = 1:N_line
            for t = 1:N_time
                @constraint(m, (fp[lines[l].head,t,s])^2 + (fq[lines[l].head,t,s])^2 <= ( sigma_l[l,t,s] * lines[lines[l].arcID].u)^2 )
                @constraint(m, (fp[lines[l].head,t,s] - a[lines[l].head,t,s]*lines[l].r)^2 + (fq[lines[l].head,t,s] - a[lines[l].head,t,s]*lines[l].x)^2 <= ( sigma_l[l,t,s] * lines[lines[l].arcID].u)^2 )
            end
        end
        for b = 1:N_bus
            for t = 1:N_time
                @constraint(m, fp[b,t,s]  - sum(fp[l,t,s] - buses[l].R*a[l,t,s] for l in buses[b].children) - pg[b,t,s] - sum(pdisch[k,b,t,s] * u[k,b,t,s] for k in 1:N_storage) + sum(pch[k,b,t,s] * u[k,b,t,s] for k in 1:N_storage) + sigma_d[b,t,s] * buses[b].Pd * load24[t] == 0)
                @constraint(m, fq[b,t,s]  - sum(fq[l,t,s] - buses[l].X*a[l,t,s] for l in buses[b].children) - qg[b,t,s] - sum(qdisch[k,b,t,s] * u[k,b,t,s] for k in 1:N_storage) + sum(qch[k,b,t,s] * u[k,b,t,s] for k in 1:N_storage) + sigma_d[b,t,s] * buses[b].Qd * load24[t] - v[b,t,s]*buses[b].B == 0)
            end
        end
    end

    ###### left just in case of PH algorithm ######
    # linkage among scenarios
    for k in 1:N_storage
        for b in 1:N_bus
            for s in 2:N_scen
            #    @constraint(m, u[k,b,1,s] == u[k,b,1,1])
            end
        end
    end
    # equation () | Stationary ESS location for normal operation
    for k in 1:N_storage
        for b in 1:N_bus
            for t in 2:N_time
            #    @constraint(m, u[k,b,t,1] == u[k,b,1,1])
            end
        end
    end
    ##############################################
    for t in 1:N_time
        for s in S_N
            for l in 1:N_line
                @constraint(m, sigma_l[l,t,s] == 1)
            end
            for b in 1:N_bus
                @constraint(m, sigma_d[b,t,s] == 1)
            end
        end
    end

    # Microgrid formation
    for t in 1:N_time
        for s in S_E
			temp_contin = scen_data.contin[:,t,s]
            N_cont = sum(temp_contin.!=0)
            @constraint(m, sum(sigma_l[l,t,s] for l in 1:N_line) >= N_line - max(N_cont,N_storage+N_gen))
            for l in 1:N_line
                if scen_data.contin[l,t,s] != 0
                    @constraint(m, sigma_l[l,t,s] == 0)
                    @constraint(m, fp[lines[l].head,t,s] == 0)
                    @constraint(m, fq[lines[l].head,t,s] == 0)
                    @constraint(m, a[lines[l].head,t,s] == 0)
                end
            end
        end
    end

    tic()
    status = solve(m)
    simulation_time = toc()
    obj_value = getobjectivevalue(m)
    result_v = getvalue(v)[:,:,:] # [b,t,s]
    result_a = getvalue(a)[:,:,:] # [b,t,s]
    result_fp = getvalue(fp)[:,:,:] # [b,t,s]
    result_fq = getvalue(fq)[:,:,:] # [b,t,s]
    result_pg = getvalue(pg)[:,:,:] # [i,t,s]
    result_qg = getvalue(qg)[:,:,:] # [i,t,s]
    result_e = getvalue(e)[:,:,:] # [k,t,s]
    result_SoC = getvalue(SoC)[:,:,:] # [k,t,s]
    result_pch = getvalue(pch)[:,:,:,:] # [k,b,t,s]
    result_pdisch = getvalue(pdisch)[:,:,:,:] # [k,b,t,s]
    result_u = getvalue(u)[:,:,:,:] # [k,b,t,s]

    result_z = getvalue(z)[:,:] # [k,b]
    result_delta = getvalue(delta)[:,:,:,:] # [k,b,t,s]
    result_mu = getvalue(mu)[:,:,:,:] # [k,l,t,s]
    result_sigma_l = getvalue(sigma_l)[:,:,:] # [l,t,s]
    result_sigma_d = getvalue(sigma_d)[:,:,:] # [b,t,s]
    result_EC = getvalue(EC)[:]

    result_IC = getvalue(IC)
    result_OC = getvalue(OC)[:]

    result_OCgen = getvalue(OCgen)[:]
    result_OCloadshed = getvalue(OCloadshed)[:]
    result_OCinter = getvalue(OCinter)[:]
    result_OCstorage = getvalue(OCstorage)[:]
    result_x = getvalue(x)[:] # [k]
    result_ESlocation = Array{Int64}(N_storage, N_time, N_scen)
    for s = 1:N_scen
        for k = 1:N_storage
            temp = result_u[k,:,:,s]
            temp[find(temp.<0.5)]=0
            temp[find(temp.>0.5)]=1
            temp_idx = ind2sub(temp,find(temp.==1))
            temp_location = Array{Int64}(N_time,1)
            temp_location[1:end] = -1
            for i = 1:length(temp_idx[1])
                temp_location[temp_idx[2][i]]  = temp_idx[1][i] - 1
            end
            result_ESlocation[k,:,s] = temp_location
        end
    end

    u = Array{Int64}(N_storage, N_bus, N_time, N_scen)
    u[:] = 0
    for s = 1:N_scen
        for i = 1:N_storage
            for j = 1:N_bus
                for k = 1:N_time
                    if (result_u[i,j,k,s] == 1)
                        u[i,j,k,s] = 1
                    elseif (result_u[i,j,k,s] == 0)
                        u[i,j,k,s] = 0
                    end
                end
            end
        end
    end
    z = Array{Int64}(N_storage, N_bus)
    z[:] = 0
    for k in 1:N_storage
        for b in 1:N_bus
            if (result_z[k,b] == 1)
                z[k,b] = 1
            elseif (result_z[k,b] == 0)
                z[k,b] = 0
            end
        end
    end

    return status, obj_value, simulation_time, result_ESlocation, result_SoC, u, z, result_v, result_a, result_pg, result_qg, result_IC, result_OC, result_EC, result_OCinter, result_OCgen, result_OCloadshed, result_OCstorage, result_x, result_sigma_d, result_sigma_l
end


text:
"""
Delegates loading of code for different interpolation aglorithms (in their own modules).
Also wraps these modules, so they're called easily by `pack_chunk_list_timeseries_to_matrix`.

Author: Eric Ford
Created: August 2020
"""

# First include modules for different interpolation algoritihms
# File called module will include whatever packages it depends on, and also
# provide at least one function with a common interface (currently named) interp_chunk_to_grid_ALGNAME!.

include("common.jl")
export numerical_deriv

include("linear/linear.jl")
import .LinearInterpolation
export interp_chunk_to_shifted_grid_linear, interp_chunk_to_shifted_grid_linear!
export interp_chunk_to_grid_linear, interp_chunk_to_grid_linear!

include("sinc/sinc.jl")
using .SincInterpolation
export interp_chunk_to_shifted_grid_sinc, interp_chunk_to_shifted_grid_sinc!
export interp_chunk_to_grid_sinc, interp_chunk_to_grid_sinc!

#=
include("gp_slow/gp_brute_force.jl")             # Basic, brute-force GPs with no dependancies, hasn't been updated recently
import .GPInterpolation
export interp_chunk_to_shifted_grid_gp_brute_force, interp_chunk_to_shifted_grid_gp_brute_force!
export interp_chunk_to_grid_gp_brute_force, interp_chunk_to_grid_gp_brute_force!
=#

include("gp_fast/temporalgps.jl")    # Dependenacies on Stheno, Temporal GPs and Static Arrays
import .TemporalGPInterpolation
export interp_chunk_to_shifted_grid_gp_temporal, interp_chunk_to_shifted_grid_gp_temporal!
export interp_chunk_to_grid_gp_temporal, interp_chunk_to_grid_gp_temporal!


text:
# Note that this test asumes TPF properly works

pm = PoolModel("ss1")
pm <= Setting(:data_vintage, "190822")
filepath = dirname(@__FILE__)
pm <= Setting(:dataroot, "$(filepath)/../reference/")
data = df_to_matrix(pm, load_data(pm))

# This commented code produces the saved output
tuning = Dict(:r_star => 2., :c_init => 0.3, :target_accept_rate => 0.4,
              :resampling_method => :systematic, :n_mh_steps => 1,
              :n_particles => 1000, :n_presample_periods => 0,
              :allout => true)
pm <= Setting(:tuning, tuning, "tuning parameters for TPF")
pm[:ρ].value = 0.8
pm[:μ].value = 0.
pm[:σ].value = 1.
Random.seed!(1793)
Φpost, Ψpost, F_ϵpost, F_upost, F_λpost = compute_system(pm)
s_init = reshape(rand(F_λpost, tuning[:n_particles]), 1, 1000)
s_init = [s_init; 1 .- s_init]
tpf_out, ~, ~ = tempered_particle_filter(data, Φpost, Ψpost, F_ϵpost, F_upost,
                                   s_init; tuning..., verbose = :none,
                                   fixed_sched = [1.], parallel = false, poolmodel = true)

@testset "Check likelihood and posterior calculations for dynamic weight" begin
    Random.seed!(1793)
    lh = likelihood(pm, data)
    @test lh == tpf_out

    Random.seed!(1793)
    x = map(α->α.value, pm.parameters)
    post_at_start = DSGE.posterior!(pm, x, data)
    @test post_at_start == tpf_out

    Random.seed!(1793)
    global y = x .+ [-.7; -300.; 5.]
    post_not_at_start = DSGE.posterior!(pm, y, data)
    ϵ = 0.0004
    @test abs(post_at_start - post_not_at_start) > ϵ
end

pm = PoolModel("ss1", weight_type = :static)
Random.seed!(1793)
static_out = sum(DSGE.filter_likelihood(pm, data))

@testset "Check likelihood and posterior calculations for static weight" begin
    Random.seed!(1793)
    lh = likelihood(pm, data)
    @test lh == static_out

    Random.seed!(1793)
    x = map(α->α.value, pm.parameters)
    post_at_start = DSGE.posterior!(pm, x, data)
    @test post_at_start == static_out

    Random.seed!(1793)
    global z = x .+ (1 .- x) ./ 2
    post_not_at_start = DSGE.posterior!(pm, z, data)
    ϵ = 1.
    @test abs(post_at_start - post_not_at_start) > ϵ
end

nothing


text:
export rule

@rule typeof(*)(:out, Marginalisation) (m_A::PointMass, m_in::PointMass) = PointMass(mean(m_A) * mean(m_in))

@rule typeof(*)(:out, Marginalisation) (m_A::PointMass, m_in::F) where { F <: NormalDistributionsFamily } = begin
    A = mean(m_A)
    return convert(promote_variate_type(F, NormalMeanVariance), A * mean(m_in), A * cov(m_in) * A')
end


text:
relu(x) = max(x, 0)

function propagate!(n::AbstractArray, k::AbstractArray, model::AbstractBinomialModel, dt)
    # propagate the hidden states n, k by one timestep
    m_out, m_in = size(n)

    @inbounds for i in 1:m_in
        @inbounds for j in 1:m_in
            # refill counts and probabilities
            n[i, j]  = relu(model.N[i] - n[i, j] + k[i, j])
            p_refill = 1 - exp(-dt / model.τ[i])

            # draw vesicles to be refilled, update n
            k[i, j] = rand(Binomial(n[i, j], p_refill))
            n[i, j] = model.N[i] - n[i, j] + k[i, j]

            # draw vesicles to be released
            k[i, j] = rand(Binomial(n[i, j], model.p[i]))
        end
    end
    return n, k
end

function propagate!(n::AnyCuArray, k::AnyCuArray, model::AbstractBinomialModel, dt)
    # propagate the hidden states n, k by one timestep

    # refill counts and probabilities
    @. n     = relu(model.N - n + k)
    p_refill = @. 1 - exp(-dt / model.τ)

    # draw vesicles to be refilled, update n
    BinomialGPU.rand_binomial!(k, count = n, prob = p_refill)
    @. n     = model.N - n + k

    # draw vesicles to be released
    BinomialGPU.rand_binomial!(k, count = n, prob = model.p)
    return n, k
end

function propagate!(state::BinomialState, model, dt)
    propagate!(state.n, state.k, model, dt)
    return state
end


text:
function print_log(res_dir)
    tree = Core.ResDirTree(res_dir);
    println(Core.load_log(tree));
end

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
export
    nthfibonacci, nfibonacci,
    allfibonacci, somefibonacci, exactfibonacci

struct FibonacciIterator{T<:Integer}
    x1::T
    x2::T
end

nthfibonacci(n::Int, S::Type = Int) = @pipe allfibonacci(S) |> drop(_, n - one(n)) |> first
nfibonacci(n::Int, S::Type = Int) = collect(exactfibonacci(n, S))

allfibonacci(T::Type = Int) = allfibonacci(one(T), one(T))
allfibonacci(x1::Integer, x2::Integer) = FibonacciIterator(promote(x1, x2)...)

somefibonacci(xmax::Integer) = somefibonacci(xmax, one(xmax), one(xmax))
somefibonacci(xmax::Integer, x1::Integer, x2::Integer) = somefibonacci(promote(xmax, x1, x2)...)
somefibonacci(xmax::T, x1::T, x2::T) where {T<:Integer} = @pipe allfibonacci(x1, x2) |> takewhile(x -> x <= xmax, _)

exactfibonacci(n::Int, S::Type = Int) = exactfibonacci(n, one(S), one(S))
exactfibonacci(n::Int, x1::Integer, x2::Integer) = exactfibonacci(n, promote(x1, x2)...)
exactfibonacci(n::Int, x1::S, x2::S) where {S<:Integer} = @pipe allfibonacci(x1, x2) |> take(_, n)

Base.start(it::FibonacciIterator) = (checked_sub(it.x2, it.x1), it.x1)
Base.next(::FibonacciIterator, state) = (state[2], (state[2], checked_add(state[1], state[2])))
Base.done(::FibonacciIterator, _) = false

Base.eltype(it::FibonacciIterator) = Base.eltype(typeof(it))
Base.eltype(::Type{FibonacciIterator{T}}) where {T} = T

Base.iteratorsize(::FibonacciIterator) = Base.IsInfinite()


text:
using GeoPhyInv # precompile it

"""
* `files` are in folder
* `gfiles` are files in the generated folder
"""
function run_test(files, folder, gfiles, nfail)
    exename = joinpath(Sys.BINDIR, Base.julia_exename())
    testdir = @__DIR__
    for file in files
        fp = joinpath(testdir, folder, string(file, ".jl"))
        printstyled("Testing GeoPhyInv file:\t", fp, "\n"; bold = true, color = :white)
        # if f ∈ excludedfiles
        # println("Test Skip:")
        # println("$f")
        # continue
        # end
        try
            run(`$exename -O3 --startup-file=no --check-bounds=no $fp`)
        catch ex
            nfail += 1
        end
    end

    for file in gfiles
        fpg = joinpath(testdir, "generated", folder, string(file, ".jl"))
        printstyled("Testing GeoPhyInv file:\t", fpg, "\n"; bold = true, color = :white)
        try
            run(`$exename -O3 --startup-file=no --check-bounds=no $fpg`)
        catch ex
            nfail += 1
        end
    end
    printstyled("Number of failed tests:\t", nfail, "\n"; bold = true, color = :white)
end

nfail = 0
run_test(["convert", "param_adj"], "media", ["doc"], nfail)
run_test([], "ageom", ["doc"], nfail)
run_test([], "srcwav", ["doc"], nfail)
run_test([], "data", ["doc"], nfail)
# run_test(["accuracy2D", "rho_projection"],"fdtd", ["doc","create_snaps","reuse_expt"], nfail)
#run_test(["accuracy2D","backprop", "rho_projection"],"fdtd", ["doc","create_snaps","reuse_expt"], nfail)
#run_test(["gradient_accuracy", "born_map"],"fwi", ["doc", "pizza", "born_tutorial"], nfail)
# run_test(["interp_tests"], "Interpolation", [], nfail)
# run_test(
# ["testscript_RandomEigenfns", "adj_state_expt", "adj_state", "testdAdx"],
# "Poisson",
# ["doc", "forw", "test_born"],
# nfail,
# )

# exit(nfail)


text:
function compare_assert(a::Int64, b::Int64)
    @assert(a == b)
end

function main()
    compare_assert(1, 1)
end