PRO yopos, timein, latitude, longitude, p0long=p0long, replot=replot, $
    giffile=giffile, zbuffer=zbuffer, outimg=outimg, qdebug=qdebug, $
        plotnight=plotnight
; Coefficients for Equation of Time
ceot = [-0.120553, -0.0323706, 0.000394917, 0.163718]
; Keplers constant for Earth
; (conversion between orbit period and orbit semi-major axis)
kepler = 0.01877
; Earth mean radius (km)
earth_radius = 6371.0

; Station Code for Kagoshima
stcode = BYTE(85)
; Latitude, longitude, of Kagoshima (degrees)
; Latitude = orbital inclination
kaglat = 31.345
kaglong = 131.

; Name of world map file (file location is version dependent)
worldname = 'worldelv.dat'
subdir=(['images',concat_dir('examples','data')])(since_version('4'))
worldfile=(concat_dir(!dir,concat_dir(subdir,worldname)))(0)
if not file_exist(worldfile) then begin
   message,/info,"Cannot find elevation file: " + worldname
   return
endif

; Common to keep track of display window
COMMON YOPOS_WINDOW, winnumber, startx, starty, newworld
; Scale factors for Yohkoh spacecraft image
yoscale = 0.5
ylatscale = 512./(512. - 53.*yoscale)
ylongscale = 512./(512. - 97.*yoscale)

; HANDLE TIME  -----------------------------------------------------------
; Is there an input time?
IF( N_ELEMENTS(timein) GT 0 ) THEN BEGIN
    timstr = ANYTIM2INTS( timein )
    IF( KEYWORD_SET(replot) EQ 0 ) THEN timestr = timein
   ENDIF ELSE BEGIN
    timstr = ANYTIM2INTS( UT_TIME() )
ENDELSE
; Look +/- 12 hours
; slf - fix problem (day boundry)

; GET YOHKOH ORBIT PARAMETERS  ---------------------------------------------
; Inclination never changes
inclin = kaglat
RD_WEEK_FILE, timegrid(timstr,hour=-12),timegrid(timstr,hour=12),'FEM', fem
IF( GET_NBYTES(fem) LT 10) THEN BEGIN
        PRINT,' No FEM data available for this time.'
        RETURN
END
nfem = N_ELEMENTS(fem) - 1

; Sidereal period in minutes
period = ( 1440. * (fem(nfem).day - fem(0).day) + $
        (fem(nfem).time - fem(0).time)/60000. ) /  $
        (fem(nfem).sc_rev - fem(0).sc_rev)
; Synodic period
persynod = 1./(1./period - 1./1440.)
perrat = period/persynod
; Orbit semimajor axis
semimajor = SQRT( period^3/kepler )

; When was Yohkoh over Kagoshima (U.T. seconds)?
nstation = N_ELEMENTS(fem(0).st_station)
timsec = fem.time/1000. + (fem.day - fem(0).day) * 86400.
timsec = TRANSPOSE( REBIN( timsec, nfem+1, nstation) )
contact = WHERE( fem.st$station EQ stcode )
starts = fem.st_station + timsec
ends   = fem.en_station + timsec
duration = ends(contact) - starts(contact)
centers  = (starts(contact) + ends(contact))/2.
ldur = MAX(duration, top)
timref = centers(top)


; Choose map central longitude
IF (n_elements(p0long) GT 0) then $
    plong = p0long $
else $
    plong = 139.                     ; ISAS longitude

; TOP OF REPEAT LOOP
REPEATLOOP:

; Force the size scale
yos = yoscale

; If time is input, translate, otherwise, use the clock
IF (N_ELEMENTS(timestr) EQ 0) THEN BEGIN
        timenow = UT_TIME()
    ENDIF ELSE BEGIN
        timenow = timestr
ENDELSE

timstr = ANYTIM2INTS( timenow )
uthours = timstr.time/3600000.
ptitle = FMT_TIM( timenow ) +  ' U.T.'

; MAKE WINDOW  ------------------------------------------------------------
; Make a 512 x 512 window so we can use YO_LOGO later
; Reuse window if it exists and is the right size
; optionally use Z-buffer (via wdef call)
zbuffer = keyword_set(zbuffer) or keyword_set(giffile)

IF( N_ELEMENTS(winnumber) GT 0 and (1-zbuffer)) THEN BEGIN
    DEVICE, WINDOW=win_arr
    IF(win_arr(winnumber)) THEN BEGIN   ; "old" window exists
            wtemp=!d.window
                wset,winnumber
                if (!d.x_size ne 512) or (!d.y_size ne 512) then wset,wtemp
    ENDIF
        if !d.window ne winnumber then delvarx,winnumber
ENDIF 

if zbuffer then set_plot,'z'
if n_elements(winnumber) eq 0 or zbuffer then WDEF, winnumber, 512, 512, zbuffer=zbuffer

; PLOT MAP  ------------------------------------------------------------
; Set the scaling
left = plong - 180.
right = plong + 180.
LOADCT, 8
TEK_COLOR
MAP_SET, 0., plong, /CYLINDRICAL,  TITLE = ptitle, $
     COLOR=1
; Get the background image
IF( N_ELEMENTS(newworld) LE 0) THEN BEGIN
    world = BYTARR(360,360)
    OPENR, lun, worldfile, /GET_LUN
    READU, lun, world
    CLOSE, lun
    FREE_LUN, lun
    sea = WHERE(world LT 125)
    world(sea) = 4
    south = world(*,0:64)
    hi = WHERE(south ge 125)
    south(hi) = 255
    world(*,0:64) = south
    north = world(250:*,295:*)
    hi = WHERE(north GE 140)
    north(hi) = 255
    world(250:*,295:*) = north
; Remap image
        if since_version('4.') then begin
           newworld = map_patch(world, $
                      xstart=startx,ystart=starty, $ 
                      lon0=0., lon1=360.) 
           message, /info, 'Inside since_version loop' 
        endif else $
          newworld = MAP_IMAGE(world, startx, starty,COMPRESS=1,LONMIN=0., $
              /WHOLE_MAP)
ENDIF
world = 0

TV, newworld, startx, starty

MAP_CONTINENTS, COLOR=9

; PLOT TERMINATOR, DAY/NIGHT  -------------------------------------------
; Call for ephemeris
ephsun = GET_SUN( timstr )
rasun = ephsun(6)
ras = rasun * !PI / 12.
; Equation of time (in hours)
eot = ceot(0)*COS(ras) + ceot(1)*SIN(ras) + ceot(2)*COS(2.*ras) $
         + ceot(3)*SIN(2.*ras)
; Longitude of noon
zenlong = (12. - (uthours + eot) ) * 15.
; Declination of Sun
decsun = ephsun(9)

; Terminator
; Initial coordinate system has Sun in +z direction,
; terminator lies in (x,y) plane.
; Longitude of noon defines zero angle in (x,y) plane
; Rotate by (90 - decsun) around the y axis
theta = FINDGEN(360.) * !DTOR
x = COS(theta)
y = SIN(theta)
phi = (90. - decsun) * !DTOR
xp = x * COS(phi)
yp = y
zp = -x * SIN(phi)
; Convert to (latitude, longitude)
latitude = ASIN(zp) / !DTOR
longitude = ATAN(yp,xp) / !DTOR 
longitude = longitude + zenlong
; Eliminate discontinuity at theta = 180.
; (Changing the (+) to a (-) fixed some day/night problems - BNH)
longitude(180:*) = longitude(180:*) + 360.

; Find plot range, fit longitudes into it
; Since longitudes are in order, just do circular shift
over = WHERE(longitude GT right, ocount)
IF ( ocount GT 0 ) THEN BEGIN
    longitude(over) = longitude(over) - 360.
    longitude = SHIFT(longitude, ocount)
    latitude = SHIFT(latitude, ocount)
ENDIF
under = WHERE(longitude LT left, ocount)
IF(ocount GT 0) THEN longitude(under) = longitude(under) + 360.
; Pad out the array to get complete plot
longt = [longitude, longitude(0)]
lat = [latitude, latitude(0)]


; Plot terminator
OPLOT, longt, lat,  COLOR=0, THICK=2

; Mark which side is night
subsl = SORT(longitude)
nlats = N_ELEMENTS(longitude)
IF(decsun GE 0.) THEN edge = -90. ELSE edge = 90.
if keyword_set(qdebug) then stop

; ------ This Just Doesn't Work. -- (BNH)
if keyword_set(plotnight) then $
  polyfill, [left,left,longitude(subsl),right,right], $
     [edge,latitude(subsl(0)),latitude(subsl),latitude(subsl(nlats-1)),edge], $
         /LINE_FILL, COLOR=0, THICK=1, /DATA

if keyword_set(qdebug) then stop

; PLOT SOUTH ATLANTIC ANOMALY  -----------------------------------------
; Position, extent are:
xsaapos = -25.
ysaapos = -30.
xsaaxtent = 75.
ysaaxtent = 25.
theta = FINDGEN(25) * 15. * !DTOR
xsaa = xsaapos + xsaaxtent * COS(theta)
ysaa = ysaapos + ysaaxtent * SIN(theta)
; Handle spanning the map edge
higher = WHERE(xsaa GE left)
lower = WHERE(xsaa LE left)
IF( (lower(0) LT 0) OR (higher(0) LT 0) ) THEN BEGIN
    POLYFILL, xsaa, ysaa, /LINE_FILL, COLOR=2, /DATA
    ENDIF ELSE BEGIN
    POLYFILL, xsaa(higher), ysaa(higher), /LINE_FILL, COLOR=2, /DATA
    POLYFILL, [left,left,MIN(xsaa(higher)),MIN(xsaa(higher))],      $
        [MAX(ysaa(higher)),MIN(ysaa(higher)),MIN(ysaa(higher)), $
        MAX(ysaa(higher))], /LINE_FILL, COLOR=2, /DATA
    POLYFILL, [xsaa(lower),right,right], [ysaa(lower),MIN(ysaa(lower)), $
        MAX(ysaa(lower))], /LINE_FILL, COLOR=2, /DATA
ENDELSE

; LONGITUDE, LATITUDE LINES
MAP_GRID

; PLOT RADIUS OF KSC CONTACTS -----------------------------------------------
; Use 10 degree elevation angle
anga = !DTOR * (90. + 10.)
angc = ASIN( earth_radius * SIN(anga)/semimajor )
angb = (!PI - anga - angc)/!DTOR
x = kaglong + angb * COS(theta)/COS(kaglat*!DTOR)
y = kaglat + angb * SIN(theta)
OPLOT, x, y, COLOR=7

; PLOT YOHKOH ORBIT AND LOCATION  ----------------------------------------

; Time difference in minutes between map time and S/C overhead at KSC
diftim = (timstr.time/1000. + (timstr.day - fem(0).day)*86400.  - timref)/60.
; Orbital longitude beyond integral number of revolutions
revs = 2. * !PI * ( (diftim / period) MOD 1. )
; Convert orbital longitude to Earth
xr = SIN(revs)
yr = COS(revs) * COS(inclin)
rr = SQRT(xr^2 + yr^2)
erevs = ACOS(yr/rr)/!DTOR
IF( xr LT 0.) THEN erevs = 360. - erevs
; Earth longitude rotation
rotearth = 360. * diftim/1440.

; Spacecraft Earth longitude
sclong = (kaglong + erevs - rotearth) MOD 360. 
IF(sclong LT left) THEN sclong = sclong + 360.
IF(sclong GT right) THEN sclong = sclong - 360.
orblong = (kaglong - rotearth) MOD 360.
IF(orblong LT left) THEN orblong = orblong + 360.
IF(orblong GT right) THEN orblong = orblong - 360.

; 360 degrees of longitude correspond to more than one orbit
; because of Earth rotation
ylat = inclin * COS( !DTOR * (longt - orblong) / perrat )

subso = SORT(longt)
OPLOT, longt(subso), ylat(subso),  THICK=2., COLOR=15
                                                             
; Find  current location
mn = MIN( ABS( longt - sclong ), imn )
sclat = ylat(imn)

; Now tailor to exact window size
; Where is center of plot
pc = CONVERT_COORD( plong, 0., /DATA, /TO_NORM )
pc = pc(0:1)
; Where should the spacecraft be?
sc = CONVERT_COORD( sclong, sclat, /DATA, /TO_NORM )
sc = sc(0:1)
; Increase the offsets because YO_LOGO restricts itself
sc(0) = pc(0) +  (sc(0) - pc(0)) * ylongscale
sc(1) = pc(1) +  (sc(1) - pc(1)) * ylatscale
YO_LOGO, sc(0), sc(1), yos

if data_chk(giffile, /string) then begin
   zbuff2file, giffile, /gif
endif

; BOTTOM OF REPEAT LOOP
IF( KEYWORD_SET(replot) EQ 0) THEN RETURN
IF( replot EQ 0) THEN RETURN

WAIT, 60*replot
GOTO, REPEATLOOP

END