!#########################################################################
!		
!    Copyright (C) 2003-2012 Department of Physics and Astronomy,
!                            University of Rochester,
!                            Rochester, NY
!
!    problem.f90 of module Binary 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 OrbitingParticles
!! @brief Contains files necessary for the Orbiting Particles test problem

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

!> @defgroup OrbitingParticles Orbiting Particle Module
!! @brief Module for setting up orbiting particles
!! @ingroup Modules

!> Orbiting Particle Module
!! @ingroup OrbitingParticles
MODULE Problem
  USE DataDeclarations
  USE ParticleDeclarations
  USE Ambients
  USE Disks
  USE Winds
  IMPLICIT NONE
  SAVE
  PUBLIC ProblemModuleInit, ProblemGridInit, &
       ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
  REAL(KIND=qPREC) :: alpha, densw, velw, t1, t2
  INTEGER :: radiusw, nWinds
  LOGICAL :: windsPresent=.false., disk=.false.
  TYPE(WindDef), POINTER :: Wind
  INTEGER :: nParticles, edge
  REAL(KIND=qPREC) :: mass=0
  REAL(KIND=qPREC) :: xloc(3)
  REAL(KIND=qPREC) :: vel(3)
  TYPE(ParticleDef), POINTER :: Primary
  TYPE(ParticleDef), POINTER :: Secondary
  TYPE(ParticleDef), POINTER :: Particle

  TYPE(ParticleListDef), POINTER :: particlelist
  INTEGER :: grav_soft_rad
  INTEGER :: ids(2)
  TYPE(AmbientDef), POINTER :: Ambient
  TYPE(DiskDef), POINTER :: mydisk
  REAL(KIND=qPREC) :: time, rhoOut, pOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut, &
       buff(17) 

  !TYPE(pWindDef), DIMENSION(:), ALLOCATABLE :: MyWinds

CONTAINS

  !> Initializes module variables
   SUBROUTINE ProblemModuleInit()      
      TYPE(InfoDef) :: Info
      INTEGER :: i
      NAMELIST /AmbientData/ rhoOut, pOut, vxOut, vyOut, vzOut, BxOut, ByOut, BzOut
      NAMELIST /ProblemData/ nParticles, densw, velw, windsPresent, disk, t1, t2
      NAMELIST /ParticleData/ mass,xloc,vel,alpha,radiusw,grav_soft_rad, buff
      NAMELIST /RestartData/ ids

      time=levels(Info%level)%tnow

      OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
      READ(PROBLEM_DATA_HANDLE,NML=ProblemData)

if (.not. lrestart) then		
      CALL CreateAmbient(Ambient)
      READ(PROBLEM_DATA_HANDLE,NML=AmbientData)
      if (disk) then
		   Ambient%density=1d0
		else ; Ambient%density=rhoOut
		end if
      Ambient%pressure=pOut
      Ambient%B(:)=(/BxOut, ByOut, BzOut/)
      Ambient%velocity(:)=(/vxOut, vyOut, vzOut/)
end if      

if (windsPresent) then
        DO i=1,nDim
                DO edge=1,2
                        IF (Gmthbc(i,edge) == 1) THEN
                                If (ASSOCIATED(Wind)) Then
                                        Print *, "Yes Wind"
                                Else
                                        Print *,"No Wind"
                        End If 
                NULLIFY(Wind)
                CALL CreateWind(Wind)
                Wind%density=rhoOut
                Wind%temperature=pOut/rhoOut
                Wind%dir=i
                Wind%edge=edge
	        Wind%type=OUTFLOW_ONLY
         END IF
      END DO
   END DO
