Revision dd8daa98413988b11c594e5a3bff8f5659d7ecf7 authored by Florent Renaud on 26 June 2014, 09:41:19 UTC, committed by Florent Renaud on 26 June 2014, 09:41:19 UTC
Originally released on 18 March 2013
Based on Nbody6 downloaded on 29 January 2013
1 parent edab1b7
Raw File
fchain.f
      SUBROUTINE FCHAIN(I,IR,XI,XIDOT,FIRR,FD)
*
*
*       Irregular force & derivative due to chain.
*       ------------------------------------------
*
      INCLUDE 'common6.h'
      COMMON/CHAINC/  XC(3,NCMAX),UC(3,NCMAX),BODYC(NCMAX),ICH,
     &                LISTC(LMAX)
      REAL*8  XI(3),XIDOT(3),DX(3),DV(3),FIRR(3),FD(3),XIS(3),VIS(3)
*
*
*       Form terms for the original chain c.m. interaction.
      DR2 = 0.0
      DRDV = 0.0
      DO 5 K = 1,3
          DX(K) = X(K,ICH) - XI(K)
          DV(K) = XDOT(K,ICH) - XIDOT(K)
          DR2 = DR2 + DX(K)**2
          DRDV = DRDV + DX(K)*DV(K)
    5 CONTINUE
      DR2I = 1.0/DR2
      DR3I = BODY(ICH)*DR2I*SQRT(DR2I)
      DRDV = 3.0*DRDV*DR2I
*
*       Subtract force and first derivative from current values.
      DO 10 K = 1,3
          FIRR(K) = FIRR(K) - DX(K)*DR3I
          FD(K) = FD(K) - (DV(K) - DX(K)*DRDV)*DR3I
   10 CONTINUE
*
*       Resolve chain coordinates & velocities using the predicted c.m.
      IF (IR.EQ.0) THEN
          CALL XCPRED(0)
      END IF
*
*       Obtain contributions from all members of the chain.
      DO 20 J = 1,NCH
          DR2 = 0.0
          DRDV = 0.0
          DO 15 L = 1,3
              DX(L) = XC(L,J) - XI(L)
              DV(L) = UC(L,J) - XIDOT(L)
	      DR2 = DR2 + DX(L)**2
	      DRDV = DRDV + DX(L)*DV(L)
   15     CONTINUE
          DR2I = 1.0/DR2
          DR3I = BODYC(J)*DR2I*SQRT(DR2I)
	  DRDV = 3.0*DRDV*DR2I
*
*       Add force & first derivative.
          DO 18 L = 1,3
              FIRR(L) = FIRR(L) + DX(L)*DR3I
              FD(L) = FD(L) + (DV(L) - DX(L)*DRDV)*DR3I
   18     CONTINUE
   20 CONTINUE
*
      RETURN
*
      END
back to top