units.f

Go to the documentation of this file.

      SUBROUTINE units
*
*
*       Initialization of units & scaling factors.
*       ------------------------------------------
*
      include 'common6.h'
*
*
*       Define GM, AU and PC in cgs units (2009 IAU values).
      gm = 1.327 124 420 99 d+26
      au = 1.49597870700 d+13
      pc = 1296000.0d0/twopi*au
*
*       Form scaling factors for binary periods A*SQRT(A/M) to yrs and days.
      yrs = (rbar*1296000.0d0/twopi)**1.5/sqrt(zmbar)
      days = 365.25*yrs
*
*       Specify conversion factors for lengths to solar radii & AU.
      rsun = 6.960d+10
      su = pc/rsun*rbar
      rau = pc/au*rbar
*
*       Copy solar mass scaling to new variable (M = BODY*<M>).
      smu = zmbar
*
*       Form time scale in seconds and velocity scale in km/sec.
      tstar = sqrt(pc/gm)*pc
      vstar = 1.0d-05*sqrt(gm/pc)
*
*       Convert time scale from units of seconds to million years.
      tstar = tstar/(3.15576d+07*1.0d+06)
*
*       Ensure ZMBAR & RBAR > 0 (=0: assume <M>/Sun = 1, RBAR = 1 pc).
      IF (zmbar.LE.0.0d0) zmbar = float(n)/zmass
      IF (rbar.LE.0.0d0) rbar = 1.0
*
*       Scale to working units of RBAR in pc & ZMBAR in solar masses.
      tstar = tstar*sqrt(rbar**3/(zmass*zmbar))
      vstar = vstar*sqrt(zmass*zmbar/rbar)
*
*       Copy TSTAR to secondary time-scale factor (just in case).
      tscale = tstar
*
*       Physical scaling: X, M, V, T from RBAR*X, ZMBAR*M, VSTAR*V, TSTAR*T.
      WRITE (6,10)  rbar, zmbar, vstar, tstar, bodym*zmbar, su
   10 FORMAT (/,12x,'PHYSICAL SCALING:    R* =',f5.2,'  M* =',f8.1,
     &              '  V* =',f6.3,'  T* =',f6.3,'  <M> =',f5.2,
     &              '  SU =',1p,e8.1)
*
*       Define relevant parameter for the GR case (RZ = 6*<m>/c^2).
      IF (kz(27).EQ.3.OR.kz(28).GT.0) THEN
          clight = 3.0d+05/vstar
          rz = 6.0*zmass/(float(n)*clight**2)
          WRITE (6,20)  vstar, clight, rz
   20     FORMAT (/,12x,'GR SCALING:    V* =',1p,e10.2,'  C =',e10.2,
     &                                  '  RZ =',e10.2)
      ELSE
          rz = 0.0
          clight = 0.0
      END IF
*
      RETURN
*
      END