      module blnclc_mod

        use define_star

        implicit none

        contains
        
        subroutine blnclc(ntoe, mmatch, difmag, diff, stars1, &
        istar1, stars2, istar2, rotate, xoff, yoff, theoff, relmag, &
        imatch1, imatch2)

!       Finds the offsets between two images using two lists of stars.

!       The way to change system1 to system2 is 
!       r1**2 = stars1%x**2 + stars2%x**2
!       tan(theta1)=stars1%y/stars1%x
!       stars2%x=r1*cos(theta1+theoff)+xoff
!       stars2%y=r1*sin(theta1+theoff)+yoff.
!       Now, this is equivalent to
!       stars2%x=stars1%x*cos(theoff)-stars1%y*sin(theoff)+xoff
!       stars2%y=stars1%x*sin(theoff)+stars1%y*cos(theoff)+yoff.
!       Note the sign change change with respect to convention (there for
!       historical compatability).

!       By Tim Naylor.

!       The maximum number of stars that will be matched.  If less than this
!       are matched it is changed on output to the actual number matched.
        integer, intent(inout) :: ntoe
!       The maximum number of stars to use in the matching process.
        integer mmatch
!       The largest difference in magnitude to be tolerated.
        real difmag
!       The largest difference in x or y position to be tolerated.
        real diff
!       Number of stars.
        integer istar1, istar2
!       The stars.
        type(a_star), dimension(:) :: stars1, stars2
!       Should the code allow frame rotation?
        logical rotate

!       Output variables.
        ! The offsets in X, Y, theta and relative magnitude.
        real, intent(out) :: xoff, yoff, theoff, relmag
        ! The index numbers of the pairs of matched stars.
        integer, optional, allocatable, dimension(:) :: imatch1, imatch2

!       Locals.
        integer match, mtoe
        real :: xoff_try, yoff_try, theoff_try, relmag_try
        integer ifr1s1, ifr2s1, ifr1s2, ifr2s2, ifr1s3, ifr2s3
        real fr1x, fr1y, r1, theta1, fr2x, fr2y, r2, theta2
        real xpos3, ypos3, r3, theta3
        real xshift, yshift, dshift
        integer, dimension(ntoe) :: jmatch1, jmatch2

        logical debug
        common /bug/ debug

        if (present(imatch1)) allocate(imatch1(ntoe))
        if (present(imatch2)) allocate(imatch2(ntoe))

        if (.not. rotate) theoff=0.0

380     mtoe=0
!       Just take the stars in pairs, starting with the brightest, and
!       look for similar separations.
400     do 440 ifr1s1=1, min(istar1, mmatch)
          !print*, ifr1s1, mtoe
          do 430 ifr1s2=1, min(istar1, mmatch)
!           Are these the same star, if so skip this time.
            if (ifr1s1 .eq. ifr1s2) goto 430
!           Find the difference in magnitude between the two stars in frame 1.
            relmag_try = stars1(ifr1s1)%col(1)%data-stars1(ifr1s2)%col(1)%data
!           Find the difference in position between stars 1 and 2 in frame 1.
            fr1x=stars1(ifr1s1)%x-stars1(ifr1s2)%x
            fr1y=stars1(ifr1s1)%y-stars1(ifr1s2)%y
            r1=sqrt(fr1y*fr1y+fr1x*fr1x)
!           Compare this with a series of star pairs in frame 2.
            do 420 ifr2s1=1, min(istar2, mmatch)
              do 410 ifr2s2=1, min(istar2, mmatch)
!               Are these the same star, if so skip this time.
                if (ifr2s1 .eq. ifr2s2) goto 410
!               Is there a large difference in mags, if so skip this time,
                if (abs(stars2(ifr2s1)%col(1)%data-stars2(ifr2s2)%col(1)%data-relmag_try) .gt. &
                difmag) goto 410
                fr2x=stars2(ifr2s1)%x-stars2(ifr2s2)%x
                fr2y=stars2(ifr2s1)%y-stars2(ifr2s2)%y
                r2=sqrt(fr2y*fr2y+fr2x*fr2x)
                if (abs(r1-r2) .lt. diff) then
!                 A fair match, so go on to find a third star.
                  if (rotate) then
!                   First find the angle between the frames.
                    if ( abs(stars1(ifr1s1)%y-stars1(ifr1s2)%y) &
                    .lt. 1.0e-30) goto 410
                    if ( abs(stars1(ifr1s1)%x-stars1(ifr1s2)%x) &
                    .lt. 1.0e-30) goto 410
                    if ( abs(stars2(ifr2s1)%x-stars2(ifr2s2)%x) &
                    .lt. 1.0e-30) goto 410
                    if ( abs(stars2(ifr2s1)%x-stars2(ifr2s2)%x) &
                    .lt. 1.0e-30) goto 410
                    theta1=atan2( (stars1(ifr1s1)%y-stars1(ifr1s2)%y), &
                                (stars1(ifr1s1)%x-stars1(ifr1s2)%x) )
                    theta2=atan2( (stars2(ifr2s1)%y-stars2(ifr2s2)%y), &
                                (stars2(ifr2s1)%x-stars2(ifr2s2)%x) )
                    theoff_try=theta2-theta1

!                   Now calculate the position angle in frame 1.
                    theta1=atan2(stars1(ifr1s1)%y,stars1(ifr1s1)%x)
                    r1=sqrt( stars1(ifr1s1)%x*stars1(ifr1s1)%x + &
                    stars1(ifr1s1)%y*stars1(ifr1s1)%y )
                    xoff_try=stars2(ifr2s1)%x-r1*cos(theta1+theoff_try)
                    yoff_try=stars2(ifr2s1)%y-r1*sin(theta1+theoff_try)
                  else
                    xoff_try=stars2(ifr2s1)%x-stars1(ifr1s1)%x
                    yoff_try=stars2(ifr2s1)%y-stars1(ifr1s1)%y
                    theoff_try=0.0
                  end if

!                 Set relmag to the difference between the frames.
                  relmag_try = stars2(ifr2s1)%col(1)%data - stars1(ifr1s1)%col(1)%data
                  if (debug) then 
                    print*, '@ Preliminary offsets are;'
                    print*, '@ X2 = X1 +', xoff_try
                    print*, '@ Y2 = Y1 +', yoff_try             
                    if (rotate) print*, '@ theta2 = theta1 + ', theoff_try
                    print*, '@ and frame mag difference ', relmag_try
                    print*, '@ From star ', stars1(ifr1s1)%id
                    print*, '@ and ', stars1(ifr1s2)%id
                  end if
!       How many stars can we match at these offsets?
500     match=0
        xshift=0.0
        yshift=0.0
        dshift=0.0
505     find_other_matches: do ifr1s3=1, min(istar1,mmatch)
          if (rotate) then
            r3 = sqrt(stars1(ifr1s3)%x*stars1(ifr1s3)%x &
            + stars1(ifr1s3)%y*stars1(ifr1s3)%y)
            if (abs(stars1(ifr1s3)%y) .lt. 1.0e-30) cycle find_other_matches
            if (abs(stars1(ifr1s3)%x) .lt. 1.0e-30) cycle find_other_matches
            theta3=atan2(stars1(ifr1s3)%y,stars1(ifr1s3)%x)
            xpos3=r3*cos(theta3+theoff_try)+xoff_try
            ypos3=r3*sin(theta3+theoff_try)+yoff_try
          else
            xpos3=stars1(ifr1s3)%x+xoff_try
            ypos3=stars1(ifr1s3)%y+yoff_try
          end if
          do 510 ifr2s3=1, min(istar2,mmatch)
            if (abs(xpos3-stars2(ifr2s3)%x) .lt. diff) then
              if (abs(ypos3-stars2(ifr2s3)%y) .lt. diff) then
                if (abs(stars2(ifr2s3)%col(1)%data-stars1(ifr1s3)%col(1)%data-relmag_try) .lt. &
                  difmag) then
                  if (abs(stars2(ifr2s3)%y) .lt. 1.0e-30) goto 510
                  if (abs(stars2(ifr2s3)%x) .lt. 1.0e-30) goto 510
                  match=match+1
                  if (rotate) then
                    xshift=xshift+stars2(ifr2s3)%x-(xpos3-xoff_try)
                    yshift=yshift+stars2(ifr2s3)%y-(ypos3-yoff_try)
                  else
                    xshift=xshift+stars2(ifr2s3)%x-stars1(ifr1s3)%x
                    yshift=yshift+stars2(ifr2s3)%y-stars1(ifr1s3)%y
                  end if
                  dshift=dshift+stars2(ifr2s3)%col(1)%data-stars1(ifr1s3)%col(1)%data
                  if (debug) then
                    print*, '@ Confirmed by star', stars1(ifr1s3)%id
                  end if
                  jmatch1(match)=ifr1s3
                  jmatch2(match)=ifr2s3
                  ! Have we reached the required number of macthing stars?
                  if (match == ntoe) exit find_other_matches
!                 Found a match for the star in frame 1, so pick another
!                 star in frame 1.
                  cycle find_other_matches
                end if
              end if
            end if
510       continue
        end do find_other_matches
        ! How many did we match?
        if (mtoe == 0) then
          print*, '  Xoff    Yoff   Theta   Delta_mag   N_stars'
        end if
        if (match > mtoe) then  
          ! Update the offsets, but only if matched more stars.  (If
          ! the same number have been matched as before the earlier
          ! result is probably better as based on brighter stars.)
          xoff=xshift/real(match)
          yoff=yshift/real(match)
          relmag=dshift/real(match)
          theoff=theoff_try
          if (present(imatch1)) imatch1(1:match)=jmatch1(1:match)
          if (present(imatch2)) imatch2(1:match)=jmatch2(1:match)
          print*, xoff, yoff, theoff, relmag, match 
        end if
        mtoe=max(mtoe, match)
        if (mtoe == ntoe) goto 600
                end if
410           continue
420         continue
430       continue  
440     continue

!       Failed to match more than mtoe stars.
        if (mtoe .gt. 2) then
          ntoe=mtoe
        else if (mtoe .eq. 2) then
!         Nothing to confirm the offsets, but let's use them.
          ntoe=mtoe
          xoff=xshift/2.0
          yoff=yshift/2.0
        else
!         MAPCCD has failed!
          ntoe=0
        end if

600   end subroutine blnclc

    end module blnclc_mod
