*     789012345678901234567890123456789012345678901234567890123456789012
*        1         2         3         4         5         6         7
      program KeplerAssim

*     ==================================================================
*     Assimillation of planetary observations to determine orbital
*     parameters.  Fortran 77.

*     Observations and other data are read-in from a file.  A run of
*     this code does one of four things:
*       1. XSEC plot cross sections through the cost function.
*       2. ANAL analyse observations for their best-fit orbital params.
*       3. RECO Turn alt/azi observations in input file into RA/dec
*               observations in a new input file (for trials only).
*       4. TRAJ Use params in input file to generate planetary
*               trajectories.
*     See documentation "www.met.rdg.ac.uk/~ross/Astronomy/Kepler.html"
*     for further info.
*     ==================================================================
*     Development history:
*     Ross Bannister, February 2001 to February 2003
*     Include modifications as recommended by reviewers of paper
*       1. Option for second derivs of forward model in Hessian calc
*          (direct inversion maintained), December 2002.
*       2. Use double precision numbers, February 2003.
*       3. Include Gaussian ellimination, February 2003
*     ==================================================================

*     Declare variables
*     =================
      implicit none

*     ----- The No of observations (#alt+#azi+#RA+#dec measurements) is
*     Nobs.  This should be <= maxobs
      integer maxobs,Nobs
      parameter (maxobs=50)

*     ----- Variables that hold information about the obs.
      integer tpe(maxobs)
      real*8 y0(maxobs),Error(maxobs),Time(maxobs),long,lat,x(6)
      real*8 y(maxobs)
      character*8 filen
      character*4 rmode
      character*12 fname
      logical verbose,CorrectRefract,FinalRun,d2Hess

*     ----- Variables that are used in the analysis
      integer i,j,k,iter,maxiter,partial
      real*8 pi,ms,crit,CostFunc,amax,amin,imax,imin,range(6),x0(6)
      real*8 H(maxobs,6),Hess(6,6),Hess1(6,6),Hess2(6,6),Hd2(maxobs,6,6)
      real*8 InvH(6,6),Grad(6)
      real*8 Temp1(maxobs),Temp2(maxobs),Temp3(maxobs),Temp4(maxobs)
      real*8 Lambda(6,6),U(6,6),deltay(maxobs),deltax(6)
      real*8 total,low,high,v(6),trace,CostJ,son,lim

*     ----- Variables associated with the forward model
      integer MonthDays(0:12),zeroYr
      real*8 RefTime,NumDays,obl,Epoch,aE,eE,iE,Omega1E,omega2E,LE





*     Initialze Data
*     ==============
      data MonthDays/0,31,59,90,120,151,181,212,243,273,304,334,365/
      pi=4.0*atan(1.0)
*     Zero of time will be 1st Jan 1990
      zeroYr=1990
*     Reference time for alt\azi conversion
      RefTime=NumDays(zeroYr,1990,3,21,0,0,MonthDays)
      obl=pi*(23.0+(26.0+21.0/60.0)/60.0)/180.0
*     Epoch is called t0 in many subroutines
      Epoch=NumDays(zeroYr,2000,1,1,0,0,MonthDays)
*     Earth parameters
      data aE,eE,iE,Omega1E,omega2E,LE/1.00000011, 0.01671022, 0.00005,
     :     -11.26064, 102.94719, 100.46435/





*     Read-in observational data, incl errs and long \ lat, etc.
*     ==========================================================
      call ReadObs (tpe,y0,Error,Time,Nobs,maxobs,zeroYr,long,lat,x,
     :              rmode,verbose,amin,amax,imin,imax,crit,maxiter,
     :              partial,CorrectRefract,d2Hess,MonthDays,fname)
*     y0 array contains the observations
*     Error contains the associated errors
*     Time contains the time of observations
*     x is the inital guess at the parameters
*     rmode is the run mode
*     amin,amax is the range of parameter a (similarly for imin,imax)

      if (verbose) then
        print*,'======================================================='
        print*,'Astronomical Data Assimmilation'
        print*,'Ross Bannister, February/April 2001 & July/August 2002'
        print*,'University of Reading, UK'
        print*,'======================================================='
        print*
        print*,'longitude =',long
        print*,'latitude =',lat
        print*,'Number of observation pairs =',Nobs/2
        print*,'Initial guess: ',(x(k),k=1,6)
        print*,'Range of parameter a: ',amin,amax
        print*,'Range of parameter i: ',imin,imax
        print*,'Run mode = ',rmode
      endif






*     'IF' STATEMENTS THAT EXECUTE THE CHOSEN RUN MODE OF THE CODE
*     ============================================================

*     ------------------------------------------------------------------
      if (rmode.eq.'XSEC') then
*     -------------------------
      if (verbose) print*,'RUN MODE: Compute cost fn cross sections'
*     Probe some of physical parameter space to find out how J behaves
*     Define limits of the varibles about the x already defined above
      data range/5.0,0.4,180.,180.,180.,180./

      do k=1,6
        x0(k)=x(k)
      enddo

      do k=1,6
*       Adjust parameter k in turn
*       Open a file for the output
        write (filen,'(A5,I1)') 'Param',k
        if (verbose) print*,filen
        open (12,file=filen)
        do i=0,150
*         Set range for parameter to vary over
          if (k.eq.1) then
            x0(k)=amin+real(i)*(amax-amin)/150.
          elseif (k.eq.3) then
            x0(k)=imin+real(i)*(imax-imin)/150.
          else
            x0(k)=x(k)+range(k)*real(i-75)/75.
          endif
          call CheckRange (x,verbose)
*         Find the projected observations with this set of parameters
          call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,
     :                  long,lat,RefTime,tpe,Nobs,maxobs,y,obl)
*         Predicted observations are in y
*         Find the cost function
          write (12,202) x0(k),CostFunc(y,y0,Error,Nobs,maxobs)
        enddo
        close (12)
        x0(k)=x(k)
      enddo










      elseif (rmode.eq.'ANAL') then
*     -----------------------------
*     ------------------------------------------------------------------
*     FinalRun allows a subset of the code within the ANAL loop to be
*     rerun after the analysis has been computed (to help compute
*     analysis errors).  Set to true when converged
      FinalRun=.false.

      if (verbose) print*,'RUN MODE: Perform least squares minimization'
*     Open a file for the output
      open (12,file='Convergence')
      write (12,203) '#============================================='//
     :               '==================================='
      write (12,203) '# No.     a       e         i       Omega '//
     :               'omega-bar   Epsilon      Inc.     Cost'
      write (12,203) '#---------------------------------------------'//
     :               '-----------------------------------'

      iter=0
*     Output the initial parameters and cost function
      call AstroFM (x,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,
     :               long,lat,RefTime,tpe,Nobs,maxobs,y,obl)
      write (12,201) iter,(x(i),i=1,6),ms,
     :               CostFunc(y,y0,Error,Nobs,maxobs)


10    continue

*     This is the start of the minimization loop
*     ==========================================
      if (verbose) print*,'--------------------------------------------'
*     (x is the linearization point)
      if (FinalRun) then
        print*,'##### PERFORMING ERROR ANALYSIS #####'
      else
        iter=iter+1
        print*,'##### STARTING ITERATION NO. ',iter,' #####'
      endif

20    continue

      if (.not.FinalRun) then
        if (iter.le.partial) then
*         Do a partial minimization over inclination (parameter 3)
          call PartialMin (3,imin,imax,x,aE,eE,iE,Omega1E,omega2E,LE,
     :        Time,Epoch,long,lat,RefTime,tpe,Nobs,y,obl,y0,Error,
     :        maxobs,verbose)
*         Do a partial minimization over semi-major axis (parameter 1)
          call PartialMin (1,amin,amax,x,aE,eE,iE,Omega1E,omega2E,LE,
     :        Time,Epoch,long,lat,RefTime,tpe,Nobs,y,obl,y0,Error,
     :        maxobs,verbose)
*         Put the parameters into the appropriate range
          call CheckRange (x,verbose)
        endif

*       Calculate the projected obs from the current estimate of x
*       Use the forward model
        if (verbose) print*,'Calculating the forward model'
        call AstroFM (x,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,
     :                long,lat,RefTime,tpe,Nobs,maxobs,y,obl)
*       Array y contains projected observations (y0 contains actual obs)

        if (verbose) then
          print*,'Here are the observations and projected observations'
          do k=1,Nobs
            print*,int(real(k+1)/2.),180.*y0(k)/pi,180.*y(k)/pi
          enddo
        endif

*       Find deltay
        do k=1,Nobs
          deltay(k)=y0(k)-y(k)
        enddo
*       Due to the periodic nature of angles, deltay calculated in this
*       way may be misleading.  Correct for this here.
        do k=1,Nobs
          if (tpe(k).eq.1) then
*           The measurement is of RA or dec (within range 0-2pi)
*           Need to possibly adjust RA
            if (mod(k,2).eq.1) call PutInRange (deltay(k),pi,0)
          else
