; pro to calculate the distance between each Aeronet site with EPA sites.
; if the distance is smaller than fdt_close, then we consider them are 
; co-located.

 function dis_2pt, lat1_d, lon1_d, lat2_d, lon2_d
   ; function to calculate the distance between two points on earth
   ; if we know lat/lon of those point based on the 'Haversine' formula
   ; http://www.movable-type.co.uk/scripts/latlong.html
   ; based on great circle of two points on a sphere.

   drad=6371.d
   deg2rad=!pi/180.d
   lat1=lat1_d*deg2rad
   lat2=lat2_d*deg2rad
   lon1=lon1_d*deg2rad
   lon2=lon2_d*deg2rad
   dlat=double(lat2-lat1)
   dlon=double(lon2-lon1)
   da=sin(dlat/2.)^2+cos(lat1)*cos(lat2)*(sin(dlon/2.)^2) 
   dc=2.*atan(sqrt(da),sqrt(1.-da))
   ddis=drad*dc

   return, ddis
 end

; read epa sites info

  ssp='PM88500'
  syear='2008'

  openr, iunit, ssp+'_site_'+syear+'.dat', /get_lun

  text=' text'
  readf, iunit, format='(a7,i)', text, ns_e
  print, text, ns_e
  lat_e=fltarr(ns_e) & lon_e=fltarr(ns_e)
  for is= 0, ns_e-1 do begin
    readf, iunit, format='(i4,x,i3,x,i4,x,i5,x,2f10.3,i5)', $
     iis, ist, ict, isd, lat, lon, iline
    ;readf, iunit, iis, ist, ict, isd, lat, lon, iline
    readf, iunit, text
    lat_e(is)=lat & lon_e(is)=lon
  endfor
    
  close, iunit
  free_lun, iunit

; read in aeronet sites info
  lat_a=fltarr(60) & lon_a=lat_a
  openr, iunit, '../aeronet/aeronet_site_latlon_'+syear+'.dat', /get_lun
  k=-1
  while not eof(iunit) do begin
   k=k+1
   readf, iunit, lat, lon
   lat_a(k)=lat
   lon_a(k)=lon
  endwhile
  close, iunit
  free_lun, iunit
  lat_a=lat_a(0:k) & lon_a=lon_a(0:k)
  ns_a=k+1

; open output file
  ;n_sel=n_elements(where(iaer2epa gt 0))
  openw, ounit, 'AERONET_EPA_'+ssp+'_closeby_sites_'+syear+'.dat', /get_lun

; check distance between sites
  ddis_close=8. ; 8km
  iaer2epa=intarr(ns_a) & iaer2epa(*)=-999
  ddis_min=dblarr(ns_a)
  for is_a=0, ns_a-1 do begin
   ddis_min(is_a)=100000.
   for is_e=0, ns_e-1 do begin
     dd=dis_2pt(lat_a(is_a), lon_a(is_a), lat_e(is_e), lon_e(is_e))
     if (dd lt ddis_min(is_a)) then begin
        ddis_min(is_a)=dd
        if (ddis_min(is_a) lt ddis_close) then begin
          printf,ounit,format='(i4, 2f12.3, i5, 2f12.3, f12.3)', $
            is_a,lat_a(is_a),lon_a(is_a), $
            is_e,lat_e(is_e),lon_e(is_e), $
            ddis_min(is_a)
          ;print,'close enough!'
          ;print, is_a,lat_a(is_a),lon_a(is_a)
          ;print, is_e,lat_e(is_e),lon_e(is_e)
          if(iaer2epa(is_a) lt 0) then begin
            iaer2epa(is_a)=is_e
          endif else begin
            print, 'closer!', is_a, lat_a(is_a),lon_a(is_a)
            ;stop
            print, 'previous',lat_e(iaer2epa(is_a)),lon_e(iaer2epa(is_a))
            print, dis_2pt(lat_a(is_a),lon_a(is_a), $
                           lat_e(iaer2epa(is_a)),lon_e(iaer2epa(is_a)))
            print, 'now',lat_e(is_e),lon_e(is_e)
            print, dis_2pt(lat_a(is_a),lon_a(is_a),lat_e(is_e),lon_e(is_e))
            iaer2epa(is_a)=is_e
          endelse
        endif
     endif
   endfor
   ;stop
  endfor

   ;for is_a =0, ns_a-1 do begin
   ; if(iaer2epa(is_a) gt 0) then begin
   ;   printf,ounit,format='(i4, 2f12.3, i5, 2f12.3, f12.3)', $
   ;         is_a,lat_a(is_a),lon_a(is_a), $
   ;         iaer2epa(is_a),lat_e(iaer2epa(is_a)),lon_e(iaer2epa(is_a)), $
   ;         ddis_min(is_a)
   ; endif
   ;endfor
  close,ounit
  free_lun,ounit

  end