program skew

c  This is a slightly modified version of a program listed the 
c  solutions manual of "An Introduction to Modern Astrophysics," 
c  Bradley W. Carroll and Dale A. Ostlie,
c  Addison-Wesley Publishing Company, copyright 1996.

       implicit none
       integer*4 j, N
       parameter(N = 2001)
       real*8 v1r, v2r, theta, dtheta, t, dt, pi, G, r, v, vr
       real*8 phi, i, sini, yr, vrcom
       real*8 m1sun, m2sun, m1, m2, P, a, aAU, Msun, AU, e, M,
     1     mu, L, dAdt
c
c Compute the radial velocities of both stars in a binary system
c as functions of time for any orientation and eccentricity of 
c the binary system. We work in cgs units!

       data pi, G, Msun, AU, yr / 3.1415926535897932d+00,
     1  6.67428d-8, 1.9891d33, 1.4959787066d13,
     2  3.15581495d7/
c
c the star's masses

       print*
       write(*,*) ' Enter M1 in solar masses'
       read(*,*) m1sun
       print*

       write(*,*) ' Enter M2 in solar masses'
       read(*,*) m2sun
       print*

       m1 = m1sun*Msun
       m2 = m2sun*Msun
       M = m1 + m2
       mu = m1*m2/M  !reduced mass of the system

c the geometrical parameters

       write(*,*) ' Enter semi-major axis in AU'
       read(*,*) aAU
       a = aAU*AU
9      print*
       write(*,*) ' Enter eccentricity, e '
       read(*,*) e
       if (e .ge. 1.) then  ! change 8 Feb. 20006; it originally tested 'i'
         write(6,*) 'orbit unbound; enter another value of e < 1'
         go to 9
       endif

       print*
       write(*,*) ' Enter radial velocity of center of mass (km/s) '
       read(*,*) vrcom

10     print*
       write(*,*) ' Enter inclination wrt los in degrees  '
       read(*,*) i
       if (i .eq. 0.) then
         write(6,*) 'orbit in plane of sky; enter another value of i'
         go to 10
       elseif(i .gt. 90.) then
         write(6,*) 'angle must be less than or equal to 90 degrees'
         go to 10
       endif
       i = pi*i/180.d0  !convert to radians
       sini = sin(i)  !take the sine

11     print*
       write(*,*) ' Enter angle: major axis wrt los in degrees '
       read(*,*) phi
       phi = pi*phi/180.d0  !convert to radians
       if(phi .gt. 90.) then
         write(6,*) 'angle must be less than or equal to 90 degrees'
         go to 11
       endif
       print*
c

c angular momentum stuff
       L = mu*sqrt(G*M*a*(1-e*e))
       dAdt = 0.5d0*L/mu

c the orbital period
       P = sqrt(4.0d0*pi*pi * a*a*a/(G*M))
c       P = sqrt(aAU*aAU*aAU/(M/Msun)) * yr
c       Pyr = sqrt(aAU*aAU*aAU/(M/Msun))
c

c       i = pi/6.0d0  !radians (inclination)
c       sini = sin(i)
c       phi = pi/2.0d0  !radian (major axis relative to line of sight)
c

       open(unit=10, file='skew.dat', form='formatted',
     1  status='unknown')
c       N = 1001
       theta = 0.0d0
       t = 0.0d0
       dt = 1.3d0*P/float(N)
c
       write(10,90) m1sun
90     format('#', 'M1 = ',1p,e9.3,' solar masses')

       write(10,91) m2sun
91     format('#', 'M2 = ',1p,e9.3,' solar masses')

       write(10,92) vrcom, e, i*180.d0/pi, phi*180.d0/pi
92    format('#','vrcom, e, i, phi = ',1p,e9.3,2x,0p,f6.4,2x,2(f7.3,2x))

       write(10,99) aAU
99     format('#', 'system a = ', 1p,e15.6,' AU')

       write(10,100) P/yr
100    format('#', 'period = ', 1p,e15.6,' yrs')

       write(10,102) 
102    format('#')

       write(10,103) 
103    format('#', 9x, 'time(yrs)     v1r(km/s)      v2r(km/s)')
c
       do j = 1, N
         r = a*(1.0d0 - e*e)/(1.0d0 + e*cos(theta))
         v = sqrt(G*M*(2.0d0/r - 1.0d0/a))
         vr = v*sini*sin(theta+phi)
         v1r = (mu/m1)*vr + vrcom
         v2r = -(mu/m2)*vr + vrcom  !the opposite phase
         write(10,110) t/yr, v1r/1.0d5, v2r/1.0d5
         dtheta = 2.0d0*dAdt/(r*r) * dt
         theta = theta + dtheta
         t = t + dt
       end do
c
       
110    format(5x, 1p,3e15.6)
       write(6,*)' Output sent to file named skew.dat'
       print *
       stop
       end