Commit d01f7271 authored by Bruce Cowan's avatar Bruce Cowan

Make use of masking in kepler

parent f8cc7af9
Pipeline #2115 passed with stage
in 1 minute and 42 seconds
......@@ -7,6 +7,8 @@ implicit none
real(dp), dimension(2) :: M
real(dp) :: ecc
real(dp), dimension(size(M)) :: E
integer :: i
M = [2.0084512_dp, 0.0_dp]
ecc = 0.00877_dp
......@@ -15,30 +17,24 @@ ecc = 0.00877_dp
! E should be 0 for M(2)
E = solve_kepler(M, ecc)
print *, E(1)
print *, E(2)
print '(f10.8)', (E(i), i=1, size(E))
contains
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
allocate(E(size(M)))
do i = 1, size(M)
E(i) = M(i)
do
delta = E(i) - ecc * sin(E(i)) - M(i)
if (abs(delta) <= 1e-9_dp) then
exit
end if
E(i) = E(i) - delta / (1 - ecc * cos(E(i)))
end do
real(dp), intent(in), dimension(:) :: M
real(dp), intent(in) :: ecc
real(dp), dimension(size(M)) :: E
real(dp), dimension(size(M)) :: delta
E = M
delta = 1
do while (any(delta > 1e-9_dp))
where (delta > 1e-9_dp)
delta = E - ecc * sin(E) - M
E = E - delta / (1 - ecc * cos(E))
end where
end do
end function solve_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