program puborb2
c
c	Program to generate a set of separations for a given set of orbital 
c	elements.  points are determined at 360 equal time steps.
c
	common p,a,xi,xnode,t,ecc,omega
	dimension theta(360),rho(360)
	character*22 fname
c
c	Call "orbit" to read elements 
c
	write(6,901)
  901	format(/' Enter orbital elements in the following order:'/
     $     ' p, a, i, ln, t, e, lp'/)
	call orbit
c
c	Read name of output file
c
	write(6,902)
  902	format(/' Enter name of output data file:')
	read(5,903) fname
  903	format(a22)
	open(13,err=999,file=fname)
c
c	Call "porbit" for 360 different dates to determine theta and rho
c
	deltap=p/150.
	date=t
	do 300 n=1,360
	theta(n)=0.
	rho(n)=0.
  300	continue
c
	do 400 n=1,360
	if(n .eq. 1) date=t
	if(n .gt. 1) date=date+deltap*sqrt(rho(n-1)/a)
	if(date .le. t+p) go to 350
	nend=n-1
	go to 500
  350	call porbit(date,theta(n),rho(n))
  400	continue
c
c	Set remaining points equal to theta(nend) and rho(nend)
c
  500	do 600 n=nend+1,360
	theta(n)=theta(nend)
	rho(n)=rho(nend)
  600	continue
c
c	Write separations into output file and exit
c
 	write(13,904) (theta(n),rho(n),n=1,360)
  904	format(8(F6.1,F8.4))
c
  999	stop
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine orbit
	common p,a,xi,xnode,t,ecc,omega
	degrad=57.2957795
	read(5,*) p,a,xi,xnode,t,ecc,omega
	ee=sqrt((1.+ecc)/(1.-ecc))
	aee=a*(1.-ecc*ecc)
	omega=omega/degrad
	cosi=cos(xi/degrad)
	xnode=xnode/degrad
	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine porbit(date,pa,sep)
	common p,a,xi,xnode,t,ecc,omega
	degrad=57.2957795
	ee=sqrt((1.+ecc)/(1.-ecc))
	aee=a*(1.-ecc*ecc)
	cosi=cos(xi/degrad)
	xmu=6.2831853/p
	em=xmu*(date-t)
c
	e0=em+ecc*sin(em)+0.5*ecc*ecc*sin(2.*em)
	do 200 l=1,10
	em0=e0-ecc*sin(e0)
	e0=e0+(em-em0)/(1.-ecc*cos(e0))
  200	continue
c
	xnu=2.*atan(ee*tan(0.5*e0))
	r=aee/(1.+ecc*cos(xnu))
	pa=xnode+atan(cosi*tan(xnu+omega))
	sep=r*cos(xnu+omega)/cos(pa-xnode)
	pa=pa*degrad
	if(sep .gt. 0.) go to 300
	sep=-sep
	pa=pa+180.
  300	if(pa .lt.   0.) pa=pa+360.
	if(pa .gt. 360.) pa=pa-360.
c
  	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc