Commit 1e869cc0 authored by Bruce Cowan's avatar Bruce Cowan

Change solve_kepler to function

parent cc7d6917
......@@ -14,20 +14,22 @@ ecc = 0.00877_dp
! E should be ~2.016365 for M(1)
! E should be 0 for M(2)
call solve_kepler(M, ecc, E)
E = solve_kepler(M, ecc)
print *, E(1)
print *, E(2)
contains
subroutine solve_kepler(M, ecc, E)
real(dp), intent(in) :: M(:)
real(dp), intent(in) :: ecc
real(dp), intent(out) :: E(:)
integer :: i
real(dp) :: delta
function solve_kepler(M, ecc) result(E)
real(dp), intent(in), dimension(:) :: M
real(dp), intent(in) :: ecc
real(dp), allocatable, dimension(:) :: E
integer :: i
real(dp) :: delta
do i=1,size(M)
allocate(E(size(M)))
do i = 1, size(M)
E(i) = M(i)
do
delta = E(i) - ecc * sin(E(i)) - M(i)
......@@ -38,6 +40,6 @@ subroutine solve_kepler(M, ecc, E)
E(i) = E(i) - delta / (1 - ecc * cos(E(i)))
end do
end do
end subroutine solve_kepler
end function solve_kepler
end program kepler
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment