	program puborbs2
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	Note: this version reads a set of elements (orb6.master format)
c	for creating many files at once
c       240815 : do in JY rather than BY
c
	common p,a,xi,xnode,t,ecc,omega
	dimension theta(360),rho(360)
	character*10 wds
	character*1 acode,pcode,tcode
	character*8 ref
	character*22 fname
c
	open(10,file='puborbs.in',status='UNKNOWN')
c
c	Call "orbit" to read elements 
c
  100	read(10,901,end=999) wds,p,pcode,a,acode,xi,xnode,t,tcode,ecc,
     $     omega,ref
  901	format(t20,a10,t82,f11.6,a1,t106,f9.5,a1,t126,f8.4,t144,f8.4,
     $     t163,f12.6,a1,t188,f8.6,t206,f7.4,t252,a8)
        if (acode .eq. 'm') a=a/1000.
        if (pcode .eq. 'c') p=p*100.
        if (pcode .eq. 'd') p=p/365.25
        if (pcode .eq. 'h') p=p/365.25/24.
        if (pcode .eq. 'm') p=p/365.25/24./60.
        if (tcode .eq. 'd') t=(t-54544.5)/365.25 +2000.
        if (tcode .eq. 'm') t=(t+0.5-54544.5)/365.25 +2000.
	write(6,991) wds,p,a,xi,xnode,t,ecc,omega
  991	format(' read elements for ',a10,': ',7f12.5)
	degrad=57.2957795
	omega=omega/degrad
	xnode=xnode/degrad
c
c	Generate name of output file
c
	fname='wds'//wds//'.'//ref
	open(13,err=999,file=fname)
	write(6,993) fname
  993	format(' output file = ',a22)
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 loop back for next set of elements
c
 	write(13,904) (theta(n),rho(n),n=1,360)
  904	format(8(F6.1,F8.4))
	close(13)
	write(6,992)
  992	format(' thetas and rhos written out')
	go to 100
c
  999	stop
	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