*           The measurement is of alt or azi (within range 0-2pi)
*           Need to possibly adjust azi
            if (mod(k,2).eq.1) call PutInRange (deltay(k),pi,0)
          endif 
        enddo
      endif

*     Linearize the forward model (for the Jacobean)
*     This is where we need real*8 instead of real variables
      if (verbose) print*,'Linearizing forward model'
      call LinearFM (x,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,long,lat,
     :               RefTime,tpe,Nobs,maxobs,H,Temp1,Temp2,obl)

      if (verbose) then
*       Output the forward model
        print*,'----- Jacobean ------------------------------------'
        do k=1,Nobs
          print 204,k,(H(k,i),i=1,6)
        enddo
        print*,'---------------------------------------------------'
      endif



*     -----Calculate the Hessian in parameter space--------------------
      if (verbose) print*,'Calculating Hessian - 1st derivs of H'

*     Contribution from the 1st derivatives of H (H(T)R(-1)H)
      do i=1,6
        do j=1,6
          if (j.ge.i) then
*           Upper triangular and diagonal terms
            total=0.
            do k=1,Nobs
              total=total+H(k,i)*H(k,j)/Error(k)
            enddo
            Hess1(i,j)=total
          else
*           Lower diagonal terms
            Hess1(i,j)=Hess1(j,i)
          endif
          Hess(i,j)=Hess1(i,j)
        enddo
        if (verbose) print 200,(Hess1(i,j),j=1,6)
      enddo

      if (d2Hess) then
*       Contribution from the 2nd derivatives of H (if required)
        if (verbose) print*,'Calculating Hessian - 2nd derivs of H'
        call SecondH (x,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,long,lat,
     :                RefTime,tpe,Nobs,maxobs,Hd2,obl,Temp1,Temp2,Temp3,
     :                Temp4)
*       Compute the rest of the Hessian
        do i=1,6
          do j=1,6
            if (j.ge.i) then
*             Upper triangular and diagonal terms
              total=0.
              do k=1,Nobs
                total=total-Hd2(k,i,j)*deltay(k)/Error(k)
              enddo
              Hess2(i,j)=total
            else
*             Lower diagonal terms
              Hess2(i,j)=Hess2(j,i)
            endif
            Hess(i,j)=Hess(i,j)+Hess2(i,j)
          enddo
          if (verbose) print 200,(Hess2(i,j),j=1,6)
        enddo

        if (verbose) then
*         Display full Hessian
          print*,'Total Hessian'
          do i=1,6
            print 200,(Hess(i,j),j=1,6)
          enddo
        endif
      endif


      if (.not.FinalRun) then
*       -----Calculate minus the gradient in parameter space------------
        do i=1,6
          total=0.
          do k=1,Nobs
            total=total+H(k,i)*deltay(k)/Error(k)
          enddo
          Grad(i)=total
        enddo

*       -----Solve equations to find increment in x-space---------------
*       Solve: Hessian deltax = -gradient
*       Use Gaussian ellimination - the increment is deltax
*       Hess1 and v in this call are now just workspace variables
        call GaussEll (Hess,Hess1,Grad,v,deltax,6,6)

*       -----Calculate updated x (and find increment size)--------------
        ms=0.
        do i=1,6
          ms=ms+deltax(i)*deltax(i)
*         Find the updated value of x
          x(i)=x(i)+deltax(i)
        enddo

*       Put the parameters into the appropriate range
        call CheckRange (x,verbose)

        CostJ=CostFunc(y,y0,Error,Nobs,maxobs)
        write (12,201) iter,(x(i),i=1,6),ms,CostJ
        print 201,iter,(x(i),i=1,6),ms,CostJ

*       Check that a and i remain in the specified ranges
        if ((x(1).lt.amin).or.(x(1).gt.amax)) then
          print*,'Fatal error from main program'
          print*,'param 1 is not in the specified range'
          print*,'value : ',x(1)
          print*,'Range : ',amin,amax
          stop
        endif
        if ((x(3).lt.imin).or.(x(3).gt.imax)) then
          print*,'Fatal error from main program'
          print*,'param 3 is not in the specified range'
          print*,'value : ',x(3)
          print*,'Range : ',imin,imax
          stop
        endif
      endif

*     -----Have we converged?-------------------------------------------
      if ((.not.FinalRun).and.(ms.gt.crit).and.(iter.lt.maxiter))
     :   goto 10
      if ((ms.gt.crit).and.verbose)
     :       print*,'Premature exit - maximum No. of iterations reached'
*     Have completed the iterations - Set FinalRun flag for last run
*     through (need to repeat some of algorithm - for error analysis)
      if (.not.FinalRun) then
        FinalRun=.true.
        goto 10
      endif
*     End of iterations
      write (12,203) '#---------------------------------------------'//
     :               '-----------------------------------'


*     -----Find the inverse Hessian (for error calculation)-------------
*     A final and partial trip round the loop (FinalRun true) above
*     calculated the Hessian only
      if (verbose) print*,'Diagonalizing'
      lim=0.00001
      call Jacobi (Hess,Lambda,U,6,6,lim,trace)
*     The eigenvalues are the diagonal elements of Lambda
*     The eigenvectors are the rows of U

      if (verbose) then
        print*,'Here are the eigenvalues of the Hessian'
        print 200,(Lambda(i,i),i=1,6)
        print*,'Here is the matrix of eigenvectors'
        do i=1,6
          print 200,(U(i,j),j=1,6)
        enddo
      endif

*     Check for eigenvalues that show Hessian is not pos. definite
      do i=1,6
        if (Lambda(i,i).le.0.)
     :    print*,'Warning - bad eigenvalue No. ',i,' (',Lambda(i,i),')'
      enddo

*     Compute the inverse Hessian matrix, U(T)lambda(-1)U
      do i=1,6
        do j=1,6
          total=0.
          do k=1,6
            if (Lambda(k,k).ne.(0.))
     :         total=total+U(k,i)*U(k,j)/Lambda(k,k)
          enddo
          InvH(i,j)=total
        enddo
      enddo


      print*,'---------------------------------------------------------'
      write (12,203) '# Analysis errors (negative numbers have negative'
     :              //' curvatures)'
      write (12,205) '#   ',(son(InvH(i,i))*sqrt(abs(InvH(i,i))),i=1,6)
      write (12,203) '#============================================='//
     :               '==================================='
      write (12,203) '# Hessian eigenvalues'
      write (12,206) '# ',(Lambda(i,i),i=1,6)
      write (12,203) '#============================================='//
     :               '==================================='
      close (12)
      print*,'Convergence and analysis details are in file Convergence'









      elseif (rmode.eq.'RECO') then
*     -----------------------------
      if (verbose) then
        print*,'RUN MODE: Reconfigure input file'
        print*,'Any alt/azi details will be transformed to RA/dec'
      endif
      call Reconfig (fname,zeroYr,RefTime,MonthDays)









      elseif (rmode.eq.'TRAJ') then
*     -----------------------------
      if (verbose) then
        print*,'RUN MODE: Generate RA/dec trajectories'
        print*,'Any alt/azi details will be transformed to RA/dec'
      endif
*     Find the min and max time
      call TimeRange (Time,Nobs,low,high,maxobs)
*     Open two files 12: model trajectory 'low' to 'high' (TrajCont)
*                    13: obs and model at obs times (TrajPoints)
      open (12,file='TrajCont')
      open (13,file='TrajPoints')
*     Deal with the continuous trajectory (250 points)
      if (verbose) print*,'Continuous trajectory calculation'
      call TrajCont (low,high,250,x,aE,eE,iE,Omega1E,omega2E,LE,Epoch,
     :               obl,12)
*     Deal with the obs and projected obs
      if (verbose) print*,'Obs and predicted obs calculation'
      call TrajPoint (Time,tpe,y0,Nobs,maxobs,reftime,long,lat,x,aE,eE,
     :                iE,Omega1E,omega2E,LE,Epoch,obl,13)
      close (12)
      close (13)
      if (verbose) print*,'Output in files TrajCont and TrajPoints'


      endif
*     -----











      print*,'Program finished'
200   format (6F14.5)
201   format (I4,F8.3,F9.4,5F10.4,F10.1)
202   format (F12.6,F16.2)
203   format ((A))
204   format (I6,F10.2,F10.3,F14.7,F15.8,F17.10,F14.7)
205   format ((A),F8.3,F9.4,4F10.4)
206   format ((A),6F13.3)
      end








*     ==================================================================
*     ==================================================================

      real*8 function CostFunc (y,y0,Error,Nobs,maxobs)
*     This routine will actually calculate the value of the cost
*     function (x-representation, non linear)
*     y0 observations
*     y projected observations
      implicit none
      integer maxobs,i,Nobs
      real*8 y(maxobs),y0(maxobs),Error(maxobs),total,diff,pi

      pi=4.*atan(1.)
      total=0.
      do i=1,Nobs
        diff=abs(y0(i)-y(i))
