pro draw_petal,l2a,filename,_extra=extra,west_ct=west_ct,alt_ct=alt_ct,nan_blue=nan_blue,min_blue=min_blue,cam=cam, $
                                bottom=bottom,top=top,linecolor=linecolor,coastcolor=coastcolor,min_lat=min_lat

  ;Set up west coast color table
  if n_elements(min_blue) eq 0 then min_blue=64
  if n_elements(min_lat) eq 0 then min_lat=55

; Create the color table according to John McNabb
  red=indgen(256)
  green=indgen(256)
  blue=indgen(256)

  for i=0,255 do begin
    blue[i]=min_blue+1.75*i
    if (blue[i] gt 255) then begin
       blue[i]=255
    endif
  endfor

;  black
  red[0]=0
  green[0]=0
  blue[0]=0

;  red
  red[250]=255
  green[250]=0
  blue[250]=0

;  green
  red[251]=0
  green[251]=255
  blue[251]=0

; blue
  red[252]=0
  green[252]=0
  blue[252]=255

;option 2
; grid
  red[253]=91
  green[253]=91
  blue[253]=239
  
; coast
  gcolor=150
  red[254]=gcolor
  green[254]=gcolor
  blue[254]=gcolor
  coast_color=254

;  white
  red[255]=255
  green[255]=255
  blue[255]=255
  
  if ~keyword_set(nan_blue) then blue[0]=0
  
  ;Load the map
  restore,!common_input_data+'/coast.sav'
  south=(l2a.hemisphere eq 'S')
  kmperpix=l2a.km_per_pixel
  plot_project_lambert,lon=coast[0,*]-l2a.center_lon,lat=coast[1,*],x=x_coast,y=y_coast,south=south,kmperpix=kmperpix

  ;Set up z buffer
  old_plot=!d.name
  set_plot,'z'
  device,set_resolution=size(*l2a.albedo,/dim)
  map_size=(size(*l2a.albedo,/dim))[0]
  if n_elements(alt_ct) eq 0 then alt_ct=1
  if ~keyword_set(west_ct) then loadct,alt_ct else tvlct,red,green,blue
  
  if n_elements(bottom) eq 0 then bottom=0
  if n_elements(top) eq 0 then top=100
  ofs=1

  ;------------------- draw the data -------------------------  
  in50=*l2a.albedo
  
  w=where(~finite(in50),nw, COMPLEMENT=finite_pixels_index)     ; index of bad pixels
  top = 150
  index = WHERE( finite(in50) AND in50 gt 100 )
  num_max = N_ELEMENTS( index )
  in50[finite_pixels_index] = ((in50[finite_pixels_index]-bottom)>0<(top-bottom))+ofs
  maxdata=max(in50,/nan)              ; Max value of valid data
  track_loc = WHERE(in50 GE 1 AND in50 LT 8 )
  track_loc2 = WHERE(in50 EQ 0 OR ~finite(in50))

  if ~keyword_set(square) then begin
    s=size(in50,/dim)
    x=double(ramp(-s[0]/2,s[0]/2,s,0))
    y=double(ramp(-s[1]/2,s[1]/2,s,1))
    w=where(sqrt(x^2+y^2) gt s[0]/2,nw)
  end

   in50_min = MIN(in50[finite_pixels_index], /NAN)
   in50_max = MAX(in50[finite_pixels_index], /NAN)

  int_image = HIST_EQUAL(in50, minv=8, maxv=in50_max, TOP = 240)
  max_index = WHERE( int_image GT 249 )
  min_index = WHERE( int_image LE 8  )

  ; This sets the original background pixels to dark blue
  IF track_loc2[0] NE -1 THEN int_image[track_loc] = 1

  tv, int_image

  ;draw the coastlines
  plots,/device,x_coast-l2a.bbox[0],y_coast-l2a.bbox[1],color=coastcolor

  ;Prepare the circle of latitude arrays
  q=findgen(629)/100.0;
  c=cos(q);
  s=sin(q);
 
  ;draw circles of latitude
  r=dblarr(91)
  
  for i=85,30,-5 do begin
    plot_project_lambert,lat=i,lon=dindgen(360)-180,x=x_lat,y=y_lat,kmperpix=kmperpix
    plots,/device,x_lat-l2a.bbox[0],y_lat-l2a.bbox[1],color=linecolor,linestyle=1
    xyouts,/device,align=1,x_lat[90]-l2a.bbox[0],y_lat[90]-l2a.bbox[1],color=linecolor,string(i,format='(%"%02d")'),charsize=2    
    xyouts,/device,align=0,x_lat[270]-l2a.bbox[0],y_lat[270]-l2a.bbox[1],color=linecolor,string(i,format='(%"%02d")'),charsize=2    
  end

  ;Draw midnight sun line  
  junk=sunvec(l2a.stack_start_time,delta=delta)
  if south then mn_sun_lat=90d +delta else mn_sun_lat=90d -delta
  plot_project_lambert,lat=mn_sun_lat,lon=dindgen(360)-180,x=x_lat,y=y_lat,kmperpix=kmperpix
  plots,/device,x_lat-l2a.bbox[0],y_lat-l2a.bbox[1],color=linecolor,linestyle=2
   
  xyouts,/device,charsize=2,align=0,x_lat[270]-l2a.bbox[0],y_lat[270]-l2a.bbox[1]+20,string(mn_sun_lat,format='(%"%4.1f")'),color=linecolor
;  if mn_sun_lat gt 90 then label='polar night' else label='midnight sun'
;  xyouts,/device,           align=1,x_pole-y_50-2,y_pole-y_50-r_mn_sun,label,color=linecolor

  ;draw lines of longitude
  top=dblarr(361)
  top[*]=85;
  top[[0,90,180,270]]=90;
  for i=0,360,15 do begin
    plot_project_lambert,lat=indgen(top[i]-30+1)+30,lon=i-l2a.center_lon,x=x_lat,y=y_lat,kmperpix=kmperpix
    if(top[i] eq 90) then begin
      if(i gt 180) then begin
        lon=i-360
        hem_ew='W' 
      end else begin
        lon=i
        hem_ew='E'
      end
      if(i eq 0) then hem_ew='';
      if(i eq 180) then hem_ew =''; 
      xyouts,/device,x_lat[56]-l2a.bbox[0],y_lat[56]-l2a.bbox[1],string(format='(%"%d%s")',abs(lon),hem_ew),charsize=2
    end
    plots,/device,x_lat-l2a.bbox[0],y_lat-l2a.bbox[1],color=linecolor,linestyle=1
  end

  xyouts,/device,charsize=2,20,20,string(format='(%"AIM Orbit %d")',l2a.aim_orbit_number)
  xyouts,/device,charsize=2,20,45,strmid(usec2vms(l2a.stack_start_time),1,16)
  xyouts,/device,align=1.0,charsize=2,map_size-24,45,'AIM PMC Presence'
  xyouts,/device,align=1.0,charsize=2,map_size-24,20,'NASA/HU/VT/CU LASP'

  saveimage,/png,filename
  set_plot,old_plot
  
end