; docformat = 'rst' ; ; NAME: ; facetProperties ; PURPOSE: ; Compute the geometric properties of each facet of a shape model ; given illumination and viweing geometry. ;+ ; :Description: ; Compute the geometric properties of each facet of a shape model ; given illumination and viweing geometry. ; ; :Categories: ; Shape Models ; ; :Params: ; model: in, required, type=fltarr or structure ; The cartesian coordinates of the shape vertices, or a model ; structure as provided by readVer. ; P: in, optional, type=lonarr() ; Connectivity array ([N_vertex, ID_vertex....]), if model are ; coordinates only. ; ; :Returns: The normal vectors to all facets ; ; :Keywords: ; period: in, optional, type=float ; Rotation period in hours (required for radial velocity) ; RV: in, optional, type=boolean, default=0 ; Computes radial velocity (km/s) of each facet ; incEm: in, optional, type=boolean, default=0 ; Only returns cosines of incidence and emission angles ; shadow: in, optional, type=boolean, default=0 ; Returns after computation of self-shadowing ; ; :Examples: ; ; :Uses: ; facetNormals, orientModel ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in April 2014, B. Carry (IMCCE) ; 2016 July - B. Carry (OCA) - Added idl2 compile option ;- function facetProperties, model, P, SEP=SEP, SSP=SSP, period=period, $ RV=RV, incEm=incEm, shadow=shadow ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; COMPILE_OPT hidden, idl2 ;--I.1-- Input Syntax Verification ---------------------------------------------------- if N_params() LT 1 then begin message, /IOERROR, 'Syntax - properties = facetProperties(model [, SEP=, SSP=])' return, -1 endif ;--I.2-- Interpretation in V input: Structure of Coordinates -------------------------- dimMod=size(model) if dimMod[dimMod[0]+1] ne 8 then begin ;--I.2.1-- Not a Model Structure -> Connectivity Required if not keyword_set(P) then begin message, /IOERROR, 'Syntax - properties = facetProperties(V, P [, SEP=, SSP=])' return, -1 endif ;--I.2.2-- Create the Model Structure V = model nbFacet = (n_elements(P))/4. nbVertex= (n_elements(V))/3. model={ vertex: {N: nbVertex, xyz: V}, $ facet: {N: nbFacet, con: P} } endif ;--I.3-- Dealing with Geometry -------------------------------------------------------- if not keyword_set(sep) then begin message, /INFO, 'Sub-Earth point not submitted, using default (0.,0.)' sep={lon:0.,lat:0.} endif if not keyword_set(ssp) then begin message, /INFO, 'Subsolar point not submitted, using sub-Earth point ('+$ strtrim(string(sep.lon,format='(F5.1)'),2)+','+$ strtrim(string(sep.lat,format='(F5.1)'),2)+')' ssp=sep endif ;--I.4-- Facet Informative Structure -------------------------------------------------- empty={visible:0, illum:0, shadow:0, $ ;-Visible and Illuminated Flag cosEm:0., cosIn:0., $ ;-Cosine of Incidence and Emission angles rv: {sun:0., obs:0.}, $ ;-Radial Velocity (towar sun and observer) area:0., $ ;-Surface area scat: {lambert: 0., $ ;-Scattering laws: Lambert LS: 0., $ ;- Lommel-Seelinger hapke: 0.} } ;- Hapke facet=replicate(empty,model.facet.N) ;--I.5-- Period Input for Radial Velocity --------------------------------------------- if keyword_set(RV) then begin if not keyword_set(period) then begin message, /IOERROR, 'Provide a period (h) if you want the radial velocity' return, -1 endif else begin spinAxis = [0,0, 2*!PI/(period*3600.)] endelse endif ;--I.6-- Optimization Parameters ------------------------------------------------------ distMargin = 1.2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Illumination and Viewing Geometries -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Direction of the Sun --------------------------------------------------------- xSun = cos(ssp.lon*!DTOR) * cos(ssp.lat*!DTOR) ySun = sin(ssp.lon*!DTOR) * cos(ssp.lat*!DTOR) zSun = sin(ssp.lat*!DTOR) vSun = [xSun,ySun,zSun] ;--II.2-- Direction of the Observer ---------------------------------------------------- xObs = cos(sep.lon*!DTOR) * cos(sep.lat*!DTOR) yObs = sin(sep.lon*!DTOR) * cos(sep.lat*!DTOR) zObs = sin(sep.lat*!DTOR) vObs = [xObs,yObs,zObs] ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Visible and Illuminated Facet -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Compute Normals to Each Facets ----------------------------------------------- facetField = tag_names( model.facet ) checkNormal = where( strCmp(facetField,'NORMAL') ) if checkNormal[0] eq -1 then verticals = facetNormals( model.vertex.xyz, model.facet.con ) $ else verticals = model.facet.normal ;--III.2-- Compute Sun & Observer Altitude for Each Facet ------------------------------- facet.cosIn = reform(verticals ## transpose(vSun)) facet.cosEm = reform(verticals ## transpose(vObs)) ;--III.3-- Search Visible and Illuminated Facets ---------------------------------------- inLight = where( facet.cosIn ge 0., nbIll ) inSight = where( facet.cosEm ge 0., nbVis ) if nbIll ne 0 then facet[inLight].illum = 1 if nbVis ne 0 then facet[inSight].visible= 1 ;--III.4-- Return if Requested ---------------------------------------------------------- if keyword_set(incEm) then return, facet ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Compute Self-Shadowing -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Orient Model with Sun at +z -------------------------------------------------- verTest = orientModel( model.vertex.xyz, ssp, 0 ) ;--IV.2-- Facet Barycenters & Inter-Facet Distances ------------------------------------ baryArr = fltarr(3,model.facet.N) interDist= fltarr(model.facet.N,model.facet.N) for kFacet=0, model.facet.N-1 do begin ;--IV.2.1-- Barycenter Location baryArr[*,kFacet] = total( verTest[*, model.facet.con[4*kFacet+[1,2,3]] ], 2 ) / 3. ;--IV.2.2-- Inter-Facet Distance for kOther=0, kFacet-1 do begin interDist[kFacet,kOther] = sqrt( (baryArr[0,kOther] - baryArr[0,kFacet] )^2. + $ (baryArr[1,kOther] - baryArr[1,kFacet] )^2. ) interDist[kFacet,kOther] = interDist[kOther,kFacet] endfor endfor ;--IV.3-- Search Overlapping Facets ---------------------------------------------------- for kIll=0, nbIll-1 do begin ;--IV.3.1-- Transpose Indice to General Facets ------------------------------------------- kFacet = inLight[kIll] ;--IV.3.2-- Local Inter-Facet Distance --------------------------------------------------- locDist = interDist[*,kFacet] ordreDist=sort(locDist) ;--IV.3.3-- Where Barycenter are Included in Current Facet ------------------------------- sideVert = verTest[*,model.facet.con[4*kFacet+[1,2,3]]] sideVert[0,*] = ( sideVert[0,*] - baryArr[0,kFacet] )^2. sideVert[1,*] = ( sideVert[1,*] - baryArr[1,kFacet] )^2. sideVert[2,*] = ( sideVert[2,*] - baryArr[2,kFacet] )^2. sideDist = sqrt( total( sideVert, 1) ) trash = where( locDist le sideDist*distMargin, nbToCheck) ;--IV.3.4-- Facets Overlap when Barycenter is Encompassed -------------------------------- for kBary=1, nbToCheck-1 do begin ;--IV.3.4/A-- Vertex of the Close-by Facet cornerArr = model.facet.con[4*ordreDist[kBary]+[1,2,3]] X1=verTest[0,cornerArr[0]] & X2=verTest[0,cornerArr[1]] & X3=verTest[0, cornerArr[2]] Y1=verTest[1,cornerArr[0]] & Y2=verTest[1,cornerArr[1]] & Y3=verTest[1, cornerArr[2]] Z1=verTest[2,cornerArr[0]] & Z2=verTest[2,cornerArr[1]] & Z3=verTest[2, cornerArr[2]] ;--IV.3.4/B-- Side Equations (lines) of the Close-by Facet slope12= ( Y2-Y1 ) / ( X2-X1 ) slope13= ( Y3-Y1 ) / ( X3-X1 ) slope23= ( Y3-Y2 ) / ( X3-X2 ) ordOrigin12= Y1 - slope12*X1 ordOrigin13= Y1 - slope13*X1 ordOrigin23= Y2 - slope23*X2 ;--IV.3.4/C-- If Barycenter is Within the Facet -> Mutual Shadowing if ( (baryArr[1,kFacet] - slope12 * baryArr[0,kFacet] - OrdOrigin12) * (Y3 - slope12 * X3 - OrdOrigin12) gt 0 ) and $ ( (baryArr[1,kFacet] - slope13 * baryArr[0,kFacet] - OrdOrigin13) * (Y2 - slope13 * X2 - OrdOrigin13) gt 0 ) and $ ( (baryArr[1,kFacet] - slope23 * baryArr[0,kFacet] - OrdOrigin23) * (Y1 - slope23 * X1 - OrdOrigin23) gt 0 ) $ then begin ;--IV.3.4/D-- Facet Shadowed is the one Behind if baryArr[2,kFacet] ge baryArr[2,ordreDist[kBary]] then $ facet[ordreDist[kBary]].shadow = 1 $ else $ facet[kFacet].shadow = 1 ;--IV.3.4/E-- Speed Up: Facet Already Shadowed goto, nextBarycenter endif endfor ;-Loop on Sorted Barycenters nextBarycenter: endfor ;-Loop on Illuminated Facets ;--IV.4-- Return if Requested ---------------------------------------------------------- if keyword_set(shadow) then return, facet ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Illumination and Viewing Geometries -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--V.1-- Scattering Laws --------------------------------------------------------------- ;--V.1.1-- Lambert facet.scat.lambert = facet.cosEm ;--V.1.2-- Lommel-Seelinger facet.scat.LS = facet.cosEm / (facet.cosEm+facet.cosIn) ;--V.1.3-- Hapke facet.scat.hapke = 0. ;--V.2-- Surface Area & Flux ---------------------------------------------------------- for kFacet=0L, model.facet.N-1L do begin ;--V.2.1-- Facet Vertex vertX = model.vertex.xyz[0, model.facet.con[4*kFacet+[1,2,3] ] ] vertY = model.vertex.xyz[1, model.facet.con[4*kFacet+[1,2,3] ] ] vertZ = model.vertex.xyz[2, model.facet.con[4*kFacet+[1,2,3] ] ] ;--V.2.2-- Facet Area facet[kFacet].area = triangleArea( vertX, vertY, vertZ ) endfor ;--V.3-- Radial Velocity -------------------------------------------------------------- if keyword_set(RV) then $ for kFacet=0L, model.facet.N-1L do begin ;--V.3.1-- Facet Vertex vertX = model.vertex.xyz[0, model.facet.con[4*kFacet+[1,2,3] ] ] vertY = model.vertex.xyz[1, model.facet.con[4*kFacet+[1,2,3] ] ] vertZ = model.vertex.xyz[2, model.facet.con[4*kFacet+[1,2,3] ] ] ;--V.3.2-- Facet Average Vector vecFacet = [ mean(vertX), mean(vertY), mean(vertZ) ] ;--V.3.3-- Facet Velocity vel = CROSSP(spinAxis, vecFacet) ;--V.3.4-- Facet Velocity Projected on Line of Sights facet[kFacet].rv.sun = TRANSPOSE(vSun) # vel facet[kFacet].rv.obs = TRANSPOSE(vObs) # vel endfor return, facet end