*       Put into correct range
        call PutInRange (diff,pi,0)
        total=total+(diff)**2./Error(i)
      enddo
      CostFunc=total/2.
      return
      end


*     ==================================================================
      subroutine ReadObs (tpe,RawObs,Error,Time,c,maxobs,zeroYr,long,
     :                    lat,x,rmode,verbose,amin,amax,imin,imax,crit,
     :                    maxiter,partial,CorrectRefract,d2Hess,
     :                    MonthDays,fname)
*     This subroutine will read a AstroObs file
*     All angles output are in radians (except long and lat - degrees)

      implicit none
      integer maxobs,tpe(maxobs),c,c1,t,d,Yr,Mo,Dy,Hr,Mn,zeroYr,s,l
      integer MonthDays(0:12),switch,maxiter,partial
      real*8 RawObs(maxobs),Time(maxobs),Error(maxobs),err,NumDays,long
      real*8 lat,pi,deg,x(6),amin,amax,imin,imax,crit
      character*75 temp
      character*12 fname
      character*4 rmode
      logical verbose,CorrectRefract,d2Hess

      pi=4.*atan(1.)

      print*,'Please enter the input filename (e.g. AstroObs)'
      read*,fname
      print*,fname

      open (12,file=fname)
12    read (12,100,err=13) temp
      if (index(temp,'END OF HEADER').eq.0) goto 12
      goto 14
13    print*,'Error: No end to input file header.'
      stop
14    read (12,100,err=13) temp
      read (12,100,err=13) temp

*     Longitude
      read (12,*,end=11) s,d,Mn
      long=real(s)*(real(d)+real(Mn)/60.)
      read (12,100) temp
*     Latitude
      read (12,*,end=11) s,d,Mn
      lat=real(s)*(real(d)+real(Mn)/60.)
      read (12,100) temp

*     First guess at solution
      read (12,*) (x(l),l=1,6)
      read (12,100) temp

*     Run mode (XSEC, ANAL, RECO, TRAJ)
      read (12,100) rmode
      read (12,100) temp

*     Verbose
      read (12,*) switch
      verbose=(switch.ne.0)
      read (12,100) temp

*     Ranges (length of semi-major axis)
      read (12,*) amin,amax
      read (12,100) temp
*     (inclination)
      read (12,*) imin,imax
      read (12,100) temp

*     Check that the first guess is within this range
      if ((x(1).lt.amin).or.(x(1).gt.amax)) then
        print*,'Fatal error from subroutine ReadObs'
        print*,'Initial guess for param 1 is not in the specified range'
        print*,'Initial guess : ',x(1)
        print*,'Range : ',amin,amax
        stop
      endif
      if ((x(3).lt.imin).or.(x(3).gt.imax)) then
        print*,'Fatal error from subroutine ReadObs'
        print*,'Initial guess for param 3 is not in the specified range'
        print*,'Initial guess : ',x(3)
        print*,'Range : ',imin,imax
        stop
      endif

*     Convergence cutoff in length of increment
      read (12,*) crit
      read (12,100) temp

*     Maximum number of iterations
      read (12,*) maxiter
      read (12,100) temp

*     Number of iterations performing partial minimization
      read (12,*) partial
      read (12,100) temp

*     Refraction correction
      read (12,*) switch
      CorrectRefract=(switch.ne.0)
      read (12,100) temp

*     Include 2nd derriv of forward model in Hessian?
      read (12,*) switch
      d2Hess=(switch.ne.0)
      read (12,100) temp

      c=-1
10    read (12,100,end=11) temp
      read (12,*,end=11) t
*       Increment the counter
        c=c+2
        c1=c+1
        if (c1.gt.maxobs) then
          print*,'TOO MANY OBSERVATIONS'
          print*,'INCREASE maxobs'
          stop
        endif
*       Store the type flag
        tpe(c)=t
        tpe(c1)=t
*       Find the time
        read (12,*,end=11) Yr,Mo,Dy,Hr,Mn
        Time(c)=NumDays(zeroYr,Yr,Mo,Dy,Hr,Mn,MonthDays)
        Time(c1)=Time(c)
        if (t.eq.1) then
*         The observation is of RA \ dec variety
          read (12,*,end=11) Hr,Mn,err
          RawObs(c)=pi*(real(Hr)+real(Mn)/60.)/12.
          Error(c)=(pi*err/720.)**2.
          read (12,*,end=11) d,Mn,err
          RawObs(c1)=pi*(real(d)+real(Mn)/60.)/180.
          Error(c1)=(pi*err/10800.)**2.

          else

*         The observation is of alt \ azi variety
          read (12,*,end=11) deg,err
          RawObs(c)=pi*deg/180.
          Error(c)=(pi*err/180.)**2.
          read (12,*,end=11) deg,err
          RawObs(c1)=pi*deg/180.
          Error(c1)=(pi*err/180.)**2.
        endif
      goto 10

11    close (12)
      c=c+1

100   format ((A))
      return
      end


*     ==================================================================
      subroutine ReconFig (fname,zeroYr,RefTime,MonthDays)
*     This subroutine will reconfigure an AstroObs file
*     (Change all alt/azi observations to RA/dec type)
*     All angles output are in radians (except long and lat - degrees)

      implicit none
      integer t,d,Yr,Mo,Dy,Hr,Mn,zeroYr,s,l
      integer MonthDays(0:12),verb,maxiter,partial,switch
      real*8 Time,err,NumDays,RefTime,long,RA,dec,crit
      real*8 lat,pi,x(6),amin,amax,imin,imax,alt,azi,erralt,errazi,son
      character*75 temp
      character*4 rmode
      character*12 fname

      pi=4.*atan(1.)

      print*,'Opening input file for reading  : ',fname
      open (12,file=fname)
      print*,'Opening output file for writing : AstroObs1'
      open (13,file='AstroObs1')
      do l=1,28
        read (12,100) temp
        write (13,100) temp
      enddo


*     Longitude
      read (12,*,end=11) s,d,Mn
      long=real(s)*(real(d)+real(Mn)/60.)
      write (13,*) s,d,Mn
      read (12,100) temp
      write (13,100) temp

*     Latitude
      read (12,*,end=11) s,d,Mn
      lat=real(s)*(real(d)+real(Mn)/60.)
      write (13,*) s,d,Mn
      read (12,100) temp
      write (13,100) temp

*     First guess at solution
      read (12,*) (x(l),l=1,6)
      write (13,*) (x(l),l=1,6)
      read (12,100) temp
      write (13,100) temp

*     Run mode (1=cost fn cross sections, 2=minimize cost fn,
*               3=reconfigure input file alt/azi -> RA/dec
*               4=generate RA/dec trajectory)
      read (12,100) rmode
      write (13,100) rmode
      read (12,100) temp
      write (13,100) temp

*     Verbose
      read (12,*) verb
      write (13,*) verb
      read (12,100) temp
      write (13,100) temp

*     Ranges (length of semi-major axis)
      read (12,*) amin,amax
      write (13,*) amin,amax
      read (12,100) temp
      write (13,100) temp
*     (inclination)
      read (12,*) imin,imax
      write (13,*) imin,imax
      read (12,100) temp
      write (13,100) temp

*     Convergence cutoff in length of increment
      read (12,*) crit
      write (13,*) crit
      read (12,100) temp
      write (13,100) temp

*     Maximum number of iterations
      read (12,*) maxiter
      write (13,*) maxiter
      read (12,100) temp
      write (13,100) temp

*     Number of iterations performing partial minimization
      read (12,*) partial
      write (13,*) partial
      read (12,100) temp
      write (13,100) temp

*     Refraction correction
      read (12,*) switch
      write (13,*) switch
      read (12,100) temp
      write (13,100) temp

*     Include 2nd derriv of forward model in Hessian?
      read (12,*) switch
      write (13,*) switch
      read (12,100) temp
      write (13,100) temp

10    read (12,100,end=11) temp
        write (13,100) temp
        read (12,*,end=11) t
        write (13,*) 1

*       Find the time
        read (12,*,end=11) Yr,Mo,Dy,Hr,Mn
        write (13,*) Yr,Mo,Dy,Hr,Mn
        Time=NumDays(zeroYr,Yr,Mo,Dy,Hr,Mn,MonthDays)

        if (t.eq.1) then
*         The observation is of RA \ dec variety - no need to change
          read (12,*,end=11) Hr,Mn,err
          write (13,*) Hr,Mn,err
          read (12,*,end=11) d,Mn,err
          write (13,*) d,Mn,err

          else

*         The observation is of alt \ azi variety - need to convert
          read (12,*,end=11) alt,erralt
          read (12,*,end=11) azi,errazi
          call localinverse (RA,dec,Time,long,lat,reftime,alt,azi)

*         NEED TO DEAL WITH THE ERROR
          err=2.

*         Deal with the RA
          s=int(son(RA))
          RA=abs(RA)
          Hr=s*int(RA)
          RA=60.*(RA-real(abs(Hr)))
          Mn=int(RA)
          write (13,*) Hr,Mn,err
*         Deal with the dec
          s=int(son(dec))
          dec=abs(dec)
          d=s*int(dec)
          dec=60.*(dec-real(abs(d)))
          Mn=int(dec)
          write (13,*) d,Mn,err

        endif
      goto 10

11    close (12)
      close (13)

100   format ((A))
      return
      end


*     ==================================================================
      subroutine LinearFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,long,
     :                      lat,reftime,tpe,Nobs,maxobs,H,Obs1,Obs2,obl)

*     Linearize the forward model about the point x0.
*     The vector x0 has components:
*     x0(1)=a, x0(2)=e, x0(3)=i, x0(4)=Omega1, x0(5)=omega2, x0(6)=L

      implicit none
      integer maxobs,tpe(maxobs),Nobs,param,obs
      real*8 x0(6),aE,eE,iE,Omega1E,omega2E,LE,Time(maxobs),t0,long,lat
      real*8 reftime,H(maxobs,6),delta,Obs1(maxobs),Obs2(maxobs),d
      real*8 obl,x0up(6),x0dn(6)

      do param=1,6
        x0up(param)=x0(param)
        x0dn(param)=x0(param)
      enddo
      delta=0.00005

      do param=1,6
*       (Vary each parameter in turn)
*       Need slightly different values for gradient calculation
        x0up(param)=x0(param)+delta
        x0dn(param)=x0(param)-delta
*       For the eccentricity, check the range
        if (param.eq.2) then
          if (x0up(param).gt.(1.)) x0up(param)=1.
          if (x0dn(param).lt.(0.)) x0dn(param)=0.
        endif
        d=x0up(param)-x0dn(param)

*       Call the forward model with a slightly modified set of params
        call AstroFM (x0up,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                long,lat,reftime,tpe,Nobs,maxobs,Obs1,obl)
        call AstroFM (x0dn,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                long,lat,reftime,tpe,Nobs,maxobs,Obs2,obl)
*       Restore the working list of parameters
        x0up(param)=x0(param)
        x0dn(param)=x0(param)
*       Calculate the components of the derivative
        do obs=1,Nobs
          H(obs,param)=(Obs1(obs)-Obs2(obs))/d
        enddo
      enddo
      return
      end


*     ==================================================================
      subroutine AstroFM (xx,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                    long,lat,reftime,tpe,Nobs,maxobs,ProObs,obl)
*     Astronomical 'forward model'.
*     Given a set of orbital parameters, the times, the location on
*     Earth and the type of observation, this subroutine will generate
*     the observtions which are expected.

      implicit none
      integer maxobs,Nobs,tpe(maxobs),loop,pairs,ind1,ind2
      real*8 xx(6),a,e,i,Omega1,omega2,L,aE,eE,iE,Omega1E,omega2E,LE,t0
      real*8 Time(maxobs),long,lat,ProObs(maxobs),RA,Dec,reftime
      real*8 alt,azi,obl,xcel,ycel,zcel,xecl,yecl,zecl

*     Longitude and latitude should be in degrees
*     Nobs is actually 2*(number of observation co-ordinates)

      a=xx(1)
      e=xx(2)
      i=xx(3)
      Omega1=xx(4)
      omega2=xx(5)
      L=xx(6)

      pairs=Nobs/2
      do loop=1,pairs
        ind1=(loop-1)*2+1
        ind2=ind1+1
        call Orbit (a,e,i,Omega1,omega2,L,aE,eE,iE,Omega1E,omega2E,
     :              LE,Time(ind1),t0,obl,xcel,ycel,zcel,xecl,yecl,zecl)
*       The used o/p is xcel,ycel,zcel (geocentric equatorial co-ords)

        if (tpe(ind1).eq.1) then
*         Observation required is RA \ dec
*         Calculate RA and dec from xcel,ycel,zcel
          call Celestial (xcel,ycel,zcel,RA,Dec)
          ProObs(ind1)=RA
          ProObs(ind2)=Dec
          else
*         Observation required is alt \ azi
          call local (xcel,ycel,zcel,Time(ind1),long,lat,reftime,alt,
     :                azi)
          ProObs(ind1)=alt
          ProObs(ind2)=azi
        endif
      enddo

      return
      end


*     ==================================================================
      subroutine local (xeq,yeq,zeq,time,long,lat,reftime,alt,azi)

*     Convert to local (alt-azimuth) co-ordinates
*     Expecting xeq,yeq,zeq (equ. co-ords), long and lat in degrees
*     alt and azi are output in radians
      implicit none
      real*8 xeq,yeq,zeq,time,long,lat,reftime,reflong,radlat
      real*8 slong,clong,slat,clat,Nhor,Ehor,Zhor,alt,azi,pi,remainder

      pi=4.0*atan(1.0)

*     Determine the effective longitude (in radians)
      reflong=long/360.0+(time-reftime)/0.997269562
      reflong=2.0*pi*remainder(reflong)

*     Workout some trig
      radlat=pi*lat/180.0
      slong=sin(reflong)
      clong=cos(reflong)
      slat=sin(radlat)
      clat=cos(radlat)
*     Convert to local (horizontal) co-ordinates
      Nhor=slat*clong*xeq+slat*slong*yeq+clat*zeq
      Ehor=slong*xeq-clong*yeq
      Zhor=-clat*clong*xeq-clat*slong*yeq+slat*zeq
*     Calculate the altitude (radians)
      alt=atan(Zhor/sqrt(Ehor*Ehor+Nhor*Nhor))
*     Calculate the azimuth (radians)
      if (Nhor.eq.0) then
        if (Ehor.gt.0) then
          azi=pi/2.
          else
          azi=3.*pi/2.
        endif
        else
        azi=atan(Ehor/Nhor)
        if (Nhor.lt.(0.0)) azi=azi-pi
      endif
      if (azi.lt.0) azi=azi+2.*pi
      return
      end


*     ==================================================================
      subroutine localinverse (RA,dec,time,long,lat,reftime,alt,azi)

*     Convert from local horizontal co-ords to geocentric equatorial
*     celestial co-ords.
*     Main inputs   : alt and azi in degrees
*     Main outputs  : RA in hours and dec in degrees
*     Other intputs : time,long,lat,reftime

      implicit none
      real*8 xeq,yeq,zeq,time,long,lat,reftime,reflong,radlat
      real*8 slong,clong,slat,clat,Nhor,Ehor,Zhor,alt,azi,pi,remainder
      real*8 RA,dec,altrad,azirad

      pi=4.0*atan(1.0)

*     Convert alt and azi to radians
      altrad=pi*alt/180.
      azirad=pi*azi/180.
*     Calculate the local Cartesian co-ordinates from alt and azi
      Nhor=cos(altrad)*cos(azirad)
      Ehor=cos(altrad)*sin(azirad)
      Zhor=sin(altrad)

*     Determine the effective longitude (in radians)
      reflong=long/360.0+(time-reftime)/0.997269562
      reflong=2.0*pi*remainder(reflong)

*     Workout some trig
      radlat=pi*lat/180.0
      slong=sin(reflong)
      clong=cos(reflong)
      slat=sin(radlat)
      clat=cos(radlat)
*     Convert from local (horizontal) co-ordinates to equatorial co-ords
      xeq=slat*clong*Nhor+slong*Ehor-clat*clong*Zhor
      yeq=slat*slong*Nhor-clong*Ehor-clat*slong*Zhor
      zeq=clat*Nhor+slat*Zhor

*     Calculate RA and dec from the equatorial co-ordinates
      call Celestial (xeq,yeq,zeq,RA,dec)

*     Convert RA into hours and dec into degrees
      RA=12.*RA/pi
      dec=180.*dec/pi

      return
      end

*     ==================================================================
      subroutine Orbit (a,e,i,Omega1,omega2,L,aE,eE,iE,Omega1E,omega2E,
     :                  LE,t,t0,obl,xcel,ycel,zcel,xecl,yecl,zecl)
*     Returns the xcel,ycel,zcel (geocentric celestial co-ordinates (can
*     compute the RA/dec from these)) and xecl,yecl,zecl (ecliptic
*     co-ordinates) given the orbital parameters of a planet and the
*     Earth at the chosen time.
*     Ross Bannister, December 2000
*     Input : a,e,i,Omega1,omega2,L (orbital parameters of planet)
*             t,t0 (time and Epoch in days)
*             aE,eE,iE,Omega1E,omega2E,LE (orbital parameters of Earth)
*     Output : Geocentric celestial co-ordinates (xcel,ycel,zcel)
*              Ecliptic co-ordinates (xecl,yecl,zecl)

      implicit none
      real*8 a,e,i,Omega1,omega2,L,t,t0,aE,eE,iE,Omega1E,omega2E,LE
      real*8 xp,yp,zp,xE,yE,zE,pi,pi2,mymod,Ldot,LEdot,obl
      real*8 Lrad,Omega1rad,omega2rad,irad,LErad,Omega1Erad,omega2Erad
      real*8 iErad,con,xcel,ycel,zcel,xecl,yecl,zecl

      pi=4.0*atan(1.)
      pi2=2.*pi
