;+
; NAME:
;       Terminator_Map
;
; PURPOSE:
;
;       This is a program for viewing a map of the Earth with a day/night terminator.
;
; AUTHOR:
;
;       FANNING SOFTWARE CONSULTING
;       David Fanning, Ph.D.
;       1645 Sheely Drive
;       Fort Collins, CO 80526 USA
;       Phone: 970-221-0438
;       E-mail: davidf@dfanning.com
;       Coyote's Guide to IDL Programming: http://www.dfanning.com
;
; CATEGORY:
;
;       Graphics
;
; CALLING SEQUENCE:
;
;       Terminator_Map, center_lon, center_lat, Map_Projection='Aitoff'
;
; INPUT ARGUMENTS:
;
;       center_lon:      The longitude, in degrees, of the center of the map, (-180 to 180). Defaults to 0.0.
;
;       center_lat:      The latitude, in degrees, of the center of the map, (-90 to 90). Defaults to 0.0.
;                        Latitudes other than zero can cause program problems. See the NOTES section for details.
;
; INPUT KEYWORDS:
;
;       FLIPDAY:         The program suffers from an inability to always correctly distinguish the
;                        "day" side of the day/night terminator (see NOTES section). If the program
;                        gets it wrong, you can use this keyword to "flip" the day/night area. For
;                        example, here is a call that gets it wrong:
;
;                            Terminator_Map, 0, 90.0, MAP_PROJECTION=4, time="Thu Apr 13 08:49:34 2006"
;
;                        Unfortunately, you can only know if it is "wrong" by visual inspection. The sun
;                        is always placed correctly. I realize this limits the usefulness of the program,
;                        generally, but it is a limitation I am not motivated enough to fix. That is to
;                        say, the program meets my limited needs perfectly. :-)
;
;       KEEP_ASPECT:     If set, the aspect ratio of the map is more or less preserved.
;
;       MAP_PROJECTION:  Either the number or the name of one of the following map projection types. Defaults to "CYLINDRICAL".
;                        0 CYLINDRICAL
;                        1 GOODESHOMOLOSINE
;                        2 HAMMER
;                        3 LAMBERT
;                        4 MERCATOR
;                        5 MILLER_CYLINDRICAL
;                        6 MOLLWEIDE
;                        7 ROBINSON
;
;      TIME:             A date/time string of nearly any format. The month, however, should be represented
;                        as a three-letter string. The time is written as hh:mm:ss. Examples of date/time strings
;                        are these:
;
;                        'Wed Apr 12 22:39:52 2006'
;                        '12 APR 2006 04:16:00'
;                        'Apr 12, 2006 04:16:00'
;                        'Nov 12, 2006'
;
;                        If undefined, the current day and time (from SYSTIME) is used.
;
; OUTPUT KEYWORDS:
;
;       IMAGE:           If passed a named variable, will return a 24-bit image of the result in the display window.
;
; DEPENDENCIES:
;
;       Requires both the Coyote Library and the JHUAPL Library to be on the !PATH, since many programs
;       from both are used.
;
;           http://www.dfanning.com/programs/coyotefiles.zip
;           http://fermi.jhuapl.edu/s1r/idl/s1rlib/local_idl.html
;
; EXAMPLE:
;
;          Window, XSize=800, YSize=400
;          Terminator_Map
;
; NOTES:
;
;       The spherical geometry problems inherent with working the map projections are not
;       completely solved in this program. The main problem seems to be that if the
;       terminator polygon (represented in the program by the variables lon_terminus
;       and lat_terminus) overlaps the pole regions, then stardard methods for determining
;       if a point is inside or outside a polygon are not reliable. (*All* longitudes are
;       inside the polygon if the polygon overlaps a pole!) In the program, this manifests
;       itself as the sun side of the terminator appearing dark and the night side appearing
;       light. Hence, the FLIPDAY keyword.
;
;       To some extent, this problem can be eliminated by setting the CENTER_LAT parameter
;       to zero, at least with some map projections. I've thought about forcing it to be
;       zero and forcing you to use well-behaved map projections, but have decided to leave
;       the program the way it is so you can discover the limitations yourself.
;
; MODIFICATION HISTORY:
;
;       Written by David W. Fanning, April 12, 2006.
;-
;
;###########################################################################
;
; LICENSE
;
; This software is OSI Certified Open Source Software.
; OSI Certified is a certification mark of the Open Source Initiative.
;
; Copyright © 2006 Fanning Software Consulting.
;
; This software is provided "as-is", without any express or
; implied warranty. In no event will the authors be held liable
; for any damages arising from the use of this software.
;
; Permission is granted to anyone to use this software for any
; purpose, including commercial applications, and to alter it and
; redistribute it freely, subject to the following restrictions:
;
; 1. The origin of this software must not be misrepresented; you must
;    not claim you wrote the original software. If you use this software
;    in a product, an acknowledgment in the product documentation
;    would be appreciated, but is not required.
;
; 2. Altered source versions must be plainly marked as such, and must
;    not be misrepresented as being the original software.
;
; 3. This notice may not be removed or altered from any source distribution.
;
; For more information on Open Source Software, visit the Open Source
; web site: http://www.opensource.org.
;
;###########################################################################


