; produces eps images of map data with many additional functions
; (e.g., for plotting vectors, lines, particles, text annotation, etc.)
;
; written by Christoph Federrath, 2010-2016
;

;
; returns a Polygon object
; with the shape of an ellipse
;
function get_ellipse_obj, PTS=pts, SIZE=size, _EXTRA=extra
    if not keyword_set(size) then size = 0.5
    oEllipse = OBJ_NEW('IDLgrModel')
    for i = 0, n_elements(pts[0,*])-1 do begin
       xp = [pts[0,i]-size,pts[0,i]+size,pts[0,i]+size,pts[0,i]-size]
       yp = [pts[1,i]+size,pts[1,i]+size,pts[1,i]-size,pts[1,i]-size]
       oPolygon = obj_new('IDLgrPolygon', xp, yp, LINESTY='0', _EXTRA=extra)
       oEllipse->Add, oPolygon
    endfor
    return, oEllipse
end

;
; returns a Polygon object
; with the shape of an arrow
;
function get_arrow_obj, ARROW=arrow, _EXTRA=extra
   if not keyword_set(arrow) then arrow = [0,0]
   lfac = 0.3
   delta = 20./180.*!pi
   sind = sin(delta)
   cosd = cos(delta)
   xs = arrow[0]
   ys = arrow[1]
   xp = xs + lfac*(-xs*cosd+ys*sind)
   yp = ys + lfac*(-xs*sind-ys*cosd)
   xq = xs + lfac*(-xs*cosd-ys*sind)
   yq = ys + lfac*( xs*sind-ys*cosd)
   xs -= xs/2.
   xp -= xs
   xq -= xs
   ys -= ys/2.
   yp -= ys
   yq -= ys
   fac1 = 0.022
   fac2 = 0.7
   x1 = -xs-fac1*ys
   y1 = -ys+fac1*xs
   x2 = fac2*xs-fac1*ys
   y2 = fac2*ys+fac1*xs
   x3 = fac2*xs+fac1*ys
   y3 = fac2*ys-fac1*xs
   x4 = -xs+fac1*ys
   y4 = -ys-fac1*xs
   x = [xs,xp,x3,x4,x1,x2,xq]
   y = [ys,yp,y3,y4,y1,y2,yq]
   oVector = obj_new('IDLgrPolygon', x, y, _EXTRA=extra)
   return, oVector
end
;
; returns a Polygon object
; with the shape of a vector
;
function get_vector_obj, VECTOR=vector, _EXTRA=extra
   if not keyword_set(vector) then vector = [0,0]
   lfac = 0.3
   delta = 20./180.*!pi
   sind = sin(delta)
   cosd = cos(delta)
   xs = vector[0]
   ys = vector[1]
   xp = xs + lfac*(-xs*cosd+ys*sind)
   yp = ys + lfac*(-xs*sind-ys*cosd)
   xq = xs + lfac*(-xs*cosd-ys*sind)
   yq = ys + lfac*( xs*sind-ys*cosd)
   xs -= xs/2.
   xp -= xs
   xq -= xs
   ys -= ys/2.
   yp -= ys
   yq -= ys
   fac1 = 0.022
   fac2 = 0.7
   x1 = -xs-fac1*ys
   y1 = -ys+fac1*xs
   x2 = fac2*xs-fac1*ys
   y2 = fac2*ys+fac1*xs
   x3 = fac2*xs+fac1*ys
   y3 = fac2*ys-fac1*xs
   x4 = -xs+fac1*ys
   y4 = -ys-fac1*xs
   x = [xs,xp,x3,x4,x1,x2,xq]
   y = [ys,yp,y3,y4,y1,y2,yq]
   oVector = obj_new('IDLgrPolygon', x, y, _EXTRA=extra)
   return, oVector
end

;
; returns a Polygon object
; with the shape of a star
;
function get_star_obj, SIZE=size, _EXTRA=extra
   if not keyword_set(size) then size = 3.
   a = 1.0*size
   b = 0.25*size
   x = [0., a, b, 0., -b, -a, -b, 0., b, a]
   y = [0., 0., b, a, b, 0., -b, -a, -b, 0.]
   oStar = obj_new('IDLgrPolygon', x, y, LINESTY=6, _EXTRA=extra)
   return, oStar
end

;
; returns a Polygon object
; with the shape of a filled circle
;
function get_filled_circle_obj, SIZE=size, POINTS=points, _EXTRA=extra
   if not keyword_set(size)   then size   = 3.
   if not keyword_set(points) then points = 50
   a = 2.*!pi*dindgen(points)/double(points-1)
   x = sin(a)*size
   y = cos(a)*size
   oCircle = obj_new('IDLgrPolygon', x, y, LINESTY=6, _EXTRA=extra)
   return, oCircle
end

;
; returns a Polyline object
; with the shape of a circle
;
function get_circle_obj, SIZE=size, POINTS=points, _EXTRA=extra
   if not keyword_set(size)   then size   = 3.
   if not keyword_set(points) then points = 50
   a = 2.*!pi*dindgen(points)/double(points-1)
   x = sin(a)*size
   y = cos(a)*size
   xin = sin(a)*size * 0.85 ; control thickness
   yin = cos(a)*size * 0.85 ; control thickness
   oCircle = OBJ_NEW('IDLgrModel')
   for i = 0, points-2 do begin
      xp = [xin[i], x[i], x[i+1], xin[i+1]]
      yp = [yin[i], y[i], y[i+1], yin[i+1]]
      oPolygon = obj_new('IDLgrPolygon', xp, yp, LINESTY='0', _EXTRA=extra)
      oCircle->Add, oPolygon
   endfor
   return, oCircle
end

;
; returns a Polygon object
; with the shape of a cross
;
function get_cross_obj, SIZE=size, _EXTRA=extra
   if not keyword_set(size) then size = 3.
   a = 1.0*size
   b = 0.15*size
   x = [0., a, b, b, -b, -b, -a, -a, -b, -b, b, b, a, a]
   y = [0., b, b, a, a, b, b, -b, -b, -a, -a, -b, -b, b]
   oCross = obj_new('IDLgrPolygon', x, y, LINESTY=6, _EXTRA=extra)
   return, oCross
end

;
; returns a Crosshair object
;
function get_crosshair_obj, SIZE=size, _EXTRA=extra
   if not keyword_set(size) then size = 3.
   oCrosshair = OBJ_NEW('IDLgrModel')
   a = 1.0*size
   b = 0.3*size
   shift = 1.0
   x = [-a, -a, a, a]
   y = [-b, b, b, -b]
   oPolygon1 = obj_new('IDLgrPolygon', x+shift*size+a, y, LINESTY=6, _EXTRA=extra)
   oPolygon2 = obj_new('IDLgrPolygon', x-shift*size-a, y, LINESTY=6, _EXTRA=extra)
   x = [-b, -b, b, b]
   y = [-a, a, a, -a]
   oPolygon3 = obj_new('IDLgrPolygon', x, y+shift*size+a, LINESTY=6, _EXTRA=extra)
   oPolygon4 = obj_new('IDLgrPolygon', x, y-shift*size-a, LINESTY=6, _EXTRA=extra)
   oCrosshair->Add, oPolygon1
   oCrosshair->Add, oPolygon2
   oCrosshair->Add, oPolygon3
   oCrosshair->Add, oPolygon4
   return, oCrosshair
end

;
; returns a color
;
function get_rgb, COLOR=color
    if not keyword_set(color) then color=0
    rgb=[0,0,0] & setcolors & tvlct, r, g, b, /GET & rgb=[r[color],g[color],b[color]]
    return, rgb
end

;
; returns a Font object
;
function get_font_obj, FONTSIZE=fontsize, TEXTTHICK=textthick
    if not keyword_set(fontsize) then fontsize=12
    if not keyword_set(textthick) then textthick=1.0
    oFont = obj_new('IDLgrFont', NAME='Hershey', SIZE=fontsize, THICK=textthick)
    return, oFont
end

;
; returns a formatted Text object
;
function get_text_obj, text, TEXTCOLOR=textcolor, FONTSIZE=fontsize, TEXTDIR=textdir, TEXTTHICK=textthick, ALIGNMENT=alignment
    if not keyword_set(textcolor) then textcolor = 0
    if not keyword_set(fontsize) then fontsize=12
    if not keyword_set(alignment) then alignment = 0.0
    if not keyword_set(textthick) then textthick = 1.
    if not keyword_set(textdir) then begin
       baseline = [1.0,0.0] & updir = [0.0,1.0]
    endif else begin
       baseline = [0.0,1.0] & updir = [-1.0,0.0]
    endelse
    oText = obj_new('IDLgrText', text, FONT=get_font_obj(FONTSIZE=fontsize, TEXTTHICK=textthick), COLOR=get_rgb(COLOR=textcolor), $
                    ALIGNMENT=alignment, BASELINE=baseline, UPDIR=updir, /ENABLE_FORMATTING)
    return, oText
end

;
; returns a palette object
;
function get_palette_obj, CT=ct, RED=r, GREEN=g, BLUE=b, _EXTRA=extra
    oPalette = OBJ_NEW('IDLgrPalette')
    if not keyword_set(ct) then ct = 0
    oPalette->LoadCT, ct
    if keyword_set(top_stretch) then oPalette->SetProperty, TOP_STRETCH=top_stretch
    if keyword_set(bottom_stretch) then oPalette->SetProperty, BOTTOM_STRETCH=bottom_stretch
    oPalette->GetProperty, RED_VALUES=r
    oPalette->GetProperty, GREEN_VALUES=g
    oPalette->GetProperty, BLUE_VALUES=b
    invert = 0
    ; fix color table 'HAZE'
    if ct eq 16 then begin
       r[0] = r[2] & r[1] = r[2]
       g[0] = g[2] & g[1] = g[2]
       b[0] = b[2] & b[1] = b[2]
       r[0] = 255
       g[0] = 255
       b[0] = 255
    endif
    ; invert table(s) for sinks
    if ct eq 65 or ct eq 55 or ct eq 20 then begin
       invert = 1
    endif
    ; cut and invert prism table
    if ct eq 6 then begin
       invert = 1
       r[0] = 255
       g[0] = 255
       b[0] = 255
       r[255] = 255
       g[255] = 255
       b[255] = 255
    endif
    ; cut blue-pastel-red table
    if ct eq 17 then begin
       r[0] = 255
       g[0] = 255
       b[0] = 255
       r[255] = 255
       g[255] = 255
       b[255] = 255
    endif
    ; cut blue-red table
    if ct eq 11 then begin
       invert = 1
       startindex = 128
       endindex = 255
       tmp = r[startindex:endindex] & r = congrid(tmp,256,/interp)
       tmp = g[startindex:endindex] & g = congrid(tmp,256,/interp)
       tmp = b[startindex:endindex] & b = congrid(tmp,256,/interp)
    endif
    if invert then begin
       tmp = reverse(r[*]) & r = tmp
       tmp = reverse(g[*]) & g = tmp
       tmp = reverse(b[*]) & b = tmp
    endif
    oPalette->SetProperty, RED_VALUES=r
    oPalette->SetProperty, GREEN_VALUES=g
    oPalette->SetProperty, BLUE_VALUES=b
    return, oPalette
end

;
; returns a colorbar object
;
function get_cbar_obj, BARDIMS=barDims, RANGE=range, PALETTE=oPalette, TITLE=title, FONTSIZE=fontsize, TICKFORMAT=tickformat, $
	 BOTTOM=bottom, RIGHT=right, LOG=log, _EXTRA=extra
    ; add a colorbar
    range = double(range)
    oCBAR = OBJ_NEW('IDLgrModel')
    oBar = OBJ_NEW('IDLgrColorbar', PALETTE=oPalette)
    oBar->SetProperty, SHOW_AXIS=0, COLOR=0, /SHOW_OUTLINE, DIMENSIONS=barDims
    oCBAR->Add, oBar
    if not keyword_set(log) then log=0
    textpos=0 & tickdir=0
    location = [0, 0]
    if keyword_set(right) then begin
      textpos=1 & tickdir=1
      location = [barDims[0], 0]
    endif
    if keyword_set(bottom) then begin
      textpos=0 & tickdir=1
    endif
    ; add the colorbar axis
    if keyword_set(bottom) then begin
      s0 = -range[0]*barDims[0]/(range[1]-range[0])
      s1 = barDims[0]/(range[1]-range[0])
      if keyword_set(log) then begin
	s0 = -alog10(range[0])*barDims[0]/alog10(range[1]/range[0])
	s1 = barDims[0]/alog10(range[1]/range[0])
      endif
      oBarAxis = OBJ_NEW('IDLgrAxis', 0, RANGE=[range[0],range[1]], LOCATION=location, /EXACT, TICKLEN=0.2*barDims[1], TICKFORMAT=tickformat, $
			  TITLE=get_text_obj(title, FONTSIZE=fontsize), XCOORD_CONV=[s0,s1], TICKDIR=tickdir, TEXTPOS=textpos, LOG=log, _EXTRA=extra)
    endif else begin
      s0 = -range[0]*barDims[1]/(range[1]-range[0])
      s1 = barDims[1]/(range[1]-range[0])
      if keyword_set(log) then begin
	s0 = -alog10(range[0])*barDims[1]/alog10(range[1]/range[0])
	s1 = barDims[1]/alog10(range[1]/range[0])
      endif
      oBarAxis = OBJ_NEW('IDLgrAxis', 1, RANGE=[range[0],range[1]], LOCATION=location, /EXACT, TICKLEN=0.2*barDims[0], TICKFORMAT=tickformat, $
			  TITLE=get_text_obj(title, FONTSIZE=fontsize), YCOORD_CONV=[s0,s1], TICKDIR=tickdir, TEXTPOS=textpos, LOG=log, _EXTRA=extra)
    endelse
    oBarAxis->GetProperty, TICKTEXT=ticktext
    ticktext->GetProperty, STRINGS=strings
    ticktext->SetProperty, FONT=get_font_obj(FONTSIZE=fontsize), STRINGS=strings
    oBarAxis->SetProperty, TICKTEXT=ticktext
    oCBAR->Add, oBarAxis
    return, oCBAR
end

;
; returns an axis object
;
function get_axis_obj, AXISDIMS=axisDims, RANGE=range, TITLE=title, FONTSIZE=fontsize, TICKFORMAT=tickformat, $
    SCALE=scale, LOC=loc, LOG=log, _EXTRA=extra

    if not keyword_set (loc) then stop
    if not keyword_set(log) then log=0

    range = double(range)
    location = [0, 0]
    ticklen = 0.01*scale

    ; take care of case where the axis is inverted
    if range[0] gt range[1] then reverse_axis = 1 else reverse_axis = 0
    if reverse_axis then range = -range

    oOut = OBJ_NEW('IDLgrModel')

    if loc eq 'bottom' or loc eq 'top' then begin

      s0 = -range[0]*axisDims/(range[1]-range[0])
      s1 = axisDims/(range[1]-range[0])
      if keyword_set(log) then begin
          s0 = -alog10(range[0])*axisDims/alog10(range[1]/range[0])
          s1 = axisDims/alog10(range[1]/range[0])
      endif

      if loc eq 'bottom' then begin
          textpos=0 & tickdir=0
          here_tickformat = tickformat
          if reverse_axis then begin
              tickfmtpos = strpos(here_tickformat,'.')
              lentickfmt = strmid(here_tickformat,tickfmtpos-1,1)
              new_lentickfmt = strcompress(string(long(lentickfmt)+1),/remove_all)
              here_tickformat = strmid(here_tickformat,0,tickfmtpos-1)+new_lentickfmt+strmid(here_tickformat,tickfmtpos)
          endif
      endif
      if loc eq 'top' then begin
          textpos=1 & tickdir=1
          here_tickformat = '(A1)'
      endif

      oAxis = OBJ_NEW('IDLgrAxis', 0, RANGE=[range[0],range[1]], LOCATION=location, /EXACT, TICKFORMAT=here_tickformat, TICKLEN=ticklen, $
          TITLE=get_text_obj(title, FONTSIZE=fontsize), XCOORD_CONV=[s0,s1], TICKDIR=tickdir, TEXTPOS=textpos, LOG=log, _EXTRA=extra)
      
    endif else begin
        
      s0 = -range[0]*axisDims/(range[1]-range[0])
      s1 = axisDims/(range[1]-range[0])
      if keyword_set(log) then begin
          s0 = -alog10(range[0])*axisDims/alog10(range[1]/range[0])
          s1 = axisDims/alog10(range[1]/range[0])
      endif

      if loc eq 'left' then begin
          textpos=0 & tickdir=0
          here_tickformat = '(A1)'
      endif
      if loc eq 'right' then begin
          textpos=1 & tickdir=1
          here_tickformat = tickformat
      endif

      oAxis = OBJ_NEW('IDLgrAxis', 1, RANGE=[range[0],range[1]], LOCATION=location, /EXACT, TICKFORMAT=here_tickformat, TICKLEN=ticklen, $
          TITLE=get_text_obj(title, FONTSIZE=fontsize), YCOORD_CONV=[s0,s1], TICKDIR=tickdir, TEXTPOS=textpos, LOG=log, _EXTRA=extra)

    endelse

    oAxis->GetProperty, TICKTEXT=ticktext
    ticktext->GetProperty, STRINGS=strings
    if reverse_axis then begin
        if here_tickformat ne '(A1)' then here_tickformat = tickformat
        strings = string(-strings,format=here_tickformat)
        range = -range
    endif
    ticktext->SetProperty, FONT=get_font_obj(FONTSIZE=fontsize), STRINGS=strings
    oAxis->SetProperty, TICKTEXT=ticktext
    oOut->Add, oAxis
    
    return, oOut
end



; =========================
; ========= MAIN ==========
; =========================
pro write_eps_wcb, image_data, imageSize, filename, RANGE=range, CTAB=ctab, TOP_STRETCH=top_stretch, BOTTOM_STRETCH=bottom_stretch, $
                   TICKFORMATMAP=tickformatmap, TITLE=title, FONTSIZE=fontsize, COLORBAR=colorbar, IMAGELOG=imagelog, $
                   TEXTFIELD=textfield, TEXTCOLOR=textcolor, TEXTPOSITION=textposition, TEXTFONTSIZE=textfontsize, TEXTDIR=textdir, $
                   TEXTTHICK=textthick, TEXTALIGN=textalign, TEXTROTATE=textrotate, $
                   PARTICLEPOSITIONS=particlepositions, PARTICLERADIUS=particleradius, PARTICLECOLOR=particlecolor, PARTICLEMASSES=particlemasses, $
                   PARTICLECT=particlect, PARTICLECOLBAR=particlecolbar, MINPARTICLEMASS=minparticlemass, MAXPARTICLEMASS=maxparticlemass, $
                   TICKFORMATPARTICLECOLBAR=tickformatparticlecolbar, PCBTITLE=pcbtitle, PARTICLELOG=particlelog, $
                   DRAWLINE=drawline, LINECOLOR=linecolor, DRAWARROW=drawarrow, ARROWCOLOR=arrowcolor, $
                   DRAWCONTOUR=drawcontour, CONTOURDATA=contourdata, CONTOURVALUES=contourvalues, CONTOURCOLOR=contourcolor, CONTOURTHICK=contourthick, $
                   DRAWVECTORFIELD=drawvectorfield, VECTORCOMPONENT1=vectorcomponent1, VECTORCOMPONENT2=vectorcomponent2, VECTORFIELDCOLOR=vectorfieldcolor, $
                   VECTORSMOOTHING=vectorsmoothing, VECTORNORM=vectornorm, VECTORLEGENDTEXT=vectorlegendtext, VECTORLEGENDCOLOR=vectorlegendcolor, $
                   VECTORLEGENDPOSITION=vectorlegendposition, ONLYCB=onlycb, $
                   DRAW_ELLIPSE_PTS=draw_ellipse_pts, ELLIPSE_COLOR=ellipse_color, ELLIPSE_THICK=ellipse_thick, $
                   DRAWAXESRANGEX=drawAxesRangeX, DRAWAXESRANGEY=drawAxesRangeY, DRAWAXESTITLEX=drawAxesTitleX, DRAWAXESTITLEY=drawAxesTitleY, $
                   DRAWAXESTICKFORMATX=drawAxesTickFormatX, DRAWAXESTICKFORMATY=drawAxesTickFormatY, DRAWAXESTICKCOLOR=drawAxesTickColor, $
                   BIGFONT=bigfont

  if not keyword_set(colorbar) then colorbar = 0
  if not keyword_set(particlecolbar) then particlecolbar = 0

  ; some switches
  bfx = 0.14

  ; set dimensions
  barCanvasDims = [0.25*imageSize[0], imageSize[1]]
  if keyword_set(bigfont) then barCanvasDims = [0.30*imageSize[0], imageSize[1]]

  if colorbar ne 0 then begin
     if particlecolbar ne 0 then begin
        windowSize = [imageSize[0]+2*barCanvasDims[0],imageSize[1]]
        imagetranslate = [barCanvasDims[0], 0.]
     endif else begin
        if colorbar eq 66 then begin
           barCanvasDims = [imageSize[0], 0.187*imageSize[1]]
           windowSize = [imageSize[0],imageSize[1]+barCanvasDims[1]]
           imagetranslate = [0., barCanvasDims[1]]
        endif else begin
           windowSize = [imageSize[0]+barCanvasDims[0],imageSize[1]]
           imagetranslate = [barCanvasDims[0], 0.]
        endelse
     endelse
  endif else begin
     if particlecolbar ne 0 then begin
        if particlecolbar eq 66 then begin
           barCanvasDims = [imageSize[0], 0.187*imageSize[1]]
           windowSize = [imageSize[0],imageSize[1]+barCanvasDims[1]]
           imagetranslate = [0., barCanvasDims[1]]
        endif else begin
           windowSize = [imageSize[0]+barCanvasDims[0],imageSize[1]]
           imagetranslate = [0., 0.]
        endelse
     endif else begin
        windowSize = [imageSize[0],imageSize[1]]
        imagetranslate = [0., 0.]
     endelse
  endelse
  
  ; make room for axes
  axesSpaceMIN = [0,0]
  axesSpaceMAX = [0,0]
  if keyword_set(drawAxesRangeX) then begin
     if colorbar ne 0 then begin
        axesSpaceMIN = [0.00*imageSize[0],0.13*imageSize[0]]
     endif else begin
        axesSpaceMIN = [0.07*imageSize[0],0.13*imageSize[0]]
     endelse
     if particlecolbar ne 0 then begin
        axesSpaceMAX = [0.13*imageSize[0],0.05*imageSize[0]]
     endif else begin
        axesSpaceMAX = [0.13*imageSize[0],0.05*imageSize[0]]
     endelse
     windowSize[0] += axesSpaceMIN[0]+axesSpaceMAX[0]
     windowSize[1] += axesSpaceMIN[1]+axesSpaceMAX[1]
     imagetranslate += axesSpaceMIN
  endif
  
  if keyword_set(onlycb) then windowSize = barCanvasDims
  
	oWindow = OBJ_NEW('IDLgrBuffer', DIMENSIONS=windowSize)
	oViewgroup = OBJ_NEW('IDLgrViewgroup')
	oWindow->SetProperty, GRAPHICS_TREE = oViewgroup
	
	; Add an IDL container to the viewgroup to hold text objects.
	oContainer = OBJ_NEW('IDL_CONTAINER')
	oViewgroup->Add, oContainer

	oView = OBJ_NEW('IDLgrView', VIEWPLANE_RECT=[0., 0., windowSize[0], windowSize[1]])
	oViewgroup->Add, oView
	oModel1 = OBJ_NEW('IDLgrModel')
	oModel2 = OBJ_NEW('IDLgrModel')
	oModel3 = OBJ_NEW('IDLgrModel')
	oView->Add, oModel1
	oView->Add, oModel2
	oView->Add, oModel3

    ; Add axes around image
    if keyword_set(drawAxesRangeX) or keyword_set(drawAxesRangeY) then begin
       oModelAxes = OBJ_NEW('IDLgrModel')
       oView->Add, oModelAxes
       if not keyword_set(drawAxesTickColor) then drawAxesTickColor = '255'
    endif
    if keyword_set(drawAxesRangeX) then begin
        if not keyword_set(drawAxesTickFormatX) then drawAxesTickFormatX = '(F3.1)'
        if not keyword_set(drawAxesTitleX) then drawAxesTitleX = ''
        ; bottom axis
        oAxis = get_axis_obj(AXISDIMS=imageSize[0], RANGE=drawAxesRangeX, TITLE=drawAxesTitleX, $
            SCALE=imageSize[0], FONTSIZE=fontsize, TICKFORMAT=drawAxesTickFormatX, LOC='bottom', COLOR=drawAxesTickColor, /USE_TEXT_COLOR)
        oAxis->Translate, imagetranslate[0], imagetranslate[1], 0
        oModelAxes->Add, oAxis
        ; top axis
        oAxis = get_axis_obj(AXISDIMS=imageSize[0], RANGE=drawAxesRangeX, TITLE='', $
            SCALE=imageSize[0], FONTSIZE=fontsize, TICKFORMAT=drawAxesTickFormatX, LOC='top', COLOR=drawAxesTickColor, /USE_TEXT_COLOR)
        oAxis->Translate, imagetranslate[0], imagetranslate[1]+imageSize[1], 0
        oModelAxes->Add, oAxis
    endif
    if keyword_set(drawAxesRangeY) then begin
        if not keyword_set(drawAxesTickFormatY) then drawAxesTickFormatY = '(A1)'
        if not keyword_set(drawAxesTitleY) then drawAxesTitleY = ''
        ; left axis
        oAxis = get_axis_obj(AXISDIMS=imageSize[1], AXESLOC=imagetranslate, RANGE=drawAxesRangeY, TITLE='', $
            SCALE=imageSize[0], FONTSIZE=fontsize, TICKFORMAT=drawAxesTickFormatY, LOC='left', COLOR=drawAxesTickColor, /USE_TEXT_COLOR)
        oAxis->Translate, imagetranslate[0], imagetranslate[1], 0
        oModelAxes->Add, oAxis
        ; right axis (use bottom axis and rotate in position)
        oAxis = get_axis_obj(AXISDIMS=imageSize[1], RANGE=drawAxesRangeY, TITLE=drawAxesTitleY, $
            SCALE=imageSize[0], FONTSIZE=fontsize, TICKFORMAT=drawAxesTickFormatY, LOC='bottom', COLOR=drawAxesTickColor, /USE_TEXT_COLOR)
        oAxis->Rotate, [0,0,1], 90
        oAxis->Translate, imagetranslate[0]+imageSize[0], imagetranslate[1], 0
        oModelAxes->Add, oAxis
    endif

	; create palette
	oPaletteCBar = get_palette_obj(CT=ctab, TOP_STRETCH=top_stretch, BOTTOM_STRETCH=bottom_stretch)
        
        ; add colorbar for the image
        if colorbar ne 0 then begin
           if colorbar eq 66 then begin
              if not keyword_set(tickformatmap) then tickformatmap = '(-F4.1)'
              barDims = [0.88*barCanvasDims[0], 0.18*barCanvasDims[1]]
              oModel2->Add, get_cbar_obj(BARDIMS=barDims, RANGE=range, PALETTE=oPaletteCBar, TITLE=title, FONTSIZE=fontsize, TICKFORMAT=tickformatmap, /BOTTOM, LOG=imagelog)
              oModel2->Translate, 0.06*barCanvasDims[0]+axesSpaceMAX[0], (1.-2.0*bfx)*barCanvasDims[1], 0
           endif else begin
              if not keyword_set(tickformatmap) then tickformatmap = '(F4.1)'
              barDims = [bfx*barCanvasDims[0], 0.96*barCanvasDims[1]]
              oModel2->Add, get_cbar_obj(BARDIMS=barDims, RANGE=range, PALETTE=oPaletteCBar, TITLE=title, FONTSIZE=fontsize, TICKFORMAT=tickformatmap, LOG=imagelog)
              oModel2->Translate, (1.-1.2*bfx)*barCanvasDims[0], 0.02*barCanvasDims[1]+axesSpaceMIN[1], 0
           endelse
        endif ;colorbar
        
        ; add colorbar for the particles
        if particlecolbar ne 0 then begin
            oPalettePCBar = get_palette_obj(CT=particlect)
            if not keyword_set(pcbtitle) then pcbtitle = textoidl('log_{10} ( !8M!X_{sink} [M!D!9n!3!N] )')
            if not keyword_set(tickformatparticlecolbar) then tickformatparticlecolbar = '(F-4.1)'
            pcb_range = [minparticlemass,maxparticlemass]
            if particlecolbar eq 66 then begin
              barDims = [0.88*barCanvasDims[0], 0.18*barCanvasDims[1]]
              oModel3->Add, get_cbar_obj(BARDIMS=barDims, RANGE=pcb_range, PALETTE=oPalettePCBar, TITLE=pcbtitle, $
                  FONTSIZE=fontsize, TICKFORMAT=tickformatparticlecolbar, /BOTTOM, LOG=particlelog)
              oModel3->Translate, 0.06*barCanvasDims[0], (1.-2.0*bfx)*barCanvasDims[1], 0
            endif else begin
              barDims = [bfx*barCanvasDims[0], 0.96*barCanvasDims[1]]
              oModel3->Add, get_cbar_obj(BARDIMS=barDims, RANGE=pcb_range, PALETTE=oPalettePCBar, TITLE=pcbtitle, $
                  FONTSIZE=fontsize, TICKFORMAT=tickformatparticlecolbar, /RIGHT, LOG=particlelog)
              if keyword_set(onlycb) then begin
                 oModel3->Translate, 0.2*bfx*barCanvasDims[0], 0.02*barCanvasDims[1], 0
              endif else begin
                 oModel3->Translate, 0.2*bfx*barCanvasDims[0]+imageSize[0]+imagetranslate[0], 0.02*barCanvasDims[1]+imagetranslate[1], 0
              endelse
            endelse
        endif ;particlecolbar

        if keyword_set(onlycb) then begin
           oClipboard = OBJ_NEW('IDLgrClipboard', DIMENSIONS=windowSize, GRAPHICS_TREE=oViewgroup)
           oClipboard->Draw, FILENAME=filename, /POSTSCRIPT, /VECTOR
           obj_destroy, oClipboard
           HEAP_GC
           return
        endif
        
        ; add image
        image = congrid(image_data,imageSize[0],imageSize[1],/center,cubic=-0.5)
        if keyword_set(imagelog) then begin
           image = bytscl(alog10(image),min=alog10(range[0]),max=alog10(range[1]))
        endif else begin
           image = bytscl(image,min=range[0],max=range[1])
        endelse
        oImage = obj_new('IDLgrImage', image, PALETTE=oPaletteCBar)
        oModel1->Add, oImage
        oModel1->Translate, imagetranslate[0], imagetranslate[1], 0

        ; add textfields
        if keyword_set(textfield) then begin
           if not keyword_set(textdir) then textdir = lonarr(n_elements(textfield))
           if not keyword_set(textthick) then begin & textthick = lonarr(n_elements(textfield)) & textthick[*] = 2.0 & endif
           if not keyword_set(textalign) then begin & textalign = lonarr(n_elements(textfield)) & textalign[*] = 0.0 & endif
           for i=0, n_elements(textfield)-1 do begin
              oModelBox1 = obj_new('IDLgrModel')
              oModelBox1->Add, get_text_obj(textfield[i], TEXTCOLOR=textcolor[i], TEXTDIR=textdir[i], FONTSIZE=textfontsize[i], TEXTTHICK=textthick[i], ALIGNMENT=textalign[i])
              oModelBox2 = obj_new('IDLgrModel')
              oModelBox2->Add, get_text_obj(textfield[i], TEXTCOLOR=0, TEXTDIR=textdir[i], FONTSIZE=textfontsize[i], TEXTTHICK=textthick[i], ALIGNMENT=textalign[i])
              ; rotate
              if keyword_set(textrotate) then begin
                  axis = [0,0,1] & angle = textrotate[i]
                  oModelBox1->Rotate, axis, angle
                  oModelBox2->Rotate, axis, angle
              endif
              ; translate
              textpos=[0.0, 0.0]
              if keyword_set(textposition) then textpos=textposition[*,i]
              oModelBox1->Translate, textpos[0]*imageSize[0], textpos[1]*imageSize[1], 0
              oModelBox2->Translate, textpos[0]*imageSize[0]+0.0015*imageSize[0], textpos[1]*imageSize[1]-0.0015*imageSize[1], 0
              oModel1->Add, oModelBox1
              oModel1->Add, oModelBox2
           endfor
        endif

        ; add line(s)
        if keyword_set(drawline) then begin
           ; usage, for example for two lines: drawline = [[xs1,ys1,xe1,ye1], [xs2,ys2,xe2,ye2]]
           sl = size(drawline)
           ; check rank
           if (sl[0] eq 1) then begin
               ; artificially add another dimension for the code to work below
               dl_loc = [[drawline],[drawline]]
               nl = 1
           endif else begin
               dl_loc = drawline
               nl = sl[2]
           endelse
           if not keyword_set(linecolor) then linecolor = lonarr(nl)
           for il = 0, nl-1 do begin
               xs = dl_loc[0,il]
               ys = dl_loc[1,il]
               xe = dl_loc[2,il]
               ye = dl_loc[3,il]
               oLine = obj_new('IDLgrPlot', [xs*imageSize[0], xe*imageSize[0]], [ys*imageSize[1], ye*imageSize[1]], COLOR=linecolor[il], THICK=2.0)
               oModelBox = obj_new('IDLgrModel')
               oModelBox->Add, oLine
               oModel1->Add, oModelBox
           endfor
        endif

        ; add arrow(s) (vector)
        if keyword_set(drawarrow) then begin
            ; usage, for example for two arrows: drawarrow = [[xs1,ys1,xe1,ye1], [xs2,ys2,xe2,ye2]]
            sa = size(drawarrow)
            ; check rank
            if (sa[0] eq 1) then begin
                ; artificially add another dimension for the code to work below
                da_loc = [[drawarrow],[drawarrow]]
                na = 1
            endif else begin
                da_loc = drawarrow
                na = sa[2]
            endelse
            if not keyword_set(arrowcolor) then arrowcolor = lonarr(na)
            internal_vectornorm = ImageSize
            for ia = 0, na-1 do begin
                ax = (da_loc[2,ia]-da_loc[0,ia])*internal_vectornorm[0]
                ay = (da_loc[3,ia]-da_loc[1,ia])*internal_vectornorm[1]
                oVector = obj_new('IDLgrModel')
                oVector->Add, get_vector_obj(VECTOR=[ax,ay], COLOR=get_rgb(COLOR=arrowcolor[ia]))
                oVector->Translate, da_loc[0,ia]*internal_vectornorm[0]+ax/2, da_loc[1,ia]*internal_vectornorm[1]+ay/2, 0
                oModel1->Add, oVector
            endfor
        endif

        ; add vector field legend
        if keyword_set(drawvectorfield) then begin
           internal_vectornorm = 0.05*ImageSize[0]
           if not keyword_set(vectorfieldcolor) then vectorfieldcolor = 2
           xsq = [0.,1.,1.,0.]*internal_vectornorm*1.25 & ysq = [0.,0.,1.,1.]*internal_vectornorm*0.4
           oVectorLegend = obj_new('IDLgrModel')
           ;oSquare = obj_new('IDLgrPolygon', xsq, ysq, COLOR=get_rgb(COLOR=vectorlegendcolor))
           oVector = obj_new('IDLgrModel')
           oVector->Add, get_vector_obj(VECTOR=[internal_vectornorm,0.], COLOR=get_rgb(COLOR=vectorfieldcolor))
           oVector->Translate, 0.5*(xsq[1]-xsq[0]), 0.5*(ysq[3]-ysq[0]), 0
           oText1 = obj_new('IDLgrModel')
           oText2 = obj_new('IDLgrModel')
           if not keyword_set(vectorlegendtext) then vectorlegendtext = string(vectornorm)
           oText1->Add, get_text_obj(vectorlegendtext, TEXTCOLOR=255, FONTSIZE=fontsize, TEXTTHICK=2.0, ALIGNMENT=0.0)
           oText2->Add, get_text_obj(vectorlegendtext, TEXTCOLOR=0, FONTSIZE=fontsize, TEXTTHICK=2.0, ALIGNMENT=0.0)
           oText1->Translate, 0.0*(xsq[1]-xsq[0]), -2.0*(ysq[3]-ysq[0]), 0
           oText2->Translate, 0.0*(xsq[1]-xsq[0])+0.0015*imageSize[0], -2.0*(ysq[3]-ysq[0])-0.0015*imageSize[1], 0
           oVectorLegend->Add, oText1
           oVectorLegend->Add, oText2
           oVectorLegend->Add, oVector
           ;oVectorLegend->Add, oSquare
           if not keyword_set(vectorlegendposition) then vectorlegendposition = [0.02, 0.07]
           oVectorLegend->Translate, vectorlegendposition[0]*ImageSize[0], vectorlegendposition[1]*ImageSize[1], 0
           oModel1->Add, oVectorLegend
        endif

        ; add ellipse; draw_ellipse_pts[n_ellipses, 2, n_points]
        if keyword_set(draw_ellipse_pts) then begin
            if not keyword_set(ellipse_color) then ellipse_color = '255'
            if not keyword_set(ellipse_thick) then ellipse_thick = 0.25
            for i = 0, n_elements(draw_ellipse_pts[*,0,0])-1 do begin
                oEllipse = obj_new('IDLgrModel')
                oEllipse->Add, get_ellipse_obj(PTS=reform(draw_ellipse_pts[i,*,*]), SIZE=ellipse_thick, COLOR=get_rgb(COLOR=ellipse_color))
                oModel1->Add, oEllipse
            endfor
        endif

        ; add particles
        if keyword_set(particlepositions) then begin
           numpart = n_elements(particlepositions[*,0])
           oParticles = objarr(numpart)
           oParticlesShadow = objarr(numpart)
           index = reverse(sort(particlemasses))
           ppos_sorted = particlepositions[index,*]
           color=[0,0,0]
           if keyword_set(particlecolor) then begin
              setcolors & TVLCT, r, g, b, /GET & color=[r[particlecolor],g[particlecolor],b[particlecolor]]
           endif
           if keyword_set(particlemasses) and not keyword_set(particlecolor) then begin
              if not keyword_set(particlect) then particlect = 0
              oPalettePCBar = get_palette_obj(CT=particlect,RED=r,GREEN=g,BLUE=b)
              if not keyword_set(particlelog) then begin
                 ; scale the ordered particle masses array
                 particlemasses_scaled = bytscl(particlemasses[index],min=minparticlemass,max=maxparticlemass)
              endif else begin
                 particlemasses_scaled = bytscl(alog10(particlemasses[index]),min=alog10(minparticlemass),max=alog10(maxparticlemass))
              endelse
           endif
           for ip=0, numpart-1 do begin
              if keyword_set(particlemasses) and not keyword_set(particlecolor) then begin
                 index = particlemasses_scaled[ip]
                 color=[r[index],g[index],b[index]]
              endif
              oParticles[ip] = obj_new('IDLgrModel')
              oParticles[ip]->Add, get_filled_circle_obj(SIZE=1.0*particleradius, COLOR=color)
              oParticles[ip]->Add, get_crosshair_obj(SIZE=5.0*particleradius, COLOR=color)
              oParticlesShadow[ip] = obj_new('IDLgrModel')
              oParticlesShadow[ip]->Add, get_filled_circle_obj(SIZE=1.0*particleradius, COLOR='0')
              oParticlesShadow[ip]->Add, get_crosshair_obj(SIZE=5.0*particleradius, COLOR='0')
              ;oParticles[ip]->Add, get_cross_obj(SIZE=0.95*particleradius, COLOR=color)
              ;oParticles[ip]->Add, get_cross_obj(SIZE=0.96*particleradius, COLOR=255)
              ;oParticles[ip]->Add, get_circle_obj(SIZE=1.0*particleradius, COLOR=color)
              ;oParticles[ip]->Add, get_circle_obj(SIZE=1.00*particleradius, COLOR=255)
              oParticles[ip]->Translate, ppos_sorted[ip,0]*imageSize[0], ppos_sorted[ip,1]*imageSize[1], 0
              oParticlesShadow[ip]->Translate, ppos_sorted[ip,0]*imageSize[0]+0.001*imageSize[0], ppos_sorted[ip,1]*imageSize[1]-0.001*imageSize[1], 0
           endfor
           oModel1->Add, oParticles
           oModel1->Add, oParticlesShadow
        endif ; particles

        ; add contour
        if keyword_set(drawcontour) then begin
           if not keyword_set(contourcolor) then contourcolor = lonarr(n_elements(contourvalues))
           color = bytarr(3,n_elements(contourcolor))
           for i=0, n_elements(contourcolor)-1 do color[*,i]=get_rgb(COLOR=contourcolor[i])
           oContour = OBJ_NEW('IDLgrContour', contourdata, C_VALUE=contourvalues, /PLANAR, GEOMZ=0, $
                              C_COLOR=color, C_LINESTYLE=0, C_THICK=contourthick)
           oModel1->Add, oContour
        endif

        ; add vectorfield
        if keyword_set(drawvectorfield) then begin
           dimvec = size(vectorcomponent1,/dim)
           dimvec2 = size(vectorcomponent2,/dim)
	   if ((dimvec[0] ne dimvec2[0]) or (dimvec[1] ne dimvec2[1])) then begin & print, 'dimensions for vector component fields must be equal.' & stop & endif
	   if not keyword_set(vectorsmoothing) then vectorsmoothing = 16
	   if not keyword_set(vectornorm) then vectornorm = 1.
	   vx = rebin(vectorcomponent1, dimvec / vectorsmoothing) / vectornorm * internal_vectornorm
	   vy = rebin(vectorcomponent2, dimvec / vectorsmoothing) / vectornorm * internal_vectornorm
	   dimvec = size(vx,/dim)
           oVectors = objarr(dimvec[0],dimvec[1])
	   for j=0, dimvec[1]-1 do begin
	     for i=0, dimvec[0]-1 do begin
	      oVectors[i,j] = obj_new('IDLgrModel')
	      oVectors[i,j]->Add, get_vector_obj(VECTOR=[vx[i,j],vy[i,j]], COLOR=get_rgb(COLOR=vectorfieldcolor))
	      oVectors[i,j]->Translate, (i+0.5)/dimvec[0]*imageSize[0], (j+0.5)/dimvec[1]*imageSize[1], 0
	     endfor
	   endfor
	   oModel1->Add, oVectors
        endif

	; write eps
	oClipboard = OBJ_NEW('IDLgrClipboard', DIMENSIONS=windowSize, GRAPHICS_TREE=oViewgroup)
	oClipboard->Draw, FILENAME=filename, /POSTSCRIPT, /VECTOR
	obj_destroy, oClipboard
        HEAP_GC
	return

end