*     con calculates Ldot from a (Ldot=con/a**1.5)
      con=129597762.
      Ldot=con/a**1.5
      LEdot=con/aE**1.5

*     Convert some parameters to radians
      Lrad=pi*(L+Ldot*(t-t0)/131490000.0)/180.0
      Lrad=mymod(Lrad,pi2)
      Omega1rad=pi*Omega1/180.0
      omega2rad=pi*omega2/180.0
      irad=pi*i/180.0

      LErad=pi*(LE+LEdot*(t-t0)/131490000.0)/180.0
      LErad=mymod(LErad,pi2)
      Omega1Erad=pi*Omega1E/180.0
      omega2Erad=pi*omega2E/180.0
      iErad=pi*iE/180.0

*     Position of planet
      xp=0
      yp=0
      zp=0
*     Modified planetary co-ordinates
      call PosPlanetary (a,omega2rad,Lrad,e,Omega1rad,xp,yp)

*     Convert to ecliptic co-ordinates
      call ecliptic (xp,yp,zp,Omega1rad,irad)
      xecl=xp
      yecl=yp
      zecl=zp

*     Position of Earth
      xE=0
      yE=0
      zE=0
*     Modified planetary co-ordinates
      call PosPlanetary (aE,omega2Erad,LErad,eE,Omega1Erad,xE,yE)

*     Convert to ecliptic co-ordinates
      call ecliptic (xE,yE,zE,Omega1Erad,iErad)

*     Geocentric co-ordinates of the planet
      call geo (xp,yp,zp,xE,yE,zE)
*     Geocentric equatorial co-ordinates
      call geoequ (yp,zp,obl)
*     xcel,ycel, and zcel are output
      xcel=xp
      ycel=yp
      zcel=zp
      return
      end

*     ==================================================================
      subroutine PosPlanetary (a,longperi,mlong,e,longan,x,y)
*     To calculate the position of a planet in the modified planetary
*     co-ordinates
      implicit none
      integer count
      real*8 a,longperi,mlong,e,longan
      real*8 x,y,M,E1,E2,error,value
      real*8 deriv,change,w,r,argperi,mymod,pi,pi2,Kepler
      real*8 Kderiv

*     Determine pi
      pi=4.0*atan(1.0)
      pi2=2.*pi
*     Calculate the mean anomaly, eq. (7)
      M=mymod(mlong-longperi,pi2)
*     Use Newton-Raphson procedure to calculate the eccentric anomaly
      error=0.00001
*     First guess at answer
      E2=M
      count=0

10    E1=E2
        count=count+1
        value=Kepler(M,E1,e)
        deriv=Kderiv(E1,e)
        change=value/deriv
*       Estimate better value do  eccentric anomaly
        E2=E1-change
      if (.not.(ABS(change).lt.error).or.(count.gt.50)) goto 10

*     Calculate the true anomaly
      w=2*atan(sqrt((1+e)/(1-e))*tan(E2/2))
*     Distance from occupied focus (AU)
      r=a*(1-e*cos(E2))
*     Calculate the argument of the perihelion
      argperi=longperi-longan
*     Evaluate x,y in the modified planetary co-ordinate system (eq. 8)
      x=r*cos(w+argperi)
      y=r*sin(w+argperi)
      return
      end

*     ==================================================================
      real*8 function Kepler(M,E,ecc)
*     Kepler's equation
      implicit none
      real*8 M,E,ecc
      Kepler=M-E+ecc*sin(E)
      return
      end


      real*8 function Kderiv(E,ecc)
*     Derivative of Kepler's equation
      implicit none
      real*8 E,ecc
      Kderiv=ecc*COS(E)-1
      return
      end

*     ==================================================================
      subroutine ecliptic(x,y,z,longan,i)
*     Convert from 'modified planetary' to 'ecliptic' co-ordinates
      implicit none
      real*8 x,y,z,xx,yy,zz,longan,i,csl
      real*8 snl,csi,sni

*     Work out the longitude of the ascending node and the inclination
*     Calculate some trig functions
      csl=cos(longan)
      snl=sin(longan)
      csi=cos(i)
      sni=sin(i)
*     Do the transformation (z=0, so it is not included in the formula)
      xx=csl*x-snl*csi*y
      yy=snl*x+csl*csi*y
      zz=sni*y
      x=xx
      y=yy
      z=zz
      return
      end


      subroutine geo(x,y,z,Earthx,Earthy,Earthz)
*     Ecliptic to geocentric co-ordinates
      implicit none
      real*8 x,y,z,Earthx,Earthy,Earthz
      x=x-Earthx
      y=y-Earthy
      z=z-Earthz
      return
      end

      subroutine geoequ(y,z,obl)
*     Geocentric to equatorial geocentric co-ordinates
*     (The x-co-ordinate is not affected)
      implicit none
      real*8 y,z,obl,cs,sn,yy,zz
      cs=cos(obl)
      sn=sin(obl)
      yy=cs*y-sn*z
      zz=sn*y+cs*z
      y=yy
      z=zz
      return
      end

*     ==================================================================
      subroutine Celestial (x,y,z,RA,dec)
*     IN : x,y,z of celestial object in geocentric equatorial co-ords
*     OUT: RA (in rads) and dec (in rads)
      implicit none
      real*8 x,y,z,RA,dec,pi

      pi=4.0*atan(1.0)
*     Calculate the RA in radians
      if (x.eq.0) then
        if (y.gt.0) then
          RA=pi/2.0
          else
          RA=3.0*pi/2.0
        endif
        else
        RA=atan(y/x)
        if (x.lt.0) RA=RA+pi
        if (RA.lt.0) RA=RA+2.0*pi
      endif

*     Calculate the dec in radians
      dec=atan(z/sqrt(x*x+y*y))

      return
      end

*     ==================================================================
*     ===== Routines associated with date conversion ===================

      real*8 function NumDays (zeroYr,Yr,Mo,Dy,Hr,Mi,MonthDays)
*     Function to calculate the number of days since the zero of time
      implicit none
      integer zeroYr,Yr,Mo,Dy,Hr,Mi,MonthDays(0:12),year
      real*8 days
      logical leapyr

*     Check that the time of interest is > zero time
      if (Yr.lt.zeroYr) then
        print 101,'YOU MUST CHOOSE A TIME WHICH IS AT OR AFTER',zeroYr
        print 101,'USING TIME 00:00 GMT Jan 01 ',zeroYr
        days=0.0
        else

*       Consider the number of days due to difference in whole years
        days=0
        if (Yr.gt.zeroYr) then
          do year=zeroYr,Yr-1
            if (leapyr(year)) then
              days=days+366.0
              else
              days=days+365.0
            endif
          enddo
        endif

*       Now increment the number of days through the year
        days=days+real(MonthDays(Mo-1)+Dy-1)
        if (leapyr(Yr).and.(Mo.gt.2)) days=days+1.0
*       Finally, the fraction of a day
        days=days+(real(Hr)+real(Mi)/60.0)/24.0
      endif

      NumDays=days
101   format ((A),I)
      return
      end


      subroutine Date(zeroYr,days,Yr,Mo,Dy,Hr,Mi,MonthDays)
*     Given a number of days (since start of zeroYr), this routine works
*     out the true date
      implicit none
      integer zeroYr,Yr,Mo,Dy,Hr,Mi,MonthDays(0:12),dyear,Mday
      real*8 days,days1,Mdayprev
      logical fin,leapyr,ly

*     Calculate the number of years since start of zeroYr
      Yr=zeroYr
      days1=days
      fin=.false.

10    continue
        if (leapyr(Yr)) then
          dyear=366
          else
          dyear=365
        endif
        if(days1.ge.real(dyear)) then
          Yr=Yr+1
          days1=days1-real(dyear)
          else
          fin=.true.
        endif
      if (.not.fin) goto 10

      ly=leapyr(Yr)
*     Calculate the month number
      Mo=0
      Mday=0

20    Mo=Mo+1
        Mdayprev=Mday
        Mday=MonthDays(Mo)
        if (ly.and.(Mo.ge.2)) Mday=Mday+1
      if (.not.(real(Mday).gt.days1)) goto 20

      days1=days1-real(Mdayprev)
*     Calculate the day number
      Dy=int(days1)+1
      days1=days1-real(Dy-1)
*     Calculate the hour
      print*,'days = ',days1
      Hr=int(24.0*days1+0.01)
*     Calculate the minute
      Mi=INT(1440.0*days1-60.0*real(Hr)+0.5)
      return
      end


      logical function leapyr (yr)
*     Returns true do  a leap year, false otherwise
      implicit none
      integer yr
      leapyr=((mod(yr,4).eq.0).and.((mod(yr,100).gt.0).or.
     :        (mod(yr,400).eq.0)))
      return
      end


      subroutine PutInRange (angle,range,type)