!   Print *, "outflowonly"
!  ALLOCATE(MyWinds(6))
!  nWinds=0
!  DO i=1,nDim
!    DO edge=1,2
!       IF (Gmthbc(i,edge) == 1) THEN 
!           nWinds=nWinds+1
!           CALL CreateWind(MyWinds(nWinds)%p)
!           MyWinds(nWinds)%p%dir=i
!           MyWinds(nWinds)%p%edge=edge
!  			MyWinds(nWinds)%p%type=OUTFLOW_ONLY
!        END IF
!     END DO
!  END DO
end if
		
     IF (lRestart) THEN
       particleList=>SinkParticles
		 i=1

       DO WHILE (ASSOCIATED(particlelist))

          particle=>particlelist%self
          READ(PROBLEM_DATA_HANDLE,NML=ParticleData)

          if (i==1) then
                Secondary=>Particle
                Particle%iAccrete=KRUMHOLZ_ACCRETION
                Particle%PointGravityObj%alpha=1
                Particle%PointGravityObj%soft_length=REAL(grav_soft_rad,qPREC)*sink_dx
                Particle%PointGravityObj%soft_function=SPLINESOFT
                Particle%PointGravityObj%Mass=Particle%Q(1)
                Particle%PointGravityObj%v0(1:nDim)=Particle%Q(imom(1:nDim))
                Particle%PointGravityObj%x0=Particle%xloc
!                Particle%PointGravityObj%luminosity=0
	  end if
      

        IF (i==2) THEN
                Primary=>Particle
                Particle%iAccrete=NOACCRETION
                Particle%PointGravityObj%alpha=0.067
                Particle%PointGravityObj%soft_length     = REAL(grav_soft_rad,qPREC)*sink_dx
                Particle%PointGravityObj%soft_function   =SPLINESOFT
                Particle%PointGravityObj%Mass=Particle%Q(1)
                Particle%PointGravityObj%v0(1:nDim)=Particle%Q(imom(1:nDim))
                Particle%PointGravityObj%x0=Particle%xloc
!                Particle%PointGravityObj%luminosity=1
                Particle%OutflowObj%duration         = 200
                Particle%OutflowObj%radius=REAL(radiusw,qPREC)*sink_dx
                Particle%OutflowObj%thickness=REAL(radiusw,qPREC)*sink_dx
                Particle%OutflowObj%open_angle                   = Pi/4
                Particle%OutflowObj%theta=Pi/2
                Particle%OutflowObj%phi=Pi/2
                Particle%OutflowObj%density=densw*2.0E4
                Particle%OutflowObj%temperature                  = 1d0
                Particle%OutflowObj%speed          = 0.01
                Particle%OutflowObj%velocity(1:nDim)   =Particle%Q(imom(1:nDim))
                Particle%OutflowObj%position= Particle%xloc
                Particle%OutflowObj%begin=-1d0
                CALL UpdateOutflow(Particle%OutflowObj)
        END IF
        particlelist=>particlelist%next
        i=i+1
       END DO


     ELSE!lrestart

         DO i=1,nParticles
            READ(PROBLEM_DATA_HANDLE,NML=ParticleData)
            NULLIFY(Particle)
	!IF (i==2) then
            CALL CreateParticle(Particle)
            
            Particle%Q(1)=mass
            Particle%xloc=xloc
            Particle%Q(imom(1:nDim))=vel(1:nDim)
            Particle%Buffer=buff
            CALL AddSinkParticle(Particle)
	!END IF
             
	 if (i.eq.1) then
                Secondary=>Particle
	        Particle%iAccrete=KRUMHOLZ_ACCRETION
                CALL CreatePointGravityObject(Particle%PointGravityObj)
	        Particle%PointGravityObj%alpha=1
                Particle%PointGravityObj%soft_length=REAL(grav_soft_rad,qPREC)*sink_dx
                Particle%PointGravityObj%soft_function=SPLINESOFT
                Particle%PointGravityObj%Mass=Particle%Q(1)
                Particle%PointGravityObj%v0(1:nDim)=Particle%Q(imom(1:nDim))
                Particle%PointGravityObj%x0=Particle%xloc
!                Particle%PointGravityObj%luminosity=1
                CALL CreateOutflow(Particle%OutflowObj)
                Particle%OutflowObj%duration=200
                Particle%OutflowObj%radius=REAL(radiusw,qPREC)*sink_dx
                Particle%OutflowObj%thickness=REAL(radiusw,qPREC)*sink_dx
                Particle%OutflowObj%open_angle=Pi/18
                Particle%OutflowObj%theta=Pi/2
                Particle%OutflowObj%phi=Pi/2
                Particle%OutflowObj%density=7.768E-7!densw*2.0E4
                Particle%OutflowObj%temperature=1d0
                Particle%OutflowObj%speed=55.92d0
                Particle%OutflowObj%velocity(1:nDim)= Particle%Q(imom(1:nDim))
                Particle%OutflowObj%position= Particle%xloc
                Particle%OutflowObj%begin=-1d0
                CALL UpdateOutflow(Particle%OutflowObj)
         end if
