; --------------------------------------------------------------------------------------------- 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 convert_obj, File, SMOOTH=SMOOTH, EXTRAPOL=EXTRAPOL, grid=grid, Print=Print, Radius=Radius, save=save, $ center=center, flip=flip, equilibrage=equilibrage, raw=raw, nutmax=nutmax ; ---------------------------------------------------------------------------------- ; Procedure de conversion d'un fichier de coordonnees cartesiennes ; des sommets d'un objet polyghedrique en un objet defini sur une grille spherique ; ; AFAIRE: ; - le raccordement en longitude entre les meridiens 360 degres et 0 degres ; - l'extrapolation correcte au voisinage des poles: effets de bord importants ; - extrapolation de la surface par analyse multiresolution ; - Problemes au niveau du raccordement des poles ; - Raccordement de l'origine au centre de gravité (option /center) ; ; EXEMPLE: ; convert_obj,'eugenia.vtx',grid=5,/print ; ; Juillet 2005 ; Février 2006 ; Octobre 2006 ; Juillet 2007 ; ---------------------------------------------------------------------------------- DEVICE, DECOMPOSED=0, RETAIN=2 Objet = file_basename(File,'.vtx') XYZVertex = read_asc2(File) If Keyword_set(flip) Then XYZVertex[2,*] = -XYZVertex[2,*] Xg = 0.0 & Yg = 0.0 & Zg = 0.0 kequil = 0 tol = 2.0 If Keyword_set(nutmax) Then tol = nutmax LECTURE: XYZVertex[0,*] -= Xg XYZVertex[1,*] -= Yg XYZVertex[2,*] -= Zg 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] ; Tri suivant les longitudes 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 = SPHVertex[*,SortLat] ; On impose de se trouver aux poles pour les points extrêmes de l'échantillonage SPHVertex[1,0] = -90.0 SPHVertex[1,Nvertex-1] = 90.0 ; Echange des colonnes latitude/longitude NewSPHVertex = SPHVertex NewSPHVertex[1,*] = SPHVertex[0,*] NewSPHVertex[0,*] = SPHVertex[1,*] SPHVertex = NewSPHVertex ; Création de sommets supplémentaires entourant la grille initiale, pour interpolation à venir ExtraVertex = FLTARR(3,Nvertex*3) ExtraVertex[0,*] = TRANSPOSE([TRANSPOSE([REVERSE(SPHVertex[0,*], 2)-180.]),TRANSPOSE(SPHVertex[0,*]),$ TRANSPOSE(REVERSE(SPHVertex[0,*], 2)+180.)]) ExtraVertex[1,*] = TRANSPOSE([TRANSPOSE([REVERSE(SPHVertex[1,*], 2)]),TRANSPOSE(SPHVertex[1,*]),$ TRANSPOSE(REVERSE(SPHVertex[1,*], 2))]) ExtraVertex[2,*] = TRANSPOSE([TRANSPOSE([REVERSE(SPHVertex[2,*], 2)]),TRANSPOSE(SPHVertex[2,*]),$ TRANSPOSE(REVERSE(SPHVertex[2,*], 2))]) ExtraVertex2 = FLTARR(3,Nvertex*3) ExtraVertex2[0,*] = TRANSPOSE([TRANSPOSE([REVERSE(SPHVertex[0,*],2)]),TRANSPOSE([REVERSE(SPHVertex[0,*],2)]), $ TRANSPOSE([REVERSE(SPHVertex[0,*],2)])]) ExtraVertex2[1,*] = TRANSPOSE([TRANSPOSE([REVERSE(SPHVertex[1,*],2)+360]),TRANSPOSE([REVERSE(SPHVertex[1,*],2)]), $ TRANSPOSE([REVERSE(SPHVertex[1,*],2)-360])]) ExtraVertex2[2,*] = TRANSPOSE([TRANSPOSE([REVERSE(SPHVertex[2,*],2)]),TRANSPOSE([REVERSE(SPHVertex[2,*],2)]), $ TRANSPOSE([REVERSE(SPHVertex[2,*],2)])]) ; Nouvel objet Erase = temporary(SPHVertex) SPHVertex = ExtraVertex2 Nvertex = Nvertex*3 ; ------------------------------------------------------------------ ; Set the window size. WINDOWSIZE = 400 ; Show the data points randomly distributed. Window, /Free, Title='XY Point Distribution', XSize=WINDOWSIZE, YSize=WINDOWSIZE, XPos=0, YPos=0 Plot, Findgen(Nvertex), /NoData, XRange=[-180, 180], YRange=[-360,700],XSTYLE=1,YSTYLE=1 PlotS, SPHVertex[0,*], SPHVertex[1,*], PSym=1 ; ------------------------------------------------------------------ ; Create the set of Delaunay triangles. The variables ; "triangles" and "boundaryPts" are output variables. Triangulate, SPHVertex[0,*], SPHVertex[1,*], triangles, boundaryPts ; ------------------------------------------------------------------ ; Set the grid spacing. If Not Keyword_set(grid) Then Grid = 5. GridSpace = [Grid, Grid] ; ------------------------------------------------------------------ ; Grid the Z data, using the triangles.The variables "xvector" ; and "yvector" are output variables. Missing = MIN(SPHVertex[2,Where(SPHVertex[2,*] NE 0.)]) ;Missing = 0 GRIDDED: griddedData = TriGrid(SPHVertex[0,*], SPHVertex[1,*], SPHVertex[2,*], triangles, gridSpace, $ XGrid=xvector, YGrid=yvector , EXTRAPOLATE=boundaryPts, MISSING=Missing) ; ------------------------------------------------------------------ ; Conversion sur grille entiere ; Ici il faut que les coordonnes de grille soient des multiples du pas de grille ; compris entre -90 et 90 pour les latitudes Xvector=Fix(Xvector) Yvector=Fix(Yvector) ; Meridien le plus au sud SOuthMer = Where(Xvector EQ MIN(Xvector)) Ilat = Where(Xvector GE -90. AND Xvector LE +90., Nlat) Ilong = Where(Yvector GE 0. AND Yvector LE +360. , Nlon) ; ------------------------------------------------------------------ ; Display the surface. Window, /Free, Title=Objet+' - Initial Surface', $ XSize=WINDOWSIZE*1.5, YSize=WINDOWSIZE, XPos=0, YPos=WINDOWSIZE + 30 ;Surface, griddedData, xvector, yvector, /NoData ;Surface, griddedData, xvector, yvector, XStyle=4, $ ; YStyle=4, ZStyle=0, /NoErase Surface, griddedData[Ilat[0]:Ilat[Nlat-1], Ilong[0]:Ilong[Nlon-1]], xvector[Ilat], yvector[Ilong], /NoData Surface, griddedData[Ilat[0]:Ilat[Nlat-1], Ilong[0]:Ilong[Nlon-1]], xvector[Ilat], yvector[Ilong], XStyle=4, $ YStyle=4, ZStyle=0, /NoErase ; ------------------------------------------------------------------ ; Display the smoothed surface. If Keyword_set(smooth) Then Begin griddedData = TriGrid(SPHVertex[0,*], SPHVertex[1,*], SPHVertex[2,*], triangles, gridSpace, /Quintic) Window, /Free, Title='Smoothed Surface', $ XSize=WINDOWSIZE*1.5, YSize=WINDOWSIZE, $ XPos=WINDOWSIZE*1.5 + 10, YPos=WINDOWSIZE + 30 Surface, griddedData[Ilat[0]:Ilat[Nlat-1], Ilong[0]:Ilong[Nlon-1]], xvector[Ilat], yvector[Ilong], /NoData Surface, griddedData[Ilat[0]:Ilat[Nlat-1], Ilong[0]:Ilong[Nlon-1]], xvector[Ilat], yvector[Ilong], $ XStyle=4, YStyle=4, ZStyle=0, /NoErase EndIf ; ------------------------------------------------------------------ ; Display the extrapolated surface. If Keyword_set(extrapol) Then Begin griddedData = TriGrid(SPHVertex[0,*], SPHVertex[1,*], SPHVertex[2,*], triangles, gridSpace,$ /Quintic, Extrapolate=boundaryPts) Window, /Free, Title='Extrapolated Surface', $ XSize=WINDOWSIZE*1.5, YSize=WINDOWSIZE, $ XPos=2*WINDOWSIZE*1.5 + 20, YPos=WINDOWSIZE + 30 Surface, griddedData[Ilat[0]:Ilat[Nlat-1], Ilong[0]:Ilong[Nlon-1]], xvector[Ilat], yvector[Ilong], /NoData Surface, griddedData[Ilat[0]:Ilat[Nlat-1], Ilong[0]:Ilong[Nlon-1]], xvector[Ilat], yvector[Ilong], XStyle=4, $ YStyle=4, ZStyle=0, /NoErase EndIf ;-------------------------------------------------- ; Mise en forme des donnees pour impression finale ;-------------------------------------------------- LatVector = xvector[Ilat] LongVector = yvector[Ilong] GrilleVal = griddedData[Ilat[0]:Ilat[Nlat-1], Ilong[0]:Ilong[Nlon-1]] Nlat = N_elements(LatVector) Nlon = N_elements(LongVector) GOTO, OBSOLETE ; Moyenisation aux poles Stat = MOMENT(GrilleVal[1,*]) RPoleN = Stat[0] Stat = MOMENT(GrilleVal[Nlat-2,*]) RPoleS = Stat[0] ;If RPoleS EQ 0.0 Then RPoleS = RPoleN If RPoleS EQ 0.0 Then Begin Stat = MOMENT(GrilleVal[Nlat-2,*]) MISSING = Stat[0] GOTO, GRIDDED EnDIf If RPoleN EQ 0.0 Then Begin Stat = MOMENT(GrilleVal[1,*]) MISSING = Stat[0] GOTO, GRIDDED EnDIf OBSOLETE: ; GrilleVal[0,*] = RPoleN ; GrilleVal[Nlat-1,*] = RPoleS ; ----------------------------------- ; Calcul du Volume Result = Centre_of_masses(TRANSPOSE(GrilleVal), 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] Factor = 1.0 If Keyword_set(Radius) Then Factor = Radius/Re ; Calcul du tenseur d'inertie ; --------------------------- Result = Inertial(TRANSPOSE(GrilleVal), 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 ; Calcul de l'aplatissement dynamique: J2 ; -------------------------------------- J2 = -(1/(2.*Volume*Re^2)) * (Ix + Iy - 2.*Iz) Lambda = (MAX([Ix, Iy, Iz])/0.4)/(Volume*Re^2) ; 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 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] In2=mom_inertie[1] In3=mom_inertie[2] axes_principaux = eigenvectors Print, 'Valeurs Propres de la Matrice d''inertie : ', In1/(Volume*Re^2), In2/(Volume*Re^2), In3/(Volume*Re^2) ; 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 kequil += 1 Print, 'ITERATION : ', kequil Xg = 0.0 & Yg = 0.0 & Zg = 0.0 goto, LECTURE EndIf SORTIE: Print, '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++' Print, 'Centre de gravite : ', Xg, Yg, Zg Print, 'Rayon equivalent : ', Re*Factor Print, 'Lambda : ', lambda Print, 'J2 = ', J2 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 If Keyword_set(center) Then Begin Print, ' ******** CENTRAGE OBJET ********' erase = temporary(center) Print goto, LECTURE EndIf ; Création du fichier output .tab ; --------------------------------------- If Keyword_set(Print) OR Keyword_set(save) Then Begin FileToWrite = file_basename(File,'.vtx') + '.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 Print, LatVector[i], LongVector[j], GrilleVal[i,j]*Factor Printf, unit, LatVector[i], LongVector[j], GrilleVal[i,j]*Factor, $ FORMAT='(F6.2, 5X, F6.2, 5X, F10.5)' ; If i eq 0 OR i eq Nlat-1 Then Print, LatVector[i], LongVector[j], GrilleVal[i,j]*Factor EndFor ENDFOR close, unit EndIf END