; docformat = 'rst' ; ; NAME: ; shapeProfile ; PURPOSE: ; Extract the profile of a 3-D shape model on the plane of the sky, based on viewving/illumination geometry. ;+ ; :Description: ; Extract the profile of a 3-D shape model on the plane of the sky, based on viewving/illumination geometry. ; ; :Categories: ; Shape Models ; ; :Params: ; model: in, required, type=fltarr or structure ; A model structure as provided by readVer. ; ; :Returns: The cartesian coordinates of the profile on the plane of the sky. ; ; :Keywords: ; SEP: in, optional, type=float, default="[0.,0.]" ; Longitude and Latitude of the sub-Earth point (deg.). ; SSP: in, optional, type=float, default="[0.,0.]" ; Longitude and Latitude of the subsolar point (deg.). ; PA: in, optional, type=float, default=0. ; Pole angle to the North (deg.). ; barycenter: out, optional, type=float ; Position of the mass barycenter on the plane of the sky. ; nBin: in, optional, type=integer ; If set, interpolate the profile on nBin points. ; Az: in, optional, type=fltarr ; If set, use Az as azimuth to interpolate the profile ; ; :Uses: ; orientModel, facetProperties ; ; :Examples: ; Read a model profile and extract its profile:: ; IDL> model = readVer( 'my_model_file' ) ; IDL> sep = {lon:15., lat:+2.0} ; IDL> ssp = {lon:35., lat:+2.0} ; IDL> profile = shapeProfile( model, sep=sep, ssp=ssp ) ; IDL> cgPlot, /Iso, profile[0,*], profile[1,*] ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in June 2012, B. Carry (ESA) ; 2013 Nov. - B.Carry (IMCCE) - Changed I/O keywords ; 2016 July - B. Carry (OCA) - Added idl2 compile option ; 2016 Aug. - B. Carry (OCA) - Improved profile detection (vs shadow) & Az keyword ;- function shapeProfile, model, SEP=SEP, SSP=SSP, PA=PA, barycenter=barycenter, nBin=nBin, Az=Az COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Input Syntax Verification ----------------------------------------------------------- if not keyword_set(model) then begin message, /IOERROR, 'Syntax : xy = shapeProfile( model, [sep=, ssp=, pa=] )' return, -1 endif ;--I.2-- Dealing with Geometry --------------------------------------------------------------- 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. ;--I.3-- Check Input Model Structure & Update ------------------------------------------------ dimMod=size(model) if dimMod[dimMod[0]+1] ne 8 then begin message, /IOERROR, 'Input variable "model" is not a structure' return, -2 endif ; tags=tag_names(model) ;-TBD shapeProfile - check that all model.tags are present ;--I.4-- Useful Constants -------------------------------------------------------------------- AUtoKM = 149597870.691D DegToMAS = 3600000.0D ArcsecToRad = !DTOR/3600.0D ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Shape Model Orientation -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Orient Shape Model with Observer Direction ----------------------------------------- objVertex = orientModel( model.vertex.xyz, sep, pa ) ;--II.2-- Barycenter Position ---------------------------------------------------------------- volume = tetra_volume(objVertex, model.facet.con, moment=moments) barycenter = moments/volume ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Visible and Illuminated Facet -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Compute Sun & Observer Altitude for Each Facet ------------------------------------ info = facetProperties( model, sep=sep, ssp=ssp, /incEm ) ;--III.2-- Valid Facets are BOTH Visible and Illuminated ------------------------------------- valid = where( info.visible eq 1 and info.illum eq 1 ) facets = intarr(model.facet.n) facets[valid] = 1 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Determination of (Limb+Terminator) Profile -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Find Edges Indexes on the Profile -------------------------------------------------- edgeStatus = intarr(model.edges.n,2) for kFacet=0, model.facet.n-1 do begin pp = where( model.edges.facet eq kFacet ) edgeStatus[pp] = facets[kFacet] endfor sideID = where( total(edgeStatus,2) eq 1, nbSide ) ;--IV.2-- List Edges from Shadows for Further Removal ---------------------------------------- toRemove = -1 for kSide=0L, nbSide-1L do begin ;--IV.2.1-- Current Side Indices of Vertex id1 = model.edges.vertex[sideID[kSide],0] id2 = model.edges.vertex[sideID[kSide],1] ;--IV.2.2-- Search Intersection on Previous Sides for kOth=0, kSide-1 do begin ;--IV.2.2/A-- Other Side Vertex Indices other1 = model.edges.vertex[sideID[kOth],0] other2 = model.edges.vertex[sideID[kOth],1] ;--IV.2.2/B-- Current Side Vertex 1 Intersection intersect1 = interLine( [0,0], objVertex[*,id1], objVertex[*,other1], objVertex[*,other2], $ flag=flag1, distance=sep1, /Line1) dist1 = sqrt(objVertex[0,id1]^2 + objVertex[1,id1]^2) ;--IV.2.2/C-- Current Side Vertex 2 Intersection intersect2 = interLine( [0,0], objVertex[*,id2], objVertex[*,other1], objVertex[*,other2], $ flag=flag2, distance=sep2, /Line1) dist2 = sqrt(objVertex[0,id2]^2 + objVertex[1,id2]^2) ;--IV.2.2/D-- Remove Side if Outer Intersection if ( (flag1 eq 1 and dist1 lt sep1) and (flag2 eq 0 or flag2 eq 2) ) or $ ;-Id1 only intersects ( (flag2 eq 1 and dist2 lt sep2) and (flag1 eq 0 or flag1 eq 2) ) or $ ;-Id2 only intersects ( (flag1 eq 1 and dist1 lt sep1) and (flag2 eq 1 and dist2 lt sep2) ) then begin ;-Both intersects toRemove = [toRemove, kSide] kOth=nbSide+1 endif endfor ;--IV.2.3-- Search Intersection on Following Sides for kOth=kSide+1, nbSide-1 do begin ;--IV.2.3/A-- Other Side Vertex Indices other1 = model.edges.vertex[sideID[kOth],0] other2 = model.edges.vertex[sideID[kOth],1] ;--IV.2.3/B-- Current Side Vertex 1 Intersection intersect1 = interLine( [0,0], objVertex[*,id1], objVertex[*,other1], objVertex[*,other2], $ flag=flag1, distance=sep1, /Line1) dist1 = sqrt(objVertex[0,id1]^2 + objVertex[1,id1]^2) ;--IV.2.3/C-- Current Side Vertex 2 Intersection intersect2 = interLine( [0,0], objVertex[*,id2], objVertex[*,other1], objVertex[*,other2], $ flag=flag2, distance=sep2, /Line1) dist2 = sqrt(objVertex[0,id2]^2 + objVertex[1,id2]^2) ;--IV.2.3/D-- Remove Side if Outer Intersection if ( (flag1 eq 1 and dist1 lt sep1) and (flag2 eq 0 or flag2 eq 2) ) or $ ;-Id1 only intersects ( (flag2 eq 1 and dist2 lt sep2) and (flag1 eq 0 or flag1 eq 2) ) or $ ;-Id2 only intersects ( (flag1 eq 1 and dist1 lt sep1) and (flag2 eq 1 and dist2 lt sep2) ) then begin ;-Both intersects toRemove = [toRemove, kSide] kOth=nbSide+1 endif endfor endfor ;--IV.2.3-- Remove Multiple Entries if n_elements(toRemove) ge 2 then toRemove=toRemove[1:n_elements(toRemove)-1] uniqRemove = uniq(toRemove,sort(toRemove)) toRemove=toRemove[uniqRemove] ;--IV.3-- Find Edges Indexes on the Profile -------------------------------------------------- ;--IV.3.1-- Remove Inner Edges from sideID List select = sideId*0 select[toRemove] = 1 good = where( select eq 0, nbGood) sideId = sideId[good] ;--IV.3.2-- Retrieve Indices vertexID = 0 for kSide=0L, nbGood-1L do begin pair = model.edges.vertex[sideID[kSide],*] vertexID = [vertexID, reform(pair)] endfor nbCorner=n_elements(vertexID) vertexID=vertexID[1:nbCorner-1] ;--IV.4-- Extract Profile Corners from Shape ------------------------------------------- profile = objVertex[0:1,vertexID] ;--IV.5-- Order Profile in Azimuth ---------------------------------------------------- sideAngle = atan( profile[1,*], profile[0,*] ) sideOrd = sort( sideAngle ) profile = profile[*,sideOrd] ;--IV.6-- Return Profile if Not Interpolation Required -------------------------------- if not keyword_set(nBin) and not keyword_set(Az) then return, profile ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Profile Smoothing to Regular Angular Bins -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--V.1-- Interpolation Grid Provided in Input ------------------------------------------------ if keyword_set(Az) then begin theta = Az nBin = n_elements(Az) ;--V.2-- Define the Angular Regular Grid ----------------------------------------------------- endif else begin theta = 2.*!DPI*indgen(nBin)/float(nBin)-!PI endelse ;--V.3-- Interpolate Profile on the Regular Grid --------------------------------------------- ;--V.3.1-- Replicate Contour for Interpolation pol = cv_coord( from_rect=profile, /to_polar ) lonRad = [ reform(pol[1,*]), reform(pol[1,*]), reform(pol[1,*]) ] lonAng = [ reform(pol[0,*])-2*!PI, reform(pol[0,*]), reform(pol[0,*])+2*!PI] ;--V.3.2-- Contour Smoothing and Conversion to Cartesian smoo = interpol( lonRad, lonAng, theta ) newXY = cv_coord( from_polar=transpose([ [theta],[smoo] ]), /to_rect ) ;--V.4-- Return Interpolated Profile -------------------------------------------------------- return, newXY end