program quadflip6
c
c	the program reads a list of binary star data (separation and 
c	position angle) and a set of orbital elements, then determines
c	O-C residuals and flips theta values when necessary
c	             [ orbits.4 --> orbits.5 ]
c
c	120709: modified for new WDS format
c
c 	210512: modified to convert WDS rho/sep measures from deg/arcmin/mas
c           and clear rcode in orbits.5
c   230503: modified to keep rcode and print rcode + rho/sep measures in 
c           deg/arcmin/mas in orbits.5 since Proxima Cen overfill format in "
c   240815: convert from BY to JY

	common pm,p,a,xi,xnode,t0,ecc,omega
	character*1  decs,tcode,rcode
	character*5  apa,asep
	character*48 hcode
	character*20 infile, outfile
	character*100 star
c
	degrad=0.0174533
	pm=0.0
	equ=2000.0
c	enter file names for i/o
c
	open(10,file='orbits.4',status='UNKNOWN')
	open(11,file='orbits.5',status='UNKNOWN')
	open(12,file='quads.flipped',status='UNKNOWN')
c
c	read star name and orbital elements
c
  100	read(10,901,end=999,err=999) star,rah,ram,decs,decd,decm,pm
  901	format(a100,t1,f2.0,f3.1,a1,2f2.0,t63,f4.3)
	write(11,908) star
	write(12,908) star
  908	format(a100)
	call orbit(rx,ry,norb)
	ra=15.*rah + ram/4.
	dec=decd + decm/60.
	if(decs .eq. '-') dec=-dec
	corr=0.
	if(pm .ne. 0.0) corr = (0.0056*sin(ra*degrad)/cos(dec*degrad))
     $                       + 0.00417*pm*sin(dec*degrad)
c
c	read data points, determine residuals
c
  200	read(10,903,end=350) bess,tcode,pa,rcode,sep,hcode,apa,asep
  903	format(t7,f10.5,t19,a1,f7.3,t32,a1,f9.5,t50,a48,t20,a5,t35,a5)
	if (bess .eq. 0.0) go to 350
	if(apa .eq. ' ' .or. asep .eq. ' ') go to 200
	npt=npt+1
	pa2=pa + corr*(equ-bess)

	if (rcode .eq. 'D') then
		sep=sep*3600.
c		rcode=' '
	else if (rcode .eq. 'M') then
		sep=sep*60.
c		rcode=''
	else if (rcode .eq. 'm') then
		sep=sep/1000.
c		rcode=''
	end if

	call porbit(bess,pa,pa2,sep,xneg,yneg,tcode,rcode,hcode)
	go to 200
c
  350	write(11,907)
  907	format(98x,' ')
  	go to 100
  999	stop 'quadrant flip program finished'
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine orbit(x,y,norb)
	character*100 data
	character*8 refcode
        character*1 pcode,acode,tcode
	common pm,p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
	read(10,901) data,p,pcode,a,acode,xi,xnode,t0,tcode,ecc,
     $     omega,refcode
  901	format(a100,t3,f12.6,a1,f10.5,a1,2f9.4,f13.6,a1,f9.6,f9.4,2x,a8)
     	write(11,902) data
  902	format(a100)
        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 (acode .eq. 'M') a=a*60.
        if (acode .eq. 'm') a=a/1000.
        if (tcode .eq. 'c') t0=t0*100.
        if (tcode .eq. 'd') t0=(t0-51544.5)/365.25+2000.
        if (tcode .eq. 'm') t0=(t0+0.5-51544.5)/365.25+2000.
	ee=sqrt((1.+ecc)/(1.-ecc))
	aee=a*(1.-ecc*ecc)
	omega=omega/degrad
	cosi=cos(xi/degrad)
	xnode=xnode/degrad
	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine porbit(date,pa,pa2,sep,xneg,yneg,tcode,rcode,hcode)
	common pm,p,a,xi,xnode,t0,ecc,omega
	character*1  tcode,rcode
	character*48 hcode
	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-t0)
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))
	theta=xnode+atan(cosi*tan(xnu+omega))
	rho=r*cos(xnu+omega)/cos(theta-xnode)
	theta=theta*degrad
	if (rho .gt. 0.) go to 300
	rho=-rho
	theta=theta+180.
  300	if (theta .lt.   0.) theta=theta+360.
	if (theta .gt. 360.) theta=theta-360.
c
	xneg=-sin((theta+180.)/degrad)*rho
	yneg= cos((theta+180.)/degrad)*rho
	xobs=-sin((pa2+180.)/degrad)*sep
	yobs= cos((pa2+180.)/degrad)*sep
	resx=xobs-xneg
	resy=yobs-yneg
c
	respa=pa2-theta
	if (respa .lt. -180.) respa=respa+360.
	if (respa .gt.  180.) respa=respa-360.
	if (respa .lt. -90.)  pa=pa+180.
	if (respa .gt.  90.)  pa=pa-180.
	if (pa .lt.   0.) pa=pa+360.
	if (pa .gt. 360.) pa=pa-360.
	ressep=sep-rho
	perres=ressep/rho
	if (sep .eq. 0.0) go to 400
c
	if (rcode .eq. 'D') sep=sep/3600.
	if (rcode .eq. 'M') sep=sep/60.
	if (rcode .eq. 'm') sep=sep*1000.
	write(11,901) date,tcode,pa,rcode,sep,hcode
	if ((respa .lt. -90.) .or. (respa .gt.  90.))
     $     write(12,901) date,tcode,pa,rcode,sep,hcode
  901	format(6x,f10.5,2x,a1,f7.3,5x,a1,f9.5,8x,a48)
  400	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc