c
      program obs2vect
c
c   calculate speed and pa from MPC formatted observations.
c   input data read from STDIN, use redirect to read from file
c      e.g.   obs2vect < mpcformatteddata.dat
c
      implicit real*8(a-h,o-z)
      real*8 pi
      parameter (pi = 3.14159265358979)
      character*1 dsign,color,discover,detector
      character*7 designation
      character*3 obscode
      integer year,month,rahour,ramin,decdeg,decmin
      parameter (nmax=10)   !number of data points per object
      real*8 t(nmax),ra(nmax),dec(nmax),sig(nmax)  ! sig dummy array
c
c     
      write (6,10)
 10   format (/,"Confirm input data:")
      do 20 idata=1,nmax
         read (5,16,end=21) designation,discover,detector,year,month,
     &        day,rahour,ramin,rasec,
     &        dsign,decdeg,decmin,decsec,xmag,
     &        color,obscode
 16      format(5x,a7,a1,1x,a1,I4,1X,I2.2,1X,f8.5,1X,I2.2,1X,I2.2,
     &        1X,f5.2,1X,A1,I2.2,1X,I2.2,1X,f4.1,
     &        9x,F5.1,1x,A1,6x,a3,1x)
c
c
c   confirm input
         write (6,16) designation,discover,detector,year,month,
     &        day,rahour,ramin,rasec,
     &        dsign,decdeg,decmin,decsec,xmag,
     &        color,obscode
c
         if (idata .eq. 1) temptime = mod(day,1.0)*24.0 
         t(idata) = mod(day,1.0)*24.0 - temptime   ! time in hours
         ra(idata) = ratodeg (rahour,ramin,rasec)
         dec(idata) = dectodeg(dsign,decdeg,decmin,decsec)
cdebug          print *, idata, day,t(idata),ra(idata),dec(idata)
         ndata = idata
 20   continue
 21   continue
c
      write (6,25) ndata
 25   format(/,i4," data points read")
c
c     if ndata > 2, use least squares to fit trajectories
      if (ndata .gt. 2) then
         call fit(t,ra,ndata,sig,0,a,vra,siga,sigb,chi2)
         call fit(t,dec,ndata,sig,0,a,vdec,siga,sigb,chi2)
      else if (ndata .eq. 2) then
         vra  = (ra(2)-ra(1))/(t(2)-t(1))
         vdec = (dec(2)-dec(1))/(t(2)-t(1))
cdebug         print *, vra,vdec
      else    
         write (6,26)      
 26      format (/,"ERROR: Not enough points to find trajectory")
         vra  = -9999999.0
         vdec = -9999999.0
      end if
c
      if (ndata .ge. 2) then 
c   convert units
         decmean=(dec(ndata)+dec(1))/2.0*pi/180.0
         vra = vra*60.0*cos(decmean)    ! " per minute
         vdec = vdec*60.0   ! " per minute
cdebug        print *, vra,vdec
c
         call cartopolar(vra,vdec,v,angle)
         write (6,28) v,angle
 28      format(/,' Speed ("/min)= ',f6.2,'   P.A. = ', f6.1,/)
      end if 
c
      stop
      end
c
c-------------------------------------------------------------------
      double precision function dectodeg (dsign,decdeg,decmin,decsec)
c
c     convert declination to decimal degrees
c
      implicit real*8 (a-h,o-z)
      character*1 dsign
      integer decdeg,decmin
c
      
      temp = dble(decdeg)+dble(decmin/60.0)+dble(decsec/3600.0)
c
      if (dsign .eq. "-") then 
         dectodeg = -1.0*temp
      else
         dectodeg = temp
      endif
c
      return
      end
c
c-------------------------------------------------------------------
      double precision function ratodeg (rahour,ramin,rasec)
c
c     convert declination to decimal degrees
c
      implicit real*8 (a-h,o-z)
      integer rahour,ramin
c
      
      ratodeg = (dble(rahour)+(dble(ramin)/60.0)+(rasec/3600.0))*15.0
c
      return
      end
c
c-------------------------------------------------------------------
      subroutine cartopolar (x,y,mag,angle)
c
c     convert ra and dec motion to speed and PA
c
      implicit real*8 (a-h,o-z)
      real*8 mag
      pi = 3.141592654
c
      mag = sqrt(x**2+y**2)      
      angle = (atan (x/y))*180.0/pi
      if (y .ge. 0.0) then 
         if (x .le. 0.0) then              
            angle = angle + 360.0
         endif
      else    ! y <0
         angle = 180.0 + angle
      endif
c
      return
      end
c
c-------------------------------------------------------------------
      subroutine fit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2)
c
c  linear least squares - adapted from Numerical Recipies
c
c fits (x,y) data to a straight line y = a + bx by minimizing chi^2
c
      integer mwt,ndata
      real*8 a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
c
c USES gammq
c
c ndata:                 number of data points
c x(1:ndata),y(1:ndata): data points
c sig(1:ndata):          
c a,siga:                y-intercept and error in y-intercept
c b,sigb:                slope and error in slope
c mwt:                   set to 0 to ingnore sigma weighting
c chi2:                  chi-squared (normalized to unit standard deviation) 
c
      integer i
      real sigdat,ss,st2,sx,sxoss,sy,t,wt,gammq
      sx=0. ! Initialize sums to zero.
      sy=0.
      st2=0.
      b=0.
c
      if(mwt.ne.0) then ! Accumulate sums ...
      ss=0
c.
         do 11 i=1,ndata ! ...with weights
            wt=1./(sig(i)**2)
            ss=ss+wt
            sx=sx+x(i)*wt
            sy=sy+y(i)*wt
11       continue 
      else
         do 12 i=1,ndata ! ...or without weights.
            sx=sx+x(i)
            sy=sy+y(i)
12       continue 
         ss=float(ndata)
      endif
c
      sxoss=sx/ss
c
      if(mwt.ne.0) then
         do 13 i=1,ndata
            t=(x(i)-sxoss)/sig(i)
            st2=st2+t*t
            b=b+t*y(i)/sig(i)
13       continue 
      else
         do 14 i=1,ndata
            t=x(i)-sxoss
            st2=st2+t*t
            b=b+t*y(i)
14       continue 
      endif
c
      b=b/st2 ! Solve for a, b, a, and b.
      a=(sy-sx*b)/ss
      siga=sqrt((1.+sx*sx/(ss*st2))/ss)
      sigb=sqrt(1./st2)
c
      chi2=0. ! Calculate chi^2 .
      q=1.
      if(mwt.eq.0) then
         do 15 i=1,ndata
            chi2=chi2+(y(i)-a-b*x(i))**2
15       continue 
         sigdat=sqrt(chi2/(ndata-2))
c 
c For unweighted data evaluate typical sig using chi2, and adjust 
c the standard deviations.
         siga=siga*sigdat
         sigb=sigb*sigdat
      else
         do 16 i=1,ndata
            chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2
16       continue 
      endif
c
      return
      end