module coords_polar_m !! Coordinate transformations Cartesian (x,y) <-> Polar (rho, theta) !! rho is normalized flux surface label !! theta is geometric poloidal angle, measured ccw from +x with !! respect to magnetic axis use precision_m, only: FP, DP use constants_m, only: PI use screen_io_m, only: get_stderr use comm_handling_m, only: is_master use equilibrium_m, only: equilibrium_t use slab_equilibrium_m, only: slab_t use root_finding_m, only: find_zero, func_1D_t implicit none public :: cart_to_polar public :: polar_to_cart public :: jacobian_polar private :: aux_ptc_zero private :: set_equilibrium type, private, extends(func_1D_t) :: func_1D_equi_t class(equilibrium_t), private, pointer :: equi => null() contains procedure :: func => aux_ptc_zero procedure, private :: set_equilibrium end type contains subroutine cart_to_polar(equi, x, y, phi, rho, theta) !! Coordinate transformation from polar to Cartesian implicit none class(equilibrium_t), intent(inout) :: equi !! Equilibrium real(FP), intent(in) :: x !! Cartesian coordinate x real(FP), intent(in) :: y !! Cartesian coordinate y real(FP), intent(in) :: phi !! Toroidal angle real(FP), intent(out) :: rho !! Polar coordinate rho (normalized flux surface label) real(FP), intent(out) :: theta !! Polar coordinate theta (geometric poloidal angle) real(FP) :: xt, yt select type(equi) type is(slab_t) rho = x theta = 2.0_FP * PI * y / (equi%ymax - equi%ymin) return end select rho = equi%rho(x, y, phi) xt = x - equi%x0 yt = y - equi%y0 if (xt > 0.0_FP) then theta = atan(yt/xt) elseif (xt < 0.0_FP) then theta = atan(yt/xt) + PI else if (yt > 0.0_FP) then theta = PI/2.0_FP elseif (yt < 0.0_FP) then theta = PI*3.0_FP/2.0_FP else theta = 0.0_8 endif endif theta = mod(theta, 2.0_FP * PI) if (theta < 0.0_FP) then theta = theta + 2.0_FP*PI endif end subroutine subroutine polar_to_cart(equi, rho, theta, phi, x, y) !! Coordinate transformation from Cartesian to polar ! Basic description of algorithm: ! - Searches along contours of given theta ! - Finds distance t to magnetic axis, where rho matches via zero point search) ! - Interval of zero point search is gradually increased (tryloop, trying, e.g. to avoid to run of bspline grid for equilibria) ! Further notes: ! If an error is being returned that no point is found, try increasing ! either maximum_number_of_tries (more expensive) or ! increase_range_factor (more risky -- might run off the spline grid for ! NUMERICAL equilibrium) ! ! If an error from NUMERICAL psi (calling off the spline grid) is raised ! from this region, try reducing the increase_range_factor implicit none class(equilibrium_t), intent(inout) :: equi !! Equilibrium real(FP), intent(in) :: rho !! Polar coordinate rho real(FP), intent(in) :: theta !! Polar coordinate theta (geometric poloidal angle) real(FP), intent(in) :: phi !! Toroidal angle real(FP), intent(out) :: x !! Cartesian coordinate x real(FP), intent(out) :: y !! Cartesian coordinate y integer, parameter :: ntry_max = 100 ! Maximum number of intervalls searched real(DP), parameter :: tolerance = 1.0E-8_DP ! Tolerance for zero search real(DP), parameter :: intervall_increase_factor = 1.1_DP ! Factor that search intervall is increased at each search integer :: iuser(1), find_zero_status, try_count real(DP) :: t, ruser(3) real(DP) :: intervall_min, intervall_max real(DP) :: axis_x, axis_y type(func_1D_equi_t) :: aux_ptc_zero select type(equi) type is(slab_t) x = rho y = (equi%ymax - equi%ymin) * theta / (2.0_FP * PI) return end select if(rho < tolerance) then ! Return the magnetic axis if rho is smaller than the tolerance call equi%mag_axis_loc(phi, x, y) return endif call aux_ptc_zero%set_equilibrium(equi) ! start to search along contours of theta ruser(1) = real(rho, kind=DP) ruser(2) = real(theta, kind=DP) ruser(3) = real(phi, kind=DP) intervall_min = 0.0_DP intervall_max = 5.0E-2_DP try_count = 0 tryloop: do call find_zero(intervall_min, intervall_max, tolerance, aux_ptc_zero, iuser, ruser, t, find_zero_status) if (find_zero_status == 0) then ! normal exit, point found call equi%mag_axis_loc(phi, axis_x, axis_y) x = axis_x + abs(real(t, kind=FP)) * cos(theta) y = axis_y + abs(real(t, kind=FP)) * sin(theta) exit tryloop elseif (try_count > ntry_max) then write(get_stderr(), *) & 'error(polar_to_cart): point not found (rho, theta):', & rho, theta write(get_stderr(), *) & '... reached interval_max= ', intervall_max error stop else try_count = try_count + 1 intervall_max = intervall_increase_factor * intervall_max endif enddo tryloop end subroutine polar_to_cart real(FP) function jacobian_polar(equi, x, y, phi) !! Jacobian of polar coordinate system !! To compute derivatives on rho, finite difference is used !! N.B. currently only axisymmetric (2D), proper inclusion of toroidal !! phi will change this function implicit none class(equilibrium_t), intent(inout) :: equi !! Equilibrium real(FP), intent(in) :: x !! Cartesian coordinate x real(FP), intent(in) :: y !! Cartesian coordinate y real(FP), intent(in) :: phi !! Toroidal angle real(FP), parameter :: eps_findiff=1.0E-6_FP !! Finite distance to comput derivatives on rho real(FP) :: xa, ya, rho1,rho2, theta, drhodx, drhody real(FP) :: grhorho, gthetatheta, gphiphi, grhotheta select type(equi) type is(slab_t) jacobian_polar = (equi%ymax-equi%ymin) / (2*PI) return end select xa = x-equi%x0 ya = y-equi%y0 ! components of metric tensor gphiphi = 1.0_FP / (equi%jacobian(x, y, phi)**2) gthetatheta=1.0_FP/(xa**2+ya**2) call cart_to_polar(equi, x+eps_findiff, y, phi, rho1, theta ) call cart_to_polar(equi, x-eps_findiff, y, phi, rho2, theta ) drhodx = (rho1-rho2) / (2.0_FP*eps_findiff) call cart_to_polar(equi, x, y+eps_findiff, phi, rho1, theta ) call cart_to_polar(equi, x, y-eps_findiff, phi, rho2, theta ) drhody = (rho1-rho2) / (2.0_8*eps_findiff) grhorho=drhodx**2+drhody**2 grhotheta=(-ya*drhodx + xa*drhody) / (xa**2+ya**2) ! polar jacobian jacobian_polar = grhorho * gthetatheta * gphiphi - gphiphi * grhotheta * grhotheta jacobian_polar = 1.0_FP/sqrt(jacobian_polar) end function subroutine set_equilibrium(this, equi) !! Sets equilibrium pointer for zero point search function class(func_1D_equi_t), intent(inout) :: this !! Instance of the type class(equilibrium_t), target, intent(in) :: equi !! Equilibrium this%equi => equi end subroutine set_equilibrium real(DP) function aux_ptc_zero(this, t, iuser, ruser) ! Auxiliary function for routine polar_to_cart, which evaluates zero if rho_target = rho for fixed theta class(func_1D_equi_t), intent(in) :: this !! Instance of the type real(DP),intent(in) :: t !! Distance to magnetic axis integer,intent(in) :: iuser(:) !! Dummy not used real(DP),intent(in) :: ruser(:) !! ruser(1) : rho_target, ruser(2) : theta, ruser(3) : phi real(FP) :: rho_target, theta, phi, x, y, rho_trial, theta_trial real(FP) :: axis_x, axis_y rho_target = real(ruser(1), kind=FP) theta = real(ruser(2), kind=FP) phi = real(ruser(3), kind=FP) call this%equi%mag_axis_loc(phi, axis_x, axis_y) x = axis_x + abs(t) * cos(theta) y = axis_y + abs(t) * sin(theta) call cart_to_polar(this%equi, x, y, phi, rho_trial, theta_trial) aux_ptc_zero = real(rho_trial-rho_target, kind=DP) end function end module