gyro_surface Subroutine

subroutine gyro_surface(fields, energy, pitch, gs)
    !+ Calculates the surface of all possible trajectories
    type(LocalEMFields), intent(in) :: fields
        !+ Electromagnetic fields at guiding center
    real(Float64), intent(in)       :: energy
        !+ Energy of particle
    real(Float64), intent(in)       :: pitch
        !+ Particle pitch w.r.t the magnetic field
    type(GyroSurface), intent(out)  :: gs
        !+ Gyro-surface

    integer :: i
    real(Float64) :: alpha, vabs, omega
    real(Float64), dimension(3,3) :: s

    vabs  = sqrt(energy/(v2_to_E_per_amu*inputs%ab))
    omega= (fields%b_abs*e0)/(inputs%ab*mass_u)
    alpha = vabs/omega

    gs%omega = omega
    gs%v = vabs
    gs%axes(1) = alpha*sqrt(1-pitch**2)
    gs%axes(2) = alpha*sqrt(1-pitch**2)
    gs%axes(3) = pitch*alpha

    s = 0.d0
    s(1,1) = gs%axes(1)**(-2)
    s(2,2) = gs%axes(2)**(-2)
    s(3,3) = -gs%axes(3)**(-2)

    gs%center = fields%pos

    gs%basis(:,1) = fields%a_norm
    gs%basis(:,2) = fields%c_norm
    gs%basis(:,3) = fields%b_norm

    gs%A = matmul(gs%basis,matmul(s,transpose(gs%basis)))

end subroutine gyro_surface