PRO Terminator_Map, center_lon, center_lat, $
   FlipDay=flipday, $
   Image=image, $
   Map_Projection=map_index, $
   Keep_Aspect=keep_aspect, $
   Time=time

   ; Error handling.
   Catch, theError
   IF theError NE 0 THEN BEGIN
      Catch, /Cancel
      void = Error_Message()
      RETURN
   ENDIF

   mapProjections = ['Cylindrical', 'Goodeshomolosine', 'Hammer', $
      'Lambert', 'Mercator', 'Miller_Cylindrical', 'Mollweide', 'Robinson']

   ; Check arguments.
   IF N_Elements(center_lon) EQ 0 THEN center_lon = 0
   IF N_Elements(center_lat) EQ 0 THEN center_lat = 0
   IF N_Elements(time) EQ 0 THEN time = Systime()
   IF N_Elements(map_index) EQ 0 THEN map_index=0 ELSE BEGIN
      IF Size(map_index, /TNAME) NE 'STRING' THEN BEGIN
         map_index = 0 > map_index < 7
      ENDIF ELSE BEGIN
         index = Where(StrMid(StrUpCase(mapProjections), 0, $
     StrLen(map_index)) EQ StrUpCase(map_index), count)

         IF count EQ 0 THEN BEGIN
            Message, 'Cannot find map projection: ' + map_index +    '. Using CYLINDRICAL.', /Informational
            map_index = 0
         ENDIF ELSE map_index = index
      ENDELSE
   ENDELSE
   IF Keyword_Set(keep_aspect) THEN BEGIN
      CASE map_index OF
         3:    position = Aspect(1.0, Margin=0.05)     ; Lambert
         4:    position = Aspect(2.0/3.0, Margin=0.05) ; Mercator
         6:    position = Aspect(2.0/3.0, Margin=0.05) ; Miller_Cyindrical
         ELSE: position = Aspect(0.5, Margin=0.05)
      ENDCASE

   ENDIF ELSE position = [0.0, 0.0, 1.0, 1.0]

   ; Set up the map projection space.
   CASE map_index OF
      0: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /CYLINDRICAL
      1: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /GOODESHOMOLOSINE
      2: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /HAMMER
      3: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /LAMBERT
      4: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /MERCATOR
      5: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /ROBINSON
      6: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /MILLER_CYLINDRICAL
      7: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /MOLLWEIDE
   ENDCASE

   ; Obtain the AVHRR image, along with its color table vectors, to display.
   filename = Filepath(SubDirectory=['examples', 'data'], 'avhrr.png')
   image = Read_PNG(filename, r, g, b)

   ; Create a sun symbol to plot.
   phi = Findgen(32) * (!PI * 2 / 32.)
   phi = [ phi, phi(0)]
   UserSym, Cos(phi)*2, Sin(phi)*2, /Fill, Color=Fsc_Color('yellow')

   ; Find the latitude and longitude of the sub-solar point at the
   ; time the user specified. First, sort out the time values.
   time_copy = time
   offset = GMT_OFFSEC()                      ; Find correction from local time to Universal Time (UT).
   DT_TM_INC, time_copy, offset               ; Convert time to UT.
   DT_TM_BRK, time_copy, date_part, time_part ; Break into a date and a time part.
   DATE2YMD, date_part, y, m, d               ; Break data into year, month, and day.
   jd = YMD2JD(y, m, d)                       ; Convert date to julian day.
   ut = SECSTR(time_part)/3600.               ; Convert time string (in seconds) to sun ephemeris time in hours.

    ; Calculate solar RA and DEC and distance (in AU).
   SUN, y, m, d, ut, $
      APP_RA=ra, $           ; Apparent right accession.
      APP_DEC=dec, $         ; Apparent declination.
      DIST=sun_distance      ; Distance to sun in AU.

   ; Calculate local mean sidereal time. LMST returns value as fraction of a day,
   lm_sidereal_time = LMST(jd, ut/24.0, 0) * 24

   ; Calculate sub-solar point.
   sun_lat = dec
   sun_lon = 15.0 * (ra - lm_sidereal_time)

   ; Calculate the terminus.
   earthRadius = 6.371009D6 ; The mean Earth radius.
   scanAngle = ASIN(earthRadius / (sun_distance * 1.4956D11))
   arc_distance = !DPI - 1.57080D - scanAngle  ; 90 degrees - scanangle (but here in radians)

   lon_terminus = FltArr(36)
   lat_Terminus = FltArr(36)
   FOR j=0,35 DO BEGIN
      azimuth = j * 10.0
      results = LL_Arc_Distance([sun_lon, sun_lat], arc_distance, azimuth, /Degrees)
      lon_terminus[j] = results[0]
      lat_terminus[j] = results[1]
   ENDFOR

   ; Set up the Z graphics buffer
   thisDevice = !D.Name
   xsize = !D.X_Size
   ysize = !D.Y_Size
   Set_Plot, 'Z'
   Device, Set_Resolution=[xsize,ysize], Z_Buffer=0
   Erase, Color=0

   ; Prevent underflow warnings.
   except = !Except
   !Except = 0
   PolyFill, lon_terminus, lat_terminus, Color=128
   void = Check_Math()
   !Except = except
   snapshot = TVRD(position[0]*xsize, position[1]*ysize, $
      (position[2] - position[0])*xsize, (position[3] - position[1])*ysize)

   ; One would think that the subsolar position would "define" the day side
   ; of the day/night terminator. But it is not as simple as that. The code
   ; commented out below, for example, does NOT work because normal methods
   ; for determining if a point is inside or outside a polygon are unreliable
   ; for spherical polygons that overlap the poles.
   ;   insideit = Inside(sun_lon, sun_lat, lon_terminus, lat_terminus)
   ;   IF insideit THEN $
   ;      dark = Where(snapshot NE 128, count) ELSE $
   ;      dark = Where(snapshot EQ 128, count)
   IF Keyword_Set(flipday) THEN $
      dark = Where(snapshot NE 128, count) ELSE $
      dark = Where(snapshot EQ 128, count)

   ; Display the warped image.
   warp = Map_Image(image, x, y, xs, ys, Compress=1, Missing=254)
   TVImage, warp, x, y, XSIZE=xs, YSIZE=ys, /TV
   img = TVRD(position[0]*xsize, position[1]*ysize, $
      (position[2] - position[0])*xsize, (position[3] - position[1])*ysize)
   red = r[img]
   grn = g[img]
   blu = b[img]
   Color_Convert, red, grn, blu, h, l, s, /RGB_HLS
   rr = red
   gg = grn
   bb = blu
   Color_Convert, h, l*0.5, s, red, grn, blu, /HLS_RGB
   rr[dark] = red[dark]
   gg[dark] = grn[dark]
   bb[dark] = blu[dark]

   Set_Plot, thisDevice

   ; Have to re-establish map projection space after changing to another device,
   ; according to MAP_SET documentation.
   CASE map_index OF
      0: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /CYLINDRICAL
      1: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /GOODESHOMOLOSINE
      2: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /HAMMER
      3: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /LAMBERT
      4: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /MERCATOR
      5: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /ROBINSON
      6: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /MILLER_CYLINDRICAL
      7: Map_Set, center_lat, center_lon, 0, Position=position, /NoErase, /MOLLWEIDE
   ENDCASE

   ; Display the image.
   TVImage, [ [[rr]], [[gg]], [[bb]] ], Position=position

   ; Plot the sun on the map.
   PLOTS, sun_lon, sun_lat, psym=8


   ; Add continental outlines and time zones.
   Map_Continents, Color=Fsc_Color('Medium Gray')
   Map_Continents, Color=Fsc_Color('Medium Gray'), /Countries
   Map_Grid, Color=Fsc_Color('Light Gray'), Lons=Scale_Vector(Indgen(25), -180, 180), $
      LonNames = Reverse(StrTrim(Indgen(25),2)), LonLabel=-5, /Label

   IF Arg_Present(image) THEN image = TVRead(True=3)
END