FUNCTION visu_altitude, targetName=targetName
  
  COMMON share, objVertex, objPolygons, barr, nbFacet
  
  if (not keyword_set(objVertex)) || (not keyword_set(objPolygons)) || (not keyword_set(nbFacet)) then begin
    
    if not keyword_set(targetName) then targetName='lutetia'
    
    verfile='/astrodata/SHAPE_MODELS/VER/' + targetName + '.ver'
    
    read_ver_shape_model, verFile, $
                          objVertex, objPolygons, $
                          nbVertex=nbVertex, $
                          nbFacet=nbFacet
    
  endif
  
  ;barycenters for each facet
  barr = fltarr(3, nbFacet)
  kPoly = 0
  
  for kPoly=0, nbFacet-1 do begin
    
    barr[*, kPoly] = TOTAL(objVertex[*, objPolygons[4 * kPoly + [1, 2, 3]]], 2) / 3.
    
  endfor
  
  
;  print, total(sqrt(objVertex[0, *]^2 + objVertex[1, *]^2 + objVertex[2, *]^2)) / nbVertex
  
;  funcObj = OBJ_NEW('altitudeFunc', objVertex)
  
;  X = [49., 49., 49., 0., 0., 0.]
  
  dis = sqrt(objVertex[0, *]^2 + objVertex[1, *]^2 + objVertex[2, *]^2)
  
  mDis = total(dis) / n_elements(dis)
  
  P0 = [0., 0., 0., mDis, mDis, mDis]
  
;  print, P0
  
  scl_axis = (max(dis) - min(dis)) / 2
  
  SCL = [90., 90., 90, scl_axis, scl_axis, scl_axis]
  
;  print, SCL
  
  xi = [[1., 0., 0., 0., 0., 0.], $
        [0., 1., 0., 0., 0., 0.], $
        [0., 0., 1., 0., 0., 0.], $
        [0., 0., 0., 1., 0., 0.], $
        [0., 0., 0., 0., 1., 0.], $
        [0., 0., 0., 0., 0., 1.]]
        
  Ftol = 1.e-12
  
  
;  POWELL, P0, xi, Ftol, Fmin, 'func'
  
  
  result = AMOEBA( 1.0e-5 , FUNCTION_NAME='func' , FUNCTION_VALUE=Fmin, P0=P0, SCALE=SCL ) 
 
  altitudeMap = bytscl(reform(alt_poly(result)))
    
  return, altitudeMap
  
;  create_ver_ellipsoid, objVertex, objPolygons, $
;                          axes= [ result[3:5], $
;                                  result[3:5] ], $
;                          nbVertex=nbVertex, $
;                          nbFacet=nbFacet
  
;  resolution = [1024, 1024]
  
;  set_plot, 'Z'
;  device, set_resolution=resolution
  
  
  
;  set_plot, 'X'
;  DEVICE, DECOMPOSED = 0  
;  LOADCT, 11
;  LOADCT, 0
;  TVLCT, 255, 255, 255, !D.TABLE_SIZE - 1  
;  WINDOW, 1, XSIZE = 512, YSIZE = 512 
  
;  rMax = max(sqrt( total(objVertex^2, 1) ))
  
;  winScale = 1.
  
;  R = [winScale*rMax, winScale*rMax]
  
;  scale3, xrange=[-R(0),R(0)], $
;          yrange=[-R(1),R(1)], zrange=[-max(R), max(R)]
  
;  r = polyshade(objVertex, objPolygons, SHADES=altitudeMap)
  
;  TVSCL, polyshade(objVertex, objPolygons, POLY_SHADES=altitudeMap)
  
  
;  image = tvrd()
;  
;  image = bytscl(image)
;  
;  writefits, 'altitude.fits', image
    
END



;-------------------------------------------------------------------------------------------------------------------
FUNCTION alt_poly, X

  COMMON share
  
  C = 0.
  
  T3D, /RESET
  T3D, rotate=[-90d0, -90d0, 0d0]
  T3D, rotate=[X[0], 0d0, 0d0]
  T3D, rotate=[0d0, X[1], 0d0]
  T3D, rotate=[0d0, 0d0, X[2]]
  
  obsVertex = Vert_T3D( barr )
      
  t0 = 1 / sqrt(obsVertex[0, *]^2 / x[3]^2 + obsVertex[1, *]^2 / x[4]^2 + obsVertex[2, *]^2 / x[5]^2)
  
  D = (1 - t0) * sqrt(obsVertex[0, *]^2 + obsVertex[1, *]^2 + obsVertex[2, *]^2)
  
  return, D
  
  
END


;-------------------------------------------------------------------------------------------------------------------
FUNCTION alt_vert, X

  COMMON share
  
  C = 0.
  
  T3D, /RESET
  T3D, rotate=[-90d0, -90d0, 0d0]
  T3D, rotate=[X[0], 0d0, 0d0]
  T3D, rotate=[0d0, X[1], 0d0]
  T3D, rotate=[0d0, 0d0, X[2]]
  
  obsVertex = Vert_T3D( objVertex )
  
  t0 = 1 / sqrt(obsVertex[0, *]^2 / x[3]^2 + obsVertex[1, *]^2 / x[4]^2 + obsVertex[2, *]^2 / x[5]^2)
  
  D = (1 - t0) * sqrt(obsVertex[0, *]^2 + obsVertex[1, *]^2 + obsVertex[2, *]^2)
  
  return, D
  
  
END


;----------------------------------------------------------------------------------------------------------------------------
FUNCTION func, X
  
  COMMON share
  
  C = 0.
  
  T3D, /RESET
  T3D, rotate=[-90d0, -90d0, 0d0]
  T3D, rotate=[X[0], 0d0, 0d0]
  T3D, rotate=[0d0, X[1], 0d0]
  T3D, rotate=[0d0, 0d0, X[2]]
  
  obsVertex = Vert_T3D( barr )
  
  t0 = 1 / sqrt(obsVertex[0, *]^2 / x[3]^2 + obsVertex[1, *]^2 / x[4]^2 + obsVertex[2, *]^2 / x[5]^2)
  
  C = TOTAL((1 - t0)^2 * (obsVertex[0, *]^2 + obsVertex[1, *]^2 + obsVertex[2, *]^2))
  
  return, C
  
END