newhut.f

Go to the documentation of this file.

    subroutine hut2(spin10,spin20,spin1,spin2,nsteps,dtau)
*
*       spin evolution of circular binary.
*       ----------------------------------

         implicit real*8 (a-h,o-z)
    real*8 u(2),udot(2)

    u(1)=spin10
    u(2)=spin20

         call deriv3(u,udot)

*       include step reduction for large derivatives.
         it = 0
     1   IF (abs(udot(1))*dtau.GT.0.01*spin10) THEN
             dtau = 0.5d0*dtau
             nsteps = 2*nsteps
             it = it + 1
             IF (it.LT.5) go to 1
         END IF

     5   IF (abs(udot(2))*dtau.GT.0.01*spin20) THEN
             dtau = 0.5d0*dtau
             nsteps = 2*nsteps
             it = it + 1
             IF (it.LT.5) go to 5
         END IF

    do i=1,nsteps
       call rk4c(dtau,u)
    enddo

    spin1=u(1)
    spin2=u(2)

    end

       subroutine rk4c(dt,u)

*   runge-kutta integrator.
*       -----------------------

*       author:  rosemary mardling(3/98).

       parameter(n=2)
       implicit real*8 (a-h,o-z)
       real*8 u0(n),ut(n),du(n),u(n)
       real*8 a(4),b(4)

       a(1)=dt/2.0d0
       a(2)=a(1)
       a(3)=dt
       a(4)=0.0d0

       b(1)=dt/6.0d0
       b(2)=dt/3.0d0
       b(3)=b(2)
       b(4)=b(1)

       do i=1,n
          u0(i)=u(i)
          ut(i)=0.0d0
       enddo

       do j=1,4
          call deriv3(u,du)

          do i=1,n
             u(i)=u0(i)+a(j)*du(i)
             ut(i)=ut(i)+b(j)*du(i)
          enddo
       enddo

       do i=1,n
          u(i)=u0(i)+ut(i)
       enddo

       end
    
    subroutine deriv3(u,udot)

         implicit real*8 (a-h,m,o-z)
    real*8 u(2),udot(2)
    common/spins/angmom0,rg2(2),m21,r21,semi0,c1,c2,c3,c4,c5,semi
         SAVE ic
         DATA ic /0/

*       assume e=0.
    spin1=u(1)
    spin2=u(2)

*       e2=e**2
*       e4=e2**2
*       e6=e4*e2
*       fac=1-e2

*       f2=1+7.5*e2+5.625*e4+0.3125*e6
*       f3=1+3.75*e2+1.875*e4+0.078125*e6
*       f4=1+1.5*e2+0.125*e4
*       f5=1+3*e2+0.375*e4


    semi=angmom0-rg2(1)*spin1-m21*r21**2*rg2(2)*spin2
*       semi=(semi*(1+m21)/m21/semi0**2)**2/fac
         semi=(semi*(1+m21)/m21/semi0**2)**2
         oa = 1.0/semi
    
*       udot(1)=(oa/fac)**6*c3*(oa**1.5*f2-fac**1.5*f5*spin1)
*       udot(2)=(oa/fac)**6*c4*(oa**1.5*f2-fac**1.5*f5*spin2)
         udot(1)=oa**6*c3*(oa**1.5-spin1) - c5*spin1
         udot(2)=oa**6*c4*(oa**1.5-spin2)

*       ic = ic + 1
*       IF (ic.LT.10) THEN
*           WRITE (6,1)  ic, spin1, spin2, udot
*   1       FORMAT (' HUT DERIV    # s1 s2 udot ',i8,1p,6e9.1)
*       END IF

    end