!#########################################################################
!		
!    Copyright (C) 2003-2012 Department of Physics and Astronomy,
!                            University of Rochester,
!                            Rochester, NY
!
!    problem.f90 of module Template is part of AstroBEAR.
!
!    AstroBEAR is free software: you can redistribute it and/or modify	  
!    it under the terms of the GNU General Public License as published by 
!    the Free Software Foundation, either version 3 of the License, or    
!    (at your option) any later version.
!
!    AstroBEAR is distributed in the hope that it will be useful, 
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with AstroBEAR.  If not, see <http://www.gnu.org/licenses/>.
!
!#########################################################################
!> @dir Template
!! @brief Contains files necessary for the Template Calculation

!> @file problem.f90
!! @brief Main file for module Problem

!> @defgroup Template Template Module
!! @brief Module for calculating collapse of a uniform cloud
!! @ingroup Modules

!> Template Module 
!! @ingroup Template
MODULE Problem
  USE GlobalDeclarations
  USE DataDeclarations
  USE Clumps
  USE Ambients
  USE Winds
  USE Projections
  USE Spectras
  USE Totals
  USE Fields
  USE Cameras
  USE RiemannSolvers
  IMPLICIT NONE
  SAVE
  PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, &
         ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
  PRIVATE 
  REAL(KIND=qPREC) :: spacing, thickness, beta, mach, screen_x, rho_wire, rho_amb, rho_wind, amplitude, omega, offset(3), thickness2  
  TYPE(WindDef), POINTER :: Wind
  TYPE(ProjectionDef), POINTER :: Projection
  TYPE(TotalDef), POINTER :: Total
  TYPE(SpectraDef), POINTER :: SpectraRho, SpectraV, SpectraB
  INTEGER :: LBox=160, BoxN ! size of box for spectra domain
  LOGICAL :: lSpectra=.true., lHistogram=.false.
CONTAINS

  SUBROUTINE ProblemModuleInit()
    TYPE(AmbientDef), POINTER :: Ambient
    TYPE(ClumpDef), POINTER :: Clump
    
    NAMELIST/ProblemData/ spacing, thickness, beta, mach, screen_x, rho_wire, rho_amb, rho_wind, amplitude, omega, offset, thickness2
    OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
    READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
    CLOSE(PROBLEM_DATA_HANDLE)

    CALL CreateAmbient(Ambient)
    Ambient%density=rho_amb
    CALL CreateWind(Wind)
    wind%density=rho_wind
    Wind%velocity=mach*sqrt(gamma)
    Wind%temperature=1d0/rho_wind
    Wind%B=sqrt(2d0/beta)*(/0d0,1d0,0d0/)
    amplitude=amplitude*Wind%velocity
    omega=omega*2d0*Pi/spacing

    IF(lSpectra) THEN
       DO BoxN=1,10 
          CALL CreateSpectra(SpectraRho)
          ALLOCATE(SpectraRho%Fields(3))
          SpectraRho%Fields(:)%id=(/Mass_Field/)
          SpectraRho%type=SCALAR_SPECT
          SpectraRho%level=0

          LBox=levels(SpectraRho%level)%mX(2) !.0/levels(Spectra%level)%dx)
          SpectraRho%mB=1
          SpectraRho%mB(2:3,2)=levels(SpectraRho%level)%mX(2:3) 

          SpectraRho%mB(1,1)=1+(boxN-1)*LBox !levels(MaxLevel)%mX(1)/2 - (LBox/2)*2**MaxLevel+1
          SpectraRho%mB(1,2)=boxN*LBox !(levels(MaxLevel)%mX(2)/2 + (LBox/2)*2**MaxLevel

          CALL CreateSpectra(SpectraV)

          ALLOCATE(SpectraV%Fields(3))
          SpectraV%Fields(:)%id=(/vx_Field, vy_Field, vz_Field/)
          SpectraV%type=VECTOR_SPECT
          SpectraV%level=0

          LBox=levels(SpectraV%level)%mX(2) !.0/levels(Spectra%level)%dx)
          SpectraV%mB=1
          SpectraV%mB(2:3,2)=levels(SpectraV%level)%mX(2:3) 

          SpectraV%mB(1,1)=1+(boxN-1)*LBox !levels(MaxLevel)%mX(1)/2 - (LBox/2)*2**MaxLevel+1
          SpectraV%mB(1,2)=boxN*LBox !(levels(MaxLevel)%mX(2)/2 + (LBox/2)*2**MaxLevel

          IF (lMHD) THEN
             CALL CreateSpectra(SpectraB)
             ALLOCATE(SpectraB%Fields(3))
             SpectraB%Fields(:)%id=(/Bx_Field, By_Field, Bz_Field /)
             SpectraB%type=VECTOR_SPECT
             SpectraB%level=0

             LBox=levels(SpectraB%level)%mX(2) !.0/levels(Spectra%level)%dx)
             SpectraB%mB=1
             SpectraB%mB(2:3,2)=levels(SpectraB%level)%mX(2:3) 

             SpectraB%mB(1,1)=1+(boxN-1)*LBox !levels(MaxLevel)%mX(1)/2 - (LBox/2)*2**MaxLevel+1
             SpectraB%mB(1,2)=boxN*LBox !(levels(MaxLevel)%mX(2)/2 + (LBox/2)*2**MaxLevel
          END IF
       END DO
    END IF

    IF(lHistogram) THEN
      CALL CreateTotal(Total)
      Total%Field%id=Mass_Field
      Total%Field%component=GASCOMP
      Total%Field%name='density'
    END IF
  END SUBROUTINE ProblemModuleInit


  SUBROUTINE PostProcessInit()

    CALL CreateProjection(Projection)
    Projection%field%id=MASS_FIELD
    Projection%dim=1
    
    CALL CreateProjection(Projection)
    Projection%field%id=MASS_FIELD
    Projection%dim=2
    
    CALL CreateProjection(Projection)
    Projection%field%id=MASS_FIELD
    Projection%dim=3


  END SUBROUTINE PostProcessInit
  
  SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
    REAL(KIND=qPREC) :: x, y, z, dx, screen_fact, speed
    REAL(KIND=qPREC), DIMENSION(8) :: wl, wr, f
    INTEGER :: n
    INTEGER :: i,j,k

    CALL PlaceGrid(Info, screen_x, thickness, spacing, rho_wire,(/0d0,0d0,0d0/))
    CALL PlaceGrid(Info, screen_x, thickness2, spacing, rho_wire*(thickness/thickness2)**2,offset)

  END SUBROUTINE ProblemGridInit

    SUBROUTINE PlaceGrid(Info, xloc, t, s, rho,o)
    TYPE(InfoDef) :: Info
    REAL(KIND=qPREC) :: x, y, z, dx, screen_fact, xloc, t, s, rho, o(3)
    INTEGER :: i,j,k
      
    dx=levels(Info%level)%dx
    DO i=1,Info%mX(1)
       x=Info%xBounds(1,1)+(real(i,8)-.5d0)*dx + s*o(1)
       IF (x < xloc-t) CYCLE
       IF (x > xloc+t) CYCLE
       DO j=1,Info%mX(2)
          y=Info%xBounds(2,1)+(real(j,8)-.5d0)*dx
          y=modulo(y+s*o(2), s)-.5d0*s
          DO k=1,Info%mX(3)
             z=Info%xBounds(3,1)+(real(k,8)-.5d0)*dx
             z=modulo(z+s*o(3),s)-.5d0*s
             screen_fact=f((x-xloc)/t,y/t,z/t)
             IF (screen_fact > 0d0) THEN
                Info%q(i,j,k,1)=Info%q(i,j,k,1)*(1d0-screen_fact)+screen_fact*rho
                Info%q(i,j,k,imom(1:nDim))=0d0
                Info%q(i,j,k,iE)=1d0
                CALL Prim_To_Cons(Info%q(i,j,k,:))
             END IF
          END DO
       END DO
    END DO
  END SUBROUTINE PlaceGrid

  FUNCTION f(x,y,z)
    REAL(KIND=qPREC) :: x, y, z, f
    f=max(0d0,max(1d0-(x**2+y**2)**4, 1d0-(x**2+z**2)**4))
  END FUNCTION f

  SUBROUTINE ProblemBeforeStep(Info)
    TYPE(InfoDef) :: Info
    INTEGER :: i, rmbc
    REAL(KIND=qPREC) :: dx, t, x, vy
    rmbc=levels(Info%level)%gmbc(levels(Info%level)%step)
    dx=levels(Info%level)%dx
    t=levels(Info%level)%tnow
    DO i=1-rmbc, 1-Info%mGlobal(1,1)
       x=Info%xBounds(1,1)+(REAL(i)-.5d0)*dx
       vy=amplitude*sin(omega*(x-Wind%velocity*t))
       Info%q(i,:,:,ivy)=vy
       Info%q(i,:,:,iE)=Info%q(i,:,:,iE)+.5*Wind%Density*vy**2
    END DO
  !  print*, "Spectra%mB=", spectra%mB

  END SUBROUTINE ProblemBeforeStep

  SUBROUTINE ProblemAfterStep(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemAfterStep

  SUBROUTINE ProblemSetErrFlag(Info)
    TYPE(InfoDef) :: Info
    REAL(KIND=qPREC) :: x, y, z, dx, screen_fact
    INTEGER :: i,j,k

    CALL SetErrFlagGrid(Info, screen_x, thickness, spacing,(/0d0, 0d0,0d0/))
    CALL SetErrFlagGrid(Info, screen_x, thickness2, spacing, offset)
  END SUBROUTINE ProblemSetErrFlag


  SUBROUTINE SetErrFlagGrid(Info, xloc, t, s, o)
    TYPE(InfoDef) :: Info
    REAL(KIND=qPREC) :: x, y, z, dx, screen_fact, sloc, t, s, o(3), xloc
    INTEGER :: i,j,k

    dx=levels(Info%level)%dx
    DO i=1,Info%mX(1)
       x=Info%xBounds(1,1)+(real(i,8)-.5d0)*dx + s*o(1)
       IF (x < xloc-t-dx) CYCLE
       IF (x > xloc+t+dx) CYCLE
       DO j=1,Info%mX(2)
          y=Info%xBounds(2,1)+(real(j,8)-.5d0)*dx
          y=modulo(y+s*o(2), s)-.5d0*s
          DO k=1,Info%mX(3)
             z=Info%xBounds(3,1)+(real(k,8)-.5d0)*dx
             z=modulo(z+s*o(3),s)-.5*s
             screen_fact=f((x-xloc)/(t+dx),y/(t+dx),z/(t+dx))
             IF (screen_fact > 0d0) THEN
                Info%ErrFlag(i,j,k)=1
             END IF
          END DO
       END DO
    END DO
  END SUBROUTINE SetErrFlagGrid


  SUBROUTINE ProblemBeforeGlobalStep(n)
    INTEGER :: n
  END SUBROUTINE ProblemBeforeGlobalStep

END MODULE
