program grader
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 [ orbgrades.5* --> orbgrades.6* ]
c
common pm,p,a,xi,xnode,t0,ecc,omega,ee,aee,cosi,xmu
character*1 decs
character*8 ref,ref2
character*10 wds,wds2
character*13 data2
character*20 infile, outfile, orbfile
character*12 star,star2
character*168 data1
dimension date(2000),theta(2000),rho(2000),wt(2000),
$ diff(2000),theta2(2000),phase(2000)
dimension gr1(2),gr2(2),gr3(2),gr4(2),gr5(2)
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
c
c write(6,931)
c 931 format(' Enter orbit file: ',$)
c read(5,932) orbfile
c 932 format(a20)
c write(6,933)
c 933 format(' Enter observations file: ',$)
c read(5,932) infile
c write(6,934)
c 934 format(' Enter output file: ',$)
c read(5,932) outfile
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
901 format(a168,f3.1,a13,t20,a10,1x,a12,t177,a8)
write(6,991) data1
991 format(' read(10) ',a168)
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,t63,f4.3)
c iflag=iflag+1
read(11,903) p,a,xi,xnode2,t0,ecc,omega2,ref,dlast
903 format(f15.6,f10.7,f8.5,f8.4,1x,f11.6,f8.6,f11.7,5x,a8,t81,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
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
npt=0
200 read(11,904,end=300) bess,tobs,robs,weight
904 format(t7,f9.4,t20,f6.2,t34,f9.6,t54,f4.1)
if (bess .eq. 0.0) go to 300
npt=npt+1
tobs=tobs + corr*(equ-bess)
theta(npt)=tobs
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.5 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.49999) grade1=5.49999
if (grade2 .gt. 5.49999) grade2=5.49999
if (grade3 .gt. 5.49999) grade3=5.49999
if (grade4 .gt. 5.49999) grade4=5.49999
if (grade5 .gt. 5.49999) grade5=5.49999
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
write(12,905) data1,grade,data2,wds
905 format(a168,f3.1,a13,' wds',a10,'.png')
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(2000),diff(2000)
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