program orbgrids
c
c	Version of orbgrid7 which reads WDS-style title line.
c	Also reads pm from title line, divides by 15 for arcsec --> sec
c
c       2004 sep 24 - removed old normal point routines, cleaned up.
c       At this point will compile propoerly and all errors are essentially
c       the same as before except error in a, which for example file is about
c       a factor of 10 too high.               
c
c       2006 jun 02 - changed input to take name of form 'wds12345+1234'
c       and append '.dat', .orb', and '.plt' for input and output filenames.
c       (version orbgrid10)
c
c       2006 jun 02 - added input list option (version orbgrids)

	integer j,nobs,it,ie,ip,ipass,nloop,idsep,nnum,iip,iit,iie,n
	dimension ss(10),xma(2000),error(9,9,9),theta0(2000)
        real*8 elerr(7),elfix(7),obs(2000,6)
	character*1  decs
	character*3  code1,xcode
        character*13 filein
	character*18 code2,ycode
	character*20 file01,file02,file03
	character*6 atheta,arho
	character*130 name
	common /comchar/code1(2000),code2(2000)
	common /numcom/qt,qth,qrho,qwt,tnmin(500),tnmax(500),nnum(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),sumwt,wbar,resx,resy,res
	degrad=0.0174532925
	equ=2000.
c
c	This program calculates Thiele-Innes elements and their errors from a
c       set of weighted visual measures.  it does a grid search over a range
c	of P, E, and T0, calculating values of (A,B,F,G) by least squares.  
c	After determining a (P,E,T0) set giving minimum errors, the grid spacing
c	is reduced and the process repeated. The first grid search ends when 
c	all grid sizes decrease below 0.01 .  Weights are then recalculated
c	based on residuals and the grid search restarted.  The final grid
c	search ends when all grid sizes decrease below 0.0001 .
c
c	Enter file names for i/o
c
        open(10,file='orbgrids.in',status='UNKNOWN')
  100   read(10,901,end=999) filein
  901	format(a13)
c
        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,t63,f4.3)
	read(1,*) p0,t00,ecc0,pstep0,tstep0,estep0
	pm=pm/15.
	write(2,'(a130)') name
	write(6,'(a130)') name
	ra=rah*15. + ram/4.
	dec=decd + decm/60.
	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
	corr = (0.0056*sin(degrad*ra)/cos(degrad*dec))
     $     +0.00417*pm*sin(degrad*dec)
	ipass=1
	nloop=0
	nloop2=0
	nobs=0
c
c	Read in data  
c
  150	read(1,906,end=400) qt,qth0,qrho,xcode,qwt,ycode,atheta,arho
  906 	format(6x,f9.4,4x,f6.2,7x,f8.4,9x,a3,f5.1,5x,a18,
     $        t20,a6,t35,a6)
	if (qt .gt. 2020. .or. qt .lt. 1500.) go to 900
	if (atheta .eq. ' ' .or. arho .eq. ' ') go to 150
c
c	Store all data into arrays for plot input file
c
	qth = qth0 + corr*(equ - qt)
	if (qth .lt.   0.)   qth=qth+360.
	if (qth .gt. 360.)   qth=qth-360.
	if (qth .lt.   0.01) qth=  0.01
	if (qth .gt. 359.99) qth=359.99
c
  200	nobs=nobs+1
	t(nobs)=qt
	theta0(nobs)=qth0
	theta(nobs)=qth
	rho(nobs)=qrho
	code1(nobs)=xcode
	code2(nobs)=ycode
	w(nobs)=qwt
	xo(nobs) = rho(nobs)*cos(degrad*theta(nobs))
	yo(nobs) = rho(nobs)*sin(degrad*theta(nobs))
	go to 150
c
  400	if (ipass .eq. 1) write(2,907) (t(j),theta(j),rho(j),xo(j),
     $     yo(j),code1(j),code2(j),w(j),j=1,nobs)
  	if (ipass .eq. 1) write(6,907) (t(j),theta(j),rho(j),xo(j),
     $     yo(j),code1(j),code2(j),w(j),j=1,nobs)
  907 	format(2x,f11.4,f8.2,f8.4,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.
	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 the
c	program hops out of the loop.  Sizmin is set to 0.01 for the first 
c	pass.  For the second pass, weights are recalculated, the initial 
c	step sizes are cut in half, and sizmin is reset to 0.0001 .  Sizmin
c	for eccentricity is 10 times more stringent than for P and T0.
c
	if (ipass .eq. 1) sizmin=0.01
	if (ipass .ge. 2) sizmin=0.0001
	if (ipass .eq. 2) pstep=pstep/2.
	if (ipass .eq. 2) tstep=tstep/2.
	if (ipass .eq. 2) estep=estep/2.
	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.
	if (eccx .ge. 0.9999) ecc02=0.9999-estep*4.	
    	do 720 ip=1,9
	p=p0+pstep*real(ip-5)
	xmu = 360.0/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(ecc,xma,nobs,a,b,f,g,ss)
  	if (estep+pstep+tstep .eq. 0.0) 
     $     call elements(a,b,f,g,p,t0,ecc,semi,xinc,asc,per,pm)
	call errcalc(nobs,a,b,f,g,sigx,sigy,sigv1,sigv2,
     $     sigs1,sigs2,errtot,1)
	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.0/p
	do 730 j=1,nobs
   	xma(j)=xmu*(t(j)-t0)
  730	continue
	call thin(ecc,xma,nobs,a,b,f,g,ss)
	call errcalc(nobs,a,b,f,g,sigx,sigy,sigv1,sigv2,
     $     sigs1,sigs2,errmin,2)
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 grid 
c	    dimensions, the same values are used for the next grid search.  
c	(2) If the same step sizes have been used 5 times, it's figured a 
c	    larger step size is needed.  All step sizes are increased by
c	    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.0
	if (abs(t0-t00)   .gt. tstep) sfact=1.0
	if (abs(ecc-ecc0) .gt. estep) sfact=1.0
	nloop2=nloop2+1
	if (sfact .ne. 1.0) nloop2=0
	if (nloop2 .eq. 5) sfact=1.5
	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,p,t0,ecc,semi,xinc,asc,per,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))
	if (w(n) .eq. 0.) obs(n,4)=dble(9.99)
	obs(n,5)=dble(dpa(n))
	obs(n,6)=dble(dr(n))
  810	continue
c
	do 820 k=1,7
  820	elfix(k)=1.
	if (pstep0 .eq. 0.) elfix(1)=0.
	if (tstep0 .eq. 0.) elfix(2)=0.
	if (estep0 .eq. 0.) elfix(3)=0.

	call toko(nobs,p,t0,ecc,semi,asc,per,xinc,elerr,elfix,obs)

	pday    = p*365.2422
	errpday = elerr(1)*365.2422
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 (code1(j) .eq. 'VIS') go to 830
	if (code1(j) .eq. 'XVI') 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=atan(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),code1(j)
	write(6,915) t(j),theta(j),rho(j),pa,sep,dpa(j),dr(j),idsep,
     $     rdt(j),w(j),code1(j)
  915	format(f10.4,f8.2,f8.4,f9.2,f8.4,f8.2,f8.4,i5,f8.3,f6.1,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),code1(n),w(n),code2(n)
  917	format(6x,f9.4,4x,f6.2,7x,f8.4,9x,a3,f5.1,5x,a18)
        go to 100
	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)
c
        go to 100
  999	stop
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine kepler1(am,ecc,ae)
	degrad = 0.0174532925
	if (ecc .lt. 0.0001) ecc=0.0001
	if (ecc .gt. 0.9999) ecc=0.9999
	eo = am + (ecc*sin(degrad*am))/degrad
  100	e1 = am + (ecc*sin(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
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine thin(ecc,xma,nobs,a,b,f,g,ss)
	character*3  code1
	character*18 code2
	integer j,nobs,i,nnum
	dimension xma(2000),ss(10),ea(2000)
	common /comchar/code1(2000),code2(2000)
	common /numcom/qt,qth,qrho,qwt,tnmin(500),tnmax(500),nnum(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),sumwt,wbar,resx,resy,res
	degrad=0.0174532925
c
	if (ecc .lt. 0.0001) ecc=0.0001
	if (ecc .gt. 0.9999) ecc=0.9999
	f1 = sqrt(1.0-(ecc*ecc))
c
	do 100 j=1,nobs
	call kepler1(xma(j),ecc,ea(j))
   	x(j) = cos(degrad*ea(j)) - ecc
	y(j) = f1*sin(degrad*ea(j))
  100	continue
c
	do 200 i=1,10
  200	ss(i)=0.0
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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine interp(error,pstep,tstep,estep,p0,t00,ecc0,
     $     p,t0,ecc)
	character*3  code1
	character*18 code2
	integer it,ie,ip,nnum,ipm,itm,iem
	dimension error(9,9,9)
	common /comchar/code1(2000),code2(2000)
	common /numcom/qt,qth,qrho,qwt,tnmin(500),tnmax(500),nnum(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),sumwt,wbar,resx,resy,res
	degrad=0.0174532925
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.
	b=y2-y1-a*(2.*tm-1.)
	tm=(-0.5)*b/a
	if (y1/y2 .gt. 10. .or. y3/y2 .gt. 10.) 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.
	b=y2-y1-a*(2.*em-1.)
	em=(-0.5)*b/a
	if (y1/y2 .gt. 10. .or. y3/y2 .gt. 10.) 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.0001) ecc=0.0001
	if (ecc .gt. 0.9999) ecc=0.9999
c
	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine errcalc(nobs,a,b,f,g,sigx,sigy,sigv1,sigv2,
     $                     sigs1,sigs2,errtot,icode)
	character*3  code1
	character*18 code2
	integer j,nobs,nnum,icode
	common /comchar/code1(2000),code2(2000)
	common /numcom/qt,qth,qrho,qwt,tnmin(500),tnmax(500),nnum(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),sumwt,wbar,resx,resy,res
	degrad=0.0174532925
c
	sigx=0.
	sigy=0.
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.
	if (icode .eq. 1) return
c
	sigv1=0.
	sigv2=0.
	sigs1=0.
	sigs2=0.
	env=0.
	ens=0.
	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=atan(yc/xc)/degrad
	if (xc .lt. 0.) pa=pa+180.
	if (xc .gt. 0. .and. yc .lt. 0.) pa=pa+360.
	dr(j)=rho(j)-sep
	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.
	rdt(j)=sep*dpa(j)*degrad
c
	if (w(j) .eq. 0.) go to 400
	if (code1(j) .eq. 'VIS') go to 300
	if (code1(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,p,t0,ecc,semi,xinc,asc,per,pm)
	character*3  code1
	character*18 code2
	integer nnum
	common /comchar/code1(2000),code2(2000)
	common /numcom/qt,qth,qrho,qwt,tnmin(500),tnmax(500),nnum(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),sumwt,wbar,resx,resy,res
	degrad=0.0174532925
c
	bmf=b-f
	bpf=b+f
	amg=a-g
	apg=a+g
	pang  = abs(atan(bmf/apg))/degrad
	xmang = abs(atan(bpf/amg))/degrad
	if (bmf .lt. 0. .and. apg .lt. 0.) pang  = 180. + pang
	if (bmf .lt. 0. .and. apg .gt. 0.) pang  = 360. - pang
	if (bmf .gt. 0. .and. apg .lt. 0.) pang  = 180. - pang
	if (bpf .lt. 0. .and. amg .lt. 0.) xmang = 180. + xmang
	if (bpf .lt. 0. .and. amg .gt. 0.) xmang = 360. - xmang
	if (bpf .gt. 0. .and. amg .lt. 0.) xmang = 180. - xmang
	asc = (pang+xmang)/2.
	if (asc .lt. 0.0) asc=asc+360.
	per = (pang-xmang)/2.
	if (per .lt. 0.0) per=per+360.
	sum1 = a*a + b*b + f*f + g*g
	sum2 = a*g - b*f
	p1 = sum1/sum2
	p2 = sqrt(p1*p1-4.)
	p3 = abs((p1+p2)/2.)
	p4 = abs((p1-p2)/2.)
	if (p3 .gt. 1.) cinc=p4
	if (p4 .gt. 1.) cinc=p3
	sinc = sqrt(1.-cinc*cinc)
	if (cinc .ne.  0.) xinc = atan(abs(sinc/cinc))/degrad
	if (cinc .eq.  0.) xinc =  90.
	if (sum2 .lt.  0.) xinc = 180. - xinc
	if (xinc .ne. 90.) semi = sqrt(abs(sum2/cos(degrad*xinc)))
	if (xinc .eq. 90.) semi = 99.99
	return
	end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc	
c
	subroutine toko(nobs,p,t0,ecc,semi,asc,per,xinc,elerr,elfix,obs)
c
c	subroutines to determine formal errors, 
c	from Andrei Tokovinin's PC orbit program 
c
	real*8 alpha(7,7),b(7),beta(7),chisq,dy,el(7),elerr(7),elfix(7),
     $     gr,nd(4),obs(2000,6),ochisq,phase,pi,pi2,res1,res2,sclel(7),
     $     sd(4),sig,sig2i,t,v
        integer lista(7)
        common /kepl/EL,B,T,V,PHASE,RES1,RES2
c
	pi=3.1415927
	pi2=pi*2.
	gr=180./pi
	mfit=7
c
  	el(1)=p
	el(2)=t0
	el(3)=ecc
	el(4)=semi
	el(5)=asc
	el(6)=per
	el(7)=xinc
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.
	ELERR(i)=0.
  	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,MFIT
        DO 301 K=1,J
        ALPHA(J,K)=0.
  301	CONTINUE
        BETA(J)=0.
  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)
        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)
	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)=(-(EL(1)**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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc	
C
        SUBROUTINE COVSUM(MFIT,ALPHA,BETA,B,LISTA,DY,CHISQ,SIG2I)
        REAL*8 ALPHA(7,7),B(7),BETA(7),CHISQ,dy,SIG2I,WT
        INTEGER LISTA(7)
c
	DO 200 J=1,7
	WT=B(LISTA(J))*SIG2I
	DO 100 K=1,J
	ALPHA(J,K)=ALPHA(J,K)+WT*B(LISTA(K))
  100	CONTINUE
	BETA(J)=BETA(J)+DY*WT
  200	CONTINUE
	CHISQ=CHISQ+DY*DY*SIG2I
	RETURN
	END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc	
c
      SUBROUTINE GAUSSJ(A,N,NP,B,M,MP)
      real*8 A(7,7),B(7,1),big,dum,pivinv
      integer IPIV(50),INDXR(50),INDXC(50)
      DO 11 J=1,N
        IPIV(J)=0
11    CONTINUE
      DO 22 I=1,N
        BIG=0.
        DO 13 J=1,N
          IF (IPIV(J).NE.1) THEN
            DO 12 K=1,N
              IF (IPIV(K).EQ.0) THEN
                IF (dABS(A(J,K)).GE.BIG) THEN
                  BIG=ABS(A(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,N
            DUM=A(IROW,L)
            A(IROW,L)=A(ICOL,L)
            A(ICOL,L)=DUM
14        CONTINUE
          DO 15 L=1,M
            DUM=B(IROW,L)
            B(IROW,L)=B(ICOL,L)
            B(ICOL,L)=DUM
15        CONTINUE
        ENDIF
        INDXR(I)=IROW
        INDXC(I)=ICOL
        IF (A(ICOL,ICOL).EQ.0.) write(2,992) icol,a(icol,icol)
        IF (A(ICOL,ICOL).EQ.0.) write(6,992) icol,a(icol,icol)
  992   format(' Singular matrix: icol, a(icol,icol) = ',i6,f10.5)
        PIVINV=1./A(ICOL,ICOL)
        A(ICOL,ICOL)=1.
        DO 16 L=1,N
          A(ICOL,L)=A(ICOL,L)*PIVINV
16      CONTINUE
        DO 17 L=1,M
          B(ICOL,L)=B(ICOL,L)*PIVINV
17      CONTINUE
        DO 21 LL=1,N
          IF (LL.NE.ICOL) THEN
            DUM=A(LL,ICOL)
            A(LL,ICOL)=0.
            DO 18 L=1,N
              A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18          CONTINUE
            DO 19 L=1,M
              B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
19          CONTINUE
          ENDIF
21      CONTINUE
22    CONTINUE
      DO 24 L=N,1,-1
        IF (INDXR(L).NE.INDXC(L)) THEN
          DO 23 K=1,N
            DUM=A(K,INDXR(L))
            A(K,INDXR(L))=A(K,INDXC(L))
            A(K,INDXC(L))=DUM
23        CONTINUE
        ENDIF
24    CONTINUE
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc	
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
        real*8 a,a1,a2,a3,a4,a5,aa,anm,b(7),bb,ce,cf,cf2,cf3,ci,ctmw,
     $     cu,cv,cw,cww,dt,e,e1,ec,el(7),ff,gg,gr,k1,k2,p,phase,pi,
     $     pi2,pmu,r,res1,res2,ro,se,sf,si,stmw,su,sv,sw,sww,t,t0,ti,
     $     tmw,u,v,v0,w,ww,x,y
        common /kepl/EL,B,T,V,PHASE,RES1,RES2
	pi=3.1415927
	pi2=pi*2.
	gr=180./pi
c
        GOTO (10,20,30,40,50,60,70,80,90,100) NVAR
c
c ----- NVAR=1 ----- Initialize elements and functions of elements
c 
10      P=EL(1)
        PMU=PI2/P
        T0=EL(2)
c
        SF=EL(3)
        CF2=1.-SF*SF
        CF=dsqrt(CF2)
        CF3=CF*CF2
        EC=dsqrt((1.+SF)/(1.-SF))
c
        A=EL(4)
c
        WW=EL(5)/GR
        CWW=dcos(WW)
        SWW=dsin(WW)
c
        W=EL(6)/GR
        CW=dcos(W)
        SW=dsin(W)
c
        SI=dsin(EL(7)/GR)
        CI=dcos(EL(7)/GR)
        TI=SI/CI
c        K1=EL(8)
c        K2=EL(9)
c        V0=EL(10)
c
c Thiele-van-den-Bos elements
c
        AA=A*(CW*CWW-SW*SWW*CI)
        BB=A*(CW*SWW+SW*CWW*CI)
        FF=A*(-SW*CWW-CW*SWW*CI)
        GG=A*(-SW*SWW+CW*CWW*CI)
        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.) PHASE=PHASE+1.
c
c ----- NVAR=3 ----- Solve Kepler equation
c 
30      ANM=PHASE*PI2
        E=ANM
31      E1=E+(ANM+SF*dsin(E)-E)/(1.0-SF*dcos(E))
        IF (dABS(E1-E) .GT. 1.E-5) THEN
          E=E1
          GOTO 31
        END IF
        V=2.*datan(EC*dtan(E1/2.))
        U=V+W
        SU=dsin(U)
        CU=dcos(U)
        RETURN
c
c ----- NVAR=4 ----- Radial velocity ephemeride
c 
40      A1=SF*CW+CU
        RES1=V0+K1*A1
        RES2=V0-K2*A1
        RETURN
c
c ----- NVAR=5 ----- calc. TETA and RO
c 
50      CE=dcos(E1)
        SE=dsin(E1)
        R=A*(1.-SF*CE)
        TMW=datan(CI*dtan(U))
        IF (CU.LT.0.) TMW=TMW+PI
        STMW=dsin(TMW)
        CTMW=dcos(TMW)
c
        RO=R*CU/CTMW
        RES1=(TMW+WW)*GR
        IF (RES1 .GT. 360.) RES1=RES1-360.
        IF (RES1 .LT.   0.) RES1=RES1+360.
        RES2=RO
        RETURN
c
c ----- NVAR=6 ----- B(j)=dTETA/dEL(j), RES1=TETA
c
60      A1=(A/RO)**2
        A2=R/A
        A3=A1*CF*CI
        A4=CTMW*STMW

        B(1)=A3*DT
        B(2)=-A3*PMU
        B(3)=A1*CI*SE*(A2+CF2)/CF
        B(4)=0.
        B(5)=1.
        B(6)=CI*(CTMW/CU)**2
        B(7)=-TI*A4
        RETURN
c
c ----- NVAR=7 ----- B(j)=1/RO*dRO/dEL(j)
c 
70      A1=(A/R)**2
        A2=R/A
        A4=CTMW*STMW
        A5=A4*TI*SI
        A3=A1*(SF*SE-A5*CF)
c
        B(1)=A3*DT
        B(2)=-A3*PMU
        B(3)=-A1*((CE-SF)*CF+A4*(A2+CF2)*SE)/CF
        B(4)=1./A
        B(5)=0.
        B(6)=-A5
        B(7)=-TI*STMW*STMW
        RETURN
c
c ----- NVAR=8 ----- B(j)=1/K1*dV1/dEL(j), RES1=V1
c
80      SV=dsin(V)
        CV=dcos(V)
        A1=(1.+SF*CV)**2*SU/CF3
c
        B(1)=(-A1)*DT
        B(2)=A1*PMU
        B(3)=CW-SU*SV*(2.+SF*CV)/CF2
        B(4)=0.
        B(5)=0.
        B(6)=-SU-SF*SW
        B(7)=0.
c        B(8)=CU+SF*CW
c        RES1=V0+K1*B(8)
c        B(8)=B(8)/K1
c        B(9)=0.
c        B(10)=1./K1
        RETURN
c
c ----- NVAR=9 ----- B(j)=1/K2*dV2/dEL(j), RES1=V2
c 
90      SV=dsin(V)
        CV=dcos(V)
        A1=(1.+SF*CV)**2*SU/CF3
c
        B(1)=A1*DT
        B(2)=-A1*PMU
        B(3)=-CW+SU*SV*(2.+SF*CV)/CF2
        B(4)=0.
        B(5)=0.
        B(6)=SU+SF*SW
        B(7)=0.
c        B(8)=-CU-SF*CW
c        RES1=V0+K2*B(8)
c        B(9)=B(8)/K2
c        B(8)=0.
c        B(10)=1./K2
        RETURN
c
c ----- NVAR=10 ----- 
c
100     CV=dcos(V)
        R=CF2/(1.+SF*CV)
        X=R*CV
        Y=R*dsin(V)
        RES1=AA*X+FF*Y
        RES2=BB*X+GG*Y
        RETURN
c
        END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc