trdot.f

Go to the documentation of this file.

      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 = 1.0d+02
         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