*     To put angle in range -range to range if type=0
*                              0 to 2*range if type=1
      implicit none
      integer type
      real*8 angle,range,range2,mymod
      range2=2.*range
      if (type.eq.0) then
        angle = mymod(angle+range,range2)-range
      else
        angle = mymod(angle,range2)
      endif
      return
      end


      real*8 function mymod (x,interval)
*     This function puts x into the interval 0<=x<interval
*     Modifed to work with negative numbers
      implicit none
      real*8 x,x1,interval
      if (x.lt.(0.)) then
        x1=x-0.99999*interval
      else
        x1=x
      endif
      mymod=x-interval*real(int(x1/interval))
      return
      end


      real*8 function remainder (num)
*     Remove integer part
      implicit none
      real*8 num,sgn,result,son
      sgn=son(num)
      num=abs(num)
      result=sgn*(num-real(int(num)))
      if (result.lt.(0.0)) result=result+1.0
      remainder=result
      return
      end


*     ==================================================================
      subroutine CheckRange (x,verbose)
*     To check whether the orbital parameters are in physical range
*     The variable is changed if it is possible.
      implicit none
      real*8 x(6),HalfCircle
      logical verbose

      HalfCircle=180.

*     ##### x(1) (a), must be positive #####
      if (x(1).lt.(0.001)) then
        if (verbose) print*,' Parameter 1 should be positive ',x(1)
        x(1)=0.001
        if (verbose) print*,' Setting to ',x(1)
      endif
*     ##### x(2) (e), must be between 0 and 1 #####
      if (x(2).lt.(0.0001)) then
        if (verbose) print*,' Parameter 2 should be positive ',x(2)
        x(2)=0.0001
        if (verbose) print*,' Setting to ',x(2)
      endif
      if (x(2).gt.(0.999)) then
        if (verbose) print*,' Parameter 2 should be < 1 ',x(2)
        x(2)=0.999
        if (verbose) print*,' Setting to ',x(2)
      endif
*     ##### x(3) (i), 0 to 180 #####
      call PutInRange (x(3),HalfCircle,0)
      if (x(3).lt.(0.)) then
        if (verbose) print*,' Parameter 3 should be positive ',x(3)
        x(3)=-1.*x(3)
        x(4)=x(4)+180.
        if (verbose) print*,' Setting param 3 to ',x(3),'and param 4',
     :                      ' to ',x(4)
      endif
*     ##### x(4) (Omega), -180 to 180 #####
      call PutInRange (x(4),HalfCircle,0)
*     ##### x(5) (omega-bar), 0 to 360 #####
      call PutInRange (x(5),HalfCircle,1)
*     ##### x(6) (epsilon), 0 to 360 #####
      call PutInRange (x(6),HalfCircle,1)

      return
      end


*     ==================================================================
      subroutine PartialMin (indx,xa,xb,x,aE,eE,iE,Omega1E,omega2E,LE,
     :                       Time,Epoch,long,lat,RefTime,tpe,Nobs,y,obl,
     :                       y0,Error,maxobs,verbose)
*     Do a patial minimization of the cost function with respect to
*     parameter indx.  The min must be within the bracket (xa,xb)
      implicit none
      integer indx,maxobs,tpe(maxobs),Nobs,nmins,l,indmins(100),mainmin
      integer point
      real*8 xa,xb,x(6),aE,eE,iE,Omega1E,omega2E,LE,Time(maxobs),Epoch
      real*8 long,lat,RefTime,y(maxobs),obl,y0(maxobs),Error(maxobs)
      real*8 J(0:100),x0(6),xxa,xxb,CostFunc
      logical go,verbose

*     Put the current x into the x0 structure
      if (verbose) print*,'The state before partial min. of ',indx
      do l=1,6
        x0(l)=x(l)
        if (verbose) print*,l,x(l)
      enddo

*     Scan the bracket (xa,xb) ...
      do l=0,100
        x0(indx)=xa+real(l)*(xb-xa)/100.
*       Find the observations
        call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,
     :                long,lat,RefTime,tpe,Nobs,maxobs,y,obl)
*       Calculate the cost function for this state
        J(l)=CostFunc(y,y0,Error,Nobs,maxobs)
      enddo

*     ... note all of the min-values
      nmins=0
      do l=1,99
        if ((J(l).lt.J(l-1)).and.(J(l).lt.(J(l+1)))) then
          nmins=nmins+1
          if (nmins.gt.100) then
            print*,'Fatal error from subroutine PartialMin'
            print*,'Too many minima inside the bracket ',xa,xb
            stop
          endif
*         indmins array points to minima in J array
*         mainmin points to el. of indmins which is globally smallest
          indmins(nmins)=l
          if (nmins.eq.1) then
            mainmin=1
          elseif (J(l).lt.J(indmins(mainmin))) then
            mainmin=nmins
          endif
        endif
      enddo

*     Choose deepest minimum and bracket it with neighbouring points
      if (nmins.eq.1) then
*       A single minimum only - use main bracket
        xxa=xa
        xxb=xb
*       Store the fact that a minimum exists
        go=.true.
      elseif (nmins.gt.1) then
*       Number of minima > 1
        if (mainmin.eq.1) then
*         Minimum is leftmost one within bracket
          xxa=xa
          xxb=xa+real(indmins(1)+indmins(2))*(xb-xa)/200.
*         Store the fact that a minimum exists
          go=.true.
        elseif (mainmin.eq.nmins) then
*         Minimum is the rightmost one within bracket
          xxa=xa+real(indmins(nmins)+indmins(nmins-1))*(xb-xa)/200.
          xxb=xb
*         Store the fact that a minimum exists
          go=.true.
        else
*         The minimum is not the first or the last one
          xxa=xa+real(indmins(mainmin)+indmins(mainmin-1))*(xb-xa)/200.
          xxb=xa+real(indmins(mainmin)+indmins(mainmin+1))*(xb-xa)/200.
*         Store the fact that a minimum exists
          go=.true.
        endif
      else
*       Store the fact that no minima exist
        go=.false.
      endif

      if (go) then
        call GoldenMin1D(xxa,xxb,x,indx,y,y0,Error,Nobs,maxobs,aE,eE,iE,
     :           Omega1E,omega2E,LE,Time,Epoch,long,lat,RefTime,tpe,obl)
        else
        if (verbose) then
          print*,'Note from subroutine PartialMin'
          print*,'No minima inside bracket ',xa,xb,' for parameter ',
     :           indx
        endif
      endif

      if (verbose) then
*       Report partial J field and indicate position of minima
        print*,'----------------------------------------------'
        print*,'Partial min over index ',indx
        point=1
        do l=0,100
          if (l.eq.indmins(point)) then
            print 101,l,J(l),'*'
            if (point.lt.100) point=point+1
            else
            print 102,l,J(l)
          endif
        enddo
        print*,'Main minimum is index ',indmins(mainmin)
        print*,'New value of parameter ',indx,' is ',x(indx)
      endif

101   format (I5,F12.4,(A))
102   format (I5,F12.4)

      return
      end


*     ==================================================================
      subroutine GoldenMin1D (ax,cx,x,indx,y,y0,Error,Nobs,maxobs,aE,eE,
     :                        iE,Omega1E,omega2E,LE,Time,Epoch,long,lat,
     :                        RefTime,tpe,obl)
*     Golden ratio minimization in 1d
*     Requires function J(x) defined elsewhere
*     ------------------------------------------------------------------
*     Parameters ax,cx (input: two x-values which bracket the min)
*                x(indx) (inout: location of minimum of J w.r.t. indx)
*                indx (input: the index of the parameter 1-6 varied)
*                y (passive: array to store predicted obs)
*                y0 (input: observations)
*                Error (input: observation errors)
*                Nobs (input: number of observations)
*                maxobs (input: maximum number of observations)
*     There must be a single minimum between ax and cx
*     ------------------------------------------------------------------
*     Ross Bannister, April 2001

      implicit none
      real*8 ratio,prec
      parameter (ratio=0.381966011,prec=0.0001)
      integer maxobs,loop,status,indx,Nobs,tpe(maxobs),l
      real*8 ax,cx,a,b,c,x(6),x0(6),CostFunc,delta,Jb,Jx,xx
      real*8 y(maxobs),y0(maxobs),Error(maxobs),RefTime,obl
      real*8 aE,eE,iE,Omega1E,omega2E,LE,Time(maxobs),Epoch,long,lat

*     Note arrangement of points (min between a and c)
*       a    b         c     (status=0)
*       a          b   c     (status=1)

*     Put the current x into the x0 structure
      do l=1,6
        x0(l)=x(l)
      enddo

      a=ax
      c=cx
      loop=0
*     Calculate intermediate point, b in first instance (put status=0)
      status=0
      b=a+(c-a)*ratio

10    continue
      loop=loop+1

*     Evaluate the cost function at point b
      x0(indx)=b
      call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,
     :              long,lat,RefTime,tpe,Nobs,maxobs,y,obl)
      Jb=CostFunc(y,y0,Error,Nobs,maxobs)

