pro polarplot3,hr,min,sec,image0

;  This procedure creates a geomagnetic coordinate grid
;  with overlain coastline shapes that will be used 
;  for viewing allsky images from a variety of polar cap
;  ground stations.  UT time is specified with variables
;  hr (hours), min (minutes) and sec (seconds).
;  The output image0 is the 800 by 800 byte array that is
;  passed to procedure plot_image.pro
;
;  RAD SRI International 8/30/94
;

   degrad = !pi/180.
   image0 = bytarr(800,800)
   grid = bytarr(800,800)

;  Load previously calculated geomagnetic MLAT and MLT lookup arrays that
;  have been calculated for 00 UT.

   loadf,'mlat0.out',730,798,0,0,mlat1
   loadf,'mlt0.out',730,798,0,0,mlt1

;  Transform the MLT array to correspond to the current UT time
   UT = float(hr) + float(min)/60. + float(sec)/3600.
   fac = !pi/12.
   sum = 0.0
   param = [0.09, 0.14,-0.09, 0.04, 0.02, 0.02,-0.02, $
            0.00, 0.01, 0.01, 0.00, 0.00, 0.00, 0.00, $
            0.00, 0.00, 0.00, 0.00, 0.00]
   for n = 1,8 do begin
     sum = sum + param((2*n)-1)*sin(n*(UT*fac)) + $
                 param((2*n))  *cos(n*(UT*fac))
   endfor
   correction = param(0) + sum
   mlt1 = mlt1 + UT + correction
   ind0 = where(mlt1 ge 24, count)
   if count gt 0 then mlt1(ind0) = mlt1(ind0) - 24.0

;  Load latitude and longitude arrays for polar coastlines
   loadf,'coastlat.out',1,19584,0,0,glat
   loadf,'coastlon.out',1,19584,0,0,glon
   mlat0 = fltarr(19584)
   mlt0 = fltarr(19584)
   for i = 0,19583 do begin
     xxxx = glon(i) * (729./360.)
     yyyy = (90. - glat(i)) * (797./47.5)
     if (xxxx le 729) and (yyyy le 797) then begin
       mlat0(i) = mlat1(xxxx,yyyy)
       mlt0(i) = mlt1(xxxx,yyyy)
     endif
   endfor

;  Convert geomagnetic continent coordinates to 800 X 800 PIXEL polar grid.
   for i = 0,19583 do begin
     radius = (90. - mlat0(i))*(799./40.)
     theta = (mlt0(i)*(360./24.)) - 90.
     x = 399. + radius * cos(theta * degrad)
     y = 399. + radius * sin(theta * degrad)
     if ((x lt 0) or (x gt 799) or (y lt 0) or (y gt 799)) then goto,$ 
     continue0
     image0(x,y) = 159 - 32 * (i mod 2)
     continue0:
   endfor

;  Create grid overlay with latitude contours at 70 and 80 degrees, 
;  and cross hairs for noon-midnight and dusk-dawn meridians.
   for colat = 20,40,20 do begin
     last = colat * 100
     arc = 360./float(last)
     rad = (float(colat) * 10.) - 1.
     for i = 0,last do begin
       theta = i * arc*degrad
       x = 399. + rad * cos(theta)
       y = 399. + rad * sin(theta)
       grid(x,y) = 255
     endfor
   endfor
   grid(399,*) = 255
   grid(*,399) = 255
   grid(0:1,*)=255 & grid(798:799,*)=255
   grid(*,0:1)=255 & grid(*,798:799)=255
   
   image0 = rebin(rebin(image0,400,400),800,800)
   image0 = image0 or grid

   end