gyro_surface_coordinates Subroutine

public subroutine gyro_surface_coordinates(gs, p, u)

Calculates the parametric coordinates, u, of point p on the gyro_surface

Arguments

Type IntentOptional AttributesName
type(GyroSurface), intent(in) :: gs

Gyro_surface

real(kind=Float64), intent(in), dimension(3):: p

Point on gyro_surface

real(kind=Float64), intent(out), dimension(2):: u

Parametric coordinates (gyro-angle, t)


Called by

proc~~gyro_surface_coordinates~~CalledByGraph proc~gyro_surface_coordinates gyro_surface_coordinates proc~gyro_range gyro_range proc~gyro_range->proc~gyro_surface_coordinates proc~npa_gyro_range npa_gyro_range proc~npa_gyro_range->proc~gyro_range proc~pnpa_f pnpa_f proc~pnpa_f->proc~npa_gyro_range proc~npa_mc npa_mc proc~npa_mc->proc~npa_gyro_range proc~pnpa_mc pnpa_mc proc~pnpa_mc->proc~npa_gyro_range proc~npa_f npa_f proc~npa_f->proc~npa_gyro_range program~fidasim fidasim program~fidasim->proc~pnpa_f program~fidasim->proc~npa_mc program~fidasim->proc~pnpa_mc program~fidasim->proc~npa_f

Contents


Source Code

subroutine gyro_surface_coordinates(gs, p, u)
    !+ Calculates the parametric coordinates, `u`, of point `p` on the gyro_surface
    type(GyroSurface), intent(in)            :: gs
        !+ Gyro_surface
    real(Float64), dimension(3), intent(in)  :: p
        !+ Point on gyro_surface
    real(Float64), dimension(2), intent(out) :: u
        !+ Parametric coordinates (gyro-angle, t)

    real(Float64), dimension(3) :: pp
    real(Float64) :: t, a, b, c, d, thm, thp, dp, dm, th
    integer :: i

    pp = matmul(transpose(gs%basis),p - gs%center)
    t = pp(3)/gs%axes(3)
    a = gs%axes(1) + gs%axes(2)*t
    b = gs%axes(2) - gs%axes(1)*t
    d = pp(1) + pp(2)
    c = max(min(d/sqrt(a**2 + b**2),1.d0),-1.d0)

    thm = -acos(c) + atan2(b,a)
    thp =  acos(c) + atan2(b,a)

    dm = norm2([gs%axes(1)*(cos(thm) - t*sin(thm)), &
                gs%axes(2)*(sin(thm) + t*cos(thm)), &
                gs%axes(3)*t ] - pp)
    dp = norm2([gs%axes(1)*(cos(thp) - t*sin(thp)), &
                gs%axes(2)*(sin(thp) + t*cos(thp)), &
                gs%axes(3)*t ] - pp)

    th = thm - pi/2
    if(dp.le.dm) th = thp - pi/2
    if(th.lt.0.0) th = th + 2*pi
    u = [th, t/gs%omega]

end subroutine gyro_surface_coordinates