!
!!!
            IF (i.eq.2) THEN
                Primary=>Particle
		Particle%iAccrete=NOACCRETION
	        CALL CreatePointGravityObject(Particle%PointGravityObj)
		Particle%PointGravityObj%alpha=1
                Particle%PointGravityObj%soft_length=REAL(grav_soft_rad,qPREC)*sink_dx
                Particle%PointGravityObj%soft_function=SPLINESOFT
                Particle%PointGravityObj%Mass=Particle%Q(1)
                Particle%PointGravityObj%v0(1:nDim)=Particle%Q(imom(1:nDim))
                Particle%PointGravityObj%x0=Particle%xloc
!                Particle%PointGravityObj%luminosity=1
                CALL CreateOutflow(Particle%OutflowObj)
                Particle%OutflowObj%duration=200
                Particle%OutflowObj%radius=REAL(radiusw,qPREC)*sink_dx
                Particle%OutflowObj%thickness=REAL(radiusw,qPREC)*sink_dx
                Particle%OutflowObj%open_angle=Pi/2
                Particle%OutflowObj%theta=Pi/2
                Particle%OutflowObj%phi=Pi/2
                Particle%OutflowObj%density=9.441E-7!densw*2.0E4
                Particle%OutflowObj%temperature=1d0
                Particle%OutflowObj%speed=6.99
                Particle%OutflowObj%velocity(1:nDim)= Particle%Q(imom(1:nDim))
                Particle%OutflowObj%position= Particle%xloc
                Particle%OutflowObj%begin=-1d0
                CALL UpdateOutflow(Particle%OutflowObj)
            END IF
         END DO
   END IF!lrestart=0/1
CLOSE(PROBLEM_DATA_HANDLE)      

   END SUBROUTINE ProblemModuleInit

  !> Applies initial conditions
  !! @param Info Info object
  SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemGridInit

  !> Applies Boundary conditions
  !! @param Info Info object
  SUBROUTINE ProblemBeforeStep(Info)
    TYPE(InfoDef) :: Info
   INTEGER :: i, nParticles, grav_soft_rad
   REAL(KIND=xprec) :: beta,time, mass, xloc(3), vel(3),buff(17)
  	TYPE(ParticleListDef), POINTER :: particlelist
     TYPE(ParticleDef), POINTER :: Particle
	  LOGICAL,SAVE :: lock=.false.
     NAMELIST /ProblemData/ nParticles, densw, velw, windsPresent, disk, t1, t2
     NAMELIST /ParticleData/ mass,xloc,vel,alpha,radiusw,grav_soft_rad, buff

  END SUBROUTINE ProblemBeforeStep

  !> Could be used to update grids pre-output
  !! @param Info Info Object
  SUBROUTINE ProblemAfterStep(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemAfterStep

  !> Could be used to set force refinement
  !! @param Info Info object
  SUBROUTINE ProblemSetErrFlag(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemSetErrFlag

  SUBROUTINE ProblemBeforeGlobalStep(n)
   REAL(KIND=qPREc) :: remainder, period
   INTEGER :: n
   INTEGER :: i, nParticles, grav_soft_rad
   time=levels(n)%tnow
   particleList=>SinkParticles
   i=1
   DO WHILE (ASSOCIATED(particlelist))
      particle=>particlelist%self
      IF (i==2) Then
        Primary=>particle
      End if
      IF (i==1) Secondary=>particle
      particlelist=>particlelist%next
      i=i+1
   END DO
        

        !Primary%PointGravityObj%alpha=1
        !period=0.03328d0
        !remainder=MOD(time,period)
        !IF (time>=2.0 .AND. time<=3.0) THEN
        !        Primary%OutflowObj%open_angle=Pi/2
        !        Primary%OutflowObj%density=densw*100/4/2
        !        Primary%OutflowObj%speed=6.2d0
        !ELSE
        !        Primary%OutflowObj%open_angle=Pi/2
        !        Primary%OutflowObj%density=densw/8
        !        Primary%OutflowObj%speed=11.0d0
        !END IF

  END SUBROUTINE ProblemBeforeGlobalStep

END MODULE Problem

