PRO PlotSolarDisk, zz, ut=ut,   $
        radius  = radius    , $
        diameter= diameter  , $
        breakval= breakval  , $
        xysize  = xysize    , $
        title   = title     , $
        upto    = upto      , $
        rotate  = rotate    , $
        earth   = earth     , $
        user_position=user_position , $
        user_align  = user_align    , $
        user_string = user_string   , $
        _extra  = _extra
    @compile_opt.pro        ; On error, return to caller


InitVar, diameter, 0.8*xysize[0]
InitVar, radius  , !sun.rau
InitVar, earth   , /key


center = [0.54*!d.x_size,0.48*!d.y_size]    ; Center of disk
rdisk  = diameter/2.0                       ; Disk radius in pixels
corner = center-rdisk                       ; Lower left pixel of square enclosing disk

sun = !sun.rau/radius                       ; Radius of sun in units of disk radius


; Set up 3D transformation for adding grids

d0 = center-rdisk                           ; Square around disk in device coordinates
d1 = center+rdisk
n0 = convert_coord(d0, /device, /to_norm)   ; Square around disk in normal coordinates
n1 = convert_coord(d1, /device, /to_norm)
n0[2] = n0[1]
n1[2] = n1[1]
d0 = -1*[1,1,1]                             ; Square around disk in data coordinates (-1 to +1)
d1 =  1*[1,1,1]

origin = [0,0,0]                            ; Vectors in data coordinates
xunit  = [1,0,0]
yunit  = [0,1,0]
zunit  = [0,0,1]

; Set up transformation matrix from heliographic to x-y-z coordinates

nrot =  n_elements(rotate)/3
rot = reform(rotate,3,nrot)

xrot = rot[0,*]
yrot = rot[1,*]
zrot = rot[2,*]

zero = make_array(type=IsType(rot), dim=[1,nrot])

rot = [zero,zero,zrot,zero,yrot,zero,xrot,zero,zero]
IF nrot GT 1 THEN rot = -reverse(rot,2)
rot = reform(rot,3,3*nrot, /overwrite)

zero = where(rot[0,*] NE 0 OR rot[1,*] NE 0 OR rot[2,*] NE 0)
IF zero[0] NE -1 THEN rot = rot[*,zero]


zpole = cvt3d(rot=rot, vector=[0,0,1])      ; Solar north in rect helio coord
xpole = vectorproduct(zunit,zpole, /unit)   ; East limb on solar equator in x-y plane
ypole = vectorproduct(zpole,xpole, /unit)   ; Pointing out of the screen


lzero = cvt3d(rot=rot, vector=[1,0,0])      ; Zero helio long in rect helio coord

plot3darc, origin, zpole,  lzero, 0, !pi, circle=meridian0
plot3darc, origin, zpole, -lzero, 0, !pi, circle=meridian180

pearth = big_eph( ut,   $
    body = jpl_body(/earth,/string) , $
    /to_heliographic    , $
    /precess            , $
    /onebody            , $
    /silent             )

pearth = pearth/sqrt(total(pearth*pearth))

pearth = cvt3d(rot=rot, vector=pearth)

; Set up a coordinate system with x and y-axis in the plane of the display
; x = horizontal to the right, y = vertical up; z = out from the display towards the eye

setup3d, d0, d1, n0, n1, /xyplane


erase


CASE IsType(breakval, /defined) OF
0: zz = bytscl(zz, /nan)
1: zz = GetColors(zz, breakval, /open, /badbackground, /legend, _extra=_extra)
ENDCASE



tv, zz, corner[0], corner[1]


; Plot circles around disk, and around solar disk

plot3darc, origin, xunit, yunit, 0, 2*!pi
plot3darc, origin, sun*xunit, sun*yunit, 0, 2*!pi, thick=1.5


; Plot semicircle along solar equator on disk and on sun
; Mark disk center with symbol

plot3darc, origin,     xpole,     ypole, 0, !pi
plot3darc, origin, sun*xpole, sun*ypole, 0, !pi, thick=2
plots, /t3d, origin, psym=1

; Mark visible solar pole

i = zpole[2] gt 0
plots, /t3d, (2*i-1)*zpole, psym=1

; Plot semicircle from pole to pole along heliographic
; longitude 0 (solid) or 180 (dashed) whichever is on the visible disk

i = where(meridian0[2,*] ge 0)
IF i[0] NE -1 THEN plots, meridian0[*,i], /t3d

i = where(meridian180[2,*] ge 0)
IF i[0] NE -1 THEN plots, meridian180[*,i], /t3d, linestyle=2

; Plot Earth symbol at location Earth

IF NOT earth AND pearth[2] GE 0 THEN BEGIN
    cc = gridgen(33, range=[0,2*!pi])
    usersym, [cos(cc),-1,0,0,0], [sin(cc),0,0,1,-1]     ; Circle with cross inside
    plots, pearth, psym=8, thick=0.5, symsize=charsize;, color=F1
ENDIF

; Label UT time

xyouts, 0.99, 0.95, /norm, align=1, _extra=_extra,  $
    TimeGet(ut, /ymd, upto=upto)+' UT'

; Label title

IF IsType(title, /defined) NE 0 THEN    $
    xyouts, 0.05, 0.01, /norm, _extra=_extra, title

; Label radius of disk

xyouts, 0.05, 0.95, /norm, 'R!Ds!N='+string(format='(F5.1)',radius/!sun.rau), _extra=_extra

IF earth THEN BEGIN
    ; Label directions

    plot3dtext,  yunit, 'N', labeloffset=[0, 0.03], _extra=_extra
    plot3dtext, -yunit, 'S', labeloffset=[0,-0.09], _extra=_extra
    plot3dtext,  xunit, 'W', labeloffset=[ 0.05,0], _extra=_extra
    plot3dtext, -xunit, 'E', labeloffset=[-0.05,0], _extra=_extra

ENDIF

PlotUserstring, user_string, user_position, $
    align   = user_align    , $
    charsize= charsize      , $
    _extra=_extra

RETURN  &  END