https://github.com/florentrenaud/nbody6tt
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
trdot.f
SUBROUTINE TRDOT(I,DTM,M1)
*
*
* Time-scale for expansion of radius.
* -----------------------------------
*
INCLUDE 'common6.h'
REAL*8 TSCLS(20),LUMS(10),GB(10),TM,TN
REAL*8 M0,M1,RM,LUM,AGE,MC,MC1,RCC,RM0,AGE0,M10
REAL*8 menv,renv,k2
REAL*8 pts1,pts2,eps,alpha2,tol
PARAMETER(pts1=0.05d0,pts2=0.02d0)
PARAMETER(eps=1.0d-06,alpha2=0.09d0,tol=1.0d-10)
*
* Obtain stellar parameters at current epoch (body #I may be ghost).
KW = KSTAR(I)
M0 = BODY0(I)*ZMBAR
IF(M1.LE.0.0) M1 = RADIUS(I)*SU
M10 = M1
MC = 0.D0
AGE = TEV0(I)*TSTAR - EPOCH(I)
CALL star(KW,M0,M1,TM,TN,TSCLS,LUMS,GB,ZPARS)
CALL hrdiag(M0,AGE,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
& RM,LUM,KW,MC,RCC,MENV,RENV,K2)
*
* Quit if there is a change of type at the current TEV.
if((kstar(i).le.6.and.kw.gt.6).or.
& (kstar(i).le.9.and.kw.gt.9))then
m1 = m10
dtm = 0.d0
goto 10
endif
***
*
* Base new time scale for changes in radius & mass on stellar type.
if(kw.lt.0)then
dtm = pts1*tscls(15)
dtr = ABS(age)
elseif(kw.le.1)then
dtm = pts1*tm
dtr = tm - age
elseif(kw.ge.10)then
dtm = sqrt(10.0/lum)
dtr = dtm
elseif(kw.eq.2)then
dtm = pts1*(tscls(1) - tm)
dtr = tscls(1) - age
elseif(kw.eq.3)then
if(age.lt.tscls(6))then
dtm = pts2*(tscls(4) - age)
else
dtm = pts2*(tscls(5) - age)
endif
dtr = MIN(tscls(2),tn) - age
elseif(kw.eq.4)then
dtm = pts1*tscls(3)
dtr = MIN(tn,tscls(2) + tscls(3)) - age
elseif(kw.eq.5)then
if(age.lt.tscls(9))then
dtm = pts2*(tscls(7) - age)
else
dtm = pts2*(tscls(8) - age)
endif
dtr = MIN(tn,tscls(13)) - age
elseif(kw.eq.6)then
if(age.lt.tscls(12))then
dtm = pts2*(tscls(10) - age)
else
dtm = pts2*(tscls(11) - age)
endif
* dtm = MIN(dtm,0.005d0)
dtr = tn - age
elseif(kw.eq.7)then
dtm = pts1*tm
dtr = tm - age
elseif(kw.eq.8.or.kw.eq.9)then
if(age.lt.tscls(6))then
dtm = pts2*(tscls(4) - age)
else
dtm = pts2*(tscls(5) - age)
endif
dtr = tn - age
endif
*
* Record radius.
*
rm0 = rm
if(kw.ge.10) goto 30
age0 = age
kw0 = kw
mc1 = mc
*
* Check for type change.
*
it = 0
if((dtr-dtm).le.tol)then
*
* Check final radius for too large a jump.
*
age = MAX(age,age*(1.d0-eps) + dtr)
CALL hrdiag(M0,AGE,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
& RM,LUM,KW,MC1,RCC,MENV,RENV,K2)
dr = rm - rm0
if(ABS(dr).gt.0.1*rm0)then
dtm = dtr - age0*eps
dtdr = dtm/ABS(dr)
dtm = alpha2*MAX(rm,rm0)*dtdr
goto 20
else
dtm = dtr
goto 30
endif
endif
*
* Limit to a 10% increase assuming no further mass loss
* and thus that the pertubation functions due to small envelope mass
* will not change the radius.
*
20 age = age0 + dtm
mc1 = mc
CALL hrdiag(M0,AGE,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
& RM,LUM,KW,MC1,RCC,MENV,RENV,K2)
dr = rm - rm0
it = it + 1
if(it.eq.20.and.kw.eq.4) goto 30
IF(IT.GT.30)THEN
WRITE (6,22) IT, KSTAR(I), M0, DR, RM
22 FORMAT (' DANGER! TRDOT: IT K* M0 DR RM ',2I4,1P,3E10.2)
goto 30
ENDIF
if(ABS(dr).gt.0.1*rm0)then
dtdr = dtm/ABS(dr)
dtm = alpha2*MAX(rm0,rm)*dtdr
if(it.ge.20) dtm = 0.5d0*dtm
goto 20
endif
*
30 continue
*
* Impose a lower limit and convert time interval to scaled units.
DTM = MAX(DTM,1.0D-04)/TSTAR
*
* Use dummy for routine CHINIT (new chain or chain collision).
10 IF(IPHASE.EQ.8.OR.IPHASE.EQ.9)THEN
M1 = RM
ENDIF
*
RETURN
END