Calculates the parametric coordinates, u
, of point p
on the gyro_surface
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
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) |
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