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 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.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 (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-33282.423)/365.2421987+1950. if (tcode .eq. 'm') t0=(t0+0.5-33282.423)/365.2421987+1950. 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