/*@@
   @file      GRHydro_RotorM.F90
   @date      Aug 15, 2011
   @author    Scott Noble, Joshua Faber, Bruno Mundim
   @desc 
   Cylindrical magnetized rotor test (see Etienne et al.).
   @enddesc 
 @@*/

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "GRHydro_Macros.h"

#define velx(i,j,k) vel(i,j,k,1)
#define vely(i,j,k) vel(i,j,k,2)
#define velz(i,j,k) vel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
#define Bvecx(i,j,k) Bvec(i,j,k,1)
#define Bvecy(i,j,k) Bvec(i,j,k,2)
#define Bvecz(i,j,k) Bvec(i,j,k,3)
#define Bconsx(i,j,k) Bcons(i,j,k,1)
#define Bconsy(i,j,k) Bcons(i,j,k,2)
#define Bconsz(i,j,k) Bcons(i,j,k,3)


 /*@@
   @routine    GRHydro_Rotor
   @date       Aug 11, 2011
   @author     Scott Noble, Joshua Faber, Bruno Mundim
   @desc 
   Initial data for magnetic rotor - parameters from Etienne et al. 2010
   @enddesc 
   @calls     
   @calledby   
   @history 
   Using GRHydro_ShockTubeM.F90 as a template.
   @endhistory 

@@*/

subroutine GRHydro_RotorM(CCTK_ARGUMENTS)

  implicit none
  
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS
  
  CCTK_INT :: i, j, k, nx, ny, nz
  CCTK_REAL :: radius, det, v_max
  CCTK_REAL :: rhoin, rhoout, pressin, pressout
  CCTK_REAL :: bvcxl,bvcyl,bvczl
  CCTK_REAL :: tmp,tmp2,gam,r_rot

  CCTK_INT :: use_smoothing
  CCTK_REAL :: rsmooth_rel,rfact

!!$Hardwired bfield and gamma values
  bvcxl = 1.0
  bvcyl = 0.0
  bvczl = 0.0
  
!!$Adiabatic index for test
  gam = (5.d0/3.d0)
  
!!$radius of rotor is 0.1
  r_rot = 0.1
  
!!$Maximum velocity is 0.995
  v_max = 0.995

!!$Inner values
  rhoin   = 10.d0
  pressin = 1.d0
  
!!$Outer values
  rhoout   = 1.d0
  pressout = 1.d0


!!$ Smooth the edge?  Define the radius in relative terms if so
  use_smoothing = 0
  rsmooth_rel = 0.05
  
  nx = cctk_lsh(1)
  ny = cctk_lsh(2)
  nz = cctk_lsh(3)
  
  if(v_max>1.0) then
     write(*,*)"No superluminal speeds!!!!  v_max=",v_max," reset to 0.995"
     v_max=0.995
  endif

  do i=1,nx
     do j=1,ny
        do k=1,nz
           
           radius = sqrt(x(i,j,k)**2+y(i,j,k)**2)
           
           if(radius.le.r_rot) then
              rho(i,j,k) = rhoin
              press(i,j,k) = pressin
              velx(i,j,k) = v_max/r_rot*y(i,j,k)
              vely(i,j,k) = -1.d0*v_max/r_rot*x(i,j,k)
              velz(i,j,k) = 0.d0
           else if((use_smoothing.eq.1) .and. &
                ((radius.gt.r_rot) .and. &
                (radius.le.((1.0+rsmooth_rel)*r_rot)))) then
              rfact = (radius/r_rot - 1.d0) / rsmooth_rel
              rho(i,j,k) = rfact*rhoout + (1.d0-rfact)*rhoin
              press(i,j,k) = rfact*pressout + (1.d0-rfact)*pressin
              velx(i,j,k) = (1.d0-rfact)*v_max * y(i,j,k) / radius
              velx(i,j,k) = -1.d0*(1.d0-rfact)*v_max * x(i,j,k) / radius
              velz(i,j,k) = 0.d0
           else
              rho(i,j,k) = rhoout
              press(i,j,k) = pressout
              velx(i,j,k) = 0.d0
              vely(i,j,k) = 0.d0
              velz(i,j,k) = 0.d0
           endif
           eps(i,j,k)=press(i,j,k)/(gam-1.d0)/rho(i,j,k)
           
           Bvecx(i,j,k)=bvcxl
           Bvecy(i,j,k)=bvcyl
           Bvecz(i,j,k)=bvczl
           
           det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
           
           if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then
              call Prim2ConPolyM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
                   gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
                   det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
                   tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),rho(i,j,k),&
                   velx(i,j,k),vely(i,j,k),velz(i,j,k),&
                   eps(i,j,k),press(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),&
                   w_lorentz(i,j,k))
           else
              call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
                   gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
                   det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
                   tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),rho(i,j,k),&
                   velx(i,j,k),vely(i,j,k),velz(i,j,k),&
                   eps(i,j,k),press(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),&
                   w_lorentz(i,j,k))
           end if
           
        enddo
     enddo
  enddo
  
  densrhs = 0.d0
  srhs = 0.d0
  taurhs = 0.d0
  Bconsrhs = 0.d0


  return
  
end subroutine GRHydro_RotorM