module testmixmain
  use testHelper
  use mixmain
  use mixtypes
  use voltapp
  use gamapp
  use uservoltic

  implicit none

  character(len=strLen) :: kaijuHome

contains

  @before
  subroutine firstSerial()
  end subroutine firstSerial

  @after
  subroutine lastSerial()
  end subroutine lastSerial

  @test
  subroutine testZeroInput()
    ! testing that an input of all zeroes gets output of all zeroes
    type(voltApp_T) :: voltronApp
    real(rp) :: curTilt
    character(len=strLen) :: caseInput = 'cmiD.xml'

    write(*,*) 'Doing testZeroInput ...'
    voltronApp%vOptions%gamUserInitFunc => initUser
    call initVoltron(voltronApp, caseInput)

    call convertGameraToRemix(voltronApp%mhd2mix, voltronApp%gApp, voltronApp%remixApp)

    call mapGameraToRemix(voltronApp%mhd2mix, voltronApp%remixApp)

    ! field aligned current in remixApp%ion(h)%St%Vars(:,:,FAC)
    voltronApp%remixApp%ion(NORTH)%St%Vars(:,:,FAC) = 0
    voltronApp%remixApp%ion(SOUTH)%St%Vars(:,:,FAC) = 0

    call voltronApp%tilt%getValue(voltronApp%time,curTilt)

    call run_mix(voltronApp%remixApp%ion, curTilt)

    ! potential in remixApp%ion(h)%St%Vars(:,:,POT)
    ! error value allows 1e-3 error per cell on average
    @assertLessThanOrEqual(sum(abs(voltronApp%remixApp%ion(NORTH)%S%Solution)), 8.1_rp, 'Non-zero remix output in North hemisphere')
    @assertLessThanOrEqual(sum(abs(voltronApp%remixApp%ion(SOUTH)%S%Solution)), 8.1_rp, 'Non-zero remix output in South hemisphere')
    
    write(*,*) 'Min/Max of potential (NORTH) = ', minval(voltronApp%remixApp%ion(NORTH)%S%Solution), maxval(voltronApp%remixApp%ion(NORTH)%S%Solution)
    write(*,*) 'Min/Max of potential (SOUTH) = ', minval(voltronApp%remixApp%ion(SOUTH)%S%Solution), maxval(voltronApp%remixApp%ion(SOUTH)%S%Solution)
    
    call CouplePotentialToMhd(voltronApp)

  end subroutine testZeroInput

  @test
  subroutine testConstantSolution()
    ! testing that an input of all zeroes gets output of all zeroes
    type(voltApp_T) :: voltronApp
    real(rp) :: curTilt

    real(rp) :: testValue = 240.19
    character(len=strLen) :: caseInput = 'cmiD.xml'

    write(*,*) 'Doing testConstantSolution ...'
    voltronApp%vOptions%gamUserInitFunc => initUser
    call initVoltron(voltronApp, caseInput)

    call convertGameraToRemix(voltronApp%mhd2mix, voltronApp%gApp, voltronApp%remixApp)

    call mapGameraToRemix(voltronApp%mhd2mix, voltronApp%remixApp)

    ! field aligned current in remixApp%ion(h)%St%Vars(:,:,FAC)
    voltronApp%remixApp%ion(NORTH)%St%Vars(:,:,FAC) = 0
    voltronApp%remixApp%ion(SOUTH)%St%Vars(:,:,FAC) = 0

    ! modify low lat boundary condition for test
    voltronApp%remixApp%ion(NORTH)%P%llbc_value = testValue
    voltronApp%remixApp%ion(SOUTH)%P%llbc_value = testValue

    call voltronApp%tilt%getValue(voltronApp%time,curTilt)

    call run_mix(voltronApp%remixApp%ion, curTilt)

    ! potential in remixApp%ion(h)%St%Vars(:,:,POT)
    ! error value allows 1e-3 error per cell on average
    @assertLessThanOrEqual(sum(abs(voltronApp%remixApp%ion(NORTH)%S%Solution - testValue)), 8.1_rp, 'Remix output in North hemisphere not specified value')
    @assertLessThanOrEqual(sum(abs(voltronApp%remixApp%ion(SOUTH)%S%Solution - testValue)), 8.1_rp, 'Remix output in South hemisphere not specified value')

  end subroutine testConstantSolution

  !@test
  subroutine testAzimuthallyDependentFAC()
    ! testing that an input of all zeroes gets output of all zeroes
    type(voltApp_T) :: voltronApp
    real(rp) :: curTilt

    character(len=strLen) :: caseInput = 'cmiD.xml'
    integer :: h,ip,it
    real(rp) :: thetaMin, thetaDelta

    write(*,*) 'Doing testAzimuthallyDependentFAC ...'
    voltronApp%vOptions%gamUserInitFunc => initUser
    call initVoltron(voltronApp, caseInput)

    call convertGameraToRemix(voltronApp%mhd2mix, voltronApp%gApp, voltronApp%remixApp)

    call mapGameraToRemix(voltronApp%mhd2mix, voltronApp%remixApp)

    ! field aligned current in remixApp%ion(h)%St%Vars(:,:,FAC)
    thetaMin   = (PI/180.0)*22.0_rp
    thetaDelta = (PI/180.0)*12.0_rp
    do h=1,size(voltronApp%remixApp%ion)
        do ip=1,voltronApp%remixApp%ion(h)%G%Np
            do it=1,voltronApp%remixApp%ion(h)%G%Nt
                if(voltronApp%remixApp%ion(h)%G%t(ip,it) .ge. thetaMin .and. voltronApp%remixApp%ion(h)%G%t(ip,it) .le. (thetaMin+thetaDelta)) then
                    voltronApp%remixApp%ion(h)%St%Vars(ip,it,FAC) = sin(voltronApp%remixApp%ion(h)%G%t(ip,it)) * sin(voltronApp%remixApp%ion(h)%G%p(ip,it))
                else
                    voltronApp%remixApp%ion(h)%St%Vars(ip,it,FAC) = 0.0_rp
                endif
            end do
        end do
    end do

    call voltronApp%tilt%getValue(voltronApp%time,curTilt)

    call run_mix(voltronApp%remixApp%ion, curTilt)

    ! potential in remixApp%ion(h)%St%Vars(:,:,POT)
    call CouplePotentialToMhd(voltronApp)

  end subroutine testAzimuthallyDependentFAC

end module testmixmain