*     Choose a new point, x and determine new bracketing and status
      if (status.eq.0) then

        xx=b+(c-b)*ratio
*       Evaluate the cost function at point xx
        x0(indx)=xx
        call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,
     :               long,lat,RefTime,tpe,Nobs,maxobs,y,obl)
        Jx=CostFunc(y,y0,Error,Nobs,maxobs)

        if (Jb.lt.Jx) then
*         New arrangement is a,b,xx
          delta=c-xx
          c=xx
          status=1
          else
*         New arrangement is b,xx,c
          delta=b-a
          a=b
          b=xx
        endif

        else

        xx=b-(b-a)*ratio
        x0(indx)=xx
        call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,Epoch,
     :               long,lat,RefTime,tpe,Nobs,maxobs,y,obl)
        Jx=CostFunc(y,y0,Error,Nobs,maxobs)

        if (Jb.lt.Jx) then
*         New arrangement is xx,b,c
          delta=xx-a
          a=xx
          status=0
          else
*         New arrangement is a,xx,b
          delta=c-b
          c=b
          b=xx
        endif

      endif
      if (delta.gt.prec) goto 10
      x(indx)=b
      return
      end


*     ==================================================================
      subroutine TimeRange (Time,Nobs,low,high,maxobs)
*     To return the lowest and highest time present in the time array
      implicit none
      integer Nobs,maxobs,loop,pairs,ind1
      real*8 Time(maxobs),low,high

      pairs=Nobs/2
      low=Time(1)
      high=low
      if (pairs.gt.1) then
        do loop=2,pairs
          ind1=(loop-1)*2+1
          if (Time(ind1).gt.high) high=Time(ind1)
          if (Time(ind1).lt.low) low=Time(ind1)
        enddo
      endif

      return
      end


*     ==================================================================
      subroutine TrajCont (t1,t2,points,x,aE,eE,iE,Omega1E,omega2E,LE,
     :                     t0,obl,un)
*     To go over times between t1 and t2 (with 'points' intervals) and
*     find the RA/dec at each time (according to the model).
      implicit none
      integer points,loop,un
      real*8 t1,t2,x(6),aE,eE,iE,Omega1E,omega2E,LE,t0,t,ival,obl
      real*8 xcel,ycel,zcel,xecl,yecl,zecl,RA,dec,pi

      write (un,100) '#======================================='//
     :               '=============================='
      write (un,100) '#                                       '//
     :               '   Ecliptic co-ordinates'
      write (un,100) '#   Time (days)  RA (model) Dec (model) '//
     :               'x (model)  y (model)  z(model)'
      write (un,100) '#---------------------------------------'//
     :               '------------------------------'
      pi=4.*atan(1.)
      ival=(t2-t1)/real(points)
      do loop=1,points
        t=t1+real(loop-1)*ival
*       Find the gerocentric equatorial co-ords at this time
        call Orbit (x(1),x(2),x(3),x(4),x(5),x(6),aE,eE,iE,Omega1E,
     :              omega2E,LE,t,t0,obl,xcel,ycel,zcel,xecl,yecl,zecl)
*       Convert to RA and Dec (radians)
        call Celestial (xcel,ycel,zcel,RA,dec)
*       Convert to hours,degrees
        RA=12.*RA/pi
        dec=180.*dec/pi
        write (un,101) t,RA,dec,xecl,yecl,zecl
      enddo
      write (un,100) '#======================================='//
     :               '=============================='
100   format ((A))
101   format (F15.6,2F12.6,3F10.3)
      return
      end


*     ==================================================================
      subroutine TrajPoint (Time,tpe,y0,Nobs,maxobs,reftime,long,lat,x,
     :                      aE,eE,iE,Omega1E,omega2E,LE,t0,obl,un)
*     Subroutine to output RA/dec obs (conversion is done if obs are of
*     type alt/azi) and projected obs at these times
      implicit none
      integer maxobs,tpe(maxobs),Nobs,loop,pairs,ind1,ind2,un
      real*8 Time(maxobs),y0(maxobs),long,lat,x(6),aE,eE,iE,Omega1E
      real*8 omega2E,LE,t0,obl,RA,dec,alt,azi,RAmodel,decmodel,altmodel
      real*8 azimodel,pi,reftime,deltaRA,deltaDec,deltaAlt,deltaAzi
      real*8 xcel,ycel,zcel,xecl,yecl,zecl,HalfCircle,HalfDay
      character*7 Label(2)

      Label(1)='RA/dec '
      Label(2)='alt/azi'
      HalfCircle=180.
      HalfDay=12.

      write (un,100)
     :    '#==========================================================='
     :  //'==========================================================='

      write (un,100)
     :    '#   Time (days)   RA(o)  Dec(o)   RA(m)  Dec(m)'
     :  //' RA(O-M)Dec(O-M)  alt(o)  azi(o)  alt(m) '
     :  //' azi(m)alt(O-M)azi(O-M) Obs type'
      write (un,100)
     :    '#-----------------------------------------------------------'
     :  //'-----------------------------------------------------------'
      pi=4.*atan(1.0)
      pairs=Nobs/2
      do loop=1,pairs
        print*,loop
        ind1=(loop-1)*2+1
        ind2=ind1+1

*       Deal with the observations
        if (tpe(ind1).eq.1) then
*         RA/dec variety of observation
          RA=12.*y0(ind1)/pi
          dec=180.*y0(ind2)/pi
          call local (cos(y0(ind2))*cos(y0(ind1)),
     :                cos(y0(ind2))*sin(y0(ind1)),
     :                sin(y0(ind2)),Time(ind1),long,lat,reftime,alt,azi)
          else
*         alt/azi variety of observation
          alt=180.*y0(ind1)/pi
          azi=180.*y0(ind2)/pi
          call localinverse (RA,dec,Time(ind1),long,lat,reftime,alt,azi)
        endif

*       Now make a prediction (model state)
        call Orbit (x(1),x(2),x(3),x(4),x(5),x(6),aE,eE,iE,Omega1E,
     :              omega2E,LE,Time(ind1),t0,obl,xcel,ycel,zcel,xecl,
     :              yecl,zecl)
*       Convert to RA and Dec (radians)
        call Celestial (xcel,ycel,zcel,RAmodel,decmodel)
*       Convert to hours,degrees
        RAmodel=12.*RAmodel/pi
        decmodel=180.*decmodel/pi
*       Convert to alt/azi
        call local (xcel,ycel,zcel,Time(ind1),long,lat,reftime,altmodel,
     :              azimodel)
        altmodel=180.*altmodel/pi
        azimodel=180.*azimodel/pi

*       Compute residuals
        deltaRA=RA-RAmodel
        call PutInRange (deltaRA,HalfDay,0)
        deltaDec=dec-decmodel
        call PutInRange (deltaDec,HalfCircle,0)
        deltaAlt=alt-altmodel
        call PutInRange (deltaAlt,HalfCircle,0)
        deltaAzi=azi-azimodel
        call PutInRange (deltaAzi,HalfCircle,0)

*       Now output the data
        write (un,101) Time(ind1),RA,dec,RAmodel,decmodel,
     :                 deltaRA,deltaDec,alt,azi,altmodel,
     :                 azimodel,deltaAlt,DeltaAzi,
     :                 '  ',Label(tpe(ind1))
      enddo
      write (un,100)
     :    '#==========================================================='
     :  //'==========================================================='
100   format ((A))
101   format (F15.6,12F8.2,A2,A7)
      return
      end


*     ==================================================================
      subroutine SecondH (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                    long,lat,reftime,tpe,Nobs,maxobs,Hd2,obl,
     :                    ObsA,ObsB,ObsC,ObsD)
*     Compute the maxobsx6x6 matrix of second derivatives of the
*     forward model
*     December 2002, Ross Bannister

      implicit none
      integer maxobs,tpe(maxobs),Nobs,i,j,ob
      real*8 x0(6),aE,eE,iE,Omega1E,omega2E,LE,Time(maxobs),t0,long,lat
      real*8 reftime,Hd2(maxobs,6,6),obl
      real*8 delta,xA(6),xB(6),xC(6),xD(6)
      real*8 ObsA(maxobs),ObsB(maxobs),ObsC(maxobs),ObsD(maxobs)

      do i=1,6
        xA(i)=x0(i)
        xB(i)=x0(i)
        xC(i)=x0(i)
        xD(i)=x0(i)
      enddo
      delta=0.05

*     Compute the second derivatives of the forward model
*     Exploit symmetry of second derivatives
      do i=1,6
        do j=1,6
          if (j.gt.i) then
*           Upper-triangular components
*           Set-up the states for computation of the second derivatives
            xA(i)=x0(i)+delta
            xA(j)=x0(j)+delta
            xB(i)=x0(i)-delta
            xB(j)=x0(j)+delta
            xC(i)=x0(i)+delta
            xC(j)=x0(j)-delta
            xD(i)=x0(i)-delta
            xD(j)=x0(j)-delta
