; docformat = 'rst' ; ; NAME: ; measureVer ; PURPOSE: ; Measure dimensions, surface, volume, and orientation of a shape model ; ;+ ; :Description: ; Measure dimensions, surface, volume, and orientation of a shape model ; ; :Categories: ; Asteroid, Shape Models ; ; :Params: ; MODEL: in, required, type=string/structure ; The model structure (see readVer) or the shape model file ; (VER format) ; ; :Returns: A structure with the geometric properties of the shape ; model. Fileds are:: ; .S: Surface area ; .V: Volume ; .Origin: Cartesian coordinates of the center of mass ; .D: Volume equivalent diameter ; .R: Volume equivalent Radius ; .A: Triaxial ellipsoid radius - long ; .B: Triaxial ellipsoid radius - medium ; .C: Triaxial ellipsoid radius - short ; .ABC: Triaxial ellipsoid radius - All three ; .Area: Surface area of each facet ; .PMI: Principal Moment of Inertia -- Structure: ; .val: Moment of inertia ; .lon: Longitude of the Moment of Inertia ; .lat: Longitude of the Moment of Inertia ; ; :Keywords: ; plot: in, optional, type=string ; If set, place several validation graphics in the "plot" directory. ; ; :Examples: ; Measure a shape model and display its volume-equivalent diameter:: ; IDL> value = measureVer( 'my_model.ver' ) ; IDL> print, value.D ; ; :Uses: ; readVer(), Coyote Graphic System ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in 2011, B. Carry (ESA) ; 2014 Apr. - B.Carry (IMCCE) - Compute facets area ; 2014-Aug. - B.Carry (IMCCE) - Removed dependence on cart2sph ; 2017 Mar. - B.Carry (OCA) - idl2 - Added validation plots ;- function measureVer, model, plot=plot COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Input Syntax Verification ---------------------------------------------------- if n_params() LT 1 then begin message, /ioError, 'Syntax - values = measureVer(model)' return, -1 endif ;--I.2-- Interpretation of Input Model ------------------------------------------------ dimMod=size(model,/Type) if dimMod ne 8 then begin ;--I.2.1-- Not a Model Structure -> Read file if not file_test(model,/READ) then begin message, /ioError, 'Cannot find the input file: '+strtrim(model,2) return, -1 endif else begin model=readVer(model) endelse endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- General Properties and Overall Dimensions -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Surface and Volume ---------------------------------------------------------- surface= mesh_surfacearea( model.vertex.xyz, model.facet.con) volume = tetra_volume(model.vertex.xyz, model.facet.con, moment=moments) origin = moments / volume ;--II.2-- Volume-equivalent Radius & Diameter ----------------------------------------- ved = (6.*volume/!PI) ^ (1./3.) ;--II.3-- Encompassing Box ------------------------------------------------------------ box = {x: [min(model.vertex.xyz[0,*]), max(model.vertex.xyz[0,*])], $ y: [min(model.vertex.xyz[1,*]), max(model.vertex.xyz[1,*])], $ z: [min(model.vertex.xyz[2,*]), max(model.vertex.xyz[2,*])] } ;--II.4-- Facet Properites ------------------------------------------------------------ area=fltarr(model.facet.N) for kFacet=0, model.facet.N-1 do begin ;--II.4.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]] ] ;--II.4.2-- Facet Area area[kFacet] = triangleArea( vertX, vertY, vertZ ) endfor ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Principal Moment of Inertia -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Tetrahedra volumes & Inertia Tensor ----------------------------------------- tetraVol= dblarr(model.facet.N) inertia = dblarr(3,3) for kFacet=0L, model.facet.N-1 do begin ;--III.2-- Tetrahedra with (0,0,0) as Fourth Summit ------------------------------------ ;--III.2.1-- Vertex Coordinates xx = [ reform(model.vertex.xyz[0,model.facet.con[4*kFacet+[1,2,3]]]), origin[0] ] yy = [ reform(model.vertex.xyz[1,model.facet.con[4*kFacet+[1,2,3]]]), origin[1] ] zz = [ reform(model.vertex.xyz[2,model.facet.con[4*kFacet+[1,2,3]]]), origin[2] ] connect = [0,1,2,3] ;--III.2.2-- Tetrahedron Volume & Barycenter tetraVol[kFacet] = tetra_volume( transpose([[xx],[yy],[zz]]), connect, moment=mmt) barycenter = mmt / tetraVol[kFacet] ;--III.3-- Tetrahedron Moment of Inertia ----------------------------------------------- ;--III.3.1-- Translate Coordinates x= xx - barycenter[0] y= yy - barycenter[1] z= zz - barycenter[2] ;--III.3.2-- Moment and Products of Inertia, cf Tonon 2004 kerX = [ [replicate(x[0],4)], $ [replicate(x[1],4)], $ [replicate(x[2],4)], $ [replicate(x[3],4)] ] setX = kerX * ( diag_matrix(replicate(1.,4)) + $ diag_matrix(replicate(1.,3), -1) + $ diag_matrix(replicate(1.,2), -2) + $ diag_matrix(replicate(1.,1), -3) ) kerY = [ [replicate(y[0],4)], $ [replicate(y[1],4)], $ [replicate(y[2],4)], $ [replicate(y[3],4)] ] setY = kerY * ( diag_matrix(replicate(1.,4)) + $ diag_matrix(replicate(1.,3), -1) + $ diag_matrix(replicate(1.,2), -2) + $ diag_matrix(replicate(1.,1), -3) ) kerZ = [ [replicate(z[0],4)], $ [replicate(z[1],4)], $ [replicate(z[2],4)], $ [replicate(z[3],4)] ] setZ = kerZ * ( diag_matrix(replicate(1.,4)) + $ diag_matrix(replicate(1.,3), -1) + $ diag_matrix(replicate(1.,2), -2) + $ diag_matrix(replicate(1.,1), -3) ) ;--III.3.2/A-- Moments a = 6.* tetraVol[kFacet] *( total(diag_matrix(kerY#setY)) + total(diag_matrix(kerZ#setZ)) )/60.D b = 6.* tetraVol[kFacet] *( total(diag_matrix(kerX#setX)) + total(diag_matrix(kerZ#setZ)) )/60.D c = 6.* tetraVol[kFacet] *( total(diag_matrix(kerX#setX)) + total(diag_matrix(kerY#setY)) )/60.D ;--III.3.2/B-- Products ap= 6.* tetraVol[kFacet] *( total( diag_matrix(kerZ#(kerY+ kerY#identity(4) )) ) ) / 120.D bp= 6.* tetraVol[kFacet] *( total( diag_matrix(kerZ#(kerX+ kerX#identity(4) )) ) ) / 120.D cp= 6.* tetraVol[kFacet] *( total( diag_matrix(kerY#(kerX+ kerX#identity(4) )) ) ) / 120.D ;--III.3.2/C-- Tensor of Inertia in the Local Coordinates tensor = [ [ a, -bp, -cp], $ [-bp, b, -ap], $ [-cp, -ap, c] ] ;--III.3.3-- Tensor of Inertia at the Origin trans = (barycenter-origin) inertia += (tensor + tetraVol[kFacet]*( diag_matrix(replicate(transpose(trans)#trans,3)) - trans#trans )) endfor ;--III.4-- Reduction of the Tensor of Inertia to Principal Axes ------------------------ PMI=replicate({val:0., lon:0., lat:0.},3) ;--III.4.1-- Eigen Vectors and Values of the Tensor of Inertia eVal = hqr(elmhes(inertia), /Double) eVec = eigenvec(inertia, eVal, residual=residual) for i=0,2 do eVec[*,i] *= ABS(eVec[0,i])/eVec[0,i] ;-Normalization (for Unicity of eVec) ;--III.4.2-- Principal Moments of Inertia PMI.val = float(eVal)/volume ;-Normalization by the Mass order=sort(PMI.val) PMI.val = PMI[order].val ;--III.4.3-- Ellipsoid Dimensions along PMI a2 = 2.5*(PMI[2].val - PMI[0].val + PMI[1].val) c2 = 5*PMI[1].val - a2 b2 = 5*PMI[0].val - c2 abcPMI = sqrt( [a2, b2, c2] ) ;--III.5-- Orientation of the Principal Axes of Inertia -------------------------------- ;--III.5.1-- Principal Axes of Inertia - Cartesian Coordinates eVec=float(eVec[*,order]) v1= eVec[*,0] v2= eVec[*,1] v3= eVec[*,2] ;--III.5.2-- Principal Axes of Inertia - Spherical Coordinates eVecSPH = cv_coord( from_rect=eVec, /to_sphere, /deg ) PMI[0].lon = eVecSPH[0,0] PMI[0].lat = eVecSPH[1,0] PMI[1].lon = eVecSPH[0,1] PMI[1].lat = eVecSPH[1,1] PMI[2].lon = eVecSPH[0,2] PMI[2].lat = eVecSPH[1,2] ; v1sp = cart2sph( v1[0], v1[1], v1[2] ) ; v2sp = cart2sph( v2[0], v2[1], v2[2] ) ; v3sp = cart2sph( v3[0], v3[1], v3[2] ) ; v1sp(1:2) /= !DTOR ; v2sp(1:2) /= !DTOR ; v3sp(1:2) /= !DTOR ; ; if v1sp[1] lt 0. then v1sp[1]+=180. ; if v2sp[1] lt 0. then v2sp[1]+=180. ; if v3sp[1] lt 0. then v3sp[1]+=180. ; ; PMI[0].lon = v1sp[1] ; PMI[0].lat = v1sp[2] ; PMI[1].lon = v2sp[1] ; PMI[1].lat = v2sp[2] ; PMI[2].lon = v3sp[1] ; PMI[2].lat = v3sp[2] ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Validation Graphics -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if keyword_set(plot) then begin ;--IV.1-- Check Directory ------------------------------------------------------------------; if not file_test(plot,/Directory) then begin file_mkdir, plot endif winPS = [20, 20.] ;--IV.2-- Model VS Ellipsoid ---------------------------------------------------------------; nameArr=['XY','XZ','YZ'] xTitle=['X','X','Y'] yTitle=['Y','Z','Z'] xArr=[0,0,1] yArr=[1,2,2] rotArr=[ pmi[0].lon, pmi[0].lat, pmi[1].lat ] *!DTOR maxR = sqrt( max( mean( (model.vertex.xyz[0,*] - origin[0])^2 + $ (model.vertex.xyz[1,*] - origin[1])^2 + $ (model.vertex.xyz[2,*] - origin[2])^2 ) ) ) xRange = [-1,1] * maxR * 1.2 yRange = xRange ;--IV.2-- Model VS Ellipsoid ---------------------------------------------------------------; for k=0, 2 do begin cgPS_open, Filename=plot+'/'+nameArr[k]+'.eps', /metric, /decomposed, /portrait, $ xsize=winPS[0], ysize=winPS[1], language_level=2, /quiet cgPlot, model.vertex.xyz[xArr[k],*] - origin[xArr[k]], $ model.vertex.xyz[yArr[k],*] - origin[yArr[k]], $ /Iso, psym='Filled circle', $ xRange=xRange, $ yRange=yRange, $ position=[0.1,0.1,0.985,0.985] cgText, 0.025, 0.05, yTitle[k], /Normal cgText, 0.05, 0.025, xTitle[k], /Normal theta = 2*!PI*findgen(360)/360 xEll = abcPMI[xArr[k]]*cos(theta) yEll = abcPMI[yArr[k]]*sin(theta) xRot = xEll*cos(rotArr[k]) - yEll*sin(rotArr[k]) yRot =+xEll*sin(rotArr[k]) + yEll*cos(rotArr[k]) cgPlot, /OverPlot, xRot, yRot, color='Gray' cgPS_close, /png, Delete_PS=0 endfor endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Summary Structure -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; return, { S:surface, $ V:volume, $ origin: origin, $ D: ved, $ R: ved/2. , $ a: abcPMI[0], $ b: abcPMI[1], $ c: abcPMI[2], $ abc: abcPMI, $ area:area, $ pmi:pmi } end ;- ;- ;- ;--III.7-- Dimensions along PMIs ;- ;--III.7.1-- Shift Model to New Origin (Barycenter) ;- model.vertex.xyz[0,*] -= origin[0] ;- model.vertex.xyz[1,*] -= origin[1] ;- model.vertex.xyz[2,*] -= origin[2] ;- ;- ;--III.7.2-- Orient PMIs along XYZ ;- T3D, /RESET ;- T3D, rotate=[0.,0., -v1sp[1]] ;- T3D, rotate=[0.,v1sp[2],0.] ;- T3D, rotate=[v2sp[2],0.,0.] ;- vert = vert_t3d( model.vertex.xyz ) ;- ;- ;--III.7.3-- Dimensions along PMIs ;- ;--III.7.3.1-- Overall Dimensions ;- PMI_OD=[ max(vert[0,*])-min(vert[0,*]), $ ;- max(vert[1,*])-min(vert[1,*]), $ ;- max(vert[2,*])-min(vert[2,*]) ]/2. ;- ;- ;- volume = tetra_volume(model.vertex.xyz, model.facet.con, moment=moments) ;- rMean = ( (3.*volume)/(4.*!PI) )^(1./3.) ;-; rMean = product( PMI_OD )^(1./3.) ;- ;- ;- ;--III.8-- Return Values to Higher Level ;- results={name: ['Radius', 'Diameter', $ ;- 'PMI OD a', $ ;- 'PMI OD b', $ ;- 'PMI OD c', $ ;- 'PMI Ell a', $ ;- 'PMI Ell b', $ ;- 'PMI Ell c', $ ;- 'OD a lon', $ ;- 'OD a lat', $ ;- 'OD b lon', $ ;- 'OD b lat', $ ;- 'OD c lon', $ ;- 'OD c lat', $ ;- 'PMI a', $ ;- 'PMI b', $ ;- 'PMI c'], $ ;- value:[rMean, 2*rMean, $ ;- PMI_OD[0], $ ;- PMI_OD[1], $ ;- PMI_OD[2], $ ;- PMI_ell[0], $ ;- PMI_ell[1], $ ;- PMI_ell[2], $ ;- v1sp[1], $ ;- v1sp[2], $ ;- v2sp[1], $ ;- v2sp[2], $ ;- v3sp[1], $ ;- v3sp[2], $ ;- PMI[0], $ ;- PMI[1], $ ;- PMI[2]], $ ;- unit: ['km', 'km', $ ;- 'km', 'km', 'km', $ ;- 'km', 'km', 'km', $ ;- 'deg', 'deg', 'deg', 'deg', 'deg', 'deg', $ ;- 'km^2', 'km^2', 'km^2'], $ ;- format: ['F7.3', 'F7.3', $ ;- 'F7.3', 'F7.3', 'F7.3', $ ;- 'F7.3', 'F7.3', 'F7.3', $ ;- 'F5.1', 'F5.1', 'F5.1', 'F5.1', 'F5.1', 'F5.1', $ ;- 'E11.3', 'E11.3', 'E11.3']} ;- return, results ;- ;-; ;-; print, '' ;-; print, '----- Principal Moments of Inertia ---' ;-; print, ' I1/M = '+string(PMI[0],format='(F8.2)')+' km^2' ;-; print, ' I2/M = '+string(PMI[1],format='(F8.2)')+' km^2' ;-; print, ' I3/M = '+string(PMI[2],format='(F8.2)')+' km^2' ;-; ;-; ;-;; print, '' ;-;; print, '----- Dimensions along Principal Axes of Inertia ---' ;-; print, ' dim1 = '+string(PMI_OD[0],format='(F8.3)')+' km @ ('+string(v1sp[1],format='(F5.1)')+','+string(v1sp[2],format='(F5.1)')+')' ;-; print, ' dim2 = '+string(PMI_OD[1],format='(F8.3)')+' km @ ('+string(v2sp[1],format='(F5.1)')+','+string(v2sp[2],format='(F5.1)')+')' ;-; print, ' dim3 = '+string(PMI_OD[2],format='(F8.3)')+' km @ ('+string(v3sp[1],format='(F5.1)')+','+string(v3sp[2],format='(F5.1)')+')' ;- ;- ;- ;-endif ;--IV--End of PMI ;- ;-