program grader4
c
c	The program reads a list of binary star data (separation   
c	and position angle) and a set of orbital elements (from 
c	neworb4.new), then determines O-C residuals and other 
c	statistics. Purpose is to determine orbit grade criterion.
c	             [ orbits.5 --> orbits.6 ]
c
c	120709: modified for new WDS format
c 
c 230504: modified by Bill H to account for separations other
c         than arcseconds following edits to quadflip5
c 240815: convert code from BY to JY
c

	common       pm,p,a,xi,xnode,t0,ecc,omega,ee,aee,cosi,xmu
	character*1  decs,pcode,acode,tcode,rcode,exten(100)
	character*8  ref,ref2
	character*10 wds,wds2,wdsold
	character*13 data2
        character*18 png
	character*20 infile, outfile, orbfile
	character*12 star,star2
	character*243 data1
	dimension date(3000),theta(3000),rho(3000),wt(3000),
     $            diff(3000),theta2(3000),phase(3000)
	dimension gr1(2),gr2(2),gr3(2),gr4(2),gr5(2)
	data exten/'a','b','c','d','e','f','g','h','i','j',
     $             'k','l','m','n','o','p','q','r','s','t',
     $             'u','v','w','x','y','z','A','B','C','D',
     $             'E','F','G','H','I','J','K','L','M','N',
     $             'O','P','Q','R','S','T','U','V','W','X',
     $             'Y','Z','1','2','3','4','5','6','7','8',
     $             '9','0','0','0','0','0','0','0','0','0',
     $             '0','0','0','0','0','0','0','0','0','0',
     $             '0','0','0','0','0','0','0','0','0','0',
     $             '0','0','0','0','0','0','0','0','0','0'/
	data gr1/3.3376434, -3.2616904/
	data gr2/0.9354897,  0.0164405/
	data gr3/0.6749010,  0.1105850/
	data gr4/0.7006989,  3122.5643/
	data gr5/0.6797400,    22943.2/
c
	degrad=0.0174533
	equ=2000.0
	wdsold='9999999999'
	iexten=1
c
	open(10,file='orbits.1',status='UNKNOWN')
	open(11,file='orbits.5',status='UNKNOWN')
	open(12,file='orbits.6',status='UNKNOWN')
c
c	read line from orbit file
c
  100	read(10,901,end=999,err=970) data1,gold,data2,wds2,star2,ref2,
     $     png
  901	format(a243,f3.1,a13,t20,a10,1x,a12,t252,a8,1x,a18)
	if (wds2 .ne. wdsold) iexten=1
	if (wds2 .eq. wdsold) iexten=iexten+1
	wdsold=wds2
	write(6,991) data1
  991	format(' read(10) ',a243)
c
c	read star name and orbital elements
c
  	read(11,902,end=999,err=980) wds,star,rah,ram,decs,decd,decm,pm
  902	format(a10,a12,t1,f2.0,f3.1,a1,2f2.0,t81,f4.3)
c	iflag=iflag+1
	read(11,903) p,pcode,a,acode,xi,xnode2,t0,tcode,ecc,omega2,
     $     ref,dlast
  903	format(f14.6,a1,f10.5,a1,2f9.4,f13.6,a1,f9.6,f9.4,
     $     2x,a8,t82,f4.0)
c
	if (wds  .ne. wds2)  go to 950
	if (star .ne. star2) go to 950
	if (ref  .ne. ref2)  go to 950
c	
        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/1000.
        if (acode .eq. 'M') a=a*60.
	if (acode .eq. 'D') a=a*3600.                                   
c *** added but doubt you will ever have a semimajor axis in degrees	
        if (tcode .eq. 'd') t0=(t0-51544.5)/365.25+2000.
        if (tcode .eq. 'm') t0=(t0+0.5-51544.5)/365.25+2000.
c
	ee=sqrt((1.+ecc)/(1.-ecc))
	aee=a*(1.-ecc*ecc)
	omega=omega2*degrad
	cosi=cos(xi*degrad)
	xnode=xnode2*degrad
	xmu=6.2831853/p
