PRO vu_planarcut, hdr, ff, ihdr=ihdr, ut0=ut0, $
radius = radius , $
zoffset = zoffset , $
type = type , $
printer = printer , $
silent = silent , $
module = module , $
body = body , $
abg_euler = abg_euler , $
degrees = degrees , $
edge = edge , $
circle = circle , $
_extra = _extra
@compile_opt.pro ; On error, return to caller
@vu_fnc_include.pro ;
InitVar, silent , 0
InitVar, degrees , /key
InitVar, printer , /key
InitVar, fill , /key
InitVar, earth , /key
InitVar, circle , /key
ihdr_set = IsType(ihdr, /defined)
utt = vu_getdata(hdr, ff, ut0=ut0, ihdr=ihdr, _extra=_extra, count=count, silent=silent+1)
IF count EQ 0 THEN RETURN
; utt determines the data matrix and is used to set the relation between
; the coordinate systems.
utt = Carrington(utt, /get_time)
say, TimeGet(utt,/ymd,upto=TimeUnit(/minute))+strcompress(ihdr), threshold=0, silent=silent
CASE ihdr_set OF
0: hdri = vu_gettime(hdr, ff, ut=utt, ff_ut=ffi) ; Interpolate at time utt
1: hdri = vu_gettime(hdr, ff, ut=vu_get(hdr[ihdr], /uttime), ff_ut=ffi, /fixgrid); ihdr set explicitly
ENDCASE
IF IsType(edge,/defined) THEN BEGIN
;iedge = round( (edge-hdri.distance_start)/hdri.distance_step )
;ffi[*,*,iedge+1:*,*] = BadValue(ffi)
iedge = round( (edge-vu_get(hdri,/start_distance))/vu_get(hdri,/distance_step) )+1
IF 0 LE iedge AND iedge LT vu_get(hdri,/nrad) THEN ffi[*,*,iedge:*,*] = BadValue(ffi)
ENDIF
; Set up three unit vectors along x,y,z axes in the output frame
; (with z the normal to the plane, and x,y axes in the plane)
; The default output frame is the ecliptic seen from the north
; Then rotate these to the heliographic coordinate frame
; (z to heliographic north; x to heliographic longitude zero)
runit = [[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]]
; If abg_euler exists then apply that rotation first
IF IsType(abg_euler,/string) THEN BEGIN
; This is a fudge to allow making a planar cut through the orbital
; plane of a body whose orbital elements are stored in
; $EPHEM/orbits/orbital_elements.txt (mostly comets, I think).
tmp_body = abg_euler
tmp = strpos(tmp_body,'&')
IF tmp NE -1 THEN tmp_body = strmid(tmp_body,tmp)
tmp = filepath(root=getenv('EPHEM'),subdir='orbits','orbital_elements.txt')
IF flt_read(tmp,mask=strjoin(replicate('1',20)),data,crumbs=crumbs,/double) THEN BEGIN
crumbs = strtrim(reform(crumbs[0,*])) ; Body name
tmp = (where( strlowcase(tmp_body) EQ strlowcase(crumbs) ))[0]
IF tmp NE -1 THEN BEGIN
boost, body, abg_euler
print, data[5:6,tmp] ; Asc node and inclination
abg_euler = [-90,-data[6,tmp],90-data[5,tmp]]/ToDegrees(degrees=degrees)
ENDIF
ENDIF
IF IsType(abg_euler,/string) THEN destroyvar, abg_euler
ENDIF
IF IsType(abg_euler,/defined) THEN $
runit = EulerRotate(abg_euler, runit, degrees=degrees, /rectangular)
; Rotate to heliographic coordinates
runit = CvSky(utt, runit, _extra=_extra, /to_heliographic, $
euler_angles=euler_angles, /rectangular, degrees=degrees, euler_info=euler_info, silent=silent)
; If abg_euler exists combine it with the rotation to heliographic coordinates
IF IsType(abg_euler,/defined) THEN $
euler_angles = CombineRotations(abg_euler, euler_angles, degrees=degrees)
type = vu_type_insitu('planarcut', _extra=_extra, hdr=hdri, type=type)
FOR itype=0,n_elements(type)-1 DO BEGIN
id = (where( type[itype].data ))[0]
; Extract fields from 'type' structure
ctable = type[itype].ctable
Fmin = type[itype].Fmin
Fmax = type[itype].Fmax
display = type[itype].display
label = type[itype].label
format = type[itype].format
charsize= type[itype].charsize
xysize = type[itype].xysize
breakval= type[itype].breakval
i = where(finite(breakval), nbreak)
IF nbreak GT 0 THEN BEGIN
breakval = type[itype].breakval[i]
IF format EQ '' THEN breakval = round(breakval)
ENDIF
; Done with 'type' structure
;===========================
InitVar, radius , vu_get(hdr,/stop_distance)
InitVar, zoffset, 0.0
diameter = round(0.8*xysize[0]) ; Diameter of cross section in pixels
diameter = (diameter-1)/2*2+1
; Rectangular 2D-vectors in plane of cross section
rr = gridgen(diameter*[1,1], range=[-1.0,1.0,-1.0,1.0])
; Rectangular 3D-vectors in coordinate frame of data array
rr = (zoffset*runit[*,2])#replicate(1,diameter*diameter)+radius*runit[*,0:1]#rr
IF circle THEN outside = where(total(rr*rr,1) GT radius*radius)
rr = cv_coord(from_rect=rr, /to_sphere) ; Convert to spherical coordinates
rr[0,*] = AngleRange(rr[0,*]) ; Keep inside [0,2*!pi]
xci = vu_get(hdri, /array_range )
rri = vu_get(hdri, /start_distance)
dri = vu_get(hdri, /distance_step )
tti = vu_get(hdri, /uttime )
zz = vu_atlocation(xci,rri,dri, vu_fnc(id,ffi,/t3d_array), $
sequence=tti, location=rr, times=tti, /cv2carrington, /periodic)
zz = reform(zz, diameter, diameter, /overwrite)
IF circle THEN IF outside[0] NE -1 THEN zz[outside] = BadValue(zz)
old_device = !D.name ; Store current plot device
set_page, $
zbuffer = 1-printer , $
printer = printer , $
/silent , $
xysize = xysize , $
fullsize = printer , $
old_device = old_device , $
ctable = ctable , $
_extra = _extra
title = vu_fnc(id,/symbol,/norm)
title += ' ('+vu_fnc(id,/unit)+')'
PlotPlanarCut, zz, ut=utt , $
radius = radius , $
euler_angles= euler_angles , $
degrees = degrees , $
euler_info = euler_info , $
earth = earth , $
breakval = breakval , $
charsize = charsize , $
title = title , $
body = body , $
_extra = _extra
if (istype(module, /undefined) or (n_elements(type) eq 2)) then begin
if (istype(module, /defined) and (n_elements(type) eq 2) and (itype eq 0)) then destroyvar, module
boost, module, display + label
end
vu_get_page, hdri, title='', module[itype], ut=utt, device=old_device, $
silent=silent, new=itype ne 0, _extra=_extra
ENDFOR
RETURN & END