get_plasma_extrema Subroutine

subroutine get_plasma_extrema(r0, v0, extrema, x0, y0)
    !+ Returns extrema points where line(s) parametrized by `r0` and `v0` intersect the plasma boudnary
    real(Float64), dimension(:,:), intent(in)            :: r0, v0
    !+ Arrays the define line(s) in machine coordinates
    real(Float64), dimension(2,3), intent(out)           :: extrema
    !+ Minimum and maximumm R, Z, and Phi points
    real(Float64), dimension(:), intent(in), optional  :: x0, y0
    !+ Additional x and y points to consider

    real(Float64), dimension(:,:), allocatable :: xy_in, xy_out
    real(Float64), dimension(:,:), allocatable :: cyl_in, cyl_out
    real(Float64), dimension(:), allocatable  :: x, y, xlo, xhi
    real(Float64), dimension(:), allocatable  :: r, z, phi
    logical, dimension(:), allocatable :: skip

    real(Float64), dimension(3) :: ri, vi
    real(Float64) :: max_length, dlength
    integer :: i, iin, iout, ilo, ihi, dlo, dhi, dim, nlines, d
    logical :: inp

    nlines = size(r0(1,:))
    allocate(xy_in(2,nlines), xy_out(2,nlines))
    allocate(cyl_in(3,nlines), cyl_out(3,nlines))
    allocate(skip(nlines))

    dlength = 3.0 !cm
    skip = .False.

    loop_over_channels: do i=1, nlines
        ri = r0(:,i)
        vi = v0(:,i)
        vi = vi/norm2(vi)

        ! Find the position that the los first intersects the plasma
        call in_plasma(ri, inp, input_coords=1)
        max_length=0.0
        do while (.not.inp)
            ri = ri + vi*dlength ! move dlength
            call in_plasma(ri, inp, input_coords=1)
            max_length = max_length + dlength
            if(max_length.gt.1d9) then
                skip(i) = .True. !used below to skip los that do not intersect the plasma
                cycle loop_over_channels
            endif
        enddo

        xy_in(1,i) = ri(1) ; xy_in(2,i) = ri(2)
        call uvw_to_cyl(ri, cyl_in(:,i))

        ! Find the position that the los intersects upon exiting the plasma
        do while (inp)
            ri = ri + vi*dlength
            call in_plasma(ri, inp, input_coords=1)
        enddo

        xy_out(1,i) = ri(1) ; xy_out(2,i) = ri(2)
        call uvw_to_cyl(ri, cyl_out(:,i))

    enddo loop_over_channels

    dim = 2*count(.not.skip) ! 2 for enter and exit

    d = 0
    if (present(x0).and.present(y0)) then
        d = size(x0)
        if (size(x0).ne.size(y0)) then
            if(inputs%verbose.ge.0) then
                write(*,'("GET_PLASMA_EXTREMA: Sizes of X and Y input arrays are not identical")')
                stop
            endif
        endif
    endif

    allocate(x(dim+d), y(dim+d), r(dim), z(dim), phi(dim))

    iin = 1 ; iout = 2
    do i=1, nlines
        if (.not.skip(i)) then
            x(iin) = xy_in(1,i)  ; x(iout) = xy_out(1,i)
            y(iin) = xy_in(2,i)  ; y(iout) = xy_out(2,i)
            r(iin) = cyl_in(1,i) ; r(iout) = cyl_out(1,i)
            z(iin) = cyl_in(2,i) ; z(iout) = cyl_out(2,i)
            iin = iin + 2 ; iout = iout + 2
        endif
    enddo

    extrema(1,1) = minval(r) ; extrema(2,1) = maxval(r)
    extrema(1,2) = minval(z) ; extrema(2,2) = maxval(z)

    ! Append extra x and y points
    if(d.gt.0) then
        x(dim+1:dim+d) = x0
        y(dim+1:dim+d) = y0
    endif

    !! Domain is between 0 and 2 pi if all x points are left of the line x=0
    !! Else domain is between -pi and pi
    dlo = count(y.le.0) ; dhi = count(y.gt.0)
    allocate(xlo(dlo), xhi(dhi))

    ilo = 0 ; ihi = 0
    do i=1, size(x)
        if (y(i).le.0.d10) then
            ilo = ilo + 1
            xlo(ilo) = x(i)
        else
            ihi = ihi + 1
            xhi(ihi) = x(i)
        endif
    enddo

    phi = atan2(y, x)
    if (all(x.le.0) & !none in quadrant 1 and 4
        .or.(all(xlo.le.0).and.(size(x).ne.size(xhi))) & !quadrant 3
        .or.(all(xhi.le.0).and.(size(x).ne.size(xlo)))) then !quadrant 2
        phi = modulo(phi, 2*pi)
    endif

    extrema(1,3) = minval(phi) ; extrema(2,3) = maxval(phi)

end subroutine get_plasma_extrema