pro calculate_sm_coord, nadir

common oa_coord, gci_satellite_position, gci_satellite_velocity, $
                 gci_sat_spin, geo_satellite_position, $
                 gsm_satellite_position, gsm_satellite_velocity, $
                 sm_satellite_position, sm_sat_pos_curr, sm_satellite_velocity, $
                 sm_satellite_spin, geo_n_sm, sat_re, $
                 l_shell, mlt, mlat, invlat
common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
sm_im_x = fltarr(3)
sm_im_y = fltarr(3)
sm_im_z = fltarr(3)
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
rspin   = sqrt(sm_satellite_spin(0)*sm_satellite_spin(0)+ $
               sm_satellite_spin(1)*sm_satellite_spin(1)+ $
               sm_satellite_spin(2)*sm_satellite_spin(2))
rs2 = sm_satellite_position(0)^2+sm_satellite_position(1)^2+ $
      sm_satellite_position(2)^2
rs=sqrt(rs2)
;
;------------------------------------------------------------------------------
; angle subtened by the Earth
;------------------------------------------------------------------------------
;
e_rad = asin(1.0/rs) * 180.0 / !pi
;
;------------------------------------------------------------------------------
; calculate unit vector
;------------------------------------------------------------------------------
;
sm_im_z(*) = nadir*sm_satellite_position(*)/rs

sm_im_y(0) = sm_im_z(1)*sm_satellite_spin(2)-sm_im_z(2)*sm_satellite_spin(1)
sm_im_y(1) = sm_im_z(2)*sm_satellite_spin(0)-sm_im_z(0)*sm_satellite_spin(2)
sm_im_y(2) = sm_im_z(0)*sm_satellite_spin(1)-sm_im_z(1)*sm_satellite_spin(0)

sm_im_y(*) = sm_im_y(*)/rspin

sm_im_x(0) = sm_im_y(1)*sm_im_z(2)-sm_im_y(2)*sm_im_z(1)
sm_im_x(1) = sm_im_y(2)*sm_im_z(0)-sm_im_y(0)*sm_im_z(2)
sm_im_x(2) = sm_im_y(0)*sm_im_z(1)-sm_im_y(1)*sm_im_z(0)
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro coord_transf_erad, nadir

common oa_coord, gci_satellite_position, gci_satellite_velocity, $
                 gci_sat_spin, geo_satellite_position, $
                 gsm_satellite_position, gsm_satellite_velocity, $
                 sm_satellite_position, sm_sat_pos_curr, sm_satellite_velocity, $
                 sm_satellite_spin, geo_n_sm, sat_re, $
                 l_shell, mlt, mlat, invlat
common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
;
;------------------------------------------------------------------------------
; Find the SM components of image frame in which:
;     z-axis is along satellite position
;     y-axis (spin angle) = z-axis x sm_satillite_spin
;     x-axis completes the right-hand rule
; For image looking outward from the Earth, z-axis is pointed from the
; satellite to the Earth.
;------------------------------------------------------------------------------
;
sm_im_x = fltarr(3)
sm_im_y = fltarr(3)
sm_im_z = fltarr(3)
rspin   = sqrt(sm_satellite_spin(0)*sm_satellite_spin(0)+ $
               sm_satellite_spin(1)*sm_satellite_spin(1)+ $
               sm_satellite_spin(2)*sm_satellite_spin(2))
rs2 = sm_satellite_position(0)^2+sm_satellite_position(1)^2+ $
      sm_satellite_position(2)^2
rs=sqrt(rs2)
;
;------------------------------------------------------------------------------
; angle subtened by the Earth
;------------------------------------------------------------------------------
;
e_rad = asin(1.0/rs) * 180.0 / !pi
;
;------------------------------------------------------------------------------
; component of a unit vector
;------------------------------------------------------------------------------
;
sm_im_z(*) = nadir*sm_satellite_position(*)/rs
sm_im_y(0) = sm_im_z(1)*sm_satellite_spin(2)-sm_im_z(2)*sm_satellite_spin(1)
sm_im_y(1) = sm_im_z(2)*sm_satellite_spin(0)-sm_im_z(0)*sm_satellite_spin(2)
sm_im_y(2) = sm_im_z(0)*sm_satellite_spin(1)-sm_im_z(1)*sm_satellite_spin(0)
;
;------------------------------------------------------------------------------
; make it a unit vector
;------------------------------------------------------------------------------
;
sm_im_y(*) = sm_im_y(*)/rspin
;
;------------------------------------------------------------------------------
; x = y x z
;------------------------------------------------------------------------------
;
sm_im_x(0) = sm_im_y(1)*sm_im_z(2)-sm_im_y(2)*sm_im_z(1)
sm_im_x(1) = sm_im_y(2)*sm_im_z(0)-sm_im_y(0)*sm_im_z(2)
sm_im_x(2) = sm_im_y(0)*sm_im_z(1)-sm_im_y(1)*sm_im_z(0)
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro calculate_auroral, nadir, rs_tst

common oa_coord, gci_satellite_position, gci_satellite_velocity, $
                 gci_sat_spin, geo_satellite_position, $
                 gsm_satellite_position, gsm_satellite_velocity, $
                 sm_satellite_position, sm_sat_pos_curr, sm_satellite_velocity, $
                 sm_satellite_spin, geo_n_sm, sat_re, $
                 l_shell, mlt, mlat, invlat
common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
common earth, npt, earth_x, earth_y
common auroral, auroral_x, auroral_y, auroral_pts
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
r = fltarr(3)
rp = fltarr(3)
auroral_x = fltarr(4,npt+1)
auroral_y = fltarr(4,npt+1)
auroral_pts = lonarr(4)
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
ang0 = atan(sm_satellite_position(1),sm_satellite_position(0))+!pi
del = 2.0 * !pi / npt
for i=1,4 do begin
   i1 = i - 1
   colat=float(i)
   if (i eq 1) then colat=15.
   if (i eq 2) then colat=25.
   if (i eq 3) then colat=155.
   if (i eq 4) then colat=165.
   colatr=colat*!pi/180.       ; magnetic colatitude in radian
   ic=-1
   for ii = 0,npt do begin
      ibehind=nadir
      ang = ang0+ii*del
      r(0) = sin(colatr)*cos(ang)
      r(1) = sin(colatr)*sin(ang)
      r(2) = cos(colatr)
      rp(*)=r(*)-sm_satellite_position(*)
      rp1=sqrt(total(rp*rp))    ; dist. ./. satellite.&auroral lat.
      xp=total(rp*sm_im_x)
      yp=total(rp*sm_im_y)
      zp=total(rp*sm_im_z)
      polar=asin(xp/rp1)*180./!pi
      azimuth=atan(yp,-zp)*180./!pi
      angl=sqrt(polar*polar+azimuth*azimuth) ;angle subtended in deg
      if (angl gt e_rad or rp1 lt rs_tst or nadir eq -1) then begin
         ibehind=-1                          ; this point can be seen
         ic=ic+1
         auroral_x(i1,ic)=polar
         auroral_y(i1,ic)=azimuth
      endif
   endfor   ; for ii = 0,npt do begin
   auroral_pts(i1) = ic
endfor   ; for i=1,4 do begin
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro draw_continents, x1, x0, x_max, x_min, y1, y0, y_max, y_min, white

common oa_coord, gci_satellite_position, gci_satellite_velocity, $
                 gci_sat_spin, geo_satellite_position, $
                 gsm_satellite_position, gsm_satellite_velocity, $
                 sm_satellite_position, sm_sat_pos_curr, sm_satellite_velocity, $
                 sm_satellite_spin, geo_n_sm, sat_re, $
                 l_shell, mlt, mlat, invlat
common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
xm=(x1+x0)/2.
e_radx=e_rad*(x1-x0)/(x_max-x_min)             ; e_rad scaled to 0 to 1
;
; Find ym such that ym corresponds to zero spin angle
;
ym = (y0*y_max - y1*y_min)/(y_max - y_min)
e_rady  = e_rad*(y1-y0)/(y_max-y_min)
;
; geo latitude in deg
;
geo_lat = asin(geo_satellite_position(2)/rs)*180./!pi
geo_lon = atan(geo_satellite_position(1), $
          geo_satellite_position(0))*180./!pi
Zs = (gci_satellite_position[0] * gci_sat_spin[1]) - $
     (gci_satellite_position[1] * gci_sat_spin[0])
xN = geo_n_sm[0]*sm_im_x[0]+geo_n_sm[1]*sm_im_x[1]+geo_n_sm[2]*sm_im_x[2]
yN = geo_n_sm[0]*sm_im_y[0]+geo_n_sm[1]*sm_im_y[1]+geo_n_sm[2]*sm_im_y[2]
gamma = 90.-(atan(yN,xN))*180./!pi  ; rotation(deg), clockwise from N

map_set, geo_lat, geo_lon, $
         position=[xm-e_radx,ym-e_rady,xm+e_radx,ym+e_rady], $
         /satellite, sat_p=[rs,0.,gamma], $
         /continents, con_color=white, /noborder, /noerase, /clip
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro calculate_dipole, rs_tst, nadir

common oa_coord, gci_satellite_position, gci_satellite_velocity, $
                 gci_sat_spin, geo_satellite_position, $
                 gsm_satellite_position, gsm_satellite_velocity, $
                 sm_satellite_position, sm_sat_pos_curr, sm_satellite_velocity, $
                 sm_satellite_spin, geo_n_sm, sat_re, $
                 l_shell, mlt, mlat, invlat
common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
common earth, npt, earth_x, earth_y
common dipole, dipole_x, dipole_y, dipole_pts
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
r = fltarr(3)
rp = fltarr(3)
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
ang0 = fltarr(2)
del = fltarr(2)
rc = fltarr(2)
rc(0) = 3.0
sinang0 = sqrt(1./rc(0))
ang0(0) = asin(sinang0)
del(0) = (!pi-2.0*ang0(0))/npt
rc(1) = 6.6
sinang0 = sqrt(1./rc(1))
ang0(1) = asin(sinang0)
del(1) = (!pi-2.0*ang0(1))/npt
pi_d_180 = !pi / 180.0
dipole_x = fltarr(2,4,2,npt+1)
dipole_y = fltarr(2,4,2,npt+1)
dipole_pts = lonarr(2,4,2)
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
for ii1=0,1 do begin
   jj2 = 0
   for jj1=0,270,90 do begin
      phi=jj1*pi_d_180          ; azimuthal angle in radian
      cphi=cos(phi)
      sphi=sin(phi)
      for nn1=0,1 do begin      ; do 2 hemispheres separately
         ic=-1
         for ii2=0,npt do begin
            ang=ang0(ii1)+ii2*del(ii1)
            if (nn1 eq 1) then ang = !pi - ang
            sang=sin(ang)
            r1=rc(ii1)*sang*sang       ; dipole fieldline equation
            r2=r1*r1
            r(0)=r1*sang*cphi
            r(1)=r1*sang*sphi
            r(2)=r1*cos(ang)
            rp(*)=r(*)-sm_satellite_position(*)
            rp1=sqrt(total(rp*rp))    ; dist. ./. satellite.&fieldlines
            xp=total(rp*sm_im_x)
            yp=total(rp*sm_im_y)
            zp=total(rp*sm_im_z)
            polar=asin(xp/rp1)*180./!pi
            azimuth=atan(yp,-zp)*180./!pi
            angl=sqrt(polar*polar+azimuth*azimuth)  ;angle subtend in deg
            if (angl gt e_rad or rp1 lt rs_tst or nadir eq -1) then begin
               ic=ic+1                    ; fieldline ponit can be seen
               dipole_x(ii1,jj2,nn1,ic)=polar
               dipole_y(ii1,jj2,nn1,ic)=azimuth
               if (ic gt 0) then begin
                  ad=abs(dipole_y(ii1,jj2,nn1,ic)-dipole_y(ii1,jj2,nn1,ic-1))  ; eliminate discont. in az.
                  if(ad gt 180. and dipole_y(ii1,jj2,nn1,ic-1) gt 0.)then $
                     dipole_y(ii1,jj2,nn1,ic)=dipole_y(ii1,jj2,nn1,ic)+360.
                  if(ad gt 180. and dipole_y(ii1,jj2,nn1, ic-1) lt 0.)then $
                     dipole_y(ii1,jj2,nn1,ic)=dipole_y(ii1,jj2,nn1,ic)-360.
               endif
            endif
         endfor   ; for ii2=0,npt do begin
         dipole_pts(ii1,jj2,nn1) = ic
      endfor    ; for nn1=0,1 do begin
      jj2 = jj2 + 1
   endfor   ; for jj1=0,270,90 do begin
endfor   ; for ii1=0,1 do begin
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro calculate_earth

common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
common earth, npt, earth_x, earth_y
;
;------------------------------------------------------------------------------
; no. of points in draw fieldline, circle..
; quadruple no. of point when near earth
;------------------------------------------------------------------------------
;
npt = 480
if (rs lt 1.5) then npt=npt*4 
earth_x = fltarr(npt+1)
earth_y = fltarr(npt+1)
;
;------------------------------------------------------------------------------
; calculate the earth circle then plot it
;------------------------------------------------------------------------------
;
del = 2.0 * !pi / npt
for ii = 0,npt do begin
   ang = float(ii) * del
   earth_x(ii) = e_rad * cos(ang)
   earth_y(ii) = e_rad * sin(ang)
endfor
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro calculate_re_circles, nadir, rs_tst

common oa_coord, gci_satellite_position, gci_satellite_velocity, $
                 gci_sat_spin, geo_satellite_position, $
                 gsm_satellite_position, gsm_satellite_velocity, $
                 sm_satellite_position, sm_sat_pos_curr, sm_satellite_velocity, $
                 sm_satellite_spin, geo_n_sm, sat_re, $
                 l_shell, mlt, mlat, invlat
common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
common earth, npt, earth_x, earth_y
common re_circles, re_circles_x, re_circles_y, x_connect, y_connect, $
                   circle_pts, lt_lab

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
r=fltarr(3)  &   rp=r
three_hr=!pi/4.                          ; 3 hours in radian
ang0=atan(sm_satellite_position(1),sm_satellite_position(0))+!pi
del = 2.0 * !pi / npt
rs_tst=sqrt(rs2-1.)
re_circles_x = fltarr(2,2,npt+1)
re_circles_y = fltarr(2,2,npt+1)
circle_pts = lonarr(2,2)
ncpt=16                                  ; 8 connecting lines (16 points)
x_connect=fltarr(ncpt)
y_connect=fltarr(ncpt)
lt_lab=strarr(ncpt)
;
;------------------------------------------------------------------------------
; calculate and draw circles at 3 and 6.6 Re at the equator
;
; start drawing from opposite LT
;------------------------------------------------------------------------------
;
for i = 1,2 do begin
   i1 = i - 1
   localtime0=0                         ; initial value
   if (i eq 1) then rc=3.
   if (i eq 2) then rc=6.6
   for n = 0,1 do begin
      icn = i-1
      ic  = -1
      for ii = 0,npt do begin
          ibehind = nadir
          ang = ang0 + ii * del
          if (n eq 0) then ang = ang0 - ii * del
          r(0) = rc * cos(ang)
          r(1) = rc * sin(ang)
          r(2) = 0.
          rp(*)=r(*)-sm_satellite_position(*)
          rp1=sqrt(total(rp*rp))            ; dist. ./. satellite. & equator
          xp=total(rp*sm_im_x)
          yp=total(rp*sm_im_y)
          zp=total(rp*sm_im_z)
          polar=asin(xp/rp1)*180./!pi
          azimuth=atan(yp,-zp)*180./!pi
          angl=sqrt(polar*polar+azimuth*azimuth) ;angle subtended in deg
          if (angl gt e_rad or rp1 lt rs_tst or nadir eq -1) then begin
             ibehind=-1                          ; this point can be seen
             ic=ic+1
             re_circles_x(i1,n,ic)=polar
             re_circles_y(i1,n,ic)=azimuth
             if (ic gt 0) then begin
                ad = abs(re_circles_y(i1,n,ic) - re_circles_y(i1,n,ic-1)) ; eliminate discont. in az
                if (ad gt 180.0 and re_circles_y(i1,n,ic-1) gt 0.0) then $
                   re_circles_y(i1,n,ic) = re_circles_y(i1,n,ic) + 360.0
                if (ad gt 180.0 and re_circles_y(i1,n,ic-1) lt 0.0) then $
                   re_circles_y(i1,n,ic) = re_circles_y(i1,n,ic) - 360.0
             endif
          endif
          ang_o_3=fix(ang/three_hr)
          test=three_hr*ang_o_3
          if ((ang-test) lt del and icn le (ncpt-1)) then begin ;LabelEvery3hr
             localtime=ang_o_3*3-12          ; -12 because 180 degree -> 00 LT
             if (localtime lt 0)  then localtime=localtime+24
             if (localtime gt 24) then localtime=localtime-24
             if (icn eq (i-1) or localtime ne localtime0) then begin;no repeat
                lt_lab(icn)=string(localtime,'(i2.2)')
                if (ibehind eq -1) then begin
                   x_connect(icn) = re_circles_x(i1,n,ic)
                   y_connect(icn) = re_circles_y(i1,n,ic)
                endif else begin
                   x_connect(icn)=e_rad/angl*polar   ; if behind the earth, put point
                   y_connect(icn)=e_rad/angl*azimuth   ;    on eath's rim
                endelse
                icn=icn+2
             endif
             localtime0=localtime
          endif
      endfor   ; for ii = 0,npt do begin
      circle_pts(i1,n) = ic
   endfor   ; for n = 0,1 do begin
endfor   ; for i = 1,2 do begin
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro get_scale_minmax, l_setup, image0, image1, image2

common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
common scale_minmax, data_min, data_max, del_data_minmax

;
;------------------------------------------------------------------------------
; use the minimum and maximum values set in the input parameter file
;------------------------------------------------------------------------------
;
if (l_setup.l_color_min lt l_setup.l_color_max) then begin 
   data_min = l_setup.l_color_min
   data_max = l_setup.l_color_max
   del_data_minmax = data_max - data_min
;
;------------------------------------------------------------------------------
; get the minimum and maximum values from the 3 images to be plotted
;------------------------------------------------------------------------------
;
endif else begin
   if (l_setup.l_output_type eq 2) then begin
      temp = long((3.0 * e_rad) / 8.0)
      xbin1 = 5 - temp
      xbin2 = 6 + temp
      xbin1 = xbin1 > 0
      xbin2 = xbin2 < 11
      ybin1 = 23 - temp
      ybin2 = 22 + temp
      ybin1 = ybin1 > 0
      ybin2 = ybin2 < 44
   endif else begin
      xbin1 = 0
      ybin1 = 0
      size_info = size(image0)
      xbin2 = size_info(1) - 1
      ybin2 = size_info(2) - 1
   endelse

;   get the minimum value greater than zero from the 3 image arrays

   gt0 = where(image0(xbin1:xbin2,ybin1:ybin2) gt 0, gt0cnt)
   if (gt0cnt gt 0) then min0 = min(image0(gt0)) $
   else                  min0 = 0.01
   gt0 = where(image1(xbin1:xbin2,ybin1:ybin2) gt 0, gt0cnt)
   if (gt0cnt gt 0) then min1 = min(image1(gt0)) $
   else                  min1 = 0.01
   gt0 = where(image2(xbin1:xbin2,ybin1:ybin2) gt 0, gt0cnt)
   if (gt0cnt gt 0) then min2 = min(image2(gt0)) $
   else                  min2 = 0.01
   flx_min = min0 < min1 < min2
   if (flx_min le 0.0) then flx_min = 0.01

;  get the maximum value from the 3 image arrays

   gt0 = where(image0(xbin1:xbin2,ybin1:ybin2) gt 0, gt0cnt)
   if (gt0cnt gt 0) then max0 = max(image0(gt0)) $
   else                  max0 = 0.01
   gt0 = where(image1(xbin1:xbin2,ybin1:ybin2) gt 0, gt0cnt)
   if (gt0cnt gt 0) then max1 = max(image1(gt0)) $
   else                  max1 = 0.01
   gt0 = where(image2(xbin1:xbin2,ybin1:ybin2) gt 0, gt0cnt)
   if (gt0cnt gt 0) then max2 = max(image2(gt0)) $
   else                  max2 = 0.01
   flx_max = max0 > max1 > max2

   if (l_setup.l_linear_log eq 1) then begin
      data_min = flx_min
      data_max = flx_max
   endif else begin
      data_min = alog10(flx_min)
      data_max = alog10(flx_max)
   endelse

;  round off minimum and maximum values to nearest .5

   sign = 1.0
   if (data_min lt 0.0) then sign = -1.0
   temp = abs(data_min)
   temp1 = long(temp)
   temp2 = temp - temp1
   if (temp2 gt 0.5) then begin
      if (temp2 ge 0.75) then temp2 = 1.0 else temp2 = 0.5
   endif else begin
      if (temp2 ge 0.25) then temp2 = 0.5 else temp2 = 0.0
   endelse
   data_min = sign * (float(temp1) + temp2)

   sign = 1.0
   if (data_max lt 0.0) then sign = -1.0
   temp = abs(data_max)
   temp1 = long(temp)
   temp2 = temp - temp1
   if (temp2 gt 0.5) then begin
      if (temp2 ge 0.75) then temp2 = 1.0 else temp2 = 0.5
   endif else begin
      if (temp2 ge 0.25) then temp2 = 0.5 else temp2 = 0.0
   endelse
   data_max = sign * (float(temp1) + temp2)
   if (data_max le data_min) then data_max = data_min + 0.5

   del_data_minmax = data_max - data_min

endelse
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro load_color_table

common color_table, black, white, magenta, nc_color, md_gray, $
                    color_min, color_max, del_color_minmax
;
;------------------------------------------------------------------------------
; load TIDE color table
;------------------------------------------------------------------------------
;
rgb_1   = fltarr (3, 226)
rgb   = fltarr (3, 228)
red   = fltarr (228)
blue  = fltarr (228)
green = fltarr (228)

openr, inunit, 'color.tbl', /get_lun
  readf, inunit, rgb_1
close, inunit
free_lun, inunit

for ii = 0,4 do begin
   rgb(0,ii) = rgb_1(0,ii)
   rgb(1,ii) = rgb_1(1,ii)
   rgb(2,ii) = rgb_1(2,ii)
endfor
rgb(0,5) = 0
rgb(1,5) = 0
rgb(2,5) = 175
for ii = 5,225 do begin
   rgb(0,ii+1) = rgb_1(0,ii)
   rgb(1,ii+1) = rgb_1(1,ii)
   rgb(2,ii+1) = rgb_1(2,ii)
endfor
rgb(0,227) = 255
rgb(1,227) = 0
rgb(2,227) = 255

red   = rgb(0,*)
green = rgb(1, *)
blue  = rgb(2, *)
tvlct, red, green, blue

!p.color = 1
!p.background = 0

ncolor = !D.TABLE_SIZE
color_table_size =  !D.TABLE_SIZE

white = 0
black = 1
lt_gray = 2
md_gray = 3
dk_gray = 4
nc_color = 5     ; the color to be used for no counts
color_min = 6
color_max = 226
magenta = 227
del_color_minmax = color_max - color_min
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro plot_3constant, l_setup, image0, image1, image2, nadir

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

common rlw_plot_parm, char_size, char_size_1, $
                      window_width, window_height, $
                      x_inches, y_inches
common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
common earth, npt, earth_x, earth_y
common re_circles, re_circles_x, re_circles_y, x_connect, y_connect, $
                   circle_pts, lt_lab
common dipole, dipole_x, dipole_y, dipole_pts
common auroral, auroral_x, auroral_y, auroral_pts
common scale_minmax, data_min, data_max, del_data_minmax
common color_table, black, white, magenta, nc_color, md_gray, $
                    color_min, color_max, del_color_minmax

;
;------------------------------------------------------------------------------
; define plotting setup values
;------------------------------------------------------------------------------
;
load_color_table

x_min = -48.0
x_max = 48.0
l_clip = x_min
r_clip = x_max
y_min = -180
y_max = 180
image_x_pos_start = 0.10
image_y_pos_start = 0.47
delx_image = 2.8 / x_inches   ; image size in normalized coord
dely_image = 2.8 / y_inches   ; image size in normalized coord
x_newimg_size = window_width * delx_image   ; image size in pixels
y_newimg_size = window_height * dely_image   ; image size in pixels
x_newimg_size = round(x_newimg_size)
y_newimg_size = round(y_newimg_size)
image_x_space = delx_image + 0.04
image_y_space = dely_image + 0.04
bar_x = image_x_pos_start + (2 * image_x_space) ; color bar x location
bar_x = bar_x + 0.02

;
;------------------------------------------------------------------------------
; loop through the image frames (0=H, 1=O, 2=TOTAL)
;------------------------------------------------------------------------------
;

for image_num = 0,2 do begin

;
;------------------------------------------------------------------------------
; take the log of the image and calculate image placement
;------------------------------------------------------------------------------
;

if (image_num eq 0) then begin
   mass_label = 'H'
   original_image = image0
   log_image = image0
   x_pos = image_x_pos_start
   y_pos = image_y_pos_start
endif else if (image_num eq 1) then begin
   mass_label = 'O'
   original_image = image1
   log_image = image1
   x_pos = x_pos + image_x_space
   y_pos = image_y_pos_start
endif else begin
   mass_label = 'All Masses'
   original_image = image2
   log_image = image2
   x_pos = image_x_pos_start
   y_pos = image_y_pos_start - image_y_space
endelse
no_data = where (original_image lt 0.0, no_data_counts)
no_counts = where (original_image eq 0.0, no_counts_counts)
have_data = where (original_image gt 0.0, have_data_count)
if (have_data_count gt 0) then begin
   if (l_setup.l_linear_log eq 1) then begin
      log_image[have_data] = original_image[have_data]
   endif else begin
      log_image[have_data] = alog10(original_image[have_data])
   endelse
endif

;------------------------------------------------------------------
; scale the image data to the color range
;------------------------------------------------------------------

new_image = log_image
temp = log_image

if (have_data_count gt 0) then begin
   temp[have_data] = (new_image[have_data] - data_min) / del_data_minmax
   new_image[have_data] = round (temp[have_data] * del_color_minmax)
   new_image[have_data] = new_image[have_data] + color_min
   lt_dmin = where (log_image lt data_min, lt_dmin_count)
   if (lt_dmin_count gt 0) then new_image(lt_dmin) = color_min
   gt_dmax = where (log_image gt data_max, gt_dmax_count)
   if (gt_dmax_count gt 0) then new_image(gt_dmax) = color_max
endif

if (no_counts_counts gt 0) then $
   new_image[no_counts] = nc_color     ; give no counts this color

if (no_data_counts gt 0) then $
   new_image[no_data] = md_gray     ; give no data this color

;
;------------------------------------------------------------------------------
; define some plotting variables
;------------------------------------------------------------------------------
;

new_image = congrid(new_image, 960, 3600)   ; create one degree bins

; save old plotting min,max values in degrees

x_min_old = x_min
x_max_old = x_max
y_min_old = y_min
y_max_old = y_max

; calculate new plotting min,max values in degrees based on earth size

x_max = 3.0 * e_rad
x_min = -x_max
y_max = x_max
y_min = x_min

; define clipping values when plotting earth and field lines

l_clip = x_min
r_clip = x_max
b_clip = y_min
t_clip = y_max

; determine the size of the plot image array

temp = round(x_max * 10)
plot_image_x = 2 * temp
plot_image_y = 2 * temp

; determine bin ranges to extract from data image (+/-48, +/-180)

bin_xmin = round((x_min + 48.0) * 10)
bin_xmax = round((x_max + 48.0) * 10)
bin_ymin = 1840 + round(y_min * 10.0)
bin_ymax = 1839 + round(y_max * 10.0)
bin_xmin = bin_xmin > 0
bin_xmax = bin_xmax < 959
bin_ymin = bin_ymin > 0
bin_ymax = bin_ymax < 3599
nx_bins = (bin_xmax - bin_xmin) + 1
ny_bins = (bin_ymax - bin_ymin) + 1

; fill plotting array with extracted data array

if (x_min lt -48.0) then begin   ; plotting array larger than data array
   plot_image = bytarr(plot_image_x, plot_image_y)
   plot_image(*,*)=black                  ; black background
   start_bin = long((abs(x_min) - 48.0) * 10.0) + 1
   stop_bin = start_bin + (nx_bins -1)
   plot_image(start_bin:stop_bin, 0:plot_image_x-1) = $
                                  new_image(bin_xmin:bin_xmax,bin_ymin:bin_ymax)
endif else begin     ; plotting array equal or smaller than data array
   plot_image = bytarr(nx_bins, ny_bins)
   plot_image = new_image(bin_xmin:bin_xmax,bin_ymin:bin_ymax)
endelse

; congrid plotting to proper plotting size

x_degrees = (x_max - x_min)
y_degrees = (y_max - y_min)
new_x_size = (x_degrees / (x_max-x_min)) * x_newimg_size
new_y_size = (y_degrees / (y_max-y_min)) * y_newimg_size
plot_image = congrid(plot_image, x_newimg_size, y_newimg_size)

; smooth the image, put back in set values

mask = plot_image
plot_image = smooth(plot_image, 11, /EDGE_TRUNCATE)
nc_mask = where(mask eq nc_color, nc_mask_cnts)
if (nc_mask_cnts gt 0) then plot_image(nc_mask) = mask(nc_mask)
nd_mask = where(mask eq md_gray, nd_mask_cnts)
if (nd_mask_cnts gt 0) then plot_image(nd_mask) = mask(nd_mask)
bl_mask = where(mask eq black, bl_mask_cnts)
if (bl_mask_cnts gt 0) then plot_image(bl_mask) = mask(bl_mask)

; plot the image

tv, plot_image, x_pos, y_pos, xsize=delx_image, ysize=dely_image, /normal
str_x = x_pos + (delx_image / 2.0)
str_y = (y_pos + dely_image) + 0.0175
xyouts, str_x, str_y, mass_label, charsize = char_size, align = 0.5, /normal

;
;------------------------------------------------------------------------------
; scale the image size to 0 - 1
;------------------------------------------------------------------------------
;

x0 = x_pos
y0 = y_pos
x1 = x_pos + delx_image
y1 = y_pos + dely_image

;
;------------------------------------------------------------------------------
; Add tickmarks to images
;------------------------------------------------------------------------------
;
if (nadir eq 1) then begin

   ; draw box around image

   plot, [x_min, x_max, x_max, x_min, x_min], $
         [y_min, y_min, y_max, y_max, y_min], $
         position=[x0,y0,x1,y1], $
         xrange=[x_min,x_max], yrange=[y_min,y_max], $
         xstyle=5, ystyle=5, /noclip, /noerase

   ; draw top x tickmarks

   deltck = (x_max - x_min) / 4.0
   y_tklng = (y_max - y_min) * 0.030
   for itk = 0,4 do begin
      tkx = x_min + (itk * deltck)
      tky1 = y_max
      tky2 = y_max + y_tklng
      oplot, [tkx, tkx], [tky1, tky2], /noclip
   endfor

   ; draw and label bottom x tickmarks

   deltck = (x_max - x_min) / 4.0
   for itk = 0,4 do begin
      tkx = x_min + (itk * deltck)
      tky1 = y_min
      tky2 = y_min - y_tklng
      oplot, [tkx, tkx], [tky1, tky2], /noclip
      stry = (tky2 - y_tklng) - 1.5
      if (image_num ne 0) then begin
         label = string(format = '(f6.1)', tkx)
         xyouts, tkx, stry, label, charsize = char_size, align = 0.5
      endif
   endfor
   if (image_num ne 0) then begin
      str_x = x_pos + (delx_image / 2.0)
      str_y = y_pos - 0.045
      xyouts, str_x, str_y, 'degrees', charsize = char_size, align = 0.5, $
              /normal
   endif

   ; draw and label left y tickmarks

   deltck = (y_max - y_min) / 4.0
   x_tklng = (x_max - x_min) * 0.05
   for itk = 0,4 do begin
      tky = y_min + (itk * deltck)
      tkx1 = x_min
      tkx2 = x_min - x_tklng
      oplot, [tkx1, tkx2], [tky, tky], /noclip
      if (image_num ne 1) then begin
         label = string(format = '(f6.1)', tky)
         xyouts, tkx2, tky-0.01, label, charsize = char_size, align = 1.0
      endif
   endfor
   if (image_num ne 1) then begin
      str_x = x_pos - 0.05
      str_y = y_pos + (dely_image / 2.0)
      xyouts, str_x, str_y, 'degrees', charsize = char_size, align = 0.5, $
              orientation = 90, /normal
   endif

   ; draw right y tickmarks

   deltck = (y_max - y_min) / 4.0
   for itk = 0,4 do begin
      tky = y_min + (itk * deltck)
      tkx1 = x_max
      tkx2 = x_max + x_tklng
      oplot, [tkx1, tkx2], [tky, tky], /noclip
   endfor
endif

;
;------------------------------------------------------------------------------
; draw the earth
;------------------------------------------------------------------------------
;

if (nadir eq 1) then oplot, earth_x, earth_y, color = white

;
;------------------------------------------------------------------------------
; draw circles at 3 and 6.6 Re at the equator
; start drawing from opposite LT
;------------------------------------------------------------------------------
;

for ii1 = 0,1 do begin
   for ii2 = 0,1 do begin
      n_pts = circle_pts(ii1,ii2)
      oplot, re_circles_x(ii1,ii2,0:n_pts), re_circles_y(ii1,ii2,0:n_pts), $
             color=white, clip=[l_clip, b_clip, r_clip, t_clip]
   endfor
endfor
;  draw connections between circles and label them
for i=0,14,2 do begin
   if (abs(y_connect(i)) le 90.0 and abs(y_connect(i+1)) le 90.0) then begin
      ; connection ./. L=3.3 & 6
      oplot, x_connect(i:i+1), y_connect(i:i+1), color=white
      labx=1.1*x_connect(i+1)
      laby=1.1*y_connect(i+1)
      if (labx ge x_min and labx le x_max) then  $
         if (laby ge y_min and laby le y_max) then  $
            xyouts, labx, laby, lt_lab(i+1), color=white, charsize=char_size
   endif
endfor
;
;------------------------------------------------------------------------------
; draw dipole fieldlines at L=3, 6.6, MLT=0, 6, 12, 18
;------------------------------------------------------------------------------
;

for ii1=0,1 do begin
   for jj2=0,3 do begin
      for nn1=0,1 do begin
        n_pts = dipole_pts(ii1,jj2,nn1)
        if (n_pts gt 0) then begin
           oplot, dipole_x(ii1,jj2,nn1,0:n_pts), dipole_y(ii1,jj2,nn1,0:n_pts), $
                  color=white
        endif
      endfor
   endfor
endfor

;
;------------------------------------------------------------------------------
; draw the -75, -65, 65 and 75 deg latitudes on the surface
; This section of stuff draws the auroral lines.
; start drawing from opp. LT
;------------------------------------------------------------------------------
;

for ii1 = 0,3 do begin
   n_pts = auroral_pts(ii1)
   if (n_pts gt 0) then oplot, auroral_x(ii1,0:n_pts), auroral_y(ii1,0:n_pts), color=magenta
endfor   ; for i=1,4 do begin

;
;------------------------------------------------------------------------------
; Add the continents when nadir = 1
;------------------------------------------------------------------------------
;

if (nadir eq 1) then begin
   draw_continents, x1, x0, x_max, x_min, y1, y0, y_max, y_min, white
endif


;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

endfor     ; for image_num = 0,2 do begin

;
;------------------------------------------------------------------------------
; put on color bar, no counts bar, no data bar, and labels
;------------------------------------------------------------------------------
;

put_on_color_bar, l_setup, window_width, window_height, bar_x

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

end

;==============================================================================
pro plot_3full, l_setup, image0, image1, image2, nadir


;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

common rlw_plot_parm, char_size, char_size_1, $
                      window_width, window_height, $
                      x_inches, y_inches
common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
common earth, npt, earth_x, earth_y
common re_circles, re_circles_x, re_circles_y, x_connect, y_connect, $
                   circle_pts, lt_lab
common dipole, dipole_x, dipole_y, dipole_pts
common auroral, auroral_x, auroral_y, auroral_pts
common scale_minmax, data_min, data_max, del_data_minmax
common color_table, black, white, magenta, nc_color, md_gray, $
                    color_min, color_max, del_color_minmax

;
;------------------------------------------------------------------------------
; define plotting setup values
;------------------------------------------------------------------------------
;
load_color_table

x_min = -48.0
x_max = 48.0
l_clip = x_min
r_clip = x_max
y_min = -90
y_max = 270
b_clip = -90.0
t_clip = 90.0
image_x_pos_start = 0.10 ; 48 degrees
image_y_pos_start = 0.05 ; 48 degrees
delx_image = 1.68 / x_inches   ; image size in normalized coord 48 degrees
dely_image = 6.300 / y_inches   ; image size in normalized coord
x_newimg_size = window_width * delx_image   ; image size in pixels
y_newimg_size = window_height * dely_image   ; image size in pixels
x_newimg_size = round(x_newimg_size)
y_newimg_size = round(y_newimg_size)
image_x_space = delx_image + 0.05
bar_x = image_x_pos_start + (3 * image_x_space) ; color bar x location

;
;------------------------------------------------------------------------------
; loop through the image frames (0=H, 1=O, 2=TOTAL)
;------------------------------------------------------------------------------
;

for image_num = 0,2 do begin

;
;------------------------------------------------------------------------------
; take the log of the image and calculate image placement
;------------------------------------------------------------------------------
;

if (image_num eq 0) then begin
   mass_label = 'H'
   original_image = image0
   log_image = image0
   x_pos = image_x_pos_start
   y_pos = image_y_pos_start
endif else if (image_num eq 1) then begin
   mass_label = 'O'
   original_image = image1
   log_image = image1
   x_pos = x_pos + image_x_space
   y_pos = image_y_pos_start
endif else begin
   mass_label = 'All Masses'
   original_image = image2
   log_image = image2
   x_pos = x_pos + image_x_space
   y_pos = image_y_pos_start
endelse
no_data = where (original_image lt 0.0, no_data_counts)
no_counts = where (original_image eq 0.0, no_counts_counts)
have_data = where (original_image gt 0.0, have_data_count)
if (have_data_count gt 0) then begin
   if (l_setup.l_linear_log eq 1) then begin
      log_image[have_data] = original_image[have_data]
   endif else begin
      log_image[have_data] = alog10(original_image[have_data])
   endelse
endif

;
;------------------------------------------------------------------------------
; scale the image data to the color range
;------------------------------------------------------------------------------
;

new_image = log_image
temp = log_image

if (have_data_count gt 0) then begin
   temp[have_data] = (new_image[have_data] - data_min) / del_data_minmax
   new_image[have_data] = round (temp[have_data] * del_color_minmax)
   new_image[have_data] = new_image[have_data] + color_min
   lt_dmin = where (log_image lt data_min, lt_dmin_count)
   if (lt_dmin_count gt 0) then new_image(lt_dmin) = color_min
   gt_dmax = where (log_image gt data_max, gt_dmax_count)
   if (gt_dmax_count gt 0) then new_image(gt_dmax) = color_max
endif

if (no_counts_counts gt 0) then $
   new_image[no_counts] = nc_color     ; give no counts this color

if (no_data_counts gt 0) then $
   new_image[no_data] = md_gray     ; give no data this color

;
;------------------------------------------------------------------------------
; plot the image
;------------------------------------------------------------------------------
;

new_image_1 = new_image
new_image_1_min = min(new_image_1) > color_min
le0 = where(original_image le 0.0, le0_counts)
if (le0_counts gt 0) then new_image_1(le0) = new_image_1_min
c_image=congrid(new_image,x_newimg_size,y_newimg_size)
c_image = smooth(c_image, 11, /EDGE_TRUNCATE)
mask=congrid(new_image,x_newimg_size,y_newimg_size)
nc_mask = where(mask eq nc_color, nc_mask_cnts)
if (nc_mask_cnts gt 0) then c_image(nc_mask) = mask(nc_mask)
nd_mask = where(mask eq md_gray, nd_mask_cnts)
if (nd_mask_cnts gt 0) then c_image(nd_mask) = mask(nd_mask)

tv, c_image, x_pos, y_pos, xsize=delx_image, ysize=dely_image, /normal

str_x = x_pos + (delx_image / 2.0)
str_y = (y_pos + dely_image) + 0.0175
xyouts, str_x, str_y, mass_label, charsize = char_size, align = 0.5, /normal

;
;------------------------------------------------------------------------------
; scale the image size to 0 - 1
;------------------------------------------------------------------------------
;

x0 = x_pos
y0 = y_pos
x1 = x_pos + delx_image
y1 = y_pos + dely_image

;
;------------------------------------------------------------------------------
; add tickmarks to images
;------------------------------------------------------------------------------
;

if (nadir eq 1) then begin

   ; draw box around image

   plot, [x_min, x_max, x_max, x_min, x_min], $
         [y_min, y_min, y_max, y_max, y_min], $
         position=[x0,y0,x1,y1], $
         xrange=[x_min,x_max], yrange=[y_min,y_max], $
         xstyle=5, ystyle=5, $
         /noclip, /noerase

   ; draw top x tickmarks

   deltck = (x_max - x_min) / 4.0
   y_tklng = (y_max - y_min) * 0.015
   for itk = 0,4 do begin
      tkx = x_min + (itk * deltck)
      tky1 = y_max
      tky2 = y_max + y_tklng
      oplot, [tkx, tkx], [tky1, tky2], $
             /noclip
   endfor

   ; draw and label bottom x tickmarks

   deltck = (x_max - x_min) / 4.0
   for itk = 0,4 do begin
      tkx = x_min + (itk * deltck)
      tky1 = y_min
      tky2 = y_min - y_tklng
      oplot, [tkx, tkx], [tky1, tky2], $
             /noclip
      stry = (tky2 - y_tklng) - 1.5
      if ((itk mod 2) eq 0) then begin
         label = string(format = '(f5.1)', tkx)
         xyouts, tkx, stry, label, charsize = char_size, align = 0.5
      endif
   endfor
   str_x = x_pos + (delx_image / 2.0)
   str_y = y_pos - 0.045
   xyouts, str_x, str_y, 'degrees', charsize = char_size, align = 0.5, $
           /normal

   ; draw and label left y tickmarks

   deltck = (y_max - y_min) / 4.0
   x_tklng = (x_max - x_min) * 0.075
   for itk = 0,4 do begin
      tky = y_min + (itk * deltck)
      tkx1 = x_min
      tkx2 = x_min - x_tklng
      oplot, [tkx1, tkx2], [tky, tky], /noclip
      if (image_num eq 0) then begin
         label = string(format = '(f5.1)', tky)
         xyouts, tkx2, tky-1.5, label, charsize = char_size, align = 1.0
      endif
   endfor
   if (image_num eq 0) then begin
      str_x = x_pos - 0.05
      str_y = y_pos + (dely_image / 2.0)
      xyouts, str_x, str_y, 'degrees', $
              charsize = char_size, $
              align = 0.5, $
              orientation = 90, $
              /normal
   endif

   ; draw right y tickmarks

   deltck = (y_max - y_min) / 4.0
   for itk = 0,4 do begin
      tky = y_min + (itk * deltck)
      tkx1 = x_max
      tkx2 = x_max + x_tklng
      oplot, [tkx1, tkx2], [tky, tky], $
             /noclip
   endfor
endif

;
;------------------------------------------------------------------------------
; draw the earth
;------------------------------------------------------------------------------
;

if (nadir eq 1) then oplot, earth_x, earth_y, color = white

;
;------------------------------------------------------------------------------
; calculate and draw circles at 3 and 6.6 Re at the equator
;
; start drawing from opposite LT
;------------------------------------------------------------------------------
;

for ii1 = 0,1 do begin
   for ii2 = 0,1 do begin
      n_pts = circle_pts(ii1,ii2)
      oplot, re_circles_x(ii1,ii2,0:n_pts), re_circles_y(ii1,ii2,0:n_pts), $
             color=white, clip=[l_clip, b_clip, r_clip, t_clip]
   endfor
endfor
;  draw connections between circles and label them
for i=0,14,2 do begin
   if (abs(y_connect(i)) le 90.0 and abs(y_connect(i+1)) le 90.0) then begin
      ; connection ./. L=3.3 & 6
      oplot, x_connect(i:i+1), y_connect(i:i+1), color=white
      labx=1.1*x_connect(i+1)
      laby=1.1*y_connect(i+1)
      if (labx ge x_min and labx le x_max) then  $
         if (laby ge y_min and laby le y_max) then  $
            xyouts, labx, laby, lt_lab(i+1), color=white, charsize=char_size
   endif
endfor

;
;------------------------------------------------------------------------------
; draw dipole fieldlines at L=3, 6.6, MLT=0, 6, 12, 18
;------------------------------------------------------------------------------
;

for ii1=0,1 do begin
   for jj2=0,3 do begin
      for nn1=0,1 do begin
        n_pts = dipole_pts(ii1,jj2,nn1)
        if (n_pts gt 0) then begin
           oplot, dipole_x(ii1,jj2,nn1,0:n_pts), dipole_y(ii1,jj2,nn1,0:n_pts), $
                  color=white
        endif
      endfor
   endfor
endfor

;
;------------------------------------------------------------------------------
; draw the -75, -65, 65 and 75 deg latitudes on the surface
; This section of stuff draws the auroral lines.
; start drawing from opp. LT
;------------------------------------------------------------------------------
;

for ii1 = 0,3 do begin
   n_pts = auroral_pts(ii1)
   if (n_pts gt 0) then oplot, auroral_x(ii1,0:n_pts), auroral_y(ii1,0:n_pts), color=magenta
endfor

;
;------------------------------------------------------------------------------
; Add the continents when nadir = 1
;------------------------------------------------------------------------------
;

if (nadir eq 1) then begin
   draw_continents, x1, x0, x_max, x_min, y1, y0, y_max, y_min, white
endif

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

endfor     ; for image_num = 0,2 do begin

;
;------------------------------------------------------------------------------
; put on color bar, no counts bar, no data bar, and labels
;------------------------------------------------------------------------------
;

put_on_color_bar, l_setup, window_width, window_height, bar_x

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro plot_3square, l_setup, image0, image1, image2, nadir

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

common rlw_plot_parm, char_size, char_size_1, $
                      window_width, window_height, $
                      x_inches, y_inches
common sm_coord, rs, rs2, e_rad, sm_im_x, sm_im_y, sm_im_z
common earth, npt, earth_x, earth_y
common re_circles, re_circles_x, re_circles_y, x_connect, y_connect, $
                   circle_pts, lt_lab
common dipole, dipole_x, dipole_y, dipole_pts
common auroral, auroral_x, auroral_y, auroral_pts
common scale_minmax, data_min, data_max, del_data_minmax
common color_table, black, white, magenta, nc_color, md_gray, $
                    color_min, color_max, del_color_minmax

;
;------------------------------------------------------------------------------
; define plotting setup values
;------------------------------------------------------------------------------
;
load_color_table

x_min = -48.0
x_max = 48.0
l_clip = x_min
r_clip = x_max
y_min = x_min
y_max = x_max
b_clip = y_min
t_clip = y_max
image_x_pos_start = 0.10
image_y_pos_start = 0.47
delx_image = 2.8 / x_inches   ; image size in normalized coord
dely_image = 2.8 / y_inches   ; image size in normalized coord
x_newimg_size = window_width * delx_image   ; image size in pixels
y_newimg_size = window_height * dely_image   ; image size in pixels
x_newimg_size = round(x_newimg_size)
y_newimg_size = round(y_newimg_size)
image_x_space = delx_image + 0.04
image_y_space = dely_image + 0.04
bar_x = image_x_pos_start + (2 * image_x_space) ; color bar x location
bar_x = bar_x + 0.02

;
;------------------------------------------------------------------------------
; loop through the image frames (0=H, 1=O, 2=TOTAL)
;------------------------------------------------------------------------------
;

for image_num = 0,2 do begin

;
;------------------------------------------------------------------------------
; take the log of the image and calculate image placement
;------------------------------------------------------------------------------
;

if (image_num eq 0) then begin
   mass_label = 'H'
   original_image = image0
   log_image = image0
   x_pos = image_x_pos_start
   y_pos = image_y_pos_start
endif else if (image_num eq 1) then begin
   mass_label = 'O'
   original_image = image1
   log_image = image1
   x_pos = x_pos + image_x_space
   y_pos = image_y_pos_start
endif else begin
   mass_label = 'All Masses'
   original_image = image2
   log_image = image2
   x_pos = image_x_pos_start
   y_pos = image_y_pos_start - image_y_space
endelse
no_data = where (original_image lt 0.0, no_data_counts)
no_counts = where (original_image eq 0.0, no_counts_counts)
have_data = where (original_image gt 0.0, have_data_count)
if (have_data_count gt 0) then begin
   if (l_setup.l_linear_log eq 1) then begin
      log_image[have_data] = original_image[have_data]
   endif else begin
      log_image[have_data] = alog10(original_image[have_data])
   endelse
endif

;
;------------------------------------------------------------------------------
; scale the image data to the color range
;------------------------------------------------------------------------------
;

new_image = log_image
temp = log_image

if (have_data_count gt 0) then begin
   temp[have_data] = (new_image[have_data] - data_min) / del_data_minmax
   new_image[have_data] = round (temp[have_data] * del_color_minmax)
   new_image[have_data] = new_image[have_data] + color_min
   lt_dmin = where (log_image lt data_min, lt_dmin_count)
   if (lt_dmin_count gt 0) then new_image(lt_dmin) = color_min
   gt_dmax = where (log_image gt data_max, gt_dmax_count)
   if (gt_dmax_count gt 0) then new_image(gt_dmax) = color_max
endif

if (no_counts_counts gt 0) then $
   new_image[no_counts] = nc_color     ; give no counts this color

if (no_data_counts gt 0) then $
   new_image[no_data] = md_gray     ; give no data this color

;
;------------------------------------------------------------------------------
; plot the image
;------------------------------------------------------------------------------
;

new_image_1 = new_image
new_image_1_min = min(new_image_1) > color_min
le0 = where(original_image le 0.0, le0_counts)
if (le0_counts gt 0) then new_image_1(le0) = new_image_1_min
c_image=congrid(new_image,x_newimg_size,y_newimg_size)
c_image = smooth(c_image, 11, /EDGE_TRUNCATE)
mask=congrid(new_image,x_newimg_size,y_newimg_size)
nc_mask = where(mask eq nc_color, nc_mask_cnts)
if (nc_mask_cnts gt 0) then c_image(nc_mask) = mask(nc_mask)
nd_mask = where(mask eq md_gray, nd_mask_cnts)
if (nd_mask_cnts gt 0) then c_image(nd_mask) = mask(nd_mask)

tv, c_image, x_pos, y_pos, xsize=delx_image, ysize=dely_image, /normal

str_x = x_pos + (delx_image / 2.0)
str_y = (y_pos + dely_image) + 0.0175
xyouts, str_x, str_y, mass_label, charsize = char_size, align = 0.5, /normal

;
;------------------------------------------------------------------------------
; scale the image size to 0 - 1
;------------------------------------------------------------------------------
;

x0 = x_pos
y0 = y_pos
x1 = x_pos + delx_image
y1 = y_pos + dely_image

;
;------------------------------------------------------------------------------
; Add tickmarks to images
;------------------------------------------------------------------------------
;

if (nadir eq 1) then begin

   ; draw box around image

   plot, [x_min, x_max, x_max, x_min, x_min], $
          [y_min, y_min, y_max, y_max, y_min], $
         position=[x0,y0,x1,y1], $
         xrange=[x_min,x_max], yrange=[y_min,y_max], $
         xstyle=5, ystyle=5, $
          /noclip, /noerase

   ; draw top x tickmarks

   deltck = (x_max - x_min) / 4.0
   y_tklng = (y_max - y_min) * 0.030
   for itk = 0,4 do begin
      tkx = x_min + (itk * deltck)
      tky1 = y_max
      tky2 = y_max + y_tklng
      oplot, [tkx, tkx], [tky1, tky2], /noclip
   endfor

   ; draw and label bottom x tickmarks

   deltck = (x_max - x_min) / 4.0
   for itk = 0,4 do begin
      tkx = x_min + (itk * deltck)
      tky1 = y_min
      tky2 = y_min - y_tklng
      oplot, [tkx, tkx], [tky1, tky2], $
             /noclip
      stry = (tky2 - y_tklng) - 1.5
      if (image_num ne 0) then begin
         label = string(format = '(f5.1)', tkx)
         xyouts, tkx, stry, label, charsize = char_size, align = 0.5
      endif
   endfor
   if (image_num ne 0) then begin
      str_x = x_pos + (delx_image / 2.0)
      str_y = y_pos - 0.045
      xyouts, str_x, str_y, 'degrees', charsize = char_size, $
              align = 0.5, /normal
   endif

   ; draw and label left y tickmarks

   deltck = (y_max - y_min) / 4.0
   x_tklng = (x_max - x_min) * 0.05
   for itk = 0,4 do begin
      tky = y_min + (itk * deltck)
      tkx1 = x_min
      tkx2 = x_min - x_tklng
      oplot, [tkx1, tkx2], [tky, tky], /noclip
      if (image_num ne 1) then begin
         label = string(format = '(f5.1)', tky)
         xyouts, tkx2, tky-0.01, label, charsize = char_size, align = 1.0
      endif
   endfor
   if (image_num ne 1) then begin
      str_x = x_pos - 0.05
      str_y = y_pos + (dely_image / 2.0)
      xyouts, str_x, str_y, 'degrees', charsize = char_size, align = 0.5, $
              orientation = 90, /normal
   endif

   ; draw right y tickmarks

   deltck = (y_max - y_min) / 4.0
   for itk = 0,4 do begin
      tky = y_min + (itk * deltck)
      tkx1 = x_max
      tkx2 = x_max + x_tklng
      oplot, [tkx1, tkx2], [tky, tky], /noclip
   endfor
endif

;
;------------------------------------------------------------------------------
; draw the earth
;------------------------------------------------------------------------------
;

if (nadir eq 1) then oplot, earth_x, earth_y, color = white

;
;------------------------------------------------------------------------------
; calculate and draw circles at 3 and 6.6 Re at the equator
;
; start drawing from opposite LT
;------------------------------------------------------------------------------
;

for ii1 = 0,1 do begin
   for ii2 = 0,1 do begin
      n_pts = circle_pts(ii1,ii2)
      oplot, re_circles_x(ii1,ii2,0:n_pts), re_circles_y(ii1,ii2,0:n_pts), $
             color=white, clip=[l_clip, b_clip, r_clip, t_clip]
   endfor
endfor
;  draw connections between circles and label them
for i=0,14,2 do begin
   if (abs(y_connect(i)) le 90.0 and abs(y_connect(i+1)) le 90.0) then begin
      ; connection ./. L=3.3 & 6
      oplot, x_connect(i:i+1), y_connect(i:i+1), color=white
      labx=1.1*x_connect(i+1)
      laby=1.1*y_connect(i+1)
      if (labx ge x_min and labx le x_max) then  $
         if (laby ge y_min and laby le y_max) then  $
            xyouts, labx, laby, lt_lab(i+1), color=white, charsize=char_size
   endif
endfor

;
;------------------------------------------------------------------------------
; draw dipole fieldlines at L=3, 6.6, MLT=0, 6, 12, 18
;------------------------------------------------------------------------------
;

for ii1=0,1 do begin
   for jj2=0,3 do begin
      for nn1=0,1 do begin
        n_pts = dipole_pts(ii1,jj2,nn1)
        if (n_pts gt 0) then begin
           oplot, dipole_x(ii1,jj2,nn1,0:n_pts), dipole_y(ii1,jj2,nn1,0:n_pts), $
                  color=white
        endif
      endfor
   endfor
endfor

;
;------------------------------------------------------------------------------
; draw the -75, -65, 65 and 75 deg latitudes on the surface
; This section of stuff draws the auroral lines.
; start drawing from opp. LT
;------------------------------------------------------------------------------
;

for ii1 = 0,3 do begin
   n_pts = auroral_pts(ii1)
   if (n_pts gt 0) then oplot, auroral_x(ii1,0:n_pts), auroral_y(ii1,0:n_pts), color=magenta
endfor   ; for i=1,4 do begin

;
;------------------------------------------------------------------------------
; Add the continents when nadir = 1
;------------------------------------------------------------------------------
;

if (nadir eq 1) then begin
   draw_continents, x1, x0, x_max, x_min, y1, y0, y_max, y_min, white
endif

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

endfor     ; for image_num = 0,2 do begin

;
;------------------------------------------------------------------------------
; put on color bar, no counts bar, no data bar, and labels
;------------------------------------------------------------------------------
;

put_on_color_bar, l_setup, window_width, window_height, bar_x

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

end

;==============================================================================
pro put_on_color_bar, l_setup, window_width, window_height, bar_x

common scale_minmax, data_min, data_max, del_data_minmax
common color_table, black, white, magenta, nc_color, md_gray, $
                    color_min, color_max, del_color_minmax

;
;------------------------------------------------------------------------------
; draw the color bar
;------------------------------------------------------------------------------
;
temp_cb = bytarr (2, 221)
for i = 0, 220 do begin
   j = i + 5
   temp_cb (0,i) = j
   temp_cb(1,i) = j
endfor
bar_x_size = 15
bar_y_size = window_height / 3
colorbar=congrid(temp_cb, bar_x_size, bar_y_size)
bar_w = float(bar_x_size) / float(window_width)
bar_h = float(bar_y_size) / float(window_height)
bar_y = 0.50 - (bar_h / 2.0)
tv, colorbar, bar_x, bar_y, xsize=bar_w, ysize=bar_h, /normal
box_x = [bar_x, bar_x+bar_w, bar_x+bar_w, bar_x, bar_x]
box_y = [bar_y, bar_y, bar_y+bar_h, bar_y+bar_h, bar_y]
plot, [0.0,1.0,1.0,0.0,0.0], [0.0,0.0,1.0,1.0,0.0], xstyle=4, ystyle=4, $
      /noerase, color=0
oplot, box_x, box_y, /noclip
;
;------------------------------------------------------------------------------
; label the color bar
;------------------------------------------------------------------------------
;
if (l_setup.l_linear_log eq 1) then begin
   xyouts, bar_x-0.01, bar_y+(bar_h/2.0), 'linear counts/sec', $
           alignment=0.5, charsize=char_size, color=black, $
           /normal, orientation = 90
endif else begin
   xyouts, bar_x-0.01, bar_y+(bar_h/2.0), 'log counts/sec', $
           alignment=0.5, charsize=char_size, color=black, $
           /normal, orientation = 90
endelse
; tick marks & tick label
tck_x1 = bar_x + bar_w
tck_x2 = tck_x1 + 0.01
str_x = tck_x2 + 0.005
;if (l_setup.l_linear_log eq 1) then begin
;   if (data_max - data_min le 3.0) then val = 1.0 else val = 2.0
;   num_ticks = (data_max - data_min) / val
;   dely_tick = bar_h / num_ticks
;   for iticks = 0,num_ticks do begin
;      tick_value = data_min + (iticks * val)
;      tick_y = bar_y + (iticks * dely_tick)
;      oplot, [tck_x1, tck_x2], [tick_y, tick_y], /noclip
;      xyouts, str_x, tick_y-0.005, string(tick_value,'(f4.1)'), $
;              alignment=0.0, charsize=char_size, color=black, $
;              /normal
;   endfor
;endif else begin
   num_ticks = (data_max - data_min) / 0.50
   dely_tick = bar_h / num_ticks
   for iticks = 0,num_ticks do begin
      tick_value = data_min + (iticks * 0.50)
      tick_y = bar_y + (iticks * dely_tick)
      oplot, [tck_x1, tck_x2], [tick_y, tick_y], /noclip
      xyouts, str_x, tick_y-0.005, string(tick_value,'(f4.1)'), $
              alignment=0.0, charsize=char_size, color=black, $
              /normal
   endfor
;endelse
;
;------------------------------------------------------------------------------
; put on no counts bar and label
;------------------------------------------------------------------------------
;
nocounts = bytarr(15,15)
nocounts(*,*) = nc_color
bar_w = 15 / float(window_width)
bar_h = 15 / float(window_height)
bar_y = (0.50 - (bar_h / 2.0)) - 0.20
tv, nocounts, bar_x, bar_y, xsize=bar_w, ysize=bar_h, /normal
box_x = [bar_x, bar_x+bar_w, bar_x+bar_w, bar_x, bar_x]
box_y = [bar_y, bar_y, bar_y+bar_h, bar_y+bar_h, bar_y]
oplot, box_x, box_y, /noclip
xyouts, bar_x+.01, bar_y-0.015, 'no counts', $
        alignment=0.5, charsize=char_size, /normal
;
;------------------------------------------------------------------------------
; put on no data bar and label
;------------------------------------------------------------------------------
;
nodata = bytarr(15,15)
nodata(*,*) = md_gray
bar_w = 15 / float(window_width)
bar_h = 15 / float(window_height)
bar_y = (0.50 - (bar_h / 2.0)) - 0.25
tv, nodata, bar_x, bar_y, xsize=bar_w, ysize=bar_h, /normal
box_x = [bar_x, bar_x+bar_w, bar_x+bar_w, bar_x, bar_x]
box_y = [bar_y, bar_y, bar_y+bar_h, bar_y+bar_h, bar_y]
oplot, box_x, box_y, /noclip
xyouts, bar_x+.01, bar_y-0.015, 'no data', $
        alignment=0.5, charsize=char_size, /normal
;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro plot_3full_nodata, l_setup, image0, image1, image2

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

common rlw_plot_parm, char_size, char_size_1, $
                      window_width, window_height, $
                      x_inches, y_inches
common scale_minmax, data_min, data_max, del_data_minmax
common color_table, black, white, magenta, nc_color, md_gray, $
                    color_min, color_max, del_color_minmax

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
load_color_table

data_min = -1.0
data_max = 1.0
del_data_minmax = data_max - data_min
;
;------------------------------------------------------------------------------
; define plotting setup values
;------------------------------------------------------------------------------
;

x_min = -48.0
x_max = 48.0
y_min = -90
y_max = 270
image_x_pos_start = 0.10 ; 48 degrees
image_y_pos_start = 0.05 ; 48 degrees
delx_image = 1.68 / x_inches   ; image size in normalized coord 48 degrees
dely_image = 6.300 / y_inches   ; image size in normalized coord
x_newimg_size = window_width * delx_image   ; image size in pixels
y_newimg_size = window_height * dely_image   ; image size in pixels
x_newimg_size = round(x_newimg_size)
y_newimg_size = round(y_newimg_size)
image_x_space = delx_image + 0.05
bar_x = image_x_pos_start + (3 * image_x_space) ; color bar x location

;
;------------------------------------------------------------------------------
; loop through the image frames (0=H, 1=O, 2=TOTAL)
;------------------------------------------------------------------------------
;

for image_num = 0,2 do begin


;
;------------------------------------------------------------------------------
; set the image to no data color
;------------------------------------------------------------------------------
;

if (image_num eq 0) then begin
   mass_label = 'H'
   original_image = image0
   x_pos = image_x_pos_start
   y_pos = image_y_pos_start
endif else if (image_num eq 1) then begin
   mass_label = 'O'
   original_image = image1
   x_pos = x_pos + image_x_space
   y_pos = image_y_pos_start
endif else begin
   mass_label = 'All Masses'
   original_image = image2
   x_pos = x_pos + image_x_space
   y_pos = image_y_pos_start
endelse
new_image = original_image
new_image[*,*] = md_gray     ; give no data this color

;
;------------------------------------------------------------------------------
; plot the image
;------------------------------------------------------------------------------
;

c_image=congrid(new_image,x_newimg_size,y_newimg_size)

tv, c_image, x_pos, y_pos, xsize=delx_image, ysize=dely_image, /normal

str_x = x_pos + (delx_image / 2.0)
str_y = (y_pos + dely_image) + 0.0175
xyouts, str_x, str_y, mass_label, charsize = char_size, align = 0.5, /normal

;
;------------------------------------------------------------------------------
; scale the image size to 0 - 1
;------------------------------------------------------------------------------
;

x0 = x_pos
y0 = y_pos
x1 = x_pos + delx_image
y1 = y_pos + dely_image

;
;------------------------------------------------------------------------------
; add tickmarks to images
;------------------------------------------------------------------------------
;

; draw box around image

plot, [x_min, x_max, x_max, x_min, x_min], $
      [y_min, y_min, y_max, y_max, y_min], $
      position=[x0,y0,x1,y1], $
      xrange=[x_min,x_max], yrange=[y_min,y_max], $
      xstyle=5, ystyle=5, $
      /noclip, /noerase

; draw top x tickmarks

deltck = (x_max - x_min) / 4.0
y_tklng = (y_max - y_min) * 0.015
for itk = 0,4 do begin
   tkx = x_min + (itk * deltck)
   tky1 = y_max
   tky2 = y_max + y_tklng
   oplot, [tkx, tkx], [tky1, tky2], $
          /noclip
endfor

; draw and label bottom x tickmarks

deltck = (x_max - x_min) / 4.0
for itk = 0,4 do begin
   tkx = x_min + (itk * deltck)
   tky1 = y_min
   tky2 = y_min - y_tklng
   oplot, [tkx, tkx], [tky1, tky2], $
          /noclip
   stry = (tky2 - y_tklng) - 1.5
   if ((itk mod 2) eq 0) then begin
      label = string(format = '(f5.1)', tkx)
      xyouts, tkx, stry, label, charsize = char_size, align = 0.5
   endif
endfor
str_x = x_pos + (delx_image / 2.0)
str_y = y_pos - 0.045
xyouts, str_x, str_y, 'degrees', charsize = char_size, align = 0.5, $
        /normal

; draw and label left y tickmarks

deltck = (y_max - y_min) / 4.0
x_tklng = (x_max - x_min) * 0.075
for itk = 0,4 do begin
   tky = y_min + (itk * deltck)
   tkx1 = x_min
   tkx2 = x_min - x_tklng
   oplot, [tkx1, tkx2], [tky, tky], /noclip
   if (image_num eq 0) then begin
      label = string(format = '(f5.1)', tky)
      xyouts, tkx2, tky-1.5, label, charsize = char_size, align = 1.0
   endif
endfor
if (image_num eq 0) then begin
   str_x = x_pos - 0.05
   str_y = y_pos + (dely_image / 2.0)
   xyouts, str_x, str_y, 'degrees', $
           charsize = char_size, $
           align = 0.5, $
           orientation = 90, $
           /normal
endif

; draw right y tickmarks

deltck = (y_max - y_min) / 4.0
for itk = 0,4 do begin
   tky = y_min + (itk * deltck)
   tkx1 = x_max
   tkx2 = x_max + x_tklng
   oplot, [tkx1, tkx2], [tky, tky], $
          /noclip
endfor

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

endfor     ; for image_num = 0,2 do begin

;
;------------------------------------------------------------------------------
; put on color bar, no counts bar, no data bar, and labels
;------------------------------------------------------------------------------
;

put_on_color_bar, l_setup, window_width, window_height, bar_x

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end

;==============================================================================
pro plot_3consqr_nodata, l_setup, image0, image1, image2

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

common rlw_plot_parm, char_size, char_size_1, $
                      window_width, window_height, $
                      x_inches, y_inches
common scale_minmax, data_min, data_max, del_data_minmax
common color_table, black, white, magenta, nc_color, md_gray, $
                    color_min, color_max, del_color_minmax

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
load_color_table

data_min = -1.0
data_max = 1.0
del_data_minmax = data_max - data_min
;
;------------------------------------------------------------------------------
; define plotting setup values
;------------------------------------------------------------------------------
;

x_min = -48.0
x_max = 48.0
y_min = x_min
y_max = x_max
image_x_pos_start = 0.10
image_y_pos_start = 0.47
delx_image = 2.8 / x_inches   ; image size in normalized coord
dely_image = 2.8 / y_inches   ; image size in normalized coord
x_newimg_size = window_width * delx_image   ; image size in pixels
y_newimg_size = window_height * dely_image   ; image size in pixels
x_newimg_size = round(x_newimg_size)
y_newimg_size = round(y_newimg_size)
image_x_space = delx_image + 0.04
image_y_space = dely_image + 0.04
bar_x = image_x_pos_start + (2 * image_x_space) ; color bar x location
bar_x = bar_x + 0.02

;
;------------------------------------------------------------------------------
; loop through the image frames (0=H, 1=O, 2=TOTAL)
;------------------------------------------------------------------------------
;

for image_num = 0,2 do begin


;
;------------------------------------------------------------------------------
; set the image to no data color
;------------------------------------------------------------------------------
;

if (image_num eq 0) then begin
   mass_label = 'H'
   original_image = image0
   x_pos = image_x_pos_start
   y_pos = image_y_pos_start
endif else if (image_num eq 1) then begin
   mass_label = 'O'
   original_image = image1
   x_pos = x_pos + image_x_space
   y_pos = image_y_pos_start
endif else begin
   mass_label = 'All Masses'
   original_image = image2
   x_pos = image_x_pos_start
   y_pos = image_y_pos_start - image_y_space
endelse
new_image = original_image
new_image[*,*] = md_gray     ; give no data this color

;
;------------------------------------------------------------------------------
; plot the image
;------------------------------------------------------------------------------
;

c_image=congrid(new_image,x_newimg_size,y_newimg_size)

tv, c_image, x_pos, y_pos, xsize=delx_image, ysize=dely_image, /normal

str_x = x_pos + (delx_image / 2.0)
str_y = (y_pos + dely_image) + 0.0175
xyouts, str_x, str_y, mass_label, charsize = char_size, align = 0.5, /normal

;
;------------------------------------------------------------------------------
; scale the image size to 0 - 1
;------------------------------------------------------------------------------
;

x0 = x_pos
y0 = y_pos
x1 = x_pos + delx_image
y1 = y_pos + dely_image

;
;------------------------------------------------------------------------------
; add tickmarks to images
;------------------------------------------------------------------------------
;

; draw box around image

plot, [x_min, x_max, x_max, x_min, x_min], $
      [y_min, y_min, y_max, y_max, y_min], $
      position=[x0,y0,x1,y1], $
      xrange=[x_min,x_max], yrange=[y_min,y_max], $
      xstyle=5, ystyle=5, $
      /noclip, /noerase

; draw top x tickmarks

deltck = (x_max - x_min) / 4.0
y_tklng = (y_max - y_min) * 0.030
for itk = 0,4 do begin
   tkx = x_min + (itk * deltck)
   tky1 = y_max
   tky2 = y_max + y_tklng
   oplot, [tkx, tkx], [tky1, tky2], $
          /noclip
endfor

; draw and label bottom x tickmarks

deltck = (x_max - x_min) / 4.0
for itk = 0,4 do begin
   tkx = x_min + (itk * deltck)
   tky1 = y_min
   tky2 = y_min - y_tklng
   oplot, [tkx, tkx], [tky1, tky2], $
          /noclip
   stry = (tky2 - y_tklng) - 1.5
   if (image_num ne 0) then begin
      label = string(format = '(f5.1)', tkx)
      xyouts, tkx, stry, label, charsize = char_size, align = 0.5
   endif
endfor
if (image_num ne 0) then begin
   str_x = x_pos + (delx_image / 2.0)
   str_y = y_pos - 0.045
   xyouts, str_x, str_y, 'degrees', charsize = char_size, align = 0.5, $
           /normal
endif

; draw and label left y tickmarks

deltck = (y_max - y_min) / 4.0
x_tklng = (x_max - x_min) * 0.05
for itk = 0,4 do begin
   tky = y_min + (itk * deltck)
   tkx1 = x_min
   tkx2 = x_min - x_tklng
   oplot, [tkx1, tkx2], [tky, tky], /noclip
   if (image_num ne 1) then begin
      label = string(format = '(f5.1)', tky)
      xyouts, tkx2, tky-0.01, label, charsize = char_size, align = 1.0
   endif
endfor
if (image_num ne 1) then begin
   str_x = x_pos - 0.05
   str_y = y_pos + (dely_image / 2.0)
   xyouts, str_x, str_y, 'degrees', $
           charsize = char_size, $
           align = 0.5, $
           orientation = 90, $
           /normal
endif

; draw right y tickmarks

deltck = (y_max - y_min) / 4.0
for itk = 0,4 do begin
   tky = y_min + (itk * deltck)
   tkx1 = x_max
   tkx2 = x_max + x_tklng
   oplot, [tkx1, tkx2], [tky, tky], $
          /noclip
endfor

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;

endfor     ; for image_num = 0,2 do begin

;
;------------------------------------------------------------------------------
; put on color bar, no counts bar, no data bar, and labels
;------------------------------------------------------------------------------
;

put_on_color_bar, l_setup, window_width, window_height, bar_x

;
;------------------------------------------------------------------------------
;------------------------------------------------------------------------------
;
end