PRO dustem_plot_legend,ps=ps,st_model,_Extra=extra,help=help,top=top,left=left,right=right,bottom=bottom,middle=middle,xr=xr,yr=yr,xlog=xlog,ylog=ylog,mergePAH=mergePAH,polar=polar,nototal=nototal,total=total,graintype=graintype,letter=letter

letter = 0

if not keyword_set(graintype) then graintype = st_model.grain.grains.grain_type

ntype = n_elements(graintype)
colors=dustem_grains_colors(ntype,ls,ps=ps)

if keyword_set(ps) then colors(0) = 0

if keyword_set(polar) then aligned = st_model.align.grains.aligned

if keyword_set(mergePAH) then begin
	for i = 0, ntype-1 do begin
		if graintype(i) eq 'PAH0' then begin
			; Rename PAH0
			graintype(i) = 'PAH'
			; Remove PAH1 and shift grain types, linestyles and colors
			graintype(i+1:ntype-2) = graintype(i+2:ntype-1)
			aligned(i+1:ntype-2) = aligned(i+2:ntype-1)
			ls(i+2:ntype-1) = ls(i+3:ntype)
			colors(i+2:ntype-1) = colors(i+3:ntype)
			ntype -= 1
		endif
	endfor
endif

if keyword_set(letter) then begin
	if keyword_set(xlog) then x=xr(0)*(xr(1)/xr(0))^0.05 else x=xr(0)+(xr(1)-xr(0))*0.05
	if keyword_set(ylog) then y=yr(1)/(yr(1)/yr(0))^0.09 else y=yr(1)-(yr(1)-yr(0))*0.09
	print,'XXXXXXXX ',letter,'  ',x,y
	xyouts,x,y,'('+letter+')'
endif

if not keyword_set(xlog) then begin
	dx = (xr(1)-xr(0))/10.
	x = xr(0) + indgen(ntype)*0.
	if keyword_set(right) then x += 7.3*dx
	if keyword_set(middle) then x += dx*3
	x1 = x + dx*0.5
	x2 = x + dx*1.25
	x3 = x + dx*1.5
endif else begin
	dx = (xr(1)/xr(0))^0.1
	x = xr(0) + indgen(ntype)*0.
	if keyword_set(right) then x *= dx^7.3
	if keyword_set(middle) then x *= dx^3
	x1 = x * dx^0.5
	x2 = x * dx^1.25
	x3 = x * dx^1.5
endelse

if not keyword_set(ylog) then begin
	dy = (yr(1)-yr(0))/10.
	y0 = yr(1) - dy*1.2
	if keyword_set(bottom) then y0 -= 9*dy
	y = y0 - dy *(indgen(ntype)+1)
	if not keyword_set(total) or keyword_set(nototal) then y +=  dy
endif else begin
	dy = (yr(1)/yr(0))^(1./10.)
	y0 = yr(1) / dy^1.2
	if keyword_set(bottom) then y0 /= dy^9.
	y = y0 / dy^(indgen(ntype)+1) 
	if not keyword_set(total) or keyword_set(nototal) then y *= dy
endelse

;Do not write legend if only one polarized grain
count = 0
if keyword_set(polar) then begin
	j=where(aligned ne '',count)
	if count eq 1 then begin
		; Give same style like for total
		total = 0
		nototal = 1
		y(j)=y0
	endif
endif

if keyword_set(total) or (not keyword_set(nototal) and (keyword_set(polar) and count gt 1 or not keyword_set(polar))) then begin
	; Total
	oplot,[x1(0),x2(0)],y0*[1,1],color=colors(0)
	xyouts,x3(0),y0,'Total',color=colors(0)
endif

; Each grain type
FOR i=0,ntype-1 DO begin
	if keyword_set(polar) then if aligned(i) eq '' then begin
	;	; Shift labels upward
	;	if i lt ntype-1 then y(i+1:ntype-1) = y(i:ntype-2)
		continue
	endif
	oplot,[x1(i),x2(i)],y(i)*[1,1],line=ls(i+1),color=colors(i+1)
	xyouts,x3(i),y(i),graintype(i),color=colors(i+1)
ENDFOR


END