c
	ra=15.*rah + ram/4.
	dec=decd + decm/60.
	if(decs .eq. '-') dec=-dec
 	corr = (0.0056*sin(ra*degrad)/cos(dec*degrad))
     $     + 0.00417*pm*sin(dec*degrad)
c
c	read data points, determine residuals
c
c ***modified read & format to add rcode
c ***added if statements for rcode
c
	npt=0
  200	read(11,904,end=300) bess,tobs,rcode,robs,weight                
  904	format(t7,f10.5,t20,f7.3,t32,a1,f9.5,t53,f6.2)                  
	if (bess .eq. 0.0) go to 300
	npt=npt+1
	tobs=tobs + corr*(equ-bess)
	theta(npt)=tobs
	if (rcode .eq. 'm') robs=robs/1000.
	if (rcode .eq. 'M') robs=robs*60.    
	if (rcode .eq. 'D') robs=robs*3600.  
	rho(npt)=robs
	wt(npt)=weight
	date(npt)=bess
	call porbit(bess,tobs,robs,tres,rres,rresr,xyres)
	diff(npt)=xyres
	go to 200
c
c	determine weighted rms residuals, throw out outliers and recompute
c
  300	call rms(npt,wt,diff,sumwt,rmsdiff)
  	do 320 n=1,npt
  320	if (diff(n) .gt. 3.0*rmsdiff) wt(n)=0.
  	call rms(npt,wt,diff,sumwt2,rmsdiff2)
c
c	simple one - determine number of revolutions covered
c	up to date of orbit publication
c
	if (date(npt) .lt. dlast) dlast=date(npt)
	enrev=(dlast-date(1))/p
c
c	theta and phase coverage: extract theta and phase
c	values for observations up to publication date
c
	npt2=0
	do 350 n=1,npt
	if (wt(n) .eq. 0.) go to 350
	if (date(n) .gt. dlast) go to 350
	npt2=npt2+1
	theta2(npt2)=theta(n)
	phase(npt2)=(date(n)-t0)/p + 1000.
	phase(npt2)=phase(npt2)-float(int(phase(npt2)))
  350	continue
c
c	sort on theta, determine rms difference in
c	theta between successive observations
c	also determine maximum gap in theta
c
	do 400 n=1,npt2-1
	do 400 l=n+1,npt2
	if (theta2(l) .ge. theta2(n)) go to 400
	temp=theta2(l)
	theta2(l)=theta2(n)
	theta2(n)=temp
  400	continue
c		
	dth=0.
	dthmax=0.
	do 420 n=1,npt2-1
	tdiff=abs(theta2(n+1)-theta2(n))
  	dth=dth+tdiff*tdiff
	if (tdiff .gt. dthmax) dthmax=tdiff
  420	continue
	tdiff=abs(theta2(1)+360.-theta2(npt2))
	dth=dth+tdiff*tdiff
	if (tdiff .gt. dthmax) dthmax=tdiff
	dth=sqrt(dth/float(npt2))
c
c	sort on phase, determine rms difference in
c	phase between successive observations
c	also determine maximum gap in phase
c
	do 500 n=1,npt2-1
	do 500 l=n+1,npt2
	if (phase(l) .ge. phase(n)) go to 500
	temp=phase(l)
	phase(l)=phase(n)
	phase(n)=temp
  500	continue
c		
	dph=0.
	dphmax=0.
	do 520 n=1,npt2-1
	pdiff=abs(phase(n+1)-phase(n))
  	dph=dph+pdiff*pdiff
	if (pdiff .gt. dphmax) dphmax=pdiff
  520	continue
	pdiff=abs(phase(1)+1.-phase(npt2))
	if (pdiff .gt. dphmax) dphmax=pdiff
	dph=dph+pdiff*pdiff
	dph=sqrt(dph/float(npt2))
	dph=dph*360.
	dphmax=dphmax*360.
c
c	determine actual grades
c
	xx1=alog10(enrev)
	xx2=(dthmax+dphmax)/2.
	xx3=(dth+dph)/2.
	xx4=rmsdiff2/float(npt2)
	xx5=rmsdiff2/a * rmsdiff2/float(npt2)
