himod.f

Go to the documentation of this file.

      SUBROUTINE himod(evec,hvec,m3,mt,EOvec,HOvec,ICALL,dt0,dt,e,r,v)
*
*
*       Modification of hierarchical binary.
*       ------------------------------------
*
*       Author:  Rosemary Mardling (3/98).

      IMPLICIT REAL*8 (a-h,m,o-z)
      common/rksave/  coeff,hohat(3),e0,a,hh,mb
      REAL*8  evec(3),hvec(3),eovec(3),hovec(3),
     &        element(6),r(3),v(3),u(6),udot(6)
*
*
*   Given inner Runge-Lenz vector {\bf e}=evec and {\bf h}=hvec.
*   Calculate eccentricity e, magnitude of hvec and semi-major axis a.

      mb=mt-m3
      e0=sqrt(dot(evec,evec))
      hh=sqrt(dot(hvec,hvec))
      a=hh**2/mb/(1-e0**2)

*   Given OUTER {\bf E}=Eout and {\bf H}=Hout.
*   Calculate outer eccentricity and semi-major axis.

      eout=sqrt(dot(eovec,eovec))
      hhout=sqrt(dot(hovec,hovec))
      aout=hhout**2/mt/(1-eout**2)
      do i=1,3
         hohat(i)=hovec(i)/hhout
      enddo

*   Calculate the semi-minor axis Bout and the semi-latus rectum SLRout.

      bout=aout*sqrt(1-eout**2)
      slrout=aout*(1-eout**2)
      coeff=m3/(2*aout*bout*slrout)
      dtsum = 0.0

      IF (icall.LT.0) THEN
         do i=1,3
            u(i)=evec(i)
            u(i+3)=hvec(i)
         enddo
         CALL deriv(u,udot,icall)
         tau = e0/sqrt(udot(1)**2+udot(2)**2+udot(3)**2)
         WRITE (6,5)  e0,tau,a,eout,aout,(udot(k),k=1,3)
    5    FORMAT (' DERIV   e0 TAU a E1 A1 edot  ',
     &                     f8.4,1p,2e10.2,0p,f8.4,1p,4e10.2)
         CALL flush(3)
         icall = icall + 1
      END IF

*       Initialize integration variables.

   10 do i=1,3
     u(i)=evec(i)
     u(i+3)=hvec(i)
      enddo

*       Perform integration.

      call rkint(dt,u)

      if(e0.ge.1.0) goto 20
      eh=0.0
      do i=1,3
     evec(i)=u(i)
     hvec(i)=u(i+3)
         eh=eh+evec(i)*hvec(i)
      enddo

    e0=sqrt(dot(evec,evec))
    hh=sqrt(dot(hvec,hvec))

*       Apply symmetric re-orthogonalization (small eh; Seppo Mikkola 3/98).

        eps=-eh/(e0**2+hh**2) 
      eps=0.0
        do k=1,3
           evec(k)=evec(k)+eps*hvec(k)
           hvec(k)=hvec(k)+eps*evec(k)
        end do

      IF (dtsum + dt.GT.dt0) dt = dt0 - dtsum + 1.0d-10
      dtsum = dtsum + dt
      IF (dtsum.LT.dt0) go to 10

*   Calculate new orbital parameters.
 20   e2=dot(evec,evec)
      h2=dot(hvec,hvec)
      a=h2/mb/(1-e2)
      e=sqrt(e2)
      hh=sqrt(h2)

*   Calculate updated vectors r and v. 
*   Note that the vectors evec and hvec give no information about
*   the orbital phase. Thus there is a kind of degeneracy. Since
*   we are updating to a later time using an approximation, we can't
*   expect to know it anyway. Hence the phase of peri is set to zero.

      call transform4(evec,hvec,mb,element)
      call transform2(element,mb,r,v)

      RETURN
      END

    
    subroutine cross(u,v,w)

*       Vectorial cross product.
*       ------------------------
    real*8 u(3),v(3),w(3)

    w(1) = u(2)*v(3) - u(3)*v(2)
    w(2) = u(3)*v(1) - u(1)*v(3)
    w(3) = u(1)*v(2) - u(2)*v(1)

    end


        real*8 function dot(u,v)

