; docformat = 'rst' ; ; NAME: ; shapeShow3D ; PURPOSE: ; Draw a 3-D representation of a shape model on the plane of the sky ;+ ; :Description: ; Draw a 3-D representation of a shape model on the plane of the sky ; ; :Categories: ; Ephemeris, Visualization ; ; :Params: ; MODEL: in, required, type=string/structure ; A shape model structure or a path to a shape model file. ; Fields (see readVER) are:: ; .VERTEX: A structure with Information on Vertex:: ; .N = Number of Vertex ; .XYZ= Cartesian coordinates of the vertex ; .FACET: A structure with Information on Facet:: ; .N = Number of Facets ; .CON = Connectivity array ([N_vertex, ID_vertex....]) ; .NORMAL= Array of facet normals (only if keyword NORMALS is set) ; .EDGES: A structure with Information on Edges (only if keyword EDGES is set):: ; .N = Number of Edges ; .vertex= Array of vertex per edge ; .facet = Array of facets per edge ; .edge = Array of edges per facet ; ; :Returns: An image of the shape model on the plane of the sky. ; ; :Keywords: ; SEP: in, optional, type=fltarr(2) ; Longitude and Latitude of the sub-Earth point (deg.) ; SSP: in, optional, type=fltarr(2) ; Longitude and Latitude of the subsolar point (deg.) ; PA: in, optional, type=float ; Angle between celestial north and rotation axis (deg., ; positive through East) ; ; :Uses: ; orientModel, orientVector, facetProperties ; ; :Examples: ; Read target information, orientation, and then draw a 3-D view:: ; IDL> asterId=1 ; IDL> model = readVer('ceres.ver') ; IDL> spin = readSpin('ceres.spin') ; IDL> ep='2001-01-01T00:00:00.000' ; IDL> vGeo = ssoPositions(asterId,ep,type='aster', /geocentric) ; IDL> vSun = ssoPositions(asterId,ep,type='aster', /heliocentric) ; IDL> sep = subObserver_coord( -vGeo, s,ep) ; IDL> ssp = subObserver_coord( -vSun, s,ep) ; IDL> image = shapeShow3d( model, spin, sep=sep, ssp=ssp) ; IDL> write_png, 'my_view_of_Ceres.png', image ; ; Simplier call:: ; IDL> asterId=1 ; IDL> model = 'MYPATH/ceres.ver' ; IDL> spin = 'MYPATH/ceres.spin' ; IDL> sep = {lon: 123.0, lat:-4.5} ; IDL> ssp = {lon: 133.0, lat:-3.5} ; IDL> image = shapeShow3d( model, spin, sep=sep, ssp=ssp, pa=-2) ; IDL> write_png, 'my_view_of_Ceres.png', image ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in April 2014, B. Carry (IMCCE), based on VISU (Berthier, Carry, Hestroffer) ; 2016 Mar - B. Carry (OCA) - Completed information header + idl2 added ; 2017 Mar - B. Carry (OCA) - Now save current device and restore it at the end ;- function shapeShow3D, model, spin, res=res, fov=fov, coords=coords, $ SEP=SEP, SSP=SSP, PA=PA, dObs=dObs, dSun=dSun, phase=phase COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; spin=0 ;--I.1-- Input Syntax Verification ---------------------------------------------------- if n_params() lt 1 then begin message, /ioError, 'Syntax : out = shapeShow3D( model [, spin, res=, ep=, ...] )' return, -1 endif ;--I.2-- Shape Model & Spin Solution -------------------------------------------------- ;--I.2.1-- Model = Structure or Path dimMod=size(model,/type) if dimMod eq 7 then if file_test(model,/Read) then model=readVer(model) ;--I.2.2-- Spin Solution ---------------------------------------------------------------- dimSpin=size(spin,/type) if dimMod eq 7 then if file_test(spin,/Read) then spin=readSpin(spin) ;--I.3-- FOV & Resolution ------------------------------------------------------------- if not keyword_set(res) then res={x:2048,y:2048} if not keyword_set(fov) then $ fov={x:1.2, y:1.2, $ ;-Size of the Field of View along x & y axes pa: 0., $ ;-Rotation of the Field of View (from North to y, direct sense, deg) unit:'model' } ;-Unit for the size of the FoV (deg|min|sec|model|km) ;--I.4-- Ephemeris for Viewing Geometry ----------------------------------------------- if not keyword_set(coords) then begin if not keyword_set(sep) then sep={lon:0.,lat:0.} if not keyword_set(ssp) then ssp=sep if not keyword_set(pa) then pa=0. ephem={sep:sep, ssp:ssp, pa:pa, dObs:dObs, dSun:dSun} endif else begin ephem=coords if keyword_set(sep) then ephem.sep=sep if keyword_set(ssp) then ephem.ssp=ssp if keyword_set(pa) then ephem.pa =pa endelse ;--I.7-- Save Current Display Settings ------------------------------------------------ SetDecomposedState, 1, CurrentState=theState thisDevice = !D.Name ;--I.8-- Useful Constants ------------------------------------------------------------- AUtoKM = 149597870.691D DegToMAS = 3600000.0D ArcsecToRad = !DTOR/3600.0D ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Shape Model Orientation & Facet Properties -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Orient Shape Model with Observer Direction ----------------------------------- vertex = orientModel( Model.vertex.xyz, ephem.sep, ephem.pa ) ;--II.2-- Place Sun around the Shape Model --------------------------------------------- ecSun = cv_coord(from_sphere=[ephem.sep.lon,ephem.sep.lat,1], /TO_RECT, /DEGREE) vSun = orientVector( ecSun, ephem.sep, ephem.pa ) ;--II.3-- Facet Properties ------------------------------------------------------------- facet = facetProperties( model, sep=ephem.sep, ssp=ephem.ssp ) ;--III.3-- Valid Facets are BOTH Visible and Illuminated -------------------------------- ; visible = where( DotProduct_Sun gt 0 and DotProduct_Earth gt 0 ) ; facets = intarr(model.facet.n) ; facets(visible) = 1 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Field of View and Illumination Setup -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Light from the Sun ---------------------------------------------------------------- set_shading, values=[0,250], LIGHT=vSun, /reject ; print, fov.unit ; print, 'model ', max(sqrt( total(vertex^2, 1) ))*[fov.x,fov.y] ; print, 'min ', ephem.dObs*AUtoKm*tan([fov.x,fov.y]*!DTOR/60.) ; print, 'sec ', ephem.dObs*AUtoKm*tan([fov.x,fov.y]*!DTOR/3600.) ; ; print, ephem.dObs, AUtoKm , tan([fov.x,fov.y]*!DTOR/60.) ; print, [fov.x,fov.y]*!DTOR/60. ; print, [fov.x,fov.y] ;--III.2-- Size of the Field of View --------------------------------------------------------- case fov.unit of ;--III.2.1-- FoV in Model Units 'model': rWin = max(sqrt( total(vertex^2, 1) ))*[fov.x,fov.y] ;--III.2.2-- FoV in Kilometers 'km': rWin = [fov.x,fov.y] ;--III.2.3-- FoV in Degrees 'deg': rWin = ephem.dObs*AUtoKm*tan([fov.x,fov.y]*!DTOR) ;--III.2.4-- FoV in Arcminutes 'min': rWin = ephem.dObs*AUtoKm*tan([fov.x,fov.y]*!DTOR/60.) ;--III.2.4-- FoV in Arcseconds 'sec': rWin = ephem.dObs*AUtoKm*tan([fov.x,fov.y]*!DTOR/3600.) else: endcase ;--III.3-- Polyshade 3-D View of Shape Model ------------------------------------------------ set_plot, 'X' scale3, xrange=rWin[0]*[-1.0,1.0], yrange=rWin[1]*[-1.0,1.0], zrange=max(rWin)*[-2,2] T3D, /RESET image = POLYSHADE( vertex, model.facet.con, xsize=res.x, ysize=res.y, /T3D) wdelete ;--III.4-- Put Display in Original State and Return ----------------------------------------- set_plot, thisDevice SetDecomposedState, theState return, image end