c
	grade1= gr1(1) + gr1(2)*xx1
	grade2= gr2(1) + gr2(2)*xx2 
	grade3= gr3(1) + gr3(2)*xx3 
	grade4= gr4(1) + gr4(2)*xx4 
	grade5= gr5(1) + gr5(2)*xx5 
c
c	set range of possible grades from 0.5 to 5.4 so bin 
c 	sizes of integer grades (1 - 5) will be uniform
c
	if (grade1 .lt. 0.50001) grade1=0.50001
	if (grade2 .lt. 0.50001) grade2=0.50001
	if (grade3 .lt. 0.50001) grade3=0.50001
	if (grade4 .lt. 0.50001) grade4=0.50001
	if (grade5 .lt. 0.50001) grade5=0.50001
	if (grade1 .gt. 5.4)     grade1=5.4
	if (grade2 .gt. 5.4)     grade2=5.4
	if (grade3 .gt. 5.4)     grade3=5.4
	if (grade4 .gt. 5.4)     grade4=5.4
	if (grade5 .gt. 5.4)     grade5=5.4
	grade = (grade1+grade2+grade4)/3.
	newgrade=nint(grade)
c
c	downgrade orbit a bit if scatter in separation too large
c
	if ((newgrade .eq. 1) .and. (rmsdiff2/a .gt. 0.2)) 
     $     grade = grade + 1. 
	if ((newgrade .eq. 2) .and. (rmsdiff2/a .gt. 0.3)) 
     $     grade = grade + 1. 
	if ((newgrade .eq. 3) .and. (rmsdiff2/a .gt. 0.5)) 
     $     grade = grade + 1. 
	if ((newgrade .eq. 4) .and. (rmsdiff2/a .gt. 0.8)) 
     $     grade = grade + 1. 
c
        if (png .eq. '                  ') 
     $     png='wds'//wds//exten(iexten)//'.png'
c
	write(12,905) data1,grade,data2,png
  905	format(a243,f3.1,a13,1x,a18)
	if (abs(grade-gold) .gt. 0.1) write(6,907) wds,data2,gold,grade
  907	format(a10,a13,2f5.1)
c
  	go to 100
  950	write(6,906) wds,star,wds2,star2,ref,ref2
  906	format(' Orbit mismatch ? ',a10,a12,2x,a10,a12,2x,a8,2x,a8)
	go to 999
  970	write(6,992)
  992	format(' error reading from input file 10')
	go to 999
  980	write(6,993)
  993	format(' error reading from input file 11')
	go to 999
c
  999	stop
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine porbit(bess,tobs,robs,tres,rres,rresr,xyres)
	common pm,p,a,xi,xnode,t0,ecc,omega,ee,aee,cosi,xmu
	degrad=57.2957795
	em=xmu*(bess-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))
	tcalc=xnode+atan(cosi*tan(xnu+omega))
	rcalc=r*cos(xnu+omega)/cos(tcalc-xnode)
	tcalc=tcalc*degrad
	if(rcalc .gt. 0.) go to 300
	rcalc=-rcalc
	tcalc=tcalc+180.
  300	if(tcalc .lt.   0.) tcalc=tcalc+360.
	if(tcalc .gt. 360.) tcalc=tcalc-360.
c
	xcalc=-sin((tcalc+180.)/degrad)*rcalc
	ycalc= cos((tcalc+180.)/degrad)*rcalc
	xobs =-sin((tobs+180.)/degrad)*robs
	yobs = cos((tobs+180.)/degrad)*robs
	xres = xobs-xcalc
	yres = yobs-ycalc
	xyres= sqrt(xres*xres+yres*yres)
c
	tres=tobs-tcalc
	if(tres .lt. -180.) tres=tres+360.
	if(tres .gt.  180.) tres=tres-360.
	rres=robs-rcalc
	rresr=rres/rcalc
c
  400	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine rms(npt,wt,diff,sumwt,rmsdiff)
	dimension wt(3000),diff(3000)
	sumwt=0.
	sumdiff=0.
	do 100 n=1,npt
	sumwt=sumwt+wt(n)
	sumdiff=sumdiff+wt(n)*diff(n)*diff(n)
  100	continue
  	rmsdiff=sqrt(sumdiff/sumwt)
  	return
  	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc