Calculates the times of intersection of a line and a gyro-surface
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=Float64), | intent(in), | dimension(3) | :: | r0 | Point on line |
|
real(kind=Float64), | intent(in), | dimension(3) | :: | v0 | Direction of line |
|
type(GyroSurface), | intent(in) | :: | gs | Gyro-surface |
||
real(kind=Float64), | intent(out), | dimension(2) | :: | t | "time" to intersect |
subroutine line_gyro_surface_intersect(r0, v0, gs, t)
!+ Calculates the times of intersection of a line and a gyro-surface
real(Float64), dimension(3), intent(in) :: r0
!+ Point on line
real(Float64), dimension(3), intent(in) :: v0
!+ Direction of line
type(GyroSurface), intent(in) :: gs
!+ Gyro-surface
real(Float64), dimension(2), intent(out) :: t
!+ "time" to intersect
real(Float64), dimension(3) :: rr
real(Float64) :: a, b, c, d, tp, tm
rr = r0 - gs%center
a = dot_product(v0, matmul(gs%A,v0))
b = dot_product(rr, matmul(gs%A,v0)) + dot_product(v0,matmul(gs%A,rr))
c = dot_product(rr, matmul(gs%A,rr)) - 1.0
d = b**2 - 4*a*c
if(d.lt.0.0) then
t = 0.0
return
endif
t(1) = (-b - sqrt(d))/(2*a)
t(2) = (-b + sqrt(d))/(2*a)
end subroutine line_gyro_surface_intersect