program orbgrid13 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c This program calculates Thiele-Innes elements and their errors c from a set of weighted visual measures. It does a grid search c over a range of P, E, and T0, calculating values of (A,B,F,G) by c least squares. After determining a (P,E,T0) set giving minimum c errors, the grid spacing is reduced and the process repeated. c The first grid search ends when all grid sizes decrease below c 0.01. Weights are then recalculated based on residuals and the c grid search restarted. The final grid search ends when all grid c sizes decrease below 0.0001 . c c Version of orbgrid7 which reads WDS-style title line. c Also reads pm from title line c c 040924: removed old normal point routines, cleaned up. At this c point will compile properly and all errors are essentially the c same as before except in a, which for example file is about a c factor of 10 too high. c c 060602: changed input to take name of form 'wds12345+1234' and c append '.dat', .orb', and '.plt' for input and output filenames. c c 120711: corrected format so reads pm correctly from wds.summ c 150320: convert to double precision c 190830: "singular matrix" problems. Calls to some sections of c subroutine kepler defined terms, while calls to other sections c then used them in populating matrix; apparently updated compiler c no longer saved values in those terms when exiting subroutine, c so matrix terms became zero, infinity, or undefined. Portions of c subroutine no longer used were removed, Tokovinin terms renamed c to those used in main program, new common block "elem" added to c save terms. c 240812 : correct to give results in JY. Only change lines 297 & 298? c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c implicit double precision (a-h,o-z) dimension ss(10),xma(2000),error(9,9,9),theta0(2000),elerr(7), $ elfix(7),obs(2000,6),nnum(500),tnmin(500),tnmax(500),t(2000), $ theta(2000),rho(2000),xo(2000),yo(2000),w(2000),xow(2000), $ x(2000),y(2000),dx(2000),dy(2000),dr(2000),dpa(2000), $ rdt(2000),yow(2000) character*1 decs character*3 kode1(2000),xcode character*6 atheta,arho character*13 filein character*18 kode2(2000),ycode character*20 file01,file02,file03 character*130 name common /numcom/qt,qth,qrho,qwt,tnmin,tnmax,nnum,t,theta,rho,xo, $ yo,w,xow,x,y,dx,dy,dr,dpa,rdt,yow,sumwt,wbar,resx,resy,res common /elem/p,t0,ecc,semi,asc,per,xinc,casc,sasc,cper,sper, $ cxinc,sxinc,txinc degrad=0.0174532925D0 equ=2000.D0 c c Enter file names for i/o c write(6,901) 901 format(/' Enter input data file name: ',$) read(5,'(a13)') filein file01=filein//'.dat' file02=filein//'.orb' file03=filein//'.plt' write(6,'(//)') c open(1,file=file01,status='UNKNOWN') open(2,file=file02,status='UNKNOWN') open(3,file=file03,status='UNKNOWN') c c Read star name, position, and proper motion; determine proper c motion correction factor in theta. Also, read in initial c guesses at period, eccentricity, and T0, as well as step sizes. c read(1,904) name,rah,ram,decs,decd,decm,pm 904 format(a130,t1,f2.0,f3.1,a1,2f2.0,t81,f4.3) read(1,*) p0,t00,ecc0,pstep0,tstep0,estep0 write(2,'(a130)') name write(6,'(a130)') name ra=rah*15.D0 + ram/4.D0 dec=decd + decm/60.D0 if (decs .eq. '-') dec=-dec write(2,905) ra,dec,pm,equ write(6,905) ra,dec,pm,equ 905 format(/' RA,Dec,PM,Equinox = ',2f12.3,f12.6,f10.0/) c c Correct for precession and proper motion. Here pm is the RA c proper motion term, in arcsec per year corr = (0.00557D0*dsin(degrad*ra)/dcos(degrad*dec)) $ +0.000278D0*pm*dsin(degrad*dec) ipass=1 nloop=0 nloop2=0 nobs=0 c c Read in data c 100 read(1,906,end=400) qt,qth0,qrho,xcode,qwt,ycode,atheta,arho 906 format(6x,f10.5,3x,f7.3,6x,f9.5,8x,a3,f6.2,3x,a18, $ t20,a6,t35,a6) c write(6,993) qt c 993 format(f10.4) if (qt .eq. 0.0) go to 400 if (qt .gt. 2030. .or. qt .lt. 1500.) go to 900 if (atheta .eq. ' ' .or. arho .eq. ' ') go to 100 c c Store all data into arrays for plot input file c qth = qth0 + corr*(equ - qt) if (qth .lt. 0.D0) qth=qth+360.D0 if (qth .gt. 360.D0) qth=qth-360.D0 if (qth .lt. 0.01D0) qth= 0.01D0 if (qth .gt. 359.99D0) qth=359.99D0 c 200 nobs=nobs+1 t(nobs)=qt theta0(nobs)=qth0 theta(nobs)=qth rho(nobs)=qrho kode1(nobs)=xcode kode2(nobs)=ycode w(nobs)=qwt xo(nobs) = rho(nobs)*dcos(degrad*theta(nobs)) yo(nobs) = rho(nobs)*dsin(degrad*theta(nobs)) go to 100 c 400 if (ipass .eq. 1) write(2,907) (t(j),theta(j),rho(j),xo(j), $ yo(j),kode1(j),kode2(j),w(j),j=1,nobs) if (ipass .eq. 1) write(6,907) (t(j),theta(j),rho(j),xo(j), $ yo(j),kode1(j),kode2(j),w(j),j=1,nobs) 907 format(2x,f11.5,f9.3,f9.5,2f9.4,2x,a3,1x,a18,f7.2) write(2,'(//)') write(6,'(//)') c c Define sum of weights and weighted x and y values for use in grid c 420 sumwt=0.D0 do 450 j=1,nobs sumwt=sumwt+w(j) xow(j)=xo(j)*w(j) 450 yow(j)=yo(j)*w(j) wbar=sumwt/real(nobs) c pstep=pstep0 tstep=tstep0 estep=estep0 c c All step sizes must all go below minimum value "sizmin" before c the program hops out of the loop. Sizmin is set to 0.01 for the c first pass. For the second pass, weights are recalculated, the c initial step sizes are cut in half, and sizmin is reset to c 0.0001 . Sizmin for eccentricity is 10 times more stringent c than for P and T0. c if (ipass .eq. 1) sizmin=0.01D0 if (ipass .ge. 2) sizmin=0.0001D0 if (ipass .eq. 2) pstep=pstep/2.D0 if (ipass .eq. 2) tstep=tstep/2.D0 if (ipass .eq. 2) estep=estep/2.D0 ipass=ipass+1 c 500 do 550 ip=1,9 do 550 it=1,9 do 550 ie=1,9 550 error(ip,it,ie)=0. c ecc02=ecc0 eccx=ecc02+estep*4.D0 if (eccx .ge. 0.9999D0) ecc02=0.9999D0-estep*4.D0 do 720 ip=1,9 p=p0+pstep*real(ip-5) xmu = 360.0D0/p if (ip .eq. 1 .or. pstep .ne. 0.0) go to 580 do 570 iit=1,9 do 570 iie=1,9 570 error(ip,iit,iie)=error(1,iit,iie) go to 720 c 580 do 710 it=1,9 t0=t00+tstep*real(it-5) do 600 j=1,nobs xma(j)=xmu*(t(j)-t0) 600 continue if (it .eq. 1 .or. tstep .ne. 0.0) go to 650 do 620 iip=1,9 do 620 iie=1,9 620 error(iip,it,iie)=error(iip,1,iie) go to 710 c 650 do 700 ie=1,9 ecc=ecc02+estep*real(ie-5) if (ie .eq. 1 .or. estep .ne. 0.0) go to 680 do 660 iip=1,9 do 660 iit=1,9 660 error(iip,iit,ie)=error(iip,iit,1) go to 700 c c Calculate time dependent parameters and solve for c Thiele-Innes elements. Calculate errors c 680 call thin(xma,nobs,a,b,f,g,ss) if (estep+pstep+tstep .eq. 0.0) $ call elements(a,b,f,g,pm) call errcalc(nobs,a,b,f,g,sigx,sigy,sigv1,sigv2, $ sigs1,sigs2,errtot,1,kode1) error(ip,it,ie)=errtot c 700 continue 710 continue 720 continue c c Interpolate error array to determine minimum c call interp(error,pstep,tstep,estep,p0,t00,ecc0,p,t0,ecc) xmu=360.0D0/p do 730 j=1,nobs xma(j)=xmu*(t(j)-t0) 730 continue call thin(xma,nobs,a,b,f,g,ss) call errcalc(nobs,a,b,f,g,sigx,sigy,sigv1,sigv2, $ sigs1,sigs2,errmin,2,kode1) c write(2,908) p,t0,ecc,errmin write(6,908) p,t0,ecc,errmin 908 format(' P,T0,ECC,ERR = ',25x,f11.5,f12.5,f9.5,f10.3) nloop=nloop+1 if (nloop .gt. 50) write(2,909) if (nloop .gt. 50) write(6,909) 909 format(/' Number of loops exceeded maximum allowed') c c if step sizes are sufficiently small or number of loops c exceeds maximum, go calculate elements c if (nloop .gt. 50) go to 800 if (max(pstep,tstep) .lt. sizmin .and. estep .lt. sizmin/10.) $ go to 800 c c Changes in P, T0, and ECC are checked for possible modification c of step sizes, as follows: c c (1) If changes exceed the step size for any one of the three c grid dimensions, the same values are used for the next grid c search. c (2) If the same step sizes have been used 5 times, it's figured c a larger step size is needed. All step sizes are increased c by a factor of 1.5 . c (3) If none of the changes exceed the step sizes, all values are c decreased by a factor of 0.3 . c sfact=0.3 if (abs(p-p0) .gt. pstep) sfact=1.0D0 if (abs(t0-t00) .gt. tstep) sfact=1.0D0 if (abs(ecc-ecc0) .gt. estep) sfact=1.0D0 nloop2=nloop2+1 if (sfact .ne. 1.0D0) nloop2=0 if (nloop2 .eq. 5) sfact=1.5D0 if (nloop2 .eq. 5) nloop2=0 c pstep=pstep*sfact tstep=tstep*sfact estep=estep*sfact c write(2,'(" Step sizes: ",2f10.5,f8.5)') pstep,tstep,estep write(6,'(" Step sizes: ",2f10.5,f8.5)') pstep,tstep,estep call elements(a,b,f,g,pm) p0=p t00=t0 ecc0=ecc go to 500 c c Call toko to calculate errors c 800 do 810 n=1,nobs obs(n,1)=dble(t(n)) obs(n,2)=dble(theta(n)) obs(n,3)=dble(rho(n)) if (w(n) .ne. 0.) obs(n,4)=dble(1./w(n)) c 2007 09 27: changed value below from 9.99 to 9999.99 for c zero-weighted data if (w(n) .eq. 0.) obs(n,4)=dble(9999.99D0) obs(n,5)=dble(dpa(n)) obs(n,6)=dble(dr(n)) 810 continue c do 820 k=1,7 820 elfix(k)=1.D0 if (pstep0 .eq. 0.) elfix(1)=0.D0 if (tstep0 .eq. 0.) elfix(2)=0.D0 if (estep0 .eq. 0.) elfix(3)=0.D0 call toko(nobs,elerr,elfix,obs) pday = p*365.25 errpday = elerr(1)*365.25 c a3p2=semi*semi*semi/(p*p) write(2,910) resx,resy,res,sigv1,sigv2,sigs1,sigs2 write(6,910) resx,resy,res,sigv1,sigv2,sigs1,sigs2 910 format(/' Weighted RMS Residuals (X, Y, X+Y) = ',3f10.4, $ /' RMS Visual Residuals (dR, RdT) = ',2f10.4, $ /' RMS Speckle Residuals (dR, RdT) = ',2f10.4) write(2,'(//10x,"Orbital Elements for ",a23//)') name write(6,'(//10x,"Orbital Elements for ",a23//)') name write(2,911) a,b,f,g,pday,errpday,p,elerr(1),semi,elerr(4), $ xinc,elerr(7),asc,elerr(5),t0,elerr(2), $ ecc,elerr(3),per,elerr(6),a3p2 write(6,911) a,b,f,g,pday,errpday,p,elerr(1),semi,elerr(4), $ xinc,elerr(7),asc,elerr(5),t0,elerr(2), $ ecc,elerr(3),per,elerr(6),a3p2 911 format(10x,'A = ',f14.7,' B = ',f14.7, $ /10x,'F = ',f14.7,' G = ',f14.7, $ //10x,'Period (days) = ',f14.7,' +/- ',f12.7, $ //10x,'Period (years) = ',f14.7,' +/- ',f12.7, $ /10x,'Semi-major Axis = ',f14.7,' +/- ',f12.7, $ /10x,'Inclination = ',f14.7,' +/- ',f12.7, $ /10x,'Longitude Node = ',f14.7,' +/- ',f12.7, $ /10x,'T0 = ',f14.7,' +/- ',f12.7, $ /10x,'Eccentricity = ',f14.7,' +/- ',f12.7, $ /10x,'Longitude per = ',f14.7,' +/- ',f12.7/ $ /10x,'a^3/P^2 = ',e14.3,/) c if (ipass .ge. 3) go to 870 c c recalculate weights c write(2,912) write(6,912) 912 format(///' Weights recalculated from mean errors; grid search', $ ' repeated'//) do 860 j=1,nobs if (w(j) .eq. 0.0) go to 860 if (kode1(j) .eq. 'VIS') go to 830 if (kode1(j) .eq. 'EYE') go to 830 if (sigs1 .ne. 0.0) sratio1=abs(dr(j)/sigs1) if (sigs1 .eq. 0.0) sratio1=1. if (sigs2 .ne. 0.0) sratio2=abs(rdt(j)/sigs2) if (sigs2 .eq. 0.0) sratio2=1. sratio=max(sratio1,sratio2) if (sigv1 .ne. 0.0) svratio1=abs(dr(j)/sigv1) if (sigv1 .eq. 0.0) svratio1=1. if (sigv2 .ne. 0.0) svratio2=abs(rdt(j)/sigv2) if (sigv2 .eq. 0.0) svratio2=1. svratio=max(svratio1,svratio2) wfact=1. if (sratio .gt. 3.0) wfact=0.25 if (svratio .gt. 3.0) wfact=0. go to 840 830 if (sigv1 .ne. 0.0) vratio1=abs(dr(j)/sigv1) if (sigv1 .eq. 0.0) vratio1=1. if (sigv2 .ne. 0.0) vratio2=abs(rdt(j)/sigv2) if (sigv2 .eq. 0.0) vratio2=1. vratio=max(vratio1,vratio2) wfact=1. if (vratio .gt. 3.0) wfact=0. 840 w(j)=w(j)*wfact if (wfact .ne. 1.) write(2,913) t(j),w(j) if (wfact .ne. 1.) write(6,913) t(j),w(j) 913 format(' Weight for ',f11.4,' reduced to',f6.1) 860 continue c c Restart grid search using new weights and new minimum step size c nloop=0 go to 400 c 870 write(2,914) write(6,914) 914 format(//' Residuals from the final orbit are as follows:'// $ ' Date Observed Calculated dT ', $ ' dR dR RdT Wt Code'/) do 880 j=1,nobs xc=xo(j)-dx(j) yc=yo(j)-dy(j) sep=sqrt(xc*xc+yc*yc) pa=datan(yc/xc)/degrad if (xc .lt. 0.) pa=pa+180. if (xc .gt. 0. .and. yc .lt. 0.) pa=pa+360. idsep=nint(100.*dr(j)/sep) if (idsep .gt. 999) idsep= 999 if (idsep .lt. -999) idsep=-999 dpa(j)=theta(j)-pa if (dpa(j) .gt. 180.) dpa(j)=dpa(j)-360. if (dpa(j) .lt. -180.) dpa(j)=dpa(j)+360. write(2,915) t(j),theta(j),rho(j),pa,sep,dpa(j),dr(j),idsep, $ rdt(j),w(j),kode1(j) write(6,915) t(j),theta(j),rho(j),pa,sep,dpa(j),dr(j),idsep, $ rdt(j),w(j),kode1(j) 915 format(f10.4,f8.2,f8.4,f9.2,f8.4,f8.2,f8.4,i5,f8.3,f7.2,2x,a3) 880 continue c c Write input file for plotting c write(3,'(a130/"YES")') name write(3,916) p,semi,xinc,asc,t0,ecc,per 916 format(f13.6,'y',f10.5,'a',2f9.4,f13.6,'y',f9.6,f9.4,2x, $ 'new orb ') c do 890 n=1,nobs 890 write(3,917) t(n),theta0(n),rho(n),kode1(n),w(n),kode2(n) 917 format(6x,f9.4,4x,f6.2,7x,f8.4,9x,a3,f5.1,4x,a18) stop 900 write(2,'(" Input data error: qt, nobs = ",f12.4,i8)') qt,nobs write(6,'(" Input data error: qt, nobs = ",f12.4,i8)') qt,nobs close(1) close(2) close(3) stop end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine kepler1(am,ae) implicit double precision (a-h,o-z) common /elem/p,t0,ecc,semi,asc,per,xinc,casc,sasc,cper,sper, $ cxinc,sxinc,txinc degrad = 0.0174532925D0 if (ecc .lt. 0.0001D0) ecc=0.0001D0 if (ecc .gt. 0.9999D0) ecc=0.9999D0 eo = am + (ecc*dsin(degrad*am))/degrad 100 e1 = am + (ecc*dsin(degrad*eo))/degrad de = abs(e1 - eo) if (de .le. 0.01) go to 200 eo = e1 go to 100 200 ae = e1 return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine thin(xma,nobs,a,b,f,g,ss) implicit double precision (a-h,o-z) dimension xma(2000),ss(10),ea(2000),tnmin(500),tnmax(500), $ t(2000),theta(2000),rho(2000),xo(2000),yo(2000),w(2000), $ xow(2000),x(2000),y(2000),dx(2000),dy(2000),dr(2000), $ dpa(2000),rdt(2000),yow(2000),nnum(500) common /numcom/qt,qth,qrho,qwt,tnmin,tnmax,nnum,t,theta,rho,xo, $ yo,w,xow,x,y,dx,dy,dr,dpa,rdt,yow,sumwt,wbar,resx,resy,res common /elem/p,t0,ecc,semi,asc,per,xinc,casc,sasc,cper,sper, $ cxinc,sxinc,txinc degrad=0.0174532925D0 c if (ecc .lt. 0.0001D0) ecc=0.0001D0 if (ecc .gt. 0.9999D0) ecc=0.9999D0 f1 = sqrt(1.0D0-(ecc*ecc)) c do 100 j=1,nobs call kepler1(xma(j),ea(j)) x(j) = dcos(degrad*ea(j)) - ecc y(j) = f1*dsin(degrad*ea(j)) 100 continue c do 200 i=1,10 200 ss(i)=0.0D0 c do 300 j=1,nobs xw=x(j)*w(j) yw=y(j)*w(j) ss(1) = ss(1) + x(j)*xw ss(2) = ss(2) + y(j)*xw ss(3) = ss(3) + y(j)*yw ss(4) = ss(4) + x(j)*xow(j) ss(5) = ss(5) + y(j)*xow(j) ss(6) = ss(6) + x(j)*yow(j) ss(7) = ss(7) + y(j)*yow(j) ss(9) = ss(9) + xw ss(10)= ss(10)+ yw 300 continue c c Solve for Thiele-Innes elements c denom = ss(2)*ss(2) - ss(1)*ss(3) f = (ss(2)*ss(4) - ss(1)*ss(5))/denom a = (ss(2)*ss(5) - ss(3)*ss(4))/denom g = (ss(2)*ss(6) - ss(1)*ss(7))/denom b = (ss(2)*ss(7) - ss(3)*ss(6))/denom c return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine interp(error,pstep,tstep,estep,p0,t00,ecc0, $ p,t0,ecc) implicit double precision (a-h,o-z) dimension error(9,9,9),tnmin(500),tnmax(500),t(2000),xow(2000), $ rho(2000),xo(2000),yo(2000),theta(2000),rdt(2000),yow(2000), $ x(2000),y(2000),dx(2000),dy(2000),dr(2000),dpa(2000),w(2000), $ nnum(500) common /numcom/qt,qth,qrho,qwt,tnmin,tnmax,nnum,t,theta,rho,xo, $ yo,w,xow,x,y,dx,dy,dr,dpa,rdt,yow,sumwt,wbar,resx,resy,res degrad=0.0174532925D0 c errmin=error(1,1,1) ipm=1 itm=1 iem=1 c do 100 ip=1,9 do 100 it=1,9 do 100 ie=1,9 if (error(ip,it,ie) .gt. errmin) go to 100 errmin=error(ip,it,ie) ipm=ip itm=it iem=ie 100 continue c pm=real(ipm) if (ipm .eq. 1 .or. ipm .eq. 9) go to 200 y1=error(ipm-1,itm,iem) y2=error(ipm, itm,iem) y3=error(ipm+1,itm,iem) a=(y3-y2-y2+y1)/2. b=y2-y1-a*(2.*pm-1.) pm=(-0.5)*b/a if (y1/y2 .gt. 10. .or. y3/y2 .gt. 10.) pm=real(ipm) c 200 tm=real(itm) if (itm .eq. 1 .or. itm .eq. 9) go to 300 y1=error(ipm,itm-1,iem) y2=error(ipm,itm ,iem) y3=error(ipm,itm+1,iem) a=(y3-y2-y2+y1)/2.D0 b=y2-y1-a*(2.D0*tm-1.D0) tm=(-0.5D0)*b/a if (y1/y2 .gt. 10.D0 .or. y3/y2 .gt. 10.D0) tm=real(itm) c 300 em=real(iem) if (iem .eq. 1 .or. iem .eq. 9) go to 400 y1=error(ipm,itm,iem-1) y2=error(ipm,itm,iem) y3=error(ipm,itm,iem+1) a=(y3-y2-y2+y1)/2.D0 b=y2-y1-a*(2.D0*em-1.D0) em=(-0.5)*b/a if (y1/y2 .gt. 10.D0 .or. y3/y2 .gt. 10.D0) em=real(iem) c 400 p=p0+pstep*(pm-5.) t0=t00+tstep*(tm-5.) ecc=ecc0+estep*(em-5.) c if (ecc .lt. 0.0001D0) ecc=0.0001D0 if (ecc .gt. 0.9999D0) ecc=0.9999D0 c return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine errcalc(nobs,a,b,f,g,sigx,sigy,sigv1,sigv2, $ sigs1,sigs2,errtot,icode,kode1) implicit double precision (a-h,o-z) character*3 kode1(2000) dimension tnmin(500),tnmax(500),t(2000),theta(2000),rho(2000), $ xo(2000),yo(2000),w(2000),xow(2000),x(2000),y(2000),dx(2000), $ dy(2000),dr(2000),dpa(2000),rdt(2000),yow(2000),nnum(500) common /numcom/qt,qth,qrho,qwt,tnmin,tnmax,nnum,t,theta,rho,xo, $ yo,w,xow,x,y,dx,dy,dr,dpa,rdt,yow,sumwt,wbar,resx,resy,res degrad=0.0174532925D0 c sigx=0.D0 sigy=0.D0 c do 200 j=1,nobs xc=a*x(j) + f*y(j) dx(j)=xo(j) - xc sigx=sigx + dx(j)*dx(j)*w(j) yc=b*x(j) + g*y(j) dy(j)=yo(j) - yc sigy=sigy + dy(j)*dy(j)*w(j) 200 continue c resx = sqrt(sigx/sumwt) resy = sqrt(sigy/sumwt) res = sqrt((sigx+sigy)/sumwt) errtot=res*10000.D0 if (icode .eq. 1) return c sigv1=0.D0 sigv2=0.D0 sigs1=0.D0 sigs2=0.D0 env=0.D0 ens=0.D0 do 400 j=1,nobs xc=a*x(j) + f*y(j) yc=b*x(j) + g*y(j) sep=sqrt(xc*xc+yc*yc) pa=datan(yc/xc)/degrad if (xc .lt. 0.D0) pa=pa+180.D0 if (xc .gt. 0.D0 .and. yc .lt. 0.D0) pa=pa+360.D0 dr(j)=rho(j)-sep dpa(j)=theta(j)-pa if (dpa(j) .gt. 180.D0) dpa(j)=dpa(j)-360.D0 if (dpa(j) .lt. -180.D0) dpa(j)=dpa(j)+360.D0 rdt(j)=sep*dpa(j)*degrad c if (w(j) .eq. 0.) go to 400 if (kode1(j) .eq. 'VIS') go to 300 if (kode1(j) .eq. 'XVI') go to 300 ens=ens+1. sigs1=sigs1+dr(j)*dr(j) sigs2=sigs2+rdt(j)*rdt(j) go to 400 300 env=env+1. sigv1=sigv1+dr(j)*dr(j) sigv2=sigv2+rdt(j)*rdt(j) 400 continue if (env .ne. 0.) sigv1=sqrt(sigv1/env) if (env .ne. 0.) sigv2=sqrt(sigv2/env) if (ens .ne. 0.) sigs1=sqrt(sigs1/ens) if (ens .ne. 0.) sigs2=sqrt(sigs2/ens) return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine elements(a,b,f,g,pm) implicit double precision (a-h,o-z) dimension tnmin(500),tnmax(500),t(2000),theta(2000),rho(2000), $ xo(2000),yo(2000),w(2000),xow(2000),x(2000),y(2000),dx(2000), $ dy(2000),dr(2000),dpa(2000),rdt(2000),yow(2000),nnum(500) common /numcom/qt,qth,qrho,qwt,tnmin,tnmax,nnum,t,theta,rho,xo, $ yo,w,xow,x,y,dx,dy,dr,dpa,rdt,yow,sumwt,wbar,resx,resy,res common /elem/p,t0,ecc,semi,asc,per,xinc,casc,sasc,cper,sper, $ cxinc,sxinc,txinc degrad=0.0174532925D0 c bmf=b-f bpf=b+f amg=a-g apg=a+g pang = abs(datan(bmf/apg))/degrad xmang = abs(datan(bpf/amg))/degrad if (bmf .lt. 0.D0 .and. apg .lt. 0.D0) pang = 180.D0 + pang if (bmf .lt. 0.D0 .and. apg .gt. 0.D0) pang = 360.D0 - pang if (bmf .gt. 0.D0 .and. apg .lt. 0.D0) pang = 180.D0 - pang if (bpf .lt. 0.D0 .and. amg .lt. 0.D0) xmang = 180.D0 + xmang if (bpf .lt. 0.D0 .and. amg .gt. 0.D0) xmang = 360.D0 - xmang if (bpf .gt. 0.D0 .and. amg .lt. 0.D0) xmang = 180.D0 - xmang asc = (pang+xmang)/2.D0 if (asc .lt. 0.0D0) asc=asc+360.D0 per = (pang-xmang)/2.D0 if (per .lt. 0.0D0) per=per+360.D0 sum1 = a*a + b*b + f*f + g*g sum2 = a*g - b*f p1 = sum1/sum2 p2 = sqrt(p1*p1-4.D0) p3 = abs((p1+p2)/2.D0) p4 = abs((p1-p2)/2.D0) if (p3 .gt. 1.D0) cinc=p4 if (p4 .gt. 1.D0) cinc=p3 sinc = sqrt(1.D0-cinc*cinc) if (cinc .ne. 0.D0) xinc = datan(abs(sinc/cinc))/degrad if (cinc .eq. 0.D0) xinc = 90.D0 if (sum2 .lt. 0.D0) xinc = 180.D0 - xinc if (xinc .ne. 90.D0) semi = sqrt(abs(sum2/dcos(degrad*xinc))) if (xinc .eq. 90.D0) semi = 99.99D0 return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine toko(nobs,elerr,elfix,obs) c c subroutines to determine formal errors, from Andrei Tokovinin's c PC orbit program c implicit double precision (a-h,o-z) real*8 nd(4) dimension alpha(7,7),b(7),beta(7,1),elerr(7),elfix(7), $ obs(2000,6),sclel(7),sd(4),lista(7) common /kepl/B,T,V,PHASE,RES1,RES2, $ e1,ce,se,r,tmw,stmw,ctmw,ro,cf,cf2,dt,ec common /elem/p,t0,ecc,semi,asc,per,xinc,casc,sasc,cper,sper, $ cxinc,sxinc,txinc c pi=3.14159265358979D0 pi2=pi*2.D0 gr=180.D0/pi mfit=7 c c Least-squares fitting of elements c c Make LISTA - list of elements to improve c k=0 DO 300 I=1,7 SCLEL(i)=1.D0 ELERR(i)=0.D0 lista(i-k)=i if (ELFIX(i) .EQ. 0) K=K+1 if (i .ge. 5) sclel(i)=gr 300 continue mfit=7-k c c Accumulate covariance matrix c DO 302 J=1,7 DO 301 K=1,7 ALPHA(J,K)=0.d0 301 continue BETA(J,1)=0.d0 302 continue c CHISQ=0. DO 303 J=1,4 SD(J)=0. ND(J)=0. 303 continue CALL KEPLER(1) c c Data on visual binary orbit c DO 315 I=1,NOBS T=OBS(I,1) CALL KEPLER(2) CALL KEPLER(5) SIG2I=(OBS(I,3)/OBS(I,4))**2 c c Process theta c CALL KEPLER(6) c DY=OBS(I,5)/GR SD(3)=SD(3)+OBS(I,5)*OBS(I,5) CALL COVSUM(MFIT,ALPHA,BETA,B,LISTA,DY,CHISQ,SIG2I) c c Process rho c CALL KEPLER(7) c DY=OBS(I,6)/OBS(I,3) SD(4)=SD(4)+OBS(I,6)*OBS(I,6) CALL COVSUM(MFIT,ALPHA,BETA,B,LISTA,DY,CHISQ,SIG2I) 315 continue c ND(3)=dfloat(NOBS) ND(4)=dfloat(NOBS) c DO 317 J=2,MFIT DO 316 K=1,J-1 ALPHA(K,J)=ALPHA(J,K) 316 continue 317 continue c DO 318 J=1,4 IF (ND(J) .GT. 0.) SD(J)=dsqrt(SD(J)/ND(J)) 318 continue c write(2,'(/5x,"ChiSq=",f10.3,4f8.3/)') CHISQ, (SD(J), J=1,4) write(6,'(/5x,"ChiSq=",f10.3,4f8.3/)') CHISQ, (SD(J), J=1,4) CALL GAUSSJ(ALPHA,MFIT,7,BETA,1,1) SCLEL(1)=(-(p**2))/PI2 SIG=dsqrt(CHISQ/(2*NOBS-MFIT)) c DO 400 I=1,MFIT J=LISTA(I) ELERR(J)=ABS(SCLEL(J))*SIG*dsqrt(ABS(ALPHA(I,I)))*ELFIX(J) 400 continue c return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C SUBROUTINE COVSUM(MFIT,ALPHA,BETA,B,LISTA,DY,CHISQ,SIG2I) implicit double precision (a-h,o-z) dimension LISTA(7),ALPHA(7,7),B(7),BETA(7,1) c DO 200 J=1,7 WT=B(LISTA(J))*SIG2I c c write(6,991) j,lista(j),b(lista(j)),wt c 991 format(' j,lista(j),b(lista(j)),wt: ',2i3,2f12.6) c DO 100 K=1,J ALPHA(J,K)=ALPHA(J,K)+WT*B(LISTA(K)) 100 continue BETA(J,1)=BETA(J,1)+DY*WT 200 continue CHISQ=CHISQ+DY*DY*SIG2I RETURN END c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE GAUSSJ(alpha,mfit,NP,beta,M,MP) implicit double precision (a-h,o-z) dimension IPIV(50),INDXR(50),INDXC(50),alpha(7,7),beta(7,1) c DO 11 J=1,mfit IPIV(J)=0 11 continue c DO 22 I=1,mfit BIG=0.d0 DO 13 J=1,mfit IF (IPIV(J).NE.1) THEN DO 12 K=1,mfit IF (IPIV(K).EQ.0) THEN IF (ABS(alpha(J,K)).GE.BIG) THEN BIG=ABS(alpha(J,K)) IROW=J ICOL=K ENDIF ELSE IF (IPIV(K).GT.1) THEN write(2,'(" Singular matrix: k,ipiv(k) =",2i6)') k,ipiv(k) write(6,'(" Singular matrix: k,ipiv(k) =",2i6)') k,ipiv(k) ENDIF 12 continue ENDIF 13 continue IPIV(ICOL)=IPIV(ICOL)+1 IF (IROW.NE.ICOL) THEN DO 14 L=1,mfit DUM=alpha(IROW,L) alpha(IROW,L)=alpha(ICOL,L) alpha(ICOL,L)=DUM 14 continue DO 15 L=1,M DUM=beta(IROW,L) beta(IROW,L)=beta(ICOL,L) beta(ICOL,L)=DUM 15 continue ENDIF INDXR(I)=IROW INDXC(I)=ICOL IF (alpha(ICOL,ICOL).EQ.0.) write(2,992) icol,alpha(icol,icol) IF (alpha(ICOL,ICOL).EQ.0.) write(6,992) icol,alpha(icol,icol) 992 format(' Singular matrix: icol, alpha(icol,icol) = ',i6,f10.5) PIVINV=1./alpha(ICOL,ICOL) alpha(ICOL,ICOL)=1. DO 16 L=1,mfit alpha(ICOL,L)=alpha(ICOL,L)*PIVINV 16 continue DO 17 L=1,M beta(ICOL,L)=beta(ICOL,L)*PIVINV 17 continue DO 21 LL=1,mfit IF (LL.NE.ICOL) THEN DUM=alpha(LL,ICOL) alpha(LL,ICOL)=0. DO 18 L=1,mfit alpha(LL,L)=alpha(LL,L)-alpha(ICOL,L)*DUM 18 continue DO 19 L=1,M beta(LL,L)=beta(LL,L)-beta(ICOL,L)*DUM 19 continue ENDIF 21 continue 22 continue c DO 24 L=mfit,1,-1 IF (INDXR(L).NE.INDXC(L)) THEN DO 23 K=1,mfit DUM=alpha(K,INDXR(L)) alpha(K,INDXR(L))=alpha(K,INDXC(L)) alpha(K,INDXC(L))=DUM 23 continue ENDIF 24 continue c RETURN END c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE KEPLER(NVAR) c c *************************************************** c Collection of subroutines to calculate orbits c Programmed by A.A.Tokovinin, last change 22-Feb-91 c Elements are stored in the REAL*8 array EL(i): c 1 2 3 4 5 6 7 8 9 10 c P T e a node w i K1 K2 V0 c Data are transfered via common block KEPL: c /EL(j), B(j)=d/dEL(j), epoch T, PHASE, RES1, RES2/ c The control parameter NVAR has the following meaning: c 1 - initialize elements and functions of elements c 2 - calculate orbital phase and solve Kepler equation c 3 - solve Kepler equation c 4 - radial velocity ephemerid (res1,2 = Vel 1,2) c 5 - visual binary ephemerid (res1,2 = TETA,RO) c 6-9 - Derivatives of TETA,RO,V1,V2 over elements c 10 - visual binary ephemerides, RES1=x, RES2=y c ***************************************************** c implicit double precision (a-h,o-z) dimension b(7) common /kepl/B,T,V,PHASE,RES1,RES2, $ e1,ce,se,r,tmw,stmw,ctmw,ro,cf,cf2,dt,ec common /elem/p,t0,ecc,semi,asc,per,xinc,casc,sasc,cper,sper, $ cxinc,sxinc,txinc pi=3.14159265358979D0 pi2=pi*2.D0 gr=180.D0/pi c GO TO (10,20,30,40,50,60,70,80,90,100) NVAR c ----- NVAR=1 ----- Initialize elements and functions of elements c 10 PMU=PI2/P c CF2=1.D0-ecc*ecc CF=dsqrt(CF2) EC=dsqrt((1.D0+ecc)/(1.D0-ecc)) c casc=dcos(asc/GR) sasc=dsin(asc/GR) c cper=dcos(per/GR) sper=dsin(per/GR) c sxinc=dsin(xinc/GR) cxinc=dcos(xinc/GR) txinc=sxinc/cxinc c c Thiele-van-den-Bos elements c AA=semi*(cper*casc-sper*sasc*cxinc) BB=semi*(cper*sasc+sper*casc*cxinc) FF=semi*(-sper*casc-cper*sasc*cxinc) GG=semi*(-sper*sasc+cper*casc*cxinc) RETURN c c ----- NVAR=2 ----- Phase for the moment T and... c 20 DT=T-T0 PHASE=MOD((T-T0)/P,1.0) IF (PHASE .LT. 0.D0) PHASE=PHASE+1.D0 c c ----- NVAR=3 ----- Solve Kepler equation c 30 ANM=PHASE*PI2 E=ANM 31 E1=E+(ANM+ecc*SIN(E)-E)/(1.0-ecc*COS(E)) IF (ABS(E1-E).GT.1.E-5) THEN E=E1 GOTO 31 END IF V=2.*ATAN(EC*TAN(E1/2.)) U=V+W SU=SIN(U) CU=COS(U) RETURN c c ----- NVAR=4 ----- Radial velocity ephemeride c 40 RETURN c c ----- NVAR=5 ----- calc. TETA and RO c 50 CE=dcos(E1) SE=dsin(E1) R=semi*(1.D0-ecc*CE) TMW=datan(cxinc*dtan(U)) IF (CU.LT.0.) TMW=TMW+PI STMW=dsin(TMW) CTMW=dcos(TMW) c RO=R*CU/CTMW RES1=(TMW+asc/GR)*GR IF (RES1 .GT. 360.D0) RES1=RES1-360.D0 IF (RES1 .LT. 0.D0) RES1=RES1+360.D0 RES2=RO RETURN c c ----- NVAR=6 ----- B(j)=dTETA/dEL(j), RES1=TETA c 60 A1=(semi/RO)**2 A2=R/semi A3=A1*CF*cxinc A4=CTMW*STMW B(1)=A3*DT B(2)=-A3*PMU B(3)=A1*cxinc*SE*(A2+CF2)/CF B(4)=0.D0 B(5)=1.D0 B(6)=cxinc*(CTMW/CU)**2 B(7)=-txinc*A4 c c write(6,992) nvar,semi,ro,cf,cxinc,ctmw,stmw,(b(jj),jj=1,7) c 992 format(' nvar,semi,ro,cf,cxinc,ctmw,stmw,b: ',i3,6f12.6, c $ ' / ',7f12.6) c RETURN c c ----- NVAR=7 ----- B(j)=1/RO*dRO/dEL(j) c 70 A1=(semi/R)**2 A2=R/semi A4=CTMW*STMW A5=A4*txinc*sxinc A3=A1*(ecc*SE-A5*CF) c B(1)=A3*DT B(2)=-A3*PMU B(3)=-A1*((CE-ecc)*CF+A4*(A2+CF2)*SE)/CF B(4)=1.D0/semi B(5)=0.D0 B(6)=-A5 B(7)=-txinc*STMW*STMW c c write(6,993) nvar,semi,r,ctmw,stmw,ti,si,ecc,se,cf, c $ (b(jj),jj=1,7) c 993 format(' nvar,semi,r,ctmw,stmw,ti,si,ecc,se,cf,b: ',i3,9f12.6, c $ ' / ',7f12.6) c RETURN c c ----- NVAR=8 ----- B(j)=1/K1*dV1/dEL(j), RES1=V1 c 80 RETURN c c ----- NVAR=9 ----- B(j)=1/K2*dV2/dEL(j), RES1=V2 c 90 RETURN c c ----- NVAR=10 ----- c 100 RETURN c END c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc