module mod_interpolate

	implicit none
	
	!This module contains functions for interpolating between grid points
	
contains

function ly(yval) result(low)

	use mod_globvar

	implicit none

	!******* input/output variables *******
	real(8), intent(in) :: yval
	real(8) :: yval1, y_mod, y1
	integer :: low

	!******* compute low *******
	if (yval >= 0.0) then
		if (curv_ypos > 1.d0 + 1e-6) then
			yval1 = min(max_y,yval)
			y_mod = 1.d0/curv_ypos*log(yval1*(exp(curv_ypos) - 1.d0)/max_y + 1.d0)
			low   = ny_neg + 1 + floor(y_mod*real(ny_pos - 1,8))
		else
			yval1 = min(max_y,yval)
			y_mod = yval1**curv_ypos
			low   = ny_neg + 1 + floor(y_mod/dy_pos)
		endif
	else
		if (curv_yneg > 1.d0 + 1e-6) then
			yval1 = max(min_y,yval)
			y_mod = -1.d0/curv_yneg*log(-((yval1 - min_y)*(exp(curv_yneg) - 1.d0)/(-min_y) - exp(curv_yneg)))
			low   = 1 + floor((y_mod - (-1.d0))*real(ny_neg,8))
		else
			yval1 = max(min_y,yval)
			y_mod = -(-yval1)**curv_yneg
			y1    = -(-min_y)**curv_yneg
			low   = 1 + floor((y_mod - y1)/dy_neg)
		endif
	endif

end function ly

function uy(yval) result(high)

	use mod_globvar

	implicit none

	!******* input/output variables *******
	real(8), intent(in) :: yval
	real(8) :: yval1, y_mod, y1
	integer :: high

	!******* compute high *******
	if (yval >= 0.0) then
		if (curv_ypos > 1.d0 + 1e-6) then
			yval1 = min(max_y,yval)
			y_mod = 1.d0/curv_ypos*log(yval1*(exp(curv_ypos) - 1.d0)/max_y + 1.d0)
			high  = ny_neg + 1 + ceiling(y_mod*real(ny_pos - 1,8))
		else
			yval1 = min(max_y,yval)
			y_mod = yval1**curv_ypos
			high  = ny_neg + 1 + ceiling(y_mod/dy_pos)
		endif
	else
		if (curv_yneg > 1.d0 + 1e-6) then
			yval1 = max(min_y,yval)
			y_mod = -1.d0/curv_yneg*log(-((yval1 - min_y)*(exp(curv_yneg) - 1.d0)/(-min_y) - exp(curv_yneg)))
			high  = 1 + ceiling((y_mod - (-1.d0))*real(ny_neg,8))
		else
			yval1 = max(min_y,yval)
			y_mod = -(-yval1)**curv_yneg
			y1    = -(-min_y)**curv_yneg
			high  = 1 + ceiling((y_mod - y1)/dy_neg)
		endif
	endif

end function uy

function wy(yval) result(weight)

	use mod_globvar

	implicit none

	!******* input/output variables *******
	real(8), intent(in) :: yval
	real(8) :: yval1, y_mod, y1, weight
	integer :: low, high

	!******* compute weight *******
	if (yval >= 0.0) then
		if (curv_ypos > 1.d0 + 1e-6) then
			yval1 = min(max_y,yval)
			y_mod = 1.d0/curv_ypos*log(yval1*(exp(curv_ypos) - 1.d0)/max_y + 1.d0)
			low   = ny_neg + 1 + floor(y_mod*real(ny_pos - 1,8))
			high  = ny_neg + 1 + ceiling(y_mod*real(ny_pos - 1,8))
		else
			yval1 = min(max_y,yval)
			y_mod = yval1**curv_ypos
			low   = ny_neg + 1 + floor(y_mod/dy_pos)
			high  = ny_neg + 1 + ceiling(y_mod/dy_pos)
		endif
	else
		if (curv_yneg > 1.d0 + 1e-6) then
			yval1 = max(min_y,yval)
			y_mod = -1.d0/curv_yneg*log(-((yval1 - min_y)*(exp(curv_yneg) - 1.d0)/(-min_y) - exp(curv_yneg)))
			low   = 1 + floor((y_mod - (-1.d0))*real(ny_neg,8))
			high  = 1 + ceiling((y_mod - (-1.d0))*real(ny_neg,8))
		else
			yval1 = max(min_y,yval)
			y_mod = -(-yval1)**curv_yneg
			y1    = -(-min_y)**curv_yneg
			low   = 1 + floor((y_mod - y1)/dy_neg)
			high  = 1 + ceiling((y_mod - y1)/dy_neg)
		endif
	endif
	if (low == high) then
		weight = 0.0
	else
		weight = (yval1 - yy(low))/(yy(high) - yy(low))
	endif

end function wy

function lyfine(yfineval) result(low)

	use mod_globvar

	implicit none

	!******* input/output variables *******
	real(8), intent(in) :: yfineval
	real(8) :: yfineval1, yfine_mod, yfine1
	integer :: low

	!******* compute low *******
	if (yfineval >= 0.0) then
		if (curv_yfinepos > 1.d0 + 1e-6) then
			yfineval1 = min(max_yfine,yfineval)
			yfine_mod = 1.d0/curv_yfinepos*log(yfineval1*(exp(curv_yfinepos) - 1.d0)/max_yfine + 1.d0)
			low       = nyfine_neg + 1 + floor(yfine_mod*real(nyfine_pos - 1,8))
		else
			yfineval1 = min(max_yfine,yfineval)
			yfine_mod = yfineval1**curv_yfinepos
			low       = nyfine_neg + 1 + floor(yfine_mod/dyfine_pos)
		endif
	else
		if (curv_yfineneg > 1.d0 + 1e-6) then
			yfineval1 = max(min_yfine,yfineval)
			yfine_mod = -1.d0/curv_yfineneg*log(-((yfineval1 - min_yfine)*(exp(curv_yfineneg) - 1.d0)/(-min_yfine) - exp(curv_yfineneg)))
			low       = 1 + floor((yfine_mod - (-1.d0))*real(nyfine_neg,8))
		else
			yfineval1 = max(min_yfine,yfineval)
			yfine_mod = -(-yfineval1)**curv_yfineneg
			yfine1    = -(-min_yfine)**curv_yfineneg
			low       = 1 + floor((yfine_mod - yfine1)/dyfine_neg)
		endif
	endif

end function lyfine

function uyfine(yfineval) result(high)

	use mod_globvar

	implicit none

	!******* input/output variables *******
	real(8), intent(in) :: yfineval
	real(8) :: yfineval1, yfine_mod, yfine1
	integer :: high

	!******* compute high *******
	if (yfineval >= 0.0) then
		if (curv_yfinepos > 1.d0 + 1e-6) then
			yfineval1 = min(max_yfine,yfineval)
			yfine_mod = 1.d0/curv_yfinepos*log(yfineval1*(exp(curv_yfinepos) - 1.d0)/max_yfine + 1.d0)
			high      = nyfine_neg + 1 + ceiling(yfine_mod*real(nyfine_pos - 1,8))
		else
			yfineval1 = min(max_yfine,yfineval)
			yfine_mod = yfineval1**curv_yfinepos
			high      = nyfine_neg + 1 + ceiling(yfine_mod/dyfine_pos)
		endif
	else
		if (curv_yfineneg > 1.d0 + 1e-6) then
			yfineval1 = max(min_yfine,yfineval)
			yfine_mod = -1.d0/curv_yfineneg*log(-((yfineval1 - min_yfine)*(exp(curv_yfineneg) - 1.d0)/(-min_yfine) - exp(curv_yfineneg)))
			high      = 1 + ceiling((yfine_mod - (-1.d0))*real(nyfine_neg,8))
		else
			yfineval1 = max(min_yfine,yfineval)
			yfine_mod = -(-yfineval1)**curv_yfineneg
			yfine1    = -(-min_yfine)**curv_yfineneg
			high      = 1 + ceiling((yfine_mod - yfine1)/dyfine_neg)
		endif
	endif

end function uyfine

function wyfine(yfineval) result(weight)

	use mod_globvar

	implicit none

	!******* input/output variables *******
	real(8), intent(in) :: yfineval
	real(8) :: yfineval1, yfine_mod, yfine1, weight
	integer :: low, high

	!******* compute weight *******
	if (yfineval >= 0.0) then
		if (curv_yfinepos > 1.d0 + 1e-6) then
			yfineval1 = min(max_yfine,yfineval)
			yfine_mod = 1.d0/curv_yfinepos*log(yfineval1*(exp(curv_yfinepos) - 1.d0)/max_yfine + 1.d0)
			low       = nyfine_neg + 1 + floor(yfine_mod*real(nyfine_pos - 1,8))
			high      = nyfine_neg + 1 + ceiling(yfine_mod*real(nyfine_pos - 1,8))
		else
			yfineval1 = min(max_yfine,yfineval)
			yfine_mod = yfineval1**curv_yfinepos
			low       = nyfine_neg + 1 + floor(yfine_mod/dyfine_pos)
			high      = nyfine_neg + 1 + ceiling(yfine_mod/dyfine_pos)
		endif
	else
		if (curv_yfineneg > 1.d0 + 1e-6) then
			yfineval1 = max(min_yfine,yfineval)
			yfine_mod = -1.d0/curv_yfineneg*log(-((yfineval1 - min_yfine)*(exp(curv_yfineneg) - 1.d0)/(-min_yfine) - exp(curv_yfineneg)))
			low       = 1 + floor((yfine_mod - (-1.d0))*real(nyfine_neg,8))
			high      = 1 + ceiling((yfine_mod - (-1.d0))*real(nyfine_neg,8))
		else
			yfineval1 = max(min_yfine,yfineval)
			yfine_mod = -(-yfineval1)**curv_yfineneg
			yfine1    = -(-min_yfine)**curv_yfineneg
			low       = 1 + floor((yfine_mod - yfine1)/dyfine_neg)
			high      = 1 + ceiling((yfine_mod - yfine1)/dyfine_neg)
		endif
	endif
	if (low == high) then
		weight = 0.0
	else
		weight = (yfineval1 - yyfine(low))/(yyfine(high) - yyfine(low))
	endif

end function wyfine

end module mod_interpolate