; docformat = 'rst' ; ; NAME: ; spinToSpinIAU ; PURPOSE: ; Convert a spin file in LCI format into IAU format ; ;+ ; :Description: ; Convert a spin file in LCI format into IAU format ; ; :Categories: ; Shape Models ; ; :Params: ; spinFile: in, required, type=string/structure ; The spin structure (see readSpin) or the spin file (LCI format) ; iauFile: in, optional, type=string ; Path to the file to contain the spin in IAU format ; ; :Keywords: ; verbose: in, optional, type=boolean, default=0 ; If set, write IAU spin in stdout ; ; :Examples: ; Convert a spin from LCI:: ; IDL> spinToSpinIAU, 'myspin.spin', 'myspin.iau' ; ; :Uses: ; readSpin ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in July 2017, B. Carry (OCA) ;- pro spinToSpinIAU, spinFile, iauFile, verbose=verbose COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Input Syntax Verification ----------------------------------------------------------- if not keyword_set(spinFile) then begin message, /ioError, 'Syntax: spinToSpinIAU, spinFile [,iauFile, /verbose' return endif ;--I.2-- If Input is a File, is it Readable? ------------------------------------------------- if size( spinFile, /Type ) eq 7 then if not file_test( spinFile ) then begin message, /Informative, 'Spin file cannot be read: '+strTrim(spinFile,2) endif ;--I.3-- If Input is a Structure, Does it contain appropriate keywords? ---------------------- if size( spinFile, /Type ) eq 8 then begin ;to be coded endif ;--I.4-- Defaults Options -------------------------------------------------------------------- if not keyword_set(verbose) then verbose=0 ;--I.5-- Obliquity of Earth ------------------------------------------------------------------ epsilon = ( 23. + 26. / 60. + 21.4119 / 3600. ) * !DTOR jd0=0.D ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- From LCI to IAU Format -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Read Spin if not Provided as Structure --------------------------------------------- if size( spinFile, /Type ) ne 8 then spin=readSpin(spinFile) ;--II.2-- Set Up Individual Rotation Matrix -------------------------------------------------- ;--II.2.1-- Z-axis - Absolute rotation angle = -spin.phi0 R_rot = [ [ cos(angle), sin(angle), 0] ,$ [-sin(angle), cos(angle), 0] ,$ [ 0 , 0 , 1] ] ;--II.2.2-- Y-axis - Pole Declination angle = -spin.lat R_beta = [ [cos(angle), 0, -sin(angle)] ,$ [ 0 , 1, 0 ] ,$ [sin(angle), 0, cos(angle)] ] ;--II.2.3-- Z-axis - Pole Right Ascension angle = -spin.lon R_lambda = [ [ cos(angle), sin(angle), 0] ,$ [-sin(angle), cos(angle), 0] ,$ [ 0 , 0 , 1] ] ;-II.2.4-- X-axis - Earth Obliquity angle = -epsilon R_eps = [ [1, 0 , 0 ] ,$ [0, cos(angle), sin(angle)] ,$ [0, -sin(angle), cos(angle)] ] ;--II.3-- Final Rotation Matrix -------------------------------------------------------------- M = R_eps ## R_lambda ## R_beta ## R_rot ;--II.4-- Defining Directions in Body's Frame ------------------------------------------------- z_body = [0, 0, 1] x_body = [1, 0, 0] ;--II.5-- Rotate Body Defining Directions in Equatorial Frame --------------------------------- z_eq = reform( M ## z_body ) x_eq = reform( M ## x_body ) ;--II.6-- Convert Cartesian Coordinates of Pole Into Spherical -------------------------------- z_sph = cv_coord( from_rect=z_eq, /To_Sphere, /Degrees ) alpha = z_sph[1] mod 360. delta = z_sph[2] ;--II.7-- Absolute Rotation at Reference Epoch ------------------------------------------------ orient = cv_coord( from_sphere=[90. + alpha, 1e-8, 1], /To_Rect ) W = acos( x_eq##transpose(orient) ) /!DTOR if (crossp(orient, x_eq) ## transpose(z_eq)) lt 0 then $ W = 360 - W ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Outputs -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; omega = 2*!PI* (24./ spin.ph) jd0 = spin.jd0 ;-2.8 - Report results print, 'Pole solution according to IAU standards:' ; print, alpha, delta, omega/!DTOR, format='(F5.1,3x,F5.1,3x,F10.5)' ; print, Jd0, W, format='(D14.5,2x,F6.1)' ; print, '' print, frameCoord_ec2eq( spin.lon, spin.lat), omega/!DTOR, format='(F5.1,3x,F5.1,3x,F10.5)' ;--I.3-- Useful Constants -------------------------------------------------------------------- AUtoKM = 149597870.691D ;-km/AU cSI= 299792.4580D ;-km/s cKMperDay = cSI * 86400.0D ;-km/day J2000 = 2451545d iau = 3 coordGeo = -ssoPositions( iau, J2000, type='aster', obs='500' ) jdLTC = j2000 - sqrt( total( coordGeo^2., 1) )*AUtoKM/cKMperDay coordGeo = -ssoPositions( iau, jdLTC, rPlane=2, type='aster', obs='500' ) ;-ECJ2000 sep = subObserver_coord( -coordGeo, spin, jdLTC ) print, sep.lon sep = subObserver_coord( coordGeo, spin, jdLTC ) print, sep.lon coordGeo = -ssoPositions( iau, j2000, rPlane=2, type='aster', obs='500' ) ;-ECJ2000 sep = subObserver_coord( -coordGeo, spin, J2000 ) print, sep.lon sep = subObserver_coord( coordGeo, spin, J2000 ) print, sep.lon ;108 50 1198.413600 ;2451545.0 287.7 stop ; ; ; ; ;;----------------------------------------------------; ;;--- 2 --- Object pole and orientation conversion ---; ;;----------------------------------------------------; ; ; ;-2.1 - Reading the original spin file ; if file_exist( spinFile ) then begin ; ; ;-Open the file and read the information ; openr, 1, spinFile ; readf, 1, lambda, beta, period, jd0, phi0 ; close, 1 ; ; ;-Convert quantities into SPHERICAL (RADIAN) ; lambda *= !DTOR ; beta = ( 90. - beta ) *!DTOR ; omega = 2*!PI* (24./ period) ; phi0 *= !DTOR ; ; endif ; ; ; ;-2.2 - Conversion from BODY to EQUATORIAL reference frame ; ;-2.2.a - Z-axis - Absolute rotation ; angle = -phi0 ; R_rot = [ [ cos(angle), sin(angle), 0] ,$ ; [-sin(angle), cos(angle), 0] ,$ ; [ 0 , 0 , 1] ] ; ; ;-2.2.b - Y-axis - Pole Declination ; angle = -beta ; R_beta = [ [cos(angle), 0, -sin(angle)] ,$ ; [ 0 , 1, 0 ] ,$ ; [sin(angle), 0, cos(angle)] ] ; ; ;-2.2.c - Y-axis - Pole Right Ascension ; angle = -lambda ; R_lambda = [ [ cos(angle), sin(angle), 0] ,$ ; [-sin(angle), cos(angle), 0] ,$ ; [ 0 , 0 , 1] ] ; ; ;-2.2.d - Z-axis - Earth Obliquity ; angle = -epsilon ; R_eps = [ [1, 0 , 0 ] ,$ ; [0, cos(angle), sin(angle)] ,$ ; [0, -sin(angle), cos(angle)] ] ; ; ; ;-2.3 - Rotation matrix ; M = R_eps ## R_lambda ## R_beta ## R_rot ; ; ;-2.4 - Body's pole and origin IN BODY FRAME ; pole_body = [0, 0, 1] ; x_body = [1, 0, 0] ; ; ;-2.5 - Body's pole and origin IN EQUATORIAL FRAME ; pole_eq = reform( M ## pole_body ) ; x_eq = reform( M ## x_body ) ; coordEQ = reform( cart2sph( pole_eq(0), pole_eq(1), pole_eq(2) ) ) ; ; ;-2.6 - Equatorial Pole coordinates ; alpha = coordEQ( 1 ) mod 360. ; delta = coordEQ( 2 ) ; ; ;-2.7 - Absolute rotation at reference time ; orient = reform(sph2cart( 1, 90. + alpha, 0.00000000001D )) ; W = acos( x_eq##transpose(orient) ) /!DTOR ; if (crossp(orient, x_eq) ## transpose(pole_eq)) lt 0 then $ ; W = 360 - W ; ; ; ;-2.8 - Report results ; print, 'Pole solution according to IAU standards:' ; print, alpha, delta, omega/!DTOR, format='(F5.1,3x,F5.1,3x,F10.5)' ; print, Jd0, W, format='(D14.5,2x,F6.1)' ; ; end