Newer
Older
8001
8002
8003
8004
8005
8006
8007
8008
8009
8010
8011
8012
8013
8014
8015
8016
8017
8018
8019
8020
8021
8022
8023
8024
8025
8026
8027
8028
8029
8030
8031
8032
8033
8034
8035
8036
8037
8038
8039
8040
8041
8042
8043
8044
8045
8046
8047
8048
8049
8050
8051
8052
8053
8054
8055
8056
8057
8058
8059
8060
8061
8062
8063
8064
8065
8066
8067
8068
8069
8070
8071
8072
8073
8074
8075
8076
8077
8078
8079
8080
8081
8082
8083
8084
8085
8086
8087
8088
8089
8090
8091
8092
8093
8094
8095
8096
8097
8098
8099
8100
8101
8102
8103
8104
8105
8106
8107
8108
8109
8110
8111
8112
8113
8114
8115
8116
8117
8118
8119
8120
8121
8122
8123
8124
8125
8126
8127
8128
8129
8130
8131
8132
8133
8134
8135
8136
8137
8138
8139
8140
8141
8142
8143
8144
8145
8146
8147
8148
8149
8150
8151
8152
8153
*
* ALGORITHM: Uses power series for N in terms of F and Newton,s method
* REMARKS: ONLY GOOD FOR LOW VALUES OF N (N < 0.636*e -0.6)
* AUTHOR: M. Duncan
* DATE WRITTEN: May 26, 1992.
* REVISIONS:
***********************************************************************
real*8 function orbel_flon(e,capn)
include 'swift.inc'
c... Inputs Only:
real*8 e,capn
c... Internals:
integer iflag,i,IMAX
real*8 a,b,sq,biga,bigb
real*8 x,x2
real*8 f,fp,dx
real*8 diff
real*8 a0,a1,a3,a5,a7,a9,a11
real*8 b1,b3,b5,b7,b9,b11
PARAMETER (IMAX = 10)
PARAMETER (a11 = 156.d0,a9 = 17160.d0,a7 = 1235520.d0)
PARAMETER (a5 = 51891840.d0,a3 = 1037836800.d0)
PARAMETER (b11 = 11.d0*a11,b9 = 9.d0*a9,b7 = 7.d0*a7)
PARAMETER (b5 = 5.d0*a5, b3 = 3.d0*a3)
c----
c... Executable code
c Function to solve "Kepler's eqn" for F (here called
c x) for given e and CAPN. Only good for smallish CAPN
iflag = 0
if( capn .lt. 0.d0) then
iflag = 1
capn = -capn
endif
a1 = 6227020800.d0 * (1.d0 - 1.d0/e)
a0 = -6227020800.d0*capn/e
b1 = a1
c Set iflag nonzero if capn < 0., in which case solve for -capn
c and change the sign of the final answer for F.
c Begin with a reasonable guess based on solving the cubic for small F
a = 6.d0*(e-1.d0)/e
b = -6.d0*capn/e
sq = sqrt(0.25*b*b +a*a*a/27.d0)
biga = (-0.5*b + sq)**0.3333333333333333d0
bigb = -(+0.5*b + sq)**0.3333333333333333d0
x = biga + bigb
c write(6,*) 'cubic = ',x**3 +a*x +b
orbel_flon = x
c If capn is tiny (or zero) no need to go further than cubic even for
c e =1.
if( capn .lt. TINY) go to 100
do i = 1,IMAX
x2 = x*x
f = a0 +x*(a1+x2*(a3+x2*(a5+x2*(a7+x2*(a9+x2*(a11+x2))))))
fp = b1 +x2*(b3+x2*(b5+x2*(b7+x2*(b9+x2*(b11 + 13.d0*x2)))))
dx = -f/fp
c write(6,*) 'i,dx,x,f : '
c write(6,432) i,dx,x,f
432 format(1x,i3,3(2x,1p1e22.15))
orbel_flon = x + dx
c If we have converged here there's no point in going on
if(abs(dx) .le. TINY) go to 100
x = orbel_flon
enddo
c Abnormal return here - we've gone thru the loop
c IMAX times without convergence
if(iflag .eq. 1) then
orbel_flon = -orbel_flon
capn = -capn
endif
write(6,*) 'FLON : RETURNING WITHOUT COMPLETE CONVERGENCE'
diff = e*sinh(orbel_flon) - orbel_flon - capn
write(6,*) 'N, F, ecc*sinh(F) - F - N : '
write(6,*) capn,orbel_flon,diff
return
c Normal return here, but check if capn was originally negative
100 if(iflag .eq. 1) then
orbel_flon = -orbel_flon
capn = -capn
endif
return
end ! orbel_flon
c
***********************************************************************
c ORBEL_ZGET.F
***********************************************************************
* PURPOSE: Solves the equivalent of Kepler's eqn. for a parabola
* given Q (Fitz. notation.)
*
* Input:
* q ==> parabola mean anomaly. (real scalar)
* Returns:
* orbel_zget ==> eccentric anomaly. (real scalar)
*
* ALGORITHM: p. 70-72 of Fitzpatrick's book "Princ. of Cel. Mech."
* REMARKS: For a parabola we can solve analytically.
* AUTHOR: M. Duncan
* DATE WRITTEN: May 11, 1992.
* REVISIONS: May 27 - corrected it for negative Q and use power
* series for small Q.
***********************************************************************
real*8 function orbel_zget(q)
include 'swift.inc'
c... Inputs Only:
real*8 q
c... Internals:
integer iflag
real*8 x,tmp
c----
c... Executable code
iflag = 0
if(q.lt.0.d0) then
iflag = 1
q = -q
endif
if (q.lt.1.d-3) then
orbel_zget = q*(1.d0 - (q*q/3.d0)*(1.d0 -q*q))
else
x = 0.5d0*(3.d0*q + sqrt(9.d0*(q**2) +4.d0))
tmp = x**(1.d0/3.d0)
orbel_zget = tmp - 1.d0/tmp
endif
if(iflag .eq.1) then
orbel_zget = -orbel_zget
q = -q
endif
return
end ! orbel_zget
c----------------------------------------------------------------------