; --------------------------------------------------------------------------------------------- FUNCTION Inertial, Spherdat, Longitude, Latitude, N, R1, Re ;IF ndim NE 2 THEN BEGIN ; Message, 'Array must be three-dimensional. Returning...', /Informational ; RETURN, -1 ;ENDIF s = Size(array, /Dimensions) Ndim = n_elements(Longitude) Mdim = n_elements(Latitude) Result = fltarr(Ndim) ; Calcul du tenseur d'inertie ; --------------------------- for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^(N+5) * cos(Latitude)^3 /(N+5.) Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Iz = int_tabulated(Longitude,Result,/sort) for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^(N+5) * (1.- cos(Longitude[i])^2 * cos(Latitude)^2)*cos(Latitude) /(N+5.) Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Ix = int_tabulated(Longitude,Result,/sort) for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^(N+5) * (1. - sin(Longitude[i])^2 * cos(Latitude)^2)*cos(Latitude) /(N+5.) Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Iy = int_tabulated(Longitude,Result,/sort) for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^(N+5) * cos(Latitude)^3 /(N+5.) Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Ixy = int_tabulated(Longitude,Result*sin(Longitude)*cos(Longitude),/sort) for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^(N+5) * cos(Latitude)^2 * sin(Latitude) /(N+5.) Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Ixz = int_tabulated(Longitude,Result*cos(Longitude),/sort) for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^(N+5) * cos(Latitude)^2 * sin(Latitude) /(N+5.) Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Iyz = int_tabulated(Longitude,Result*sin(Longitude),/sort) IF (R1 ne 0) then begin Ix = Ix * (1/Re^N) Iy = Iy * (1/Re^N) Iz = Iz * (1/Re^N) Ixy = Ixy * (1/Re^N) Ixz = Ixz * (1/Re^N) Iyz = Iyz * (1/Re^N) ENDIF RETURN, [Ix, Iy, Iz, Ixy, Ixz, Iyz] END ; --------------------------------------------------------------------------------------------- FUNCTION Centre_of_masses, Spherdat, Longitude, Latitude Ndim = n_elements(Longitude) Mdim = n_elements(Latitude) ; Calcul du volume de l'objet defini en coordonnees spheriques ; ----------------------------------------------------------- Result = fltarr(Ndim) for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^3 * cos(Latitude)/3. Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Volume = int_tabulated(Longitude,Result,/sort) ; Calcul du centre d'inertie ; -------------------------- for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^4 * cos(Latitude)^2 /4. Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Xg = int_tabulated(Longitude,cos(Longitude)*Result,/sort)/Volume for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^4 * cos(Latitude)^2 /4. Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Yg = int_tabulated(Longitude,sin(Longitude)*Result,/sort)/Volume for i=0L,Ndim-1 do begin Fint = Spherdat[i,indgen(Mdim)]^4 * sin(2.*Latitude) /8. Result[i] = int_tabulated(Latitude,Fint,/sort) endfor Zg = int_tabulated(Longitude,Result,/sort)/Volume Return,[Volume, Xg, Yg, Zg] END ;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ PRO create_shape, File, grid=grid, Print=Print, Radius=Radius, save=save, mikko=mikko, astrodata=astrodata, $ centre=centre, equilibrage=equilibrage, raw=raw, nutmax=nutmax, novisu=novisu, shift=shift, $ OutFile=OutFile, decimation=decimation ; ---------------------------------------------------------------------------------- ; Procedure de conversion d'un fichier de coordonnees cartesiennes ; des sommets d'un objet polyghedrique en un objet defini sur une grille spherique ; ; EXEMPLE: ; create_shape,'eugenia.vtx',grid=5,/print ; IDL> create_shape,'ceres_opnav5_256.txt',grid=1,/mik,/print ; IDL> create_shape,'berbericia.ver',grid=5,/print,/mikko,/save ; ; IDL> create_shape,'shape.dat',grid=2,/print,/astrodata,/save ; ; Juillet 2005 ; Février 2006 ; Octobre 2006 ; Juillet 2007 ; Novembre 2007 ; septembre 2012 (option /shift pour deplacement du centre selon l'axe longitudinal figure ; ---------------------------------------------------------------------------------- DEVICE, DECOMPOSED=0, RETAIN=2 on_ioerror, FichierNonTrouve kequil = 0 tol = 2 tol2 = 0.005 ;; Handle error conditions gracefully catcherror = 0 catch_msg = 'in CREATE_SHAPE.PRO' ;Catch, catcherror If (catcherror NE 0) Then Begin catch, /cancel message, 'Error detected while '+catch_msg+':', /info message, !err_string, /info tol2 = 0.05 EndIf catcherror = 0 If Keyword_set(nutmax) Then tol = nutmax ; Lecture fichier If Keyword_set(Mikko) OR Keyword_set(astrodata) Then Begin openr, lun, File,/GET_LUN readf, lun, header line = STRSPLIT(header, /EXTRACT) Nvertex = LONG64(line[0]) XYZVertex = FLTARR(3, Nvertex) For i=0L, Nvertex-1 do begin readf, lun, line CoordVertex = STRSPLIT(line, /EXTRACT, COUNT=COunt) If (COUNT EQ 3) Then Begin XYZVertex[0,i] = CoordVertex[0] XYZVertex[1,i] = CoordVertex[1] XYZVertex[2,i] = CoordVertex[2] EndIf If (COUNT EQ 4) Then Begin XYZVertex[0,i] = CoordVertex[1] XYZVertex[1,i] = CoordVertex[2] XYZVertex[2,i] = CoordVertex[3] EndIf EndFor close, lun free_lun, lun EndIf Else Begin Objet = file_basename(File,'.vtx') Print, 'Lecture du fichier '+File+' ...' If not keyword_set(decimation) then decimation = 1 If decimation NE 1 then Print, 'WARNING ! : on ne conserve que 1 point sur ',STRING(decimation, format='(I3)') result = FILE_TEST('/usr/local/exelis/user_contrib/IMCCE/descamps/RotAster/AsteRot/shape_template.sav') restore, '/usr/local/exelis/user_contrib/IMCCE/descamps/RotAster/AsteRot/shape_template.sav' XYZVtx = read_ascii(File, count=count, template=sTemplate) print, 'Nbre de donnees : ',count count2 = count/decimation Ind = indgen(count2)*decimation XYZVertex = fltarr(3,LONG(count2)) XYZVertex[0,*] = XYZVtx.X[Ind] XYZVertex[1,*] = XYZVtx.Y[Ind] XYZVertex[2,*] = XYZVtx.Z[Ind] Rvecteur = SQRT(XYZVertex[0,*]^2 + XYZVertex[1,*]^2 + XYZVertex[2,*]^2) Latitude = ASIN(XYZVertex[2,*]/Rvecteur)/!dtor sortlat = SORT(Latitude) rayon_vecteur_moyen = MEAN(Rvecteur) Print, 'Rayon vecteur moyen : ',rayon_vecteur_moyen PLOT, Latitude[sortlat], Rvecteur[sortLat], yrange = [rayon_vecteur_moyen-20, rayon_vecteur_moyen+20], psym=3, $ XTITLE='Latitude', YTITLE='Rayon vecteur' wait, 3 If Keyword_set(shift) Then XYZVertex[0,*] += shift ; 0.46 par exemple pour HW1 0.495 ;XYZVertex = read_asc2(File, npoints=count) ;print, 'Nbre de donnees : ',count EndElse ;AxeX = WHERE(XYZVertex[1,*] LE 0.001 AND XYZVertex[0,*] GT 0.0 AND $ ; XYZVertex[1,*] GE 0.0 , Found) nbiter=0 POINTDENTREE: SPHVertex = CV_COORD(FROM_RECT=XYZVertex, /TO_SPHERE, /DEGREES) Nvertex = N_elements(SPHVertex[0,*]) NegLong = WHERE(SPHVertex[0,*] LT 0.0, Found_neg) If (Found_neg NE 0) Then SPHVertex[0,NegLong] = 360.+SPHVertex[0,NegLong] ;Longitude_shift = SPHVertex[0, AxeX[0]] ;SPHVertex[0,*] = SPHVertex[0,*] - Longitude_shift ; Tri suivant les longitudes rangées entre 0 et 360 degres SortLong = SORT(SPHVertex[0,*]) SPHVertex = SPHVertex[*,SortLong] ; Tri suivant les latitudes ; Les pbs se posent le plus souvent avec les latitudes et le sous-échantillonage jusqu'aux poles SortLat = SORT(SPHVertex[1,*]) SPHVertex = FLOAT(SPHVertex[*,SortLat]) Lon = DBLARR(Nvertex) Lat = DBLARR(Nvertex) Fver = DBLARR(Nvertex) Lon[*] = SPHVertex[0,*] Lat[*] = SPHVertex[1,*] Fver[*] = SPHVertex[2,*] rayon_vecteur_moyen = MEAN(Fver) Print, 'Rayon vecteur moyen : ',rayon_vecteur_moyen PLOT, Lon, Fver, yrange = [rayon_vecteur_moyen-20, rayon_vecteur_moyen+20], psym=3, $ XTITLE='Longitude', YTITLE='Rayon vecteur' wait, 1 PLOT, Lat, Fver, yrange = [rayon_vecteur_moyen-20, rayon_vecteur_moyen+20], psym=3, $ XTITLE='Latitude', YTITLE='Rayon vecteur' wait, 3 LatSmall = WHERE(ABS(Lat) LT 0.1, NSmall) If NSmall NE 0 Then Lat[LatSmall] = 0.0 ; SURFACE, Fver, Lon, Lat ; ------------------------------------------------------------------ ; Set the grid spacing. If Not Keyword_set(grid) Then Grid = 5 GridSpace = [Grid, Grid] ;griddedData = SPH_SCAT_OBJ(Lon, Lat, Fver, BOUNDS=[0, -90, 360, 90], GS=GridSpace) griddedData = SPH_SCAT(Lon, Lat, Fver, BOUNDS=[0, -90, 360, 90], GS=GridSpace) ; Traitement des points aberrants en rayon vecteur IGrid =WHERE(ABS(griddedData) GE MEAN(griddedData)+ 3.*STDDEV(griddedData), foundGrid) IGrid2 =WHERE(griddedData LT MEAN(griddedData), foundGridNegatif) IGrid2 =WHERE(griddedData LT 0, foundGridNegatif) If FoundGridNegatif NE 0 Then Begin ;griddedData[IGrid2] = MEAN(griddedData) EndIf If FoundGrid NE 0 Then Begin ;griddedData[IGrid] = MEAN(griddedData) EndIf ; Visu de la surface dépliée If Not Keyword_set(NoVisu) Then Begin If Not Keyword_set(Mikko) Then WINDOW, 0, TITLE=file_basename(File,'.vtx') If Keyword_set(Mikko) Then WINDOW, 0, TITLE=file_basename(File,'.ver') SURFACE, griddedData WSHOW, 0 EndIf ;-------------------------------------------------- ; Mise en forme des donnees pour impression finale ;-------------------------------------------------- Nlat = (180)/Grid+1 Nlon = (360)/Grid+1 LatVector = -90 + Grid*INDGEN(Nlat) LongVector = Grid*INDGEN(Nlon) ; ----------------------------------- ; Calcul du Volume Result = Centre_of_masses((griddedData), LongVector*!DTOR, LatVector*!DTOR) Volume = Result[0] ; Rayon equivalent: Re Re = ((3./4.)*Volume/!PI)^(1./3.) Xg = Result[1] Yg = Result[2] Zg = Result[3] Print, 'Centre de gravite : ', Xg, Yg, Zg Print, 'Volume : ', Volume Factor = 1.0 If Keyword_set(Radius) Then Factor = Radius/Re Print, 'Rayon equivalent : ', Re*Factor ; Calcul du tenseur d'inertie ; --------------------------- Result = Inertial((griddedData), LongVector*!DTOR, LatVector*!DTOR, 0, 2, Re) Ix = Result[0] Iy = Result[1] Iz = Result[2] Ixy = Result[3] Ixz = Result[4] Iyz = Result[5] Tenseur_inertie = fltarr(3,3) Tenseur_inertie[0,0] = Ix Tenseur_inertie[1,1] = Iy Tenseur_inertie[2,2] = Iz Tenseur_inertie[0,1] = Ixy Tenseur_inertie[0,2] = Ixz Tenseur_inertie[1,0] = Ixy Tenseur_inertie[1,2] = Iyz Tenseur_inertie[2,0] = Ixz Tenseur_inertie[2,1] = Iyz Print, 'Lambda : ', (MAX([Ix, Iy, Iz])/0.4)/(Volume*Re^2) print, 'Tenseur d''inertie: Ix = ',Ix, ' = ',Ix /(Volume*Re^2),' MRe^2' print,' Iy = ',Iy, ' = ',Iy /(Volume*Re^2),' MRe^2' print,' Iz = ',Iz, ' = ',Iz /(Volume*Re^2),' MRe^2' print,' Ixy = ',Ixy, ' = ',Ixy /(Volume*Re^2),' MRe^2' print,' Ixz = ',Ixz, ' = ',Ixz /(Volume*Re^2),' MRe^2' print,' Iyz = ',Iyz, ' = ',Iyz /(Volume*Re^2),' MRe^2' print Print, '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++' Print nbiter += 1 ;If nbiter EQ 3 Then goto, FILEWRITE ; Calcul de l'aplatissement dynamique: J2 ; -------------------------------------- J2 = -(1/(2.*Volume*Re^2)) * (Ix + Iy - 2.*Iz) Print, 'J2 = ', J2 ; Controle equilibrage: il faut que l'on ait Ix < Iy < Iz ; pour avoir une rotation stable autour de l'axe de plus grande Inertie ; qui est l'axe Oz par definition. AxeInertie = [Ix, Iy, Iz] SortAxes = SORT(AxeInertie) If ((SortAxes[2] NE SortAxes[1]+1) OR (SortAxes[1] NE SortAxes[0]+1) AND Not keyword_set(raw)) Then Begin Print, 'Moments d''inertie principaux :', Ix/(Volume*Re^2), Iy/(Volume*Re^2), Iz/(Volume*Re^2) Print, 'Rotation instable ==> Changement de l''ordre des axes', SortAxes XYZSommets = XYZVertex XYZVertex[0,*] = XYZSommets[SortAxes[0],*] XYZVertex[1,*] = XYZSommets[SortAxes[1],*] XYZVertex[2,*] = XYZSommets[SortAxes[2],*] Erase = temporary(XYZSommets) NewTenseur_Inertie = Tenseur_inertie NewTenseur_Inertie[0,0] = Tenseur_inertie[SortAxes[0],SortAxes[0]] NewTenseur_Inertie[1,1] = Tenseur_inertie[SortAxes[1],SortAxes[1]] NewTenseur_Inertie[2,2] = Tenseur_inertie[SortAxes[2],SortAxes[2]] NewTenseur_Inertie[0,1] = Tenseur_inertie[SortAxes[0],SortAxes[1]] NewTenseur_Inertie[0,2] = Tenseur_inertie[SortAxes[0],SortAxes[2]] NewTenseur_Inertie[1,0] = Tenseur_inertie[SortAxes[1],SortAxes[0]] NewTenseur_Inertie[1,2] = Tenseur_inertie[SortAxes[1],SortAxes[2]] NewTenseur_Inertie[2,0] = Tenseur_inertie[SortAxes[2],SortAxes[0]] NewTenseur_Inertie[2,1] = Tenseur_inertie[SortAxes[2],SortAxes[1]] Tenseur_inertie = NewTenseur_Inertie Xg = 0.0 & Yg = 0.0 & Zg = 0.0 goto, POINTDENTREE EndIf If Keyword_set(equilibrage) Then Begin ; Axe X: axe de figure --> Ix minimum ; Axe Y: axe moyen --> Iy ; Axe Z: axe de rotation --> Iz maximum mom_inertie = eigenql(Tenseur_inertie, eigenvectors = eigenvectors, /ascending, /double, /absolute) In1=mom_inertie[0] & I1 = In1/(Volume*Re^2) In2=mom_inertie[1] & I2 = In2/(Volume*Re^2) In3=mom_inertie[2] & I3 = In3/(Volume*Re^2) axes_principaux = eigenvectors Print, 'Valeurs Propres de la Matrice d''inertie : ', In1/(Volume*Re^2), In2/(Volume*Re^2), In3/(Volume*Re^2) Matrix = INVERT(TRANSPOSE(axes_principaux) ) ; Angles d'Euler ; A --> Phi : rotation propre : ne sera donc pas utilise ; B --> Theta : nutation ; C --> Psi : precession MTOEA, Matrix,A,B,C T3D, /RESET T3d, Rotate=[0,0,-C/!DTOR] T3d, Rotate=[B/!DTOR,0,0] T3d, Rotate=[0,0,-A/!DTOR] Print, 'Rotation par angles d''Euler pour equilibrage autour de Z (axe vertical) , Precession : ', -C/!DTOR Print, 'Rotation par angles d''Euler pour equilibrage autour de X (axe horizontal), Nutation : ', -B/!DTOR Print, 'Rotation par angles d''Euler pour equilibrage autour de Z (axe vertical) , Rotation : ', -A/!DTOR XYZVertex = Vert_T3d(XYZVertex) T3D, /RESET ;If (ABS(B/!DTOR) LT tol) Then Begin ; Erase = temporary(equilibrage) ; goto, SORTIE ; EndIf If (ABS(Ix-In1)/(Volume*Re^2) LE tol2 AND ABS(Iy-In2)/(Volume*Re^2) LE tol2 AND ABS (Iz-In3)/(Volume*Re^2) LE tol2) Then Begin Erase = temporary(equilibrage) goto, SORTIE EndIf kequil += 1 Print, 'ITERATION : ', kequil Xg = 0.0 & Yg = 0.0 & Zg = 0.0 goto, POINTDENTREE EndIf If Keyword_set(centre) Then Begin XYZVertex[0,*] = XYZVertex[0,*]-Xg XYZVertex[1,*] = XYZVertex[1,*]-Yg XYZVertex[2,*] = XYZVertex[2,*]-Zg RemoveKey = TEMPORARY(centre) GOTO, POINTDENTREE EndIf SORTIE: ; Création du fichier .tab ; --------------------------------------- If Keyword_set(Print) OR Keyword_set(save) Then Begin If Keyword_set(shift_enattente) Then Begin SPH_COORD = FLTARR(3, (Nlon)*(Nlat)) k = 0 FOR i=0, Nlat-1 Do BEGIN For j=0, Nlon-1 Do Begin SPH_COORD[0,k] = LongVector[j] SPH_COORD[1,k] = LatVector[i] SPH_COORD[2,k] = griddedData[j,i] k = k+1 EndFor EndFor XYZGridded = CV_COORD(FROM_SPHERE=SPH_COORD, /TO_RECT, /DEGREES) XYZGridded[0,*] -= shift SPH_COORD = CV_COORD(FROM_RECT=XYZGridded, /TO_SPHERE, /DEGREES) k = 0 FOR i=0, Nlat-1 Do BEGIN For j=0, Nlon-1 Do Begin griddedData[j,i] = SPH_COORD[2,k] print, SPH_COORD[1,k], LatVector[i], SPH_COORD[0,k]+180., LongVector[j] k = k + 1 EndFor EndFor wait, 1 SURFACE, griddedData WSHOW, 0 EndIf FILEWRITE: FileToWrite = file_basename(File,'.vtx') + '.tab' If Keyword_set(Mikko) Then FileToWrite = file_basename(File,'.ver') + '.tab' If Keyword_set(astrodata) Then FileToWrite = file_basename(File,'.dat') + '.tab' Test_file = FILE_TEST(FileToWrite) ;If (Test_file EQ 1) Then Result = FILE_DELETE(FileToWrite) GET_LUN, Unit openw, unit, FileToWrite, error = err if (err ne 0) then printf, -2, !error_state.msg FOR i=0, Nlat-1 Do BEGIN For j=0, Nlon-1 Do Begin If Keyword_set(Print) Then Print, LatVector[i], LongVector[j], griddedData[j,i]*Factor Printf, unit, LatVector[i], LongVector[j], griddedData[j,i]*Factor, $ FORMAT='(F6.2, 5X, F6.2, 5X, F10.5)' EndFor ENDFOR close, unit free_lun, Unit EndIf RETURN FichierNonTrouve: Print, ' !!! Fichier de forme '+FIle+' non trouve !!!' END