*           Call the forward model with these incremented state vectors
            call AstroFM (xA,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                long,lat,reftime,tpe,Nobs,maxobs,ObsA,obl)
            call AstroFM (xB,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                long,lat,reftime,tpe,Nobs,maxobs,ObsB,obl)
            call AstroFM (xC,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                long,lat,reftime,tpe,Nobs,maxobs,ObsC,obl)
            call AstroFM (xD,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                long,lat,reftime,tpe,Nobs,maxobs,ObsD,obl)

*           Compute the second derivatives
            do ob=1,Nobs
              Hd2(ob,i,j)=(ObsA(ob)-ObsB(ob)-ObsC(ob)+ObsD(ob))/
     :                    (4.*delta*delta)
            enddo
*           Return states to original values
            xA(i)=x0(i)
            xA(j)=x0(j)
            xB(i)=x0(i)
            xB(j)=x0(j)
            xC(i)=x0(i)
            xC(j)=x0(j)
            xD(i)=x0(i)
            xD(j)=x0(j)
          elseif (j.eq.i) then
*           Diagonal components
*           Set-up the states for computation of the second derivatives
            xA(i)=x0(i)+delta
            xB(i)=x0(i)-delta
*           Call the forward model with these incremented state vectors
            call AstroFM (xA,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                long,lat,reftime,tpe,Nobs,maxobs,ObsA,obl)
            call AstroFM (xB,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                long,lat,reftime,tpe,Nobs,maxobs,ObsB,obl)
            call AstroFM (x0,aE,eE,iE,Omega1E,omega2E,LE,Time,t0,
     :                long,lat,reftime,tpe,Nobs,maxobs,ObsC,obl)

*           Compute the second derivatives
            do ob=1,Nobs
              Hd2(ob,i,j)=(ObsA(ob)-2.*ObsC(ob)+ObsB(ob))/(delta*delta)
            enddo
*           Return states to original values
            xA(i)=x0(i)
            xB(i)=x0(i)
          else
*           Exploit symmetry of second derivatives
            do ob=1,Nobs
              Hd2(ob,i,j)=Hd2(ob,j,i)
            enddo
          endif
        enddo
      enddo

      return
      end


*     ==================================================================
      subroutine GaussEll (A,A1,b,b1,x,N,Nmax)
*     Solve for x in Ax = b (matrix A vector b) (Gaussian ellimination)
*     Ross Bannister, December 2002

      implicit none
      integer N,Nmax,i,ip,j
      real*8 A(Nmax,Nmax),A1(Nmax,Nmax),b(Nmax),b1(Nmax),x(Nmax),factor
      real*8 sum

*     Copy A into A1 and b into b1 (workspace arrays)
      do i=1,N
        do j=1,N
          A1(i,j)=A(i,j)
        enddo
        b1(i)=b(i)
      enddo

*     Stage 1 of the solver : preparation to upper triangular form
      do i=1,N-1
        do ip=i+1,N
          factor=A1(ip,i)/A1(i,i)
*         Do replacements
          do j=i+1,N
            A1(ip,j)=A1(ip,j)-factor*A1(i,j)
          enddo
          b1(ip)=b1(ip)-factor*b1(i)
*         Elliminate terms
          A1(ip,i)=0.
        enddo
      enddo

*     Stage 2 of the solver : find x
      x(N)=b1(N)/A1(N,N)
      do i=N-1,1,-1
        sum=0.
        do j=i+1,N
          sum=sum+A1(i,j)*x(j)
        enddo
        x(i)=(b1(i)-sum)/A1(i,i)
      enddo
*     The correct answer should now be in x

      return
      end


*     ==================================================================
      subroutine Jacobi (A,B,T,N,max,lim,trace)
*     To diagonalize a square symmetric matrix by a sequence of Jacobi
*     rotations
*     Ross Bannister, January 2001 / August 2002
*     A (input) matrix in original representation
*     B (output) matrix in diagonal representation
*     T (output) matrix of eigenvectors (rows)
*     N (input) the order of the matrix
*     max (input) the maximum order of the matrix (for array dim)
*     lim (input) convergence criteria, e.g. lim=0.01
*     trace (output) Sum of eigenvalues

      implicit none
      integer N,max,i,j,p,q,k,sweeps
*      integer lp1,lp2
      real*8 A(max,max),B(max,max),T(max,max),lim,off,diag,theta
      real*8 t1,t2,tt,cc,ss,sc,Bip,Bqj,Tpj,Bpp,Bqq,dis,s,c,norm,trace

      if (N.gt.max) then
        print*,'Report from matrix JACOBI.'
        print*,'N =',N
        print*,'max =',max
        print*,'Please increase the value of max.'
        stop
      endif

*     The matrix B is initially equal to A
      do j=1,N
        do i=1,N
          B(i,j)=A(i,j)
        enddo
      enddo

*     The matrix T is initially the identity matrix
      do j=1,N
        do i=1,N
          if (i.eq.j) then
            T(i,j)=1.0
            else
            T(i,j)=0.0
          endif
        enddo
      enddo

*     Calculate the sum of the square of off-diagonal elements of B
      off=0.0
      do j=2,N
        do i=1,j-1
          off=off+B(i,j)*B(i,j)
        enddo
      enddo
      off=2.0*off

*     Calculate the sum of the square of diagonal elements of B
      diag=0.0
      do i=1,N
        diag=diag+B(i,i)*B(i,i)
      enddo

*     The number of sweeps completed
      sweeps=0
10    continue

*     Loop round each upper-right element (a 'sweep')
      do p=2,N
        do q=1,p-1
          if (abs(B(q,p)).gt.(0.0000001)) then
*           In order to find the new representation of the matrix,
*           work out some details regarding the rotation, U to eliminate
*           element q,p of B (p>q).
            theta=(B(p,p)-B(q,q))/(2.0*B(q,p))
            dis=sqrt(theta*theta+1)
*           Two possible values of t
            t1=-theta+dis
            t2=-theta-dis
*           Choose the smaller value
            if (abs(t1).lt.abs(t2)) then
              tt=t1
              else
              tt=t2
            endif
*           Find the values of s and c
            c=1/sqrt(tt*tt+1)
            s=tt*c
*           Modify the operator UBU(trans) (piecewise)
*           Update only upper-right part of matrix
            cc=c*c
            ss=s*s
            sc=s*c
*           Columns p and q:
            do i=1,p-1
              if (i.ne.q) then
                Bip=B(i,p)
                B(i,p)=c*Bip+s*B(i,q)
                if (i.le.q) B(i,q)=c*B(i,q)-s*Bip
              endif
            enddo
*           Rows p and q:
            do j=q+1,N
              if (j.ne.p) then
                Bqj=B(q,j)
                B(q,j)=c*Bqj-s*B(p,j)
                if (j.gt.p) B(p,j)=c*B(p,j)+s*Bqj
              endif
            enddo
*           Mixed elements
*            B(q,p)=(cc-ss)*B(q,p)+sc*(B(q,q)-B(p,p))
            B(q,p)=0.0

*           New diagonal elements
            Bpp=B(p,p)
            Bqq=B(q,q)
            B(p,p)=cc*Bpp+ss*B(q,q)+2.0*sc*B(p,q)
            B(q,q)=cc*B(q,q)+ss*Bpp-2.0*sc*B(p,q)

*           Modified value of sum of square of off-diagonal elements...
            off=off-2.0*B(p,q)*B(p,q)
*           ...and diagonal elements
            diag=diag-Bpp*Bpp-Bqq*Bqq+B(p,p)*B(p,p)+B(q,q)*B(q,q)

*           Update lower-left of matrix
            do k=1,N
              if (k.lt.p) then
                B(p,k)=B(k,p)
                if (k.lt.q) B(q,k)=B(k,q)
              endif
              if (k.gt.q) then
                B(k,q)=B(q,k)
                if (k.gt.p) B(k,p)=B(p,k)
              endif
            enddo

*           Modify the transformation T(new)=U T(old)
*           Only rows p and q are affected
            do j=1,N
              Tpj=T(p,j)
              T(p,j)=c*Tpj+s*T(q,j)
              T(q,j)=c*T(q,j)-s*Tpj
            enddo
          endif
        enddo
      enddo

      sweeps=sweeps+1
      if (((off/diag).gt.lim).or.(sweeps.lt.3)) goto 10

*     Normalize the eigenvectors
      do i=1,N
        norm=0.0
        do j=1,N
          norm=norm+T(i,j)*T(i,j)
        enddo
        norm=sqrt(norm)
        do j=1,N
          T(i,j)=T(i,j)/norm
        enddo
      enddo

*     Calculate the sum of eigenvalues
      trace=0.
      do i=1,6
        trace=trace+B(i,i)
      enddo

      return
      end


*     ==================================================================
      real*8 function son (input)
*     Returns the sign of the input +1.,-1.,0.
      implicit none
      real*8 input,result
      if (input.gt.(0.)) then
        result=1.
      elseif (input.lt.(0.)) then
        result=-1.
      else
        result=0.
      endif
      son=result
      return
      end
*     ==================================================================
