FUNCTION cost_fn(ifov, isite)
      IMPLICIT NONE

      REAL*8 cost_fn
      INTEGER*4 ifov, isite

      INCLUDE "ovpindex.prm"
      INCLUDE "geom.cmn"
      INCLUDE "tomsframe.cmn"
      INCLUDE "sites.cmn"
      INCLUDE "instdata.cmn"

      REAL*8 radeg      /57.295779513082323/    !deg-rad conversion

      REAL*8  lat_s, lon_s, lat_t, lon_t, dil, d_ns, d_ew, sum
      REAL*8  box_lat(4), box_lon(4), csatlat, sinxi, cosxi
      REAL*8  tpenalty
      INTEGER*4 isza, ascdesc


c--Copy lats and lons to local variables

      lat_s= lat_site(isite)
      lon_s= lon_site(isite)

      lat_t= lat(ifov)
      lon_t= long(ifov)


c--Handle things differently near a pole.

      IF (ABS(lat_s) .GT. 85.) GOTO 100


c--Calculate vector from site (s) to toms (t) FOV centroid, in km

      dil= 110.*COS(lat_t/radeg)

      d_ew= (lon_t - lon_s)*dil
      d_ns= (lat_t - lat_s)*110.


c--Rotate the FOV box according to the satellite latitude (use lat of
c  the central FOV to approximate this.  
c  When it's changing fastest (at the equator)
c  the rotation is changing the slowest, as a function of latitude. Thus
c  we won't be far off. 
c

      ascdesc= errflg(ifov)/10          ! 1 for descending, 0 for ascending

      csatlat= COS(lat((i_nscans+1)/2)/radeg)
      IF(csatlat .GT. cosi) THEN
        sinxi= cosi/csatlat
      ELSE
        sinxi= 1.d0
      END IF

      cosxi= SQRT(1.d0 - sinxi*sinxi)
      IF(ascdesc .EQ. 1) cosxi= - cosxi

c--The following loop is unrolled for speed

c     DO icorn= 1, 4
c       box_lon(icorn)= cosxi*fov_lon(icorn,ifov) - 
c    &                  sinxi*fov_lat(icorn,ifov)
c       box_lat(icorn)= sinxi*fov_lon(icorn,ifov) + 
c    &                  cosxi*fov_lat(icorn,ifov)
c     END DO

        box_lon(    1)= cosxi*fov_lon(    1,ifov) -
     &                  sinxi*fov_lat(    1,ifov)
        box_lat(    1)= sinxi*fov_lon(    1,ifov) +
     &                  cosxi*fov_lat(    1,ifov)

        box_lon(    2)= cosxi*fov_lon(    2,ifov) -
     &                  sinxi*fov_lat(    2,ifov)
        box_lat(    2)= sinxi*fov_lon(    2,ifov) +
     &                  cosxi*fov_lat(    2,ifov)

        box_lon(    3)= cosxi*fov_lon(    3,ifov) -
     &                  sinxi*fov_lat(    3,ifov)
        box_lat(    3)= sinxi*fov_lon(    3,ifov) +
     &                  cosxi*fov_lat(    3,ifov)

        box_lon(    4)= cosxi*fov_lon(    4,ifov) -
     &                  sinxi*fov_lat(    4,ifov)
        box_lat(    4)= sinxi*fov_lon(    4,ifov) +
     &                  cosxi*fov_lat(    4,ifov)


c--Now compute the sum of the squares of the distances from (s) to each
c  of the corners

      sum=0.d0
c     DO icorn= 1, 4
c       sum= sum + (box_lon(icorn) - d_ew)**2 +
c    &             (box_lat(icorn) - d_ns)**2
c     END DO

c--unrolled loop

        sum= (box_lon(1) - d_ew)**2 +
     &       (box_lat(1) - d_ns)**2 +
     &       (box_lon(2) - d_ew)**2 +
     &       (box_lat(2) - d_ns)**2 +
     &       (box_lon(3) - d_ew)**2 +
     &       (box_lat(3) - d_ns)**2 +
     &       (box_lon(4) - d_ew)**2 +
     &       (box_lat(4) - d_ns)**2




c--Return cost function equal to the sum times a precomputed weighting
c  factor for SZA

      isza= INT(sza(ifov))
      IF (isza .GT. 90) isza= 90

      cost_fn= sum*szweight(isza)

      RETURN


c--If the site is near a pole, modify the cost function to be smaller closer
c  to 0h gmt, and smaller closer to the pole


  100 CONTINUE

        tpenalty= MOD(time+43200 ,86400) - 43200
        tpenalty= (tpenalty*tpenalty)/1.e2
        tpenalty= 1.
        cost_fn=  tpenalty*(90. - ABS(lat_t))

      RETURN
      END