*       Scalar dot product.
*       ------------------
        real*8 u(3),v(3)

    dot=u(1)*v(1)+u(2)*v(2)+u(3)*v(3)

    end


    subroutine transform2(element,mb,r,v)
*
*   Calculates vectors r and v given 6 orbital elements.
*       ----------------------------------------------------
*
        IMPLICIT REAL*8 (a-h,o-z)
    real*8 r(3),v(3),element(6),mb
    real*8 tempr(3),tempv(3)
    real*8 inc,b(3,3)


        a=element(1)
    e=element(2)
    inc=element(3)
    w=element(4)
    om=element(5)
    phi=element(6)

    cosp=cos(phi)
    sinp=sin(phi)
    cosi=cos(inc)
    sini=sin(inc)
    cosw=cos(w)
    sinw=sin(w)
    cosom=cos(om)
    sinom=sin(om)

        b(1,1) = cosw*cosom - sinw*cosi*sinom
        b(1,2) = cosw*sinom + sinw*cosi*cosom
        b(1,3) = sinw*sini
        b(2,1) = -sinw*cosom - cosw*cosi*sinom
        b(2,2) = -sinw*sinom + cosw*cosi*cosom
        b(2,3) = cosw*sini
        b(3,1) = sini*sinom
        b(3,2) = -sini*cosom
        b(3,3) = cosi

    h=sqrt(mb*a*(1-e**2))
    rr=a*(1-e**2)/(1+e*cosp)
    rd=e*h*sinp/(a*(1-e**2))
    phid=h/rr**2

    r(1)=rr*cosp
    r(2)=rr*sinp
    r(3)=0.0

    v(1)=rd*cosp-rr*phid*sinp
    v(2)=rd*sinp+rr*phid*cosp
    v(3)=0.0

    do i=1,3
       sum1=0.0
       sum2=0.0
       do j=1,3
          sum1=sum1 + b(j,i)*r(j)
          sum2=sum2 + b(j,i)*v(j)
       enddo
       tempr(i)=sum1
       tempv(i)=sum2
    enddo
    do i=1,3
       r(i)=tempr(i)
       v(i)=tempv(i)
    enddo

    end


    subroutine transform4(e,h,mb,element)
*
*   Calculates 5 orbital elements given vectors evec and h and mass mb.
*   -------------------------------------------------------------------
*
        IMPLICIT REAL*8 (a-h,o-z)
    real*8 element(6)
    real*8 h(3),e(3),n(3),hn(3),he(3)
    real*8 eh(3)
    real*8 ii(3),jj(3),kk(3),nn,mb
    data ii/1.d0,0.d0,0.d0/
    data jj/0.d0,1.d0,0.d0/
    data kk/0.d0,0.d0,1.d0/


    hh=sqrt(dot(h,h))

        do i=1,3
           h(i)=h(i)/hh
        enddo

    call cross(kk,h,n)

    nn=sqrt(dot(n,n))
    ecc=sqrt(dot(e,e))
        a=hh**2/mb/(1-ecc**2)

    do i=1,3
       n(i)=n(i)/nn
       eh(i)=e(i)/ecc
    enddo

    call cross(h,n,hn)
    call cross(h,eh,he)

    cosom=dot(ii,n)
    sinom=dot(jj,n)
    cosi=dot(kk,h)
    sini=nn
    cosw=dot(n,eh)
    sinw=dot(eh,hn)

    element(1)=a
    element(2)=ecc
    element(3)=atan2(sini,cosi)
    element(4)=atan2(sinw,cosw)
    element(5)=atan2(sinom,cosom)
    element(6)=0.0
*       Pericentre phase (can be specified).

*     ZI = 360.0*ELEMENT(3)/(8.0*ATAN(1.0D0))
*     WRITE (6,5)  ECC,ELEMENT(3),ELEMENT(4),ELEMENT(5),ZI
*   5 FORMAT (' TRANSFORM    e i w Om IN ',F8.4,F9.3,3F9.3)
    end