program quadflip3
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	             [ orbgrades.4* --> orbgrades.5* ]
c
	common pm,p,a,xi,xnode,t0,ecc,omega
	character*1  decs
	character*3  code1
	character*5  apa,asep
	character*4 code2,code3
	character*31 code4
	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,pa,sep,code1,code2,code3,
     $     code4,apa,asep
  903	format(t7,f9.4,t20,f6.2,t33,f10.6,t50,a3,1x,a4,5x,a4,
     $     t50,a31,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)
	call porbit(bess,pa,pa2,sep,xneg,yneg,code1,code2,
     $     code3,code4)
	go to 200
c
  350	write(11,907)
  907	format(80x,' ')
  	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.2421987
	if (pcode .eq. 'h') p=p/365.2421987/24.
	if (pcode .eq. 'm') p=p/365.2421987/24./60.
        if (acode .eq. 'm') a=a/1000.
        if (tcode .eq. 'd') t0=(t0-33282.423)/365.2421987+1950.
	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,pa2,sep,xneg,yneg,code1,code2,
     $     code3,code4)
	common pm,p,a,xi,xnode,t0,ecc,omega
	character*3  code1
	character*4 code2,code3
	character*31 code4
	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
	write(11,901) date,pa,sep,code4
	if (respa .lt. -90.)  write(12,901) date,pa,sep,code4
	if (respa .gt.  90.)  write(12,901) date,pa,sep,code4
  901	format(6x,f9.4,4x,f6.2,7x,f10.6,7x,a31)
  400	return
  902	format(1x,f10.4,2x,a2,1x,a12,18x,f10.2,f8.3)
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc