pro coincidence,date1,time1,lat1,lon1,date2,time2,lat2,lon2,  $
                   iopt,multi,maxtime,maxlat,maxlon,maxgcd, $
                   w_time,w_lat,w_lon,w_gcd, $
                   coin1,coin2

;----------------------------------------------------------------------
;
; Purpose: 
;   Routine takes two sets of measurements,  and finds the coincident 
;   measurements.  The user specifies coincidence criteria as either:
;
;      1) time-great circle distance separations 
;      2) time-latitude-longitude separations
;
;   The output is 2 arrays of indices,  that line up the coincident measurements.
;
;   Routine allows option to only use an observation for one coincidence,  
;   or to find multiple coincidences
;
; Input:
;   date1....date, yyyyddd, where ddd is julian day, for data set 1
;   time1....time, decimal hours,  for data set 1
;   lat1.....latitude, + is N, - is S,  "
;   lon1.....longitude, deg E            "
;
;   date2....date, yyyyddd, where ddd is julian day, for data set 2
;   time2....time, decimal hours,  for data set 2
;   lat2.....latitude, + is N, - is S,  "
;   lon2.....longitude, deg E            "
;
;   iopt.....1 = use distance/time,  2 = lat/lon/time
;   multi....1 = do not allow multiple coincidences
;            2 = find multiple coincidences (use a measurment more than once)
;   maxtime..maximum time separation, hours
;   matlat...maximum latitude separation, deg
;   maxlon...maximum longitude separation, deg
;   maxgcd...maximum great circle distance (gcd) separation,  km
;
;   w_time...weighting on importance of time, 0-1, 1 is most important
;   w_lat....     "     "     "       " lat                 "
;   w_lon....     "     "     "       " lon                 "
;   w_gcd....     "     "     "       " gcd                 "
;
; Output:
;   coin1...array indices where data set 1 matched data set 2
;   coin2...array indices where data set 2 matched data set 1
;
;   e.g.,  date1(coin1(0)) was coincident with date2(coin2(0))
;          date1(coin1(1)) was coincident with date2(coin2(1))
;          etc....
;     
; Source:  Mark Hervig 
;
;----------------------------------------------------------------------

;- initialize some variables

   nrec1  = n_elements(date1)
   nrec2  = n_elements(date2)

   nrec   = min([nrec1,nrec2])

   coin1  = lonarr(nrec)
   coin2  = lonarr(nrec)

   use2      = fix(findgen(nrec2)) ; array locations of records not used from "2"
   use2      = lonarr(nrec2) - 1   ; array locations of records used from "2"
   ncoin     = 0
   sum_dlat  = 0.
   sum_dlon  = 0.
   sum_gcd   = 0.
   sum_dt    = 0.
   dlat      = 0.
   dlon      = 0.
   gcd       = fltarr(nrec2)

;- make sure longitudes are all East (0-360)

   k = where(lon1 lt 0,n) & if (n gt 0) then lon1(k) = 360.+lon1(k)
   k = where(lon2 lt 0,n) & if (n gt 0) then lon2(k) = 360.+lon2(k)
 
;- compute elapsed time as hours since midnight on 1-jan-minyr

   minyr = long( min( [date1,date2] ) / 1000)
            
;                days in year        +       yrs. since minyr            + hrs since midnight

   dt1 = ( date1 - long(date1/1000)*1000 + 365.0 *long(date1/1000-minyr) ) * 24. + time1
   dt2 = ( date2 - long(date2/1000)*1000 + 365.0 *long(date2/1000-minyr) ) * 24. + time2

;- look for coincident measurements

   for i1 = 0,nrec1-1 do begin     ; loop over data set 1

     got1     = 0 
     c2       = -1
     deltamin = 99999.
     dt       = abs(dt1(i1) - dt2)

;-   use GCD and time

     if (iopt eq 1) then begin

       great_circle, lat1(i1),lon1(i1),lat2,lon2,gcd,head 

       delta = w_time*dt/maxtime + w_gcd*gcd/maxgcd

       ok    = where(gcd le maxgcd and dt le maxtime, nok)

       if (nok eq 1) then c2   = ok(0)
       if (nok gt 1) then k    = where( delta(ok) eq min(delta(ok)) ) 
       if (nok gt 1) then c2   = ok(k(0))
       if (multi ne 1 and c2 ne -1) then got1 = 1

       if (multi eq 1 and c2 ne -1) then begin
         nk = 0
         if (nok gt 0 and ncoin gt 0) then k = where(c2 eq use2,nk) ; see if that record has been used
         if (nk  eq 0) then got1 = 1
       endif

     endif

;-   use lat, lon, and time

     if (iopt eq 2) then begin

       dlat = abs(lat1(i1) - lat2)
       dlon = abs(lon1(i1) - lon2)

       k     = where(dlon gt 180.0,nk)             ; "wrap around" longitudes 
       if (nk gt 0) then dlon(k) = 360.0-dlon(k) 

       delta = w_lat*dlat/maxlat + w_lon*dlon/maxlon + w_time*dt/maxtime

       ok    = where(dlat le maxlat and dlon le maxlon and dt le maxtime, nok)

       if (nok eq 1) then c2   = ok(0)
       if (nok gt 1) then k    = where( delta(ok) eq min(delta(ok)) ) 
       if (nok gt 1) then c2   = ok(k(0))
       if (multi ne 1 and c2 ne -1) then got1 = 1

       if (multi eq 1 and c2 ne -1) then begin
         nk = 0
         if (nok gt 0 and ncoin gt 0) then k = where(c2 eq use2,nk) ; see if that record has been used
         if (nk  eq 0) then got1 = 1
       endif

     endif

;-   if coincidence was found,  track the records

     if (got1 eq 1) then begin

;-     add measurement "2" to list already used from data set 2

       if (multi eq 1) then begin
         use2(ncoin) = c2
       endif

       if iopt eq 1 then sum_gcd   = sum_gcd  + gcd(c2)   ; sum these for stats
       if iopt eq 2 then sum_dlat  = sum_dlat + dlat(c2)
       if iopt eq 2 then sum_dlon  = sum_dlon + dlon(c2)
                         sum_dt    = sum_dt   + dt(c2)

       coin1(ncoin) = i1
       coin2(ncoin) = c2
       ncoin        = ncoin + 1

     endif

   endfor    ; loop over data set 1

;- finish up

   print ,' number of coincidences ',ncoin

   if (ncoin gt 0) then begin

     coin1 = coin1(0:ncoin-1)
     coin2 = coin2(0:ncoin-1)

     ave_dlat = sum_dlat / float(ncoin)
     ave_dlon = sum_dlon / float(ncoin)
     ave_dt   = sum_dt   / float(ncoin)
     ave_gcd  = sum_gcd  / float(ncoin)

     print ,' average delta lat ',ave_dlat
     print ,' average delta lon ',ave_dlon
     print ,' average distance  ',ave_gcd
     print ,' average delta t   ',ave_dt

   endif

   return
   end