;+
; NAME:
;   vu_thomson_antifish
; PURPOSE:
;   Wrapper for href=PlotEarthSkymap=. Displays a view of the sky from Earth
;   using a 3D heliospheric matrix from one of the tomography programs.
; CATEGORY:
;   WWW: graphics
; CALLING SEQUENCE:
	PRO vu_EarthSkymap, hdr, ff, ut0=ut0, ihdr=ihdr, type=type,     $
		nodata=nodata, maxelo=maxelo, geolng=geolng,          $
		degrees=degrees, equator=equator, _extra=_extra,      $
		solar=solar, band=band, F_full=FF_full, RA_full=RA_full,   $
		is_sequence=is_sequence, point_sources=point_sources,   $
		hardcopy=hardcopy, count=count, nRA=nRA, nDE=nDE, nRR=nRR, $
		dRR=dRR, silent=silent, module=module, zero_shift=zero_shift
; INPUTS:
;   hdr          array[1]; type: structure
;               information about the tomographic reconstruction.
;               if 'hdr' does not exist, then user is prompted for
;               a data file using the 'pickfile' dialog.
;   ff          array[n,m,l,2]; type: float
;               3D volume data; if not present then the files
;               specified in the 'hdr' argument are read.
; OPTIONAL INPUT PARAMETERS:
;   ut0=ut0      scalar or array[1]; type: time structure, float or integer
;               Specifies the time on which the time series is centered.
;               If not specified then the time from the 'hdr' argument is used.
;               See href=arg_time= for more details.
;   maxelo         scalar; type: float
;               if set and non-zero then a fish-eye map is plotted out to
;               elongation 'maxelo'. By default a Hammer-Aitoff map is drawn
;   destination=destination
;             if set to valid directory then graphics output is saved as GIF files
;   /degrees   if set, indicates that all angles are in degrees. Default: radians.
;   /equator   if set, the 'horizontal' in the sky map is parallel to the
;               equatorial plane. By default the horizontal is along the ecliptic
;   geolng         scalar; type: float
;               geographic longitude of location on Earth (in units implied
;               by /degrees keyword). If 'geolng' is specified then a 'transit map'
;               is calculated for a one-day period centered on 'ut'
;
;   zero_shift       passed to href=EarthSky3DLoc=
;             (NOT used for sky sweeps, i.e. NOT used if 'geolng' is set).
;             By default the skymap is centered on the Sun. Setting 'zero_shift' moves
;             the center to increasing RA (or ecliptic longitude) by 'zero_shift'.
;             E.g. to produces an anti-solar skymap set zero_shift=!pi
;             (or zero_shift=180, /degrees).
;
;   point_sources=point_sources
;             array[n]; type: structure
;               information about IPS observation (location in sky, IPS speed
;               and g-level; see href=GetNagoyaSources=.
;   /names         if set, then the names of the point sources are plotted
;
;   /solar
;
;   (only if keyword 'band' is used)
;
;   band=band     scalar; type float
;               used to replace a limited strip in RA in an already existing skymap.
;               band: width of strip in hours
;               The strip is inserted into arrays specified in keywords
;               F_full and RA_full.
;   RA_full=RA_full   array(nRA); type: float
;    F_full=FF_full   array(nRA,nDE); type: float
;               skysweep arrays (RA and function value) where only a strip of width
;               band needs to be updated. Usually these are obtained from a
;               previous call to this procedure. If they don't exist yet then
;               keyword 'band' is ignored and a full skysweep map is intialized.
;
;   Several keywords are passed to href=vu_type_skymap=:
;   The most important ones set the data type to be plotted:
;     /ips_speed     IPS speed
;     /ips_glevel  IPS g-level
;     /brightness  Thomson scattering brightness
;     /pb         Thomson scattering polarization brightness
;
;   type=type     array[1]; type: structure
;
; OUTPUTS:
;   (none)
; OPTIONAL OUTPUTS:
;   (only if keyword 'band' is used)
;
;   RA_full=RA_full   array(nRA); type: float: default: none
;    F_full=FF_full   array(nRA,nDE); type: float; default: none
;               updated skysweep arrays (RA and function value)
;               with new strip of width band inserted
; INCLUDE:
    @compile_opt.pro      ; On error, return to caller
; EXTERNAL:
;   IPS_velocity, IPS_velocity_w, IPS_hlevel_w, ThomsonBrightness
; CALLS:
;   arg_time, IsType, vu_getdata, vu_gettime, vu_type_sky, set_page
;   EarthTransit3DLoc, EarthSky3DLoc, IPS_params, SuperArray
;   UBVConst, IntegrateLOS, AngleRange, PlotEarthSkymap, vu_get_page
;   reset_colors, InitVar, vu_get, vu_point_source, ThomsonRadialFilter
;   TimeString, TimeUnit
; PROCEDURE:
; > Several parameters define the grid used to calculate the sky map.
;   These are passed to href=EarthSky3DLoc= and href=EarthTransit3DLoc=:
;
;   nRA     scalar; type: int; default: 72
;          # points in skymap in longitudinal direction (either right ascensions
;          of ecliptic longitude)
;   nDE     scalar; type: int; default: 36
;          # points in skymap in latitudinal direction(either declination or
;          ecliptic latitude)
;   nRR     scalar; type: int; default: 20
;          # steps along lines of sight at each point in the sky
;   dRR     scalar; type: float; default: 0.1 AU
;          step size along line of sight (AU)
;
; > The sky map is build from nRA*nDE lines of sight integrations through the F3D array
;   The locations in the sky of the lines of sight are
;     RA/Lng    = (RA/Lng Sun) -180 + (360./nRA)*[0.5, 1.5, .. , nRA-0.5 ]
;     Decl./Lat = - 90 + (180./nDE)*[0.5, 1.5, .. , nDE-0.5 ]
;     Dist along los = dRR*[0.5,1.5,...,nRR-1]
;   Note that each line of sight is centered on a 'box' of 360/nRA by 180/nDE degrees,
;   and that all lines of sight together cover the entire sky.
;   The edges of the boxes for all lines of sight are returned in the RA and DE keywords
; MODIFICATION HISTORY:
;   MAR-2000, Paul Hick (UCSD/CASS)
;   AUG-2000, Paul Hick (UCSD/CASS)
;     Added option to accept a vu_read structure as a 1st argument
;   SEP-2001, Paul Hick (UCSD/CASS)
;     Fixed bug in vu_type_sky calls
;   JAN-2002, Paul Hick (UCSD/CASS)
;     Changed the initialization procedure for a skysweep with the 'band'
;     keyword set. The initialization now ignores the 'band' value and simply
;     calculates a full skysweep map (the earlier initialization put a new
;     strip in an otherwise empty sky).
;   JUL-2002, Paul Hick (UCSD/CASS)
;     Added interpolation between neigbouring times of ut0, instead of
;     truncating to nearest grid time.
;   AUG-2004, Paul Hick (UCSD/CASS)
;     Modifed IntegrateLOS calls. For Thomson scattering brightness maps with
;     point sources overplotted a radial filter was applied to the maps but
;     not to the point sources. Fixed.
;   OCT-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;     Added zero_shift keyword.
;-

InitVar, silent  , 0
InitVar, nodata  , /key
InitVar, equator , /key
InitVar, hardcopy, /key

ihdr_set = IsType(ihdr, /defined)

; Set maxelo to non-zero scalar value to get fisheye map (e.g. !pi/2*1.25)

hammer = IsType(maxelo, /undefined)        ; If maxelo undefined then plot Hammer-Aitoff map
IF NOT hammer THEN hammer = maxelo LE 0     ; If maxelo exists plot fisheye map only if maxelo>0

sweep = IsType(geolng, /defined)

CASE nodata OF
1: BEGIN
    utt = arg_time(ut0)

    merge_band = 0
    label = ''

    type = vu_type_sky(hammer=hammer, type=type,_extra=_extra)

END
0: BEGIN

    utt = vu_getdata(hdr, ff, ut0=ut0, ihdr=ihdr, silent=silent, count=count, _extra=_extra)
    IF count EQ 0 THEN RETURN  ; Return if no data files were read

    utt = arg_time(utt)

    ; Pick out the matrix at time utt or hdr[ihdr].time by interpolation

    CASE ihdr_set OF
    0: hdri = vu_gettime(hdr, ff, ut=utt, ff_ut=ffi)
    1: hdri = vu_gettime(hdr, ff, ut=vu_get(hdr[ihdr], /uttime), ff_ut=ffi, /fixgrid); ihdr set explicitly
    ENDCASE


    xcrange = vu_get(hdri, /array_range)

    nLng   = vu_get(hdri, /nlng)
    nLat   = vu_get(hdri, /nlat)
    nRad   = vu_get(hdri, /nrad)

    R3D       = vu_get(hdri, /start_distance)
    dR3D   = vu_get(hdri, /distance_step)

    betan  = (vu_get(hdri, /power_density ))[0]
    betar  = (vu_get(hdri, /power_distance))[0]

    n0     = vu_get(hdri, /density_at_1AU)

    type = vu_type_sky(hammer=hammer, type=type, _extra=_extra)

    InitVar, nRA, 72
    InitVar, nDE, 36
    InitVar, nRR, 20
    InitVar, dRR, 0.1

    ; For a skysweep the keyword 'band' indicates that only a strip of sky needs to be inserted
    ; into an existing skysweep map.
    ; If RA_full (and hence FF_full) doesn't exist yet, ignore the 'band' setting and initialize
    ; a full skysweep map (the band0 variable passed to EarthTransit3DLoc won't exist).

    merge_band = sweep AND IsType(RA_full, /defined) AND IsType(band, /defined)
    IF merge_band THEN band0 = band

    ; Calculate locations (in heliographic coordinates for all points along all
    ; lines of sight. Note that R[0,*,*,*] contains heliographic longitudes in [0,360]
    ;  RA =RA   array[nRA]; type: float
    ;   RA  = -pi+2*pi/nRA*[0.5,1.5,...,nRA-0.5]
    ;  DEC=DEC  array[nDE]; type: float
    ;   DEC = -pi/2+pi/nDE*[0.5,1.5,...,nDE-0.5]
    ; RA and DEC define a regular grid in sky coordinates for the centers of
    ; the boxex for the lines of sight used to build the skymap
    ; (they are output from EarthSky3DLoc or EarthTransit3DLoc)

    CASE sweep OF
    0: R = EarthSky3DLoc(utt, R3D, dR3D,             $
       nRA=nRA, nDE=nDE, nRR=nRR, dRR=dRR, RA=RA, DEC=DEC, $
       rearth=REarth, elo_earth=elo_earth, elo_sun=elo_sun,$
       equator=equator, degrees=degrees,           $
       zero_point=zero_point, zero_phase=zero_phase,     $
       zero_shift=zero_shift)
    1: R = EarthTransit3DLoc(utt, R3D, dR3D,          $
       nRA=nRA, nDE=nDE, nRR=nRR, dRR=dRR, RA=RA, DEC=DEC, $
       rearth=REarth, elo_earth=elo_earth, elo_sun=elo_sun,$
       equator=equator, degrees=degrees,           $
       zero_point=zero_point, zero_phase=zero_phase,     $
       geolng=geolng, solar=solar, band=band0)
    ENDCASE

    label = type.label

END
ENDCASE



nagoya  = 1
fillbad = 0
minRR   = 10



BreakVal = type.breakval

i = where(finite(BreakVal), nbreak)
IF nbreak GT 0 THEN BEGIN
    BreakVal = type.breakval[i]
    CASE type.format of
    ''  : BreakVal = round(BreakVal)
    ELSE: format = type.format
    ENDCASE
ENDIF


CASE label OF

''  : IF silent LE 1 THEN message, /info, 'no sky data'

'_v': BEGIN               ; IPS speed
    IF silent LE 1 THEN message, /info, 'IPS velocity at '+TimeString(utt,/ymd,upto=TimeUnit(/minute))

    D3D = ffi[*,*,*,1]

    id = where( finite(D3D) )
    IF id[0] NE -1 THEN D3D[id] = (D3D[id]/n0)^(2*betan)

    ips_info = [IPS_params(nagoya=nagoya, cambridg=cambridge),betar]

    F = IntegrateLOS(R, ffi[*,*,*,0], D3D, R3D=R3D, dR3D=dR3D,     $
       xcrange=xcrange, minRR=minRR, fillbad=fillbad, degrees=degrees, /longitudes, $
       f3dfnc='IPS_velocity', f3darg=ips_info, w3dfnc='IPS_velocity_w', w3darg=ips_info, $
       REarth=REarth, elo_earth=elo_earth, elo_sun=elo_sun, silent=silent)
END
'_g': BEGIN               ; IPS g-level
    IF silent LE 1 THEN message, /info, 'IPS g-level at '+TimeString(utt,/ymd,upto=TimeUnit(/minute))

    D3D = ffi[*,*,*,1]

    id = where( finite(D3D) )
    IF id[0] NE -1 THEN D3D[id] = (D3D[id]/n0)^(2*betan)

    ips_info = [IPS_params(nagoya=nagoya, cambridg=cambridge),betar]

    F = IntegrateLOS(R, D3D, R3D=R3D, dR3D=dR3D,      $
       xcrange=xcrange, minRR=minRR, fillbad=fillbad, degrees=degrees, /longitudes, $
       w3dfnc='IPS_hlevel_w', w3darg=ips_info,          $
       REarth=REarth, elo_earth=elo_earth, elo_sun=elo_sun, silent=silent)

    F = sqrt(F)
END
'_b': BEGIN                    ; Thomson scattering brightness
    IF silent LE 1 then message, /info, 'Thomson scattering brightness at '+TimeString(utt,/ymd,upto=TimeUnit(/minute))

    f3darg = [dRR,UBVConst(3),0]

    F = IntegrateLOS(R, ffi[*,*,*,1], R3D=R3D, dR3D=dR3D,  $
       xcrange=xcrange, minRR=minRR, fillbad=fillbad, degrees=degrees, /longitudes, $
       f3dfnc='ThomsonBrightness', f3darg=f3darg,         $
       REarth=REarth, elo_earth=elo_earth, elo_sun=elo_sun, silent=silent)

    F = F/ThomsonRadialFilter(elo_earth[*,*,0], degrees=degrees)

    IF IsType(point_sources,/defined) THEN BEGIN

       ; Apply same radial filter for point sources and skymap.
       ; Note that we modify the brightness in the 'point_sources' structures.
       ; The input values are restored before returning.

       point_brightness = vu_point_source(point_sources, /brightness)
       tmp = ThomsonRadialFilter(vu_point_source(point_sources, /elo, degrees=degrees), degrees=degrees)
       point_sources = vu_point_source(point_sources, point_brightness/tmp, /brightness)
    ENDIF

    ;log = 1
END
'_pb': BEGIN                   ; Polarization brightness
    IF silent LE 1 THEN message, /info, 'Thomson scattering pB at '+TimeString(utt,/ymd,upto=TimeUnit(/minute))

    Fb = replicate(n0, nLng,nLat,nRad) ; Normalized density in 1/r^2 density

    f3darg = [dRR,UBVConst(3),1]

    Itb = IntegrateLOS(R, Fb, R3D=R3D, dR3D=dR3D,       $
       xcrange=xcrange, minRR=minRR, fillbad=fillbad, degrees=degrees, /longitudes, $
       f3dfnc='ThomsonBrightness', f3darg=f3darg,         $
       REarth=REarth, elo_earth=elo_earth, elo_sun=elo_sun, silent=silent)

    It  = IntegrateLOS(R, ffi[*,*,*,1], R3D=R3D, dR3D=dR3D,    $
       xcrange=xcrange, minRR=minRR, fillbad=fillbad, degrees=degrees, /longitudes, $
       f3dfnc='ThomsonBrightness', f3darg=f3darg,         $
       REarth=REarth, elo_earth=elo_earth, elo_sun=elo_sun, silent=silent)

    f3darg = [dRR,UBVConst(3),2]

    Irb = IntegrateLOS(R, Fb, R3D=R3D, dR3D=dR3D,       $
       xcrange=xcrange, minRR=minRR, fillbad=fillbad, degrees=degrees, /longitudes, $
       f3dfnc='ThomsonBrightness', f3darg=f3darg,         $
       REarth=REarth, elo_earth=elo_earth, elo_sun=elo_sun, silent=silent)

    Ir  = IntegrateLOS(R, ffi[*,*,*,1], R3D=R3D, dR3D=dR3D,    $
       xcrange=xcrange, minRR=minRR, fillbad=fillbad, degrees=degrees, /longitudes, $
       f3dfnc='ThomsonBrightness', f3darg=f3darg,         $
       REarth=REarth, elo_earth=elo_earth, elo_sun=elo_sun, silent=silent)

    F = (It-Ir)/(Itb-Irb)
    ;F = ((It-Ir)/(It+Ir))/((Itb-Irb)/(Itb+Irb))

    ;log = 1
END
'_cd': BEGIN                    ; Columnar density
    IF silent LE 1 then message, /info, 'Column density at '+TimeString(utt,/ymd,upto=TimeUnit(/minute))

    f3darg = [dRR,UBVConst(3),0]

    F = IntegrateLOS(R, ffi[*,*,*,1], R3D=R3D, dR3D=dR3D,  $
       xcrange=xcrange, minRR=minRR, fillbad=fillbad, degrees=degrees, /longitudes, $
       f3dfnc='ThomsonColumnDensity', f3darg=f3darg,        $
       REarth=REarth, elo_earth=elo_earth, elo_sun=elo_sun, silent=silent)

    F = F/ThomsonRadialFilter(elo_earth[*,*,0], degrees=degrees)

    IF IsType(point_sources,/defined) THEN BEGIN

       ; Apply same radial filter for point sources and skymap.
       ; Note that we modify the brightness in the 'point_sources' structures.
       ; The input values are restored before returning.

       point_brightness = vu_point_source(point_sources, /brightness)
       tmp = ThomsonRadialFilter(vu_point_source(point_sources, /elo, degrees=degrees), degrees=degrees)
       point_sources = vu_point_source(point_sources, point_brightness/tmp, /brightness)
    ENDIF

    ;log = 1
END

ENDCASE


IF merge_band THEN BEGIN         ; Insert strip into skysweep map

    ; Keyword 'band' is set. Only a strip of band hours in RA is updated.
    ; RA (RA in band) is set in EarthTransit3DLoc using utt and band.

    nRA = n_elements(RA)

    ; Insert the strip of band hours. For each of the nRA meridian scans find the location
    ; in RA nearest to it. Insert at that position.

    np = intarr(nRA, /nozero)
    FOR i=0,nRA-1 DO BEGIN
       tmp = min(abs(RA[i]-RA_full), n)
       np[i] = n
    ENDFOR
    FF_full[np,*] = F
    RA_full[np  ] = RA

    F  = FF_full
    RA = RA_full

ENDIF


old_device = !D.name

set_page, zbuffer=1-hardcopy, printer=hardcopy, /silent, xysize=type.xysize, $
    fullsize=hardcopy, _extra=_extra, old_device=old_device, ctable=type.ctable

reset_colors, /reverse
IF NOT hardcopy THEN stretch, !d.n_colors-1, 0

PlotEarthSkymap, utt, RA, DEC, F, log=log, equator=equator, breakval=BreakVal, charsize=type.charsize,  $
    maxelo=maxelo, geolng=geolng, degrees=degrees, title=type.title,   $
    zero_point=zero_point, zero_phase=zero_phase, format=format,     $
    silent=silent, _extra=_extra,                       $
    point_sources=point_sources

CASE nodata OF
0: hdri = hdr[ihdr]
1: hdri = filepath(root=getenv('TEMP'),'source_map.tmp')
ENDCASE

module = type.display+'_'+(['ecl','equ'])[equator]+'_'+(['fisheye','hammer'])[hammer]+label
vu_get_page, hdri, title='', module, ut=utt, device=old_device, silent=silent, _extra=_extra

IF IsType(point_brightness,/defined) THEN   $
    point_sources = vu_point_source(point_sources, point_brightness, /brightness)

RETURN  &  END


;******************************************************************************************************

FUNCTION GetColors, Value, BreakVal, BreakName, open=open,   $
    legend=legend, format=format, ctable=ctable, cpart=cpart,  $
    usedcolors=usedcolors, badforeground=badforeground,         $
    badbackground=badbackground, charsize=charsize,       $
    flip=flip, noedge=noedge, ncolors=ncolors, step=step, _extra=_extra

;+
; NAME:
;   GetColors
; PURPOSE:
;   Convert array of function values to equivalent color indexes.
;   Optionally plot a legend.
; CALLING SEQUENCE:
;   Colors = GetColors( Value, BreakVal, BreakName,   $
;     /legend, format=format, charsize=charsize, ctable=ctable, cpart=cpart )
; INPUTS:
;   Value          array; type: any
;                function values
;   BreakVal       array; type: any
;                function values between adjacent colors
; OPTIONAL INPUT PARAMETERS:
;   /legend          if set, a legend is plotted at the left
;                edge of the screen
;   format=format     scalar; type: string
;                format used to label the legend (only used if
;                BreakVal is a float array)
;   ctable=ctable     scalar; type: integer
;                used to load a color table with LOADCT
;                (only if !d.n_colors=16)
;   cpart=cpart    array[2]; type: float
;                fractions of one; limits the range of color indices
;                used to cpart*(!d.n_colors-1)
;   /noedge        avoids using color indices 0 and !d.n_colors-1
;                this overrides the cpart setting.
;                (these sometimes are set to black and white and don't
;                fit in with the rest of the color table).
;
;   _extra=_extra
;   /badforeground       if set then bad values (detected with the
;                finite function) are set to the foreground
;                color (!p.color)
;   /badbackground       if set then bad values are set to the foreground
;                color (!p.background)
;   /flip          uses color table in reverse order, i.e. with increasing
;                fnc value corresponding to decreasing color index
; OUTPUTS:
;   GetColors        array; type: byte
;                    array with color indices (same dimensions as Value)
; INCLUDE:
    @compile_opt.pro               ; On error, return to caller
; CALLS:
;   IsType, gridgen
; SIDE EFFECTS:
;   The BreakVal array is sorted
; RESTRICTIONS:
;   #break values = n_elements(BreakVal) must be less than
;   #colors = !d.n_colors
; PROCEDURE:
;   The BreakVal array with N = n_elements(BreakVal) elements divides the data
;   range into N+1 interval (the 1st interval covers data values below
;   BreakVal[0], the last interval covers data value above BreakVal[N-1]. A
;   data value is matched to a color index depending on the interval it belongs to.
;
;   The # colors needed is 1+N (N=n_elements(BreakVal).
;   The 1+N color indices are calculated as
;     Col = round( (!d.n_colors-1.)/nBreak*indgen(nCols) )
;   i.e. as nearly evenly spaced over the full range 0,!d.n_colors-1
;   [Col[0]=0, Col[nBreak]=!d.n_colors-1]
;
;   The legend is a vertical color bar (color 0 to nBreak from bottom
;   to top) at the left side of the screen. The BreakVal values are
;   plotted at the borders between adjacent colors. A maximum of 8 evenly
;   spaced numbers is plotted
; MODIFICATION HISTORY:
;   SEP-1992, Paul Hick (UCSD/CASS)
;   APR-1993, Paul Hick (UCSD/CASS)
;     removed restrictions to !d.n_colors=8 and 16
;   AUG-1999, Paul Hick (UCSD/CASS)
;     added a check for bad values using the 'finite' function. Corresponding
;     entries in the Color output array are now set to -1 (these will be ignored
;     by href=ColorSkybox=).
;   SEP-1999, Paul Hick (UCSD/CASS)
;     added badforeground and badbackground options
;     added /open keyword
;   SEP-2003, Paul Hick (UCSD/CASS)
;     added cpart keyword
;   OCT-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;     added /flip and /noedge
;-

nBreak = n_elements(BreakVal)   ; # break values between colors

InitVar, badforeground  , /key
InitVar, badbackground  , /key
InitVar, open      , /key
InitVar, legend        , /key
InitVar, flip      , /key
InitVar, noedge        , /key
InitVar, ncolors       , !d.n_colors

CASE noedge OF
0: InitVar, cpart     , [0,1]
1: cpart = [1.0,ncolors-2.0]/(ncolors-1.0)
ENDCASE

IF nBreak GE ncolors THEN message,  $
    '#breaks ('+strcompress(nBreak,/remove_all)+') must be less then #colors ('$
    +strcompress(ncolors,/remove_all)+')'

InitVar, BreakName, BreakVal

nCols = nBreak+1                 ; # colors to be used
norm = cpart*(ncolors-1.0)

CASE open OF
0B: Col =  gridgen(nCols  ,norm=norm)      ; Color index array
1B: Col = (gridgen(nCols+2,norm=norm))[1:nCols] ; Color index array
ENDCASE
Col = round(Col)

IF flip THEN Col = reverse(Col)

BreakVal = BreakVal[sort(BreakVal)]

bad = where(1B-finite(Value))      ; Locate bad values

IF bad[0] NE -1 THEN BEGIN
    badv = Value[bad]            ; Save all bad values

    ; Set bad values to minimum value. If everything is bad use the lowest break value

    ValueMin = min(Value, /nan)         ; Pick up minimum value
    IF NOT finite(ValueMin) THEN ValueMin = BreakVal[0]

    Value[bad] = ValueMin          ; Replace bad values by minimum value
ENDIF

; Set up the colors array by assigning colors to all elements in between
; successive break values.

Colors = Col[0]+long(0*Value)      ; Set everything to first color

FOR I=1,nBreak-1 DO BEGIN
    A = where ( BreakVal[I-1] LT Value AND Value LE BreakVal[I] , N )
    IF N NE 0 THEN Colors[A] = Col[I]
ENDFOR

A = where ( Value GT BreakVal[nBreak-1] , N )
IF N NE 0 then Colors[A] = Col[nBreak]

IF bad[0] NE -1 THEN BEGIN
    CASE 1 OF
    badbackground: Colors[bad] = !p.background
    badforeground: Colors[bad] = !p.color
    ELSE      : Colors[bad] = -1
    ENDCASE

    Value [bad] = badv
ENDIF

; Load color table if specified

IF IsType(ctable, /defined) THEN loadct, ctable[0]

usedcolors = Colors[ uniq(Colors) ]

IF legend THEN BEGIN

    ; Draw vertical legend at left side of screen

    DyOff = .02
    Dy = (1.-2*DyOff)/nCols  &  Dx = .015  ; Size of box for each color

    ; Loop from bottom to top of screen

    Xleft = 0.01  &  Xright = Xleft+Dx  &  Ytop = DyOff

    FOR I=0,nBreak DO BEGIN              ; Loop over nCols = nBreak+1 colors
       Ybottom = Ytop  &  Ytop = Ytop+Dy
       polyfill, color=Col[I], /normal,   $
         [Xleft  ,Xleft,Xright,Xright ],   $
         [Ybottom,Ytop ,Ytop  ,Ybottom]
    ENDFOR

    ; Convert BreakName array to string array. If BreakName is of type float then
    ; the specified FORMAT is used in the conversion

    CASE IsType(BreakName, /generic_int) OF
    1: BreakLabel = strcompress(BreakName, /remove_all)
    0: BEGIN
       IF NOT keyword_set(format) THEN   $
         message, 'Keyword FORMAT must be assigned a valid float format value'
        BreakLabel = string(BreakName, format=format)
       BreakLabel = strcompress(BreakLabel, /remove_all)
    END
    ENDCASE

    Xleft = Xright+.75*Dx
    InitVar, charsize, 1.0

    InitVar, Step, 1+nCols/8
    step = step > 1

    FOR I=0,nBreak-1,Step DO BEGIN
    ;   Ybottom = DyOff+(I+1)*Dy-0.0125*charsize
        Ybottom = DyOff+(I+1)*Dy-0.0125*charsize/4.
       xyouts, Xleft, Ybottom, /normal, BreakLabel[I], charsize=charsize, _extra=_extra
    ENDFOR

ENDIF

RETURN, Colors  &  END

;**************************************************************************************************************

PRO PlotEarthSkymap, UT, RA, DEC, F, equator=equator,      $
    maxelo=maxelo, geolng=geolng, degrees=degrees,          $
    mirror=mirror, breakval=BreakVal, log=log,charsize=charsize,$
    title=title, upto=upto, scale=scale, minelo=minelo,         $
    zero_point=zero_point, zero_phase=zero_phase, dabg=dabg,   $
    format=format, fill2edge=fill2edge, _extra=_extra,        $
    galactic=galactic, silent=silent,                 $

    point_sources=point_sources, point_names=point_names,   $
    point_size=point_size, point_onesize=point_onesize,      $

    position=position, legend=legend, compass=compass,        $
    max_dec=max_dec, ra_step=ra_step, dec_step=dec_step
;+
; NAME:
;   PlotEarthSkymap
; PURPOSE:
;   Plot sky map, with locations of point sources overplotted (optional)
; CATEGORY:
;   www
; CALLING SEQUENCE:
;   PlotEarthSkyMap, UT, RA, DE, F
; INPUTS:
;   UT         scalar; type: standard time structure
;             time; determines position of Sun in the sky
;   RA         array[n]; type: float
;             array of ecliptic longitudes or right ascensions
;             (if /equator set) in the range [-180,+180]
;             (values are mapped into this range if they are not)
;   DEC      array[m]; type: float
;             array of ecliptic latitudes or declinations
;             (if /equator is set) in the range [-90,90]
;   F        array[n,m]; type: float
;             array of function values in sky map; each function value
;             refers to bin on the sky. The edges of the skybin are
;             specified in RA and DEC
; OPTIONAL INPUT PARAMETERS:
;   breakval  array[*]; type: integer or float
;             levels between colors (passed to ColorSkybox)
;             if not set, a set of break values is calculated
;             equally spaced between minimum and maximum
;   /log   if set, changes to logarithmic scale
;   /degrees  if set, indicates that all angles are in degrees
;             Default: radians.
;   maxelo       scalar; type: float
;             used to decide on map projection
;             > 0: fish-eye map out to elongation 'maxelo'
;             = 0: Mercator projection (similar to synoptic map)
;             < 0: Hammer-Aitoff projection
;             if absent then a Hammer-Aitoff map is drawn
;   minelo       scalar; type: float
;             Erases plot inside of minelo
;   /equator  if set, the 'horizontal' in the sky map is parallel to the
;             equatorial plane. By default the horizontal is along the ecliptic
;   upto=upto scalar; type: integer; default: TimeUnit(/hour)
;             controls the length of the UT string plotted
;   scale=scale   scalar; type: float; default: 1.0
;             controls the overall size of the skymap relative to the plot window
;   /galactic if set, add a line for the galactic plane
;
;   point_sources=point_sources
;          array; type: sky_point_source structure
;             contains information about point sources to be overplotted
;             on the skymap. See href=vu_point_source=.
;   /point_names
;          if set, then the names of the point sources are plotted
;
;   point_size=point_size
;          scalar or array[2]; controls the size of the circles used to plot the
;          point source position. The size is specified in radians or degrees
;          (depending on setting of /degrees).
;
;          point_size[0]: minimim size of circle.
;          point_size[1]: increment in circle size per unit of the function
;             value in skymap F.
;
;          point_size[1] is used to increase the circle size proportional to
;          the difference between the point source fnc value and the value in
;          the skymap F.
;
;          If point_size[1] is not set bad (and /point_onesize is NOT set)
;          then only sources with valid function values will be plotted.
;          (I.e. to plot only good sources with the same cirle size set
;          point_size[1] to zero.)
;
;          If point_size[1] is bad or omitted then /point_onesize is assumed set.
;
;   /point_onesize
;          if set then all sources (including those with bad fnc values)
;          are plotted with the same size circle point_size[0].
;          This is useful to override the information stored
;          in 'point_sources' or in keyword element point_size[1].

;          If only IPS sources are plotted (i.e. if no skymap F is specified) then
;          setting /point_onesize will plot all source position with the same size.
;
;   The following three keywords are clumsy to use. The idea is to control the orientation
;   of the map (i.e. direction of origin, and tilt of horizontal on sky).
;   href=vu_EarthSkymap= uses these keywords to set up sky 'snapshots' and sky 'sweeps'.
;
;   dabg   array[3]; type: float; default: [0,0,0]
;             (used by the projection functions FishEye, HammerAitoff and
;             MercatorProj, either directly in this procedure or indirectly in
;             ColorSkybox)
;
;   zero_point    scalar; type: float; default: 0.0
;             (used internally and passed to FishEye, HammerAitoff and MercatorProj)
;             Defines the longitude or RA of the center of the FishEye,HammerAitoff
;             or Mercator map.
;
;   zero_phase    scalar, or array with same structure as 'RA'; type: float; default: zero_point
;             (passed to ColorSkybox)
;             ColorSkybox uses this keyword to rearrange its first three arguments
;             (RA,DEC and fnc-values) to put zero_phase in the center of the map.
;             This modifies input arrays RA and DEC, but not F.
;
;   For sky 'snapshots' zero_phase is the same as zero_point.
;   For sky 'sweeps' zero_phase is a monotonic array with a value for each sweep, and
;   the value of 'zero_point' (almost) exactly in the center of 'zero_phase'. Note that
;   inconsistent settings of zero_point and zero_phase make an incorrect map.
; OUTPUTS:
; INCLUDE:
    @compile_opt.pro      ; On error, return to caller
; CALLS:
;   ToRadians, ToDegrees, ColorSkybox, GetColors, NewcombSun, CvSky, FishEye, HammerAitoff
;   PlotCurve, AngleRange, gridgen, GeographicInfo, TimeString, TimeUnit, ANiceNum
;   IsType, EulerRotate, MercatorProj, InitVar, vu_point_source
; SEE AlSO:
;   ips_point_source__define
; PROCEDURE:
; > Time UT is needed to calculate the location of equator and ecliptic
; > The ecliptic is drawn on an equatorial skymap and v.v
; > Longitude/RA is plotted increasing right-to-left across the map
;   (as it is for a viewer at Earth looking up at the sky). For sun-centered maps
;   this means that east is left and west is right. A mirror image of the map
;   can be made by setting keyword /mirror.
; > The data range for the horizontal axis is -180,+180 degrees relative to
;   'zero_point' as needed for the Hammer-Aitoff and fisheye projections. I.e.
;   longitude/RA 'zero_point' will appear at the center of the map.
; MODIFICATION HISTORY:
;   SEP-1999, Paul Hick (UCSD/CASS)
;   FEB-2002, Paul Hick (UCSD/CASS)
;     Added option to plot a sky map in Mercator projection by setting maxelo=0.
;   APR-2002, Paul Hick (UCSD/CASS)
;     Added /fill2edge keyword (passed to ColorSkybox)
;   SEP-2003, Paul Hick (UCSD/CASS)
;     Minor tweaking of labels to get Mercator projection to look better
;     Added /galactic keyword.
;   AUG-2004, Paul Hick (UCSD/CASS)
;     Started adding some code to deal with Thomson scattering brightness
;     (mainly for overplotting point sources).
;     The time plotted at the top is now rounded instead of truncating
;     to timeunit 'upto'.
;   SEP-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;     For fisheye plots point sources with outside maxelo are not plotted anymore
;     Reworked the determination of the circle sizes for the point sources, and
;     added some keywords to control them.
;-

rpm     = ToRadians(degrees=Degrees)
dpm     = ToDegrees(degrees=Degrees)
pi   = !dpi/rpm
pi2     = pi/2
pi4     = pi/4

skymap  = IsType(F, /defined)         ; Check whether skymap is specified
pnt_n   = n_elements(point_sources)     ; Check whether any point sources are specified
InitVar, point_names  , /key
InitVar, point_onesize, /key
InitVar, format       , '(F4.1)'

InitVar, galactic   , /key
InitVar, zero_point , 0d0
InitVar, silent     , 0
InitVar, charsize   , 1

IF pnt_n NE 0 THEN BEGIN

    CASE 1 OF
    IsType(title, /undefined): BEGIN
       pos_fnc = (where( tag_names( point_sources[0] ) EQ 'VV'   ))[0]
       pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZVV' ))[0]
    END
    strmid(title,0,9) EQ 'IPS speed': BEGIN
       pos_fnc = (where( tag_names( point_sources[0] ) EQ 'VV' ))[0]
       pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZVV' ))[0]
    END
    strmid(title,0,11) EQ 'IPS g-level': BEGIN
       pos_fnc = (where( tag_names( point_sources[0] ) EQ 'GG' ))[0]
       pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZGG' ))[0]
    END
    strmid(title,0,11) EQ 'Thomson sca': BEGIN
       pos_fnc = (where( tag_names( point_sources[0] ) EQ 'BB' ))[0]
       pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZBB' ))[0]
    END
    strmid(title,0,10) EQ 'Thomson pB': BEGIN
       pos_fnc = (where( tag_names( point_sources[0] ) EQ 'PB' ))[0]
       pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZPB' ))[0]
    END
    ELSE: BEGIN              ; Safety belt; probably redundant
       pos_fnc = (where( tag_names( point_sources[0] ) EQ 'VV' ))[0]
       pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZVV' ))[0]
    END
    ENDCASE

    IF IsType(point_size,/defined) THEN BEGIN

       ; If point_size[1] is not present or set bad then assume point_onesize=1
       ; i.e. plot all sources (including bad ones) with the same circle size.

       point_sz = point_size
       IF n_elements(point_sz) EQ 1 THEN point_sz = [point_sz, BadValue(point_sz)]
       point_onesize = point_onesize OR point_sz[1] EQ BadValue(point_sz)

       point_sz0 = replicate(point_sz[0], pnt_n)
       point_sz  = replicate(point_sz[1], pnt_n)

    ENDIF

    pnt_n_ = pnt_n
    CASE point_onesize OF
    0: good_points = where(finite(vu_point_source(point_sources,field=pos_fnc)), pnt_n)
    1: good_points = indgen(pnt_n)          ; Plot all sources
    ENDCASE

    IF pnt_n NE 0 THEN BEGIN

       tmp = point_sources[good_points]

       point_nam = strcompress(vu_point_source(tmp,/name), /remove_all)
       point_pos = vu_point_source(tmp, /location, degrees=degrees)

       point_fnc = vu_point_source(tmp, field=pos_fnc)
       point_tim = vu_point_source(tmp, /time)

       IF IsType(point_sz, /undefined) THEN BEGIN
         point_sz0 = vu_point_source(tmp, /size_base, degrees=degrees)
         IF total(point_sz0) LE 0.0 THEN point_sz0 = replicate(2.5/dpm,pnt_n)
         point_sz  = vu_point_source(tmp, field=pos_sz, degrees=degrees, /size_ratio)
       ENDIF

    ENDIF

ENDIF


IF NOT skymap AND pnt_n EQ 0 THEN message, 'no sky map or point source info specified'

sweep = IsType(geolng, /defined)

CASE IsType(maxelo, /defined) OF
0: BEGIN                    ; Hammer-Aitoff if maxelo undefined
    hammer  = 1
    fish    = 0
    mercator = 0
    skyedge  = -1
END
1: BEGIN                    ; 'maxelo' is defined
    hammer   = maxelo LT 0.0
    mercator = maxelo EQ 0.0
    fish     = maxelo GT 0.0
    skyedge  = maxelo
END
ENDCASE

InitVar, scale  , 1
InitVar, log    , /key
InitVar, equator, /key
InitVar, mirror , /key

mirror  = 2*mirror-1


; Calculate position of Sun in angular coordinates at time UT

SunEcl = NewcombSun(UT, /longitude, /latitude, degrees=degrees)
SunEqu = CvSky(UT, fromecliptic=SunEcl, /toequator, degrees=degrees)

CASE sweep OF

0: BEGIN

    ; A 'snapshot' skymap is always centered on the Sun (indicated by zero_point)

    Sun = SunEqu*equator+SunEcl*(1-equator)    ; Snapshot

END

1: BEGIN                        ; Sweep of sky

    ; Sweep of the sky centered on the local meridian at time UT, or on the Sun
    ; utnoon is nearest noon to UT

    GeographicInfo, UT, location=geolng, utnoon=utnoon, degrees=degrees

    ; Why utnoon??

    Sun = NewcombSun(utnoon, /longitude, /latitude, degrees=degrees)
    Sun = CvSky(utnoon, fromecliptic=Sun, /toequator, degrees=degrees)

    Sun[1] = equator*Sun[1]

END

ENDCASE



; Range of Hammer-Aitoff projection is -2*sqrt(2),2*sqrt(2) in horizontal
; and -1*sqrt(2),1*sqrt(2) in vertical direction.

CASE 1 OF
fish    : projsz = FishEye    ( [pi2,0], maxelo=pi2, degrees=degrees )
hammer  : projsz = HammerAitoff( [pi ,0], degrees=degrees )
mercator: projsz = MercatorProj( [pi ,0], degrees=degrees )
ENDCASE

projsz = projsz[0]

; 'xrange' is the horizontal range in data units. The window size exists already at this
; point (at the latest set by the 'erase' command above).

xrange = scale*mirror*projsz*[-1,1]
InitVar, position, [0.55,0.50]#[1,1]+0.4*[1,1]#[-1,1]

erase

plot, xrange, xrange*(float(!d.y_size)/!d.x_size),   $
    xrange=xrange  ,    $       ; Needed to flip the x-axis if mirror is set
    /nodata      ,    $
    /noerase     ,    $
    position=position, $
    xstyle=5     ,    $
    ystyle=5

!p.clip[2:3] = !p.clip[2:3]+1   ; Fixes bug with IDL clip window

IF IsType(BreakVal, /undefined) THEN BEGIN

    CASE 1 OF
    skymap     : fmin = min(F, /nan, max=fmax)
    pnt_n ne 0 : fmin = min(point_fnc, /nan, max=fmax)
    ENDCASE

    IF log THEN BEGIN
       fmin = alog10(fmin)
       fmax = alog10(fmax)
    ENDIF
    fmin = ANiceNum(fmin)
    fmax = ANiceNum(fmax, /upper)

    N = 40
    BreakVal = gridgen(N, /open, norm=[fmin,fmax])

    IF log THEN BreakVal = 10^BreakVal

ENDIF


nRA = n_elements(RA )
nDE = n_elements(DEC)

CASE 1 OF
skymap: BEGIN               ; ColorSkybox may change skyedge if it's too big
    InitVar, zero_phase, zero_point
    ColorSkybox, RA, DEC, GetColors(F, BreakVal, legend=legend, format=format,     $
       /open, charsize=charsize, _extra=_extra), skyedge=skyedge,              $
       zero_phase=zero_phase, dabg=dabg, degrees=degrees, /black, fill2edge=fill2edge
END
pnt_n NE 0: good = GetColors (point_fnc, BreakVal, legend=legend, format='(F6.2)',  $
       /open, charsize=charsize)
ENDCASE

nangle  = 61
angle   = gridgen(nangle, norm=[0,2*pi])
xcircle = cos(angle*rpm)
ycircle = sin(angle*rpm)

IF IsType(minelo, /defined) THEN BEGIN

    ; Erase innermost minelo elongations from plot.
    ; Set up a cone with opening angle minelo around the Sun in a coordinate system
    ; with the z-axis pointing to the Sun, and the x-axis in the ecliptic/equatorial plane.
    ; Then convert to ecliptic/equatorial coordinates.

    pos = transpose([[angle], [replicate(pi2-minelo, nangle)]])
    pos = EulerRotate([-pi2,-(pi2-Sun[1]),-Sun[0]], pos, degrees=degrees)

    CASE 1 OF
    fish    : pos = FishEye       (pos, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)
    hammer  : pos = HammerAitoff(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees)
    mercator: pos = MercatorProj(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees)
    ENDCASE

    polyfill, pos[0,*], pos[1,*], color=!p.background

ENDIF

CASE 1 OF
fish    : Sun = FishEye     (Sun, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)
hammer  : Sun = HammerAitoff(Sun, zero_phase=zero_point, dabg=dabg, degrees=degrees)
mercator: Sun = MercatorProj(Sun, zero_phase=zero_point, dabg=dabg, degrees=degrees)
ENDCASE

rr = 0.015*projsz             ; Mark location Sun
polyfill, Sun[0]+rr*xcircle, Sun[1]+rr*ycircle

; Overplot point sources over skymap

IF pnt_n NE 0 THEN BEGIN
    pos = point_pos              ; RA and declination of sources

    IF NOT equator THEN BEGIN

       CASE sweep OF
       0: pos = CvSky(UT, fromequatorial=pos, /toecliptic, degrees=degrees)
       1: BEGIN            ; This is the same trick as used in EarthTransit3Dloc
         B   = 23.439291/!radeg   ; Angle equator/ecliptic
         AB  = reform(pos[0,*])*rpm
         tmp = sin( atan( sin(AB)/cos(B), cos(AB) ) )
         tmp = atan( tmp*sin(B), tmp*cos(B)/sin(AB) )/rpm
         pos[1,*] = pos[1,*]-tmp  ; Subtract decl of ecliptic at RA pos[0,*]
       END
       ENDCASE

    ENDIF

    CASE skymap OF
    0: BEGIN

       ; If /point_onesize is set then point_fnc[i] may be bad
       ; If NOT set then all point_fnc values will be valid.

       CASE point_onesize OF
       0: dG = point_sz0+point_sz*point_fnc
       1: dG = point_sz0
       ENDCASE

       dG = transpose([[dG],[replicate(0,pnt_n)]])

       CASE 1 OF
       fish  : dG = FishEye      ( dG, maxelo=pi2, degrees=degrees )
       hammer    : dG = HammerAitoff( dG, degrees=degrees )
       mercator: dG = MercatorProj( dG, degrees=degrees )
       ENDCASE

       dG = reform(dG[0,*])

       color = replicate(!p.color, pnt_n)

    END
    1: BEGIN

       lng = reform(pos[0,*])
       lat = reform(pos[1,*])

       lng = AngleRange(lng, degrees=degrees)

       RAs = AngleRange(RA, degrees=degrees)
       tmp = max(RAs, izero)
       RAs = shift(RAs,-izero-1)

       i = round( interpol(indgen(nRA), RAs, lng) )
       j = round( interpol(indgen(nDE), DEC, lat) )

       a = where(0 LE i AND i LE nRA-1 AND 0 LE j AND j LE nDE-1)

       ; Background values in skymap at grid points nearest to the
       ; point source location. Note that some of these may be bad.

       dF = replicate(BadValue(0.0),pnt_n)
       IF a[0] NE -1 THEN dF[a] = F[(i[a]+izero+1) mod nRA,j[a]]
       dF = point_fnc-dF     ; Difference with nearest grid point

       ; Set up size of circle in degrees. Then convert to data units.
       ; Don't merge this to dG = point_sz0+(1-point_onesize)*point_sz*abs(dF)
       ; If point_onesize is set we don't want bad values in dG.

       CASE point_onesize OF
       0: dG = point_sz0+point_sz*abs(dF)
       1: dG = point_sz0
       ENDCASE

       dG = transpose([[dG],[replicate(0,pnt_n)]])

       CASE 1 OF
       fish  : dG = FishEye      ( dG, maxelo=pi2, degrees=degrees )
       hammer    : dG = HammerAitoff( dG, degrees=degrees )
       mercator: dG = MercatorProj( dG, degrees=degrees )
       ENDCASE

       dG = reform(dG[0,*])

       ; Select colors for point sources from the same color table as the skymap

       good = GetColors(point_fnc, BreakVal, /open)

       ; The circle edge is drawn using the foreground color for positive,
       ; and background color for negative differences.

       color = (dF GE 0 OR 1-finite(dF))*!p.color+(dF LT 0)*!p.background
    END
    ENDCASE

    ;print, projsz, dg[0], dg[0]/projsz
    CASE 1 OF
    fish    : pos = FishEye       (pos, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)
    hammer  : pos = HammerAitoff(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees)
    mercator: pos = MercatorProj(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees)
    ENDCASE

    ; For fish-eye maps plot only point sources inside skyedge.
    ; Note that the elongation is calculated for the time of the skymap (NOT for
    ; the time of observation of the point source).

    point_elo = vu_point_source( vu_point_source(point_sources[good_points],UT,/time), /elong, degrees=degrees)

    n = 0
    FOR i=0,pnt_n-1 DO BEGIN
       IF finite(dG[i]) AND (NOT fish OR point_elo[i] LT skyedge) THEN BEGIN
         polyfill, pos[0,i]+dG[i]*xcircle, pos[1,i]+dG[i]*ycircle,color=good [i]
         oplot   , pos[0,i]+dG[i]*xcircle, pos[1,i]+dG[i]*ycircle,color=color[i]
         n = n+1
         IF point_names THEN BEGIN
          sgn = pos[0,i] LT 0
          xyouts, pos[0,i]+(1-2*sgn)*dG[i], pos[1,i], point_nam[i], align=-.1*sgn + 1.1*(1-sgn)
         ENDIF
       ENDIF
    ENDFOR

    IF silent LE 0 THEN message, /info, 'plotting'+strcompress(n)+'/'+strcompress(pnt_n_,/rem)+' sources'

    pos = convert_coord(pos, /data, /to_device)
    point_sources[good_points] = vu_point_source(point_sources[good_points], round(pos[0:1,*]), /xy_pixel )
ENDIF

goodcolor = !d.n_colors-1
badcolor  = !p.color

CASE 1 OF
hammer OR mercator: BEGIN

    IF hammer THEN max_dec = 0.5*pi ELSE InitVar, max_dec, 0.5*pi
    max_dec = max_dec < 0.5*pi

    n   = 37
    lng = replicate(1,n)
    lat = gridgen(n, norm=[-1,1])

    lng = [-lng, lng,-lng[0]]       ; Longitudes relative to zero_phase
    lat = [ lat,-lat, lat[0]]

    n = n_elements(lng)

    ; Plot outer boundary: +/- 180 deg meridian and horizontal lines at poles
    ; for the mercator projection.

    pos = transpose( [[pi*lng],[max_dec*lat]] )
    IF hammer THEN pos = HammerAitoff(pos, degrees=degrees)
    oplot, pos[0,*], pos[1,*]

    ; Meridians at steps of ra_step

    InitVar, ra_step , 0.5d0*pi
    InitVar, dec_step, max_dec ;0.5d0*pi

    this = fix(pi/ra_step)
    this = -this*ra_step

    n   = 37
    lng = replicate(1,n)
    lat = max_dec*gridgen(n, norm=[-1,1])

    WHILE this LE pi DO BEGIN

       pos_this = this
       IF mercator THEN pos_this = AngleRange(this-zero_point,degrees=degrees,/pi)

       IF pos_this NE 0 THEN BEGIN

         pos = transpose( [[pos_this*lng],[lat]] )

         IF skymap THEN BEGIN
          RAs = AngleRange(RA, degrees=degrees)
          tmp = max(RAs, izero)
          RAs = shift(RAs,-izero-1)

          i = round( interpol(indgen(nRA), RAs, AngleRange(reform(pos[0,*])+zero_point, degrees=degrees)) )
          j = round( interpol(indgen(nDE), DEC, reform(pos[1,*])) )

          a = where(0 LE i AND i LE nRA-1 AND 0 LE j AND j LE nDE-1)

          good = bytarr(n)
          IF a[0] NE -1 THEN good[a] = finite(F[(i[a]+izero+1) mod nRA,j[a]])
         ENDIF

         IF hammer THEN pos = HammerAitoff(pos, degrees=degrees)

         IF mercator OR abs(this) NE pi THEN BEGIN
          CASE goodcolor EQ badcolor OF
          0: PlotCurve, /oplot, pos[0,*], pos[1,*], good, color=[goodcolor,badcolor], /changecolor, /silent
          1: PlotCurve, /oplot, pos[0,*], pos[1,*], color=goodcolor, /silent
          ENDCASE
         ENDIF

         IF mercator THEN BEGIN
          xyouts, pos_this, -max_dec-4*0.0177*projsz*charsize, charsize=charsize, $
              strcompress(round(AngleRange(this*dpm,/degrees)),/rem), alignment=0.5
          IF pos_this EQ pi THEN BEGIN
              xyouts, -pi, -max_dec-4*0.0177*projsz*charsize, charsize=charsize, $
                 strcompress(round(AngleRange(this*dpm,/degrees)),/rem), alignment=0.5
          ENDIF
         ENDIF

       ENDIF

       this = this+ra_step

    ENDWHILE

    n   = 37
    lng = pi*gridgen(n, norm=[-1,1])
    lat = replicate(1,n)

    this = -fix(max_dec/dec_step)*dec_step

    WHILE this LE max_dec DO BEGIN

       IF this NE 0 THEN BEGIN

         pos = transpose( [[lng],[this*lat]] )

         IF skymap THEN BEGIN
          RAs = AngleRange(RA, degrees=degrees)
          tmp = max(RAs, izero)
          RAs = shift(RAs,-izero-1)

          i = round( interpol(indgen(nRA), RAs, AngleRange(reform(pos[0,*])+zero_point, degrees=degrees)) )
          j = round( interpol(indgen(nDE), DEC, reform(pos[1,*])) )

          a = where(0 LE i AND i LE nRA-1 AND 0 LE j AND j LE nDE-1)

          good = bytarr(n)
          IF a[0] NE -1 THEN good[a] = finite(F[(i[a]+izero+1) mod nRA,j[a]])
         ENDIF

         IF hammer THEN pos = HammerAitoff(pos, degrees=degrees)

         CASE goodcolor EQ badcolor OF
         0: PlotCurve, /oplot, pos[0,*], pos[1,*], good, color=[goodcolor,badcolor], /changecolor, /silent
         1: PlotCurve, /oplot, pos[0,*], pos[1,*], color=goodcolor, /silent
         ENDCASE

         IF mercator THEN BEGIN
          xyouts, !x.crange+0.0566*projsz*mirror*[-1,1]*charsize, [this,this]-0.0177*projsz*charsize, $
              strcompress(round(this*dpm),/rem), alignment=0.5, charsize=charsize
         ENDIF

       ENDIF

       this = this+dec_step

    ENDWHILE


    FOR plane=0,galactic DO BEGIN

       n = 120                    ; Plot ecliptic/equator
       pos = transpose( [[gridgen(n, /open, norm=pi*[-1,1])],[fltarr(n)]] )

       CASE plane OF
       0: BEGIN
         CASE equator OF
         0: pos = CvSky(UT,fromequatorial=pos,/toecliptic  ,degrees=degrees)
         1: pos = CvSky(UT,fromecliptic  =pos,/toequatorial,degrees=degrees)
         ENDCASE
       END
       1: BEGIN
         CASE equator OF
         0: pos = CvSky(UT,fromgalactic=pos,/toecliptic  ,degrees=degrees)
         1: pos = CvSky(UT,fromgalactic=pos,/toequatorial,degrees=degrees)
         ENDCASE
       END
       ENDCASE

       IF skymap THEN BEGIN
         RAs = AngleRange(RA, degrees=degrees)
         tmp = max(RAs, izero)
         RAs = shift(RAs,-izero-1)

         i = round( interpol(indgen(nRA), RAs, AngleRange(reform(pos[0,*]), degrees=degrees)) )
         j = round( interpol(indgen(nDE), DEC, reform(pos[1,*])) )

         a = where(0 LE i AND i LE nRA-1 AND 0 LE j AND j LE nDE-1)
         good = bytarr(n)
         IF a[0] NE -1 THEN good[a] = finite(F[(i[a]+izero+1) mod nRA,j[a]])
       ENDIF

       pos[0,*] = pos[0,*]-zero_point
       pos[0,*] = AngleRange(pos[0,*], /pi, degrees=degrees)

       ;a    = max(pos[0,*],i)          ; Make sure lng decreases
       ;pos = shift(pos,0,-i-1)
       ;IF skymap THEN good = shift(good,-i-1)

       a   = sort(pos[0,*])          ; This does the same
       pos = pos[*,a]
       IF skymap THEN good = good[a]

       IF hammer THEN pos = HammerAitoff(pos, dabg=dabg, degrees=degrees)

       CASE goodcolor EQ badcolor OF
       0: PlotCurve, /oplot, pos[0,*], pos[1,*], good, color=[goodcolor,badcolor], /changecolor, /silent
       1: PlotCurve, /oplot, pos[0,*], pos[1,*], color=goodcolor, /silent, _extra=_extra
       ENDCASE

    ENDFOR

END
fish: BEGIN
    ;origin = [440,399.5]
    ;xy = gridgen([!d.x_size,!d.y_size], origin=origin)
    ;pp = cv_coord(from_rect=xy, /to_polar)
    ;pp = where( 310 le pp[1,*] and pp[1,*] le 330 )
    ;for i=0L,n_elements(pp)-1 do tv, [!p.background], origin[0]+xy[0,pp[i]], origin[1]+xy[1,pp[i]]

    f = projsz/skyedge
;    da = 5/dpm
    da = 7/dpm

    FOR a=pi4,pi2,pi4 DO oplot, f*a*xcircle,f*a*ycircle, _extra=_extra    ; draw 45 and 90deg circle
    IF skyedge GE pi2 THEN FOR a=-1, 1,2 DO $
    xyouts, a*(pi2+da)*f,   -da*f, '90', alignment=.5, _extra=_extra, charsize=charsize
    IF skyedge GE pi4 THEN FOR a=-1,-1,2 DO $
    xyouts,   -da*f, a*(pi4+da)*f, '45', alignment=.5, _extra=_extra, charsize=charsize
END
ENDCASE


CASE 1 OF
fish    : tmp = FishEye        ( [[-pi2,0],[pi2,0],[0,-pi2],[0,pi2]], maxelo=pi2, degrees=degrees )
hammer  : tmp = HammerAitoff( [[-pi ,0],[pi ,0],[0,-max_dec],[0,max_dec]], degrees=degrees )
mercator: tmp = [[-pi ,0],[pi ,0],[0,-max_dec],[0,max_dec]]
ENDCASE

oplot, [tmp[0,0],tmp[0,1]],[tmp[1,0],tmp[1,1]], color=badcolor, _extra=_extra; Draw horizontal axis
oplot, [tmp[0,2],tmp[0,3]],[tmp[1,2],tmp[1,3]], color=badcolor, _extra=_extra  ; Draw vertical axis

; Label with directions (E,W optional)

xyouts, tmp[0,2:3], tmp[1,2:3]+0.005*projsz*[-4,2]*charsize,   $
       (['Ecl ','Eq '])[equator]+['S','N'], alignment=0.5, charsize=charsize, _extra=_extra

CASE 1 OF
IsType(compass,/string ):   $
;    xyouts, tmp[0,0:1]+0.0566*projsz*[-1,1]*charsize, tmp[1,0:1]-0.0177*projsz*charsize,   $
;      compass, alignment=0.5, charsize=charsize
    xyouts, tmp[0,0:1]+0.015*projsz*[-1,1]*charsize, tmp[1,0:1]-0.005*projsz*charsize,   $
       compass, alignment=0.5, charsize=charsize, _extra=_extra
;IsType(compass,/defined):  $
;   xyouts, tmp[0,0:1]+0.0566*projsz*[-1,1]*charsize, tmp[1,0:1]-0.0177*projsz*charsize,  $
;     ['W','E'], alignment=0.5, charsize=charsize
IsType(compass,/defined):   $
    xyouts, tmp[0,0:1]+0.015*projsz*[-1,1]*charsize, tmp[1,0:1]-0.005*projsz*charsize,   $
       ['W','E'], alignment=0.5, charsize=charsize, _extra=_extra
ELSE:
ENDCASE

InitVar, upto, TimeUnit(/min)

xyouts, 0.99, 0.95, /norm, align=1, charsize=charsize,  $
    TimeString(UT, /ymd, upto=upto, /roundt)+' UT', _extra=_extra

IF IsType(title, /defined) THEN $
    xyouts, 0.05, 0.01, /norm, charsize=charsize, title, _extra=_extra

RETURN  &  END


;*********************************************************************************************

pro vu_thomson_antifish


    infile = dialog_pickfile(path='E:\temp\time_series',filter='*.*',/read)

    ; Read input file, stop on error

    status = flt_read(infile, gg, fmt=fmt, error=error, /modify_format)
    whatis, gg

    if not status then message, error

    n_times = gg[5,0]            ; # times is assumed to be the same for all sky bins

    n_bins  = (size(gg,/dim))[1]          ; # sky bins

    ic  = reform(gg[0,*])            ; Camera number
    ra  = reform(gg[1,*])            ; Right ascensions
    dec = reform(gg[2,*])            ; Declinations
    elo = reform(gg[3,*])            ; Elongation at some specified orbit
    pa  = reform(gg[4,*])            ; Position angle at some specified orbit

    ;t_one = reform(gg[5,*])          ; First time for each sky bin

    i = 2*indgen(n_times)
    tt = gg[6+i,*]               ; Times
    af = gg[7+i,*]               ; Intensities

    ;tt = 144. + tt/86400.           ; place all times in doy for May event
    tt = 296. + tt/86400.            ; place all times in doy for October event

    ascii_bad = min(af)              ; Bad data value in ascii file
    real_bad  = BadValue(af)          ; Bad data value used for plotting

;   for i=0,n_times-1 do begin          ; Cycle on number of times
;       for j=0,n_bins-1 do begin          ; Cycle on number of locations
;         if af[i,j] le ascii_bad+1 then af[i,j] = real_bad
;       endfor
;   endfor

    ;i = where(af le ascii_bad+1)          ; Replace ascii bad by real bad
    ;if i[0] ne -1 then af[i] = real_bad

    ut=vu_getdata(hdr,ff,path='E:\temp\nv3f_files',filter='nv3f*.*',count=count)

    fil = getfilespec(hdr.file,part='name')+'_'+strmid(getfilespec(hdr.file,part='type'),1,4)
    ;print, fil
    ttmap = TimeGet(hdr.uttime, /doy)

    ff[*,*,*,1,*] = ff[*,*,*,1,*] - (2./3.)*5.

    ; Create a new point source structure (1 element for each line of sight)
    ; Set sky location for all lines of sight

    ;whatis, ra, dec, transpose([[ra],[dec]])

    ;points = vu_point_source(points, transpose([[ra],[dec]]), /deg, /location, create=n_bins)

    afdd = elo
    tttt = elo
    raaa = elo
    deee = elo
    for i=0,count-1 do begin          ; Cycle on number of nv3 files selected

    print, 'The time of the sky map in DOY is', ttmap[i]

        icount = 0
        iii = 0
        for j=0,n_bins-1 do begin          ; Cycle on the number of locations in the sky

       for ii=0,n_times-1 do begin
           if abs(tt[ii,j]-ttmap[i]) lt 51./(60.*24.) then begin

         if af[ii,j] gt ascii_bad+1 then begin
             afdd[icount] = af[ii,j]/0.55
             tttt[icount] = tt[ii,j]
             raaa[icount] = ra[j]
             deee[icount] = dec[j]
             icount = icount + 1
             iii = iii + 1
             if abs(raaa[icount]-ra[j]) lt 0.01 and abs(deee[icount]-dec[j]) lt 0.01 then icount = icount-1
         endif
           endif
       endfor

        endfor

        if iii eq 0 then print, 'Error imin = -1 There are no times in the timeseries close to the nv3 file.'

        print, 'There are', icount, ' good values in this map and', n_bins-icount, ' bad values'


        if icount ne 0 then begin

            afd = fltarr(icount)
            ttt = fltarr(icount)
            raa = fltarr(icount)
            dee = fltarr(icount)
            for ii=0,icount-1 do begin
           afd(ii) = afdd[ii]
           ttt(ii) = tttt[ii]
           raa(ii) = raaa[ii]
           dee(ii) = deee[ii]
            endfor

       points = vu_point_source(points, transpose([[raa],[dee]]), /deg, /location, create=icount)
       points = vu_point_source(points, afd, /brightness)
       points = vu_point_source(points, timeset(yr=2003,doy=ttt), /time)

       endif

       if icount eq 0 then begin
       destroyvar, points
       endif

        ; for a fisheye map
        vu_earthskymap,hdr[i],ff[*,*,*,*,i],/brightness,/degrees,xysize=[1600,1600], $
        /bmp,break=gridgen(6,norm=[0.0,1.5]),upto=timeunit(/hour), charsize=4, thick=5, charthick=5, $
        nra=360,nde=180,nrr=100,drr=0.03, /legend, $
        ;dest='E:\temp\skyview\'+fil(i), /compass, minelo=18,maxelo=110,
        dest='E:\temp\skyview\'+fil(i), zero_shift=180, compass=["E","W"], maxelo=120
            ;point_sources=points, point_size=2

        ; for a Hammer-Aitoff map
        ;vu_earthskymap,hdr[i],ff[*,*,*,*,i],minelo=18,/degrees,/brightness,xysize=[1800,1200], $
        ;/bmp,break=gridgen(6,norm=[0.0,1.5]),upto=timeunit(/hour), charsize=3, thick=5, /compass, charthick=5, $
        ;nra=360,nde=180,nrr=100,drr=0.03, /legend, $
        ;dest='E:\temp\skyview\'+fil(i), $
        ;point_sources=points, point_size=2

    endfor

return & end