;+ ; NAME: ; EulerRotate ; PURPOSE: ; Calculate spherical angles in rotated coordinate system ; CATEGORY: ; gen/idl/toolbox/math ; CALLING SEQUENCE: FUNCTION EulerRotate, ABG, In, In1, In2, degrees=Degrees, rectangular=Rectangular ; INPUTS: ; ABG array[3] or multi-dimensional array[3,..] with first dimension of 3 ; elements. If multi-dimensional, then: ; 1. either the array structure, except for the 1st dimension, is the same ; as for In. Each position in In is rotated using the matching ABG set. ; 2. or In is only 1-dimensional i.e. contains only one position. Then ; each ABG combination is applied to this one position. ; Euler angles [A, B, G] ; A phase angle of the pole Z(new) in old coordinate system ; B polar angle of the pole Z(new) in old coordinate system ; G phase angle of X(new) in coordinate system (2) ; ; /rectangular NOT SET: ; In any multi-dimensional array with first dimension of 2 or 3 elements. ; array[0,..] is the longitude and array[1,..] the latitude; ; array[3,..] (usually the radial distance) is not modified. ; /rectangular SET: ; In array[3,..], any multi-dimensional array with array[i,..] (i=0,2) ; the x,y and z-components ; OPTIONAL INPUTS: ; /rectangular if set, the input array In is interpreted as vectors in rectangular ; coordinates. ; if not set, the input array In is assumed to be in spherical coordinates ; (longitude, latitude, radial distance). ; /degrees if set, input and output is in degrees ; OUTPUTS: ; /rectangular NOT SET: ; double array with same structure as input array In ; array[0:1,..] are longitude and latitude in new coordinate system ; (0<=Lng<360, -90<=Lat<=90) ; /rectangular SET: ; array[3,..] with same structure as input array In ; x,y,z coordinates in new coordinate system ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; cvt3d, SuperArray, SubArray, SyncDims, ToRadians ; PROCEDURE: ; > /rectangular NOT SET: ; Given longitude (phase angle) and latitude (90-polar angle) of a point ; in some coordinate system [axes: X(old),Y(old),Z(old)], angles are ; calculated in the new coordinate system [axes: X(new),Y(new),Z(new)] ; rotated with respect to the original one over the Euler angles Alfa, ; Beta, Gamma. ; The rotation from old to new system is realized by ; (1) rotation around Z(old) over Alfa ---> X(1),Y(1),Z(1)=Z(old) ; (2) rotation around Y(1) over Beta ---> X(2),Y(2)=Y(1),Z(2)=Z(new) ; (3) rotation around Z(new) over Gamma ---> X(new),Y(new),Z(new). ; > /rectangular SET: ; The 3D rotation !p.t is used do do the transformation. ; MODIFICATION HISTORY: ; 1989, Paul Hick (MPAE,UCSD) ; FEB-1998, Paul Hick (UCSD/CASS), added call to cvt3d to handle conversion ; of rectangular coordinates ; JUN-2000, Paul Hick (UCSD/CASS), added option to rotate ; single position using multiple Euler angle combinations. ; SEP-2001, Paul Hick (UCSD/CASS), ; removed a restriction on ABG for a transformation of rectangular ; coordinates. Only one set of Euler angles was allowed. Now separate ; Euler angles can be specified for each vector in 'In'; these are ; processed using an explicit do-loop). ; AUG-2004, Paul Hick (UCSD/CASS) ; For rotations of spherical coordinates (/Rectangular NOT set) a safety ; check cos(angle) > (-1) made bad coordinate values disappear ; (!values.f_nan > 1 returns 1, while !values.f_nan < 1 returns !values.f_nan). ; The safety check now explicitly checks for bad values, so they don't ; get lost. ; DEC-2005, Paul Hick (UCSD/CASS, pphick@ucsd.edu) ; Bug fix. When 1-element arrays are input at In and In1, the output would ; be a scalar. The output should always have the same structure as the input. ;- nparams = n_params() CASE nparams OF 3: BEGIN ; Longitude and latitude as separate arguments sz0 = size(In ) sz1 = size(In1) In = SuperArray(In,2,/lead) In = SubArray(In,element=1,/clear,add=In1) END 4: BEGIN ; X,Y,Z components as separate arguments In = SuperArray(In,3,/lead) In = SubArray(In,element=1,/clear,add=In1) In = SubArray(In,element=2,/clear,add=In2) END ELSE: ENDCASE Out = In szOut = size(Out) ; Save input array structure (szIn[1]=2 or 3) szABG = size(ABG) ; KLUDGE: Check for special case: ; single set of coordinates with multiple Euler angle combinations nABG = szABG[szABG[0]+2]/szABG[1] ; # Euler angle sets (szABG[1]=3) IF nABG GT 1 THEN ABG = reform(ABG, szABG[1], nABG, /overwrite) ; Change to 2D array [3,nABG] IF szOut[szOut[0]+2] LE 3 AND nABG GT 1 THEN BEGIN Out = SuperArray(Out, nABG, /trail) Out = reform(Out, [szOut[1:szOut[0]], szABG[2:szABG[0]]], /overwrite) szOut = size(Out) ; Update array structure ENDIF ; Change to 2D-array Out = reform(Out, szOut[1], szOut[szOut[0]+2]/szOut[1], /overwrite) IF keyword_set(Rectangular) THEN BEGIN ; If there is more than one set of Euler angles we need an explicit do-loop ; We assume that there is a separate set of Euler angles ABG for each vector in Out. CASE nABG OF 1: Out = cvt3d(mat=cvt3d(rot=[0,0,ABG[0]],degrees=Degrees), $ rot=[0,ABG[1],ABG[2]],degrees=Degrees,vector=Out) ELSE: $ FOR i=0L,nABG-1 DO $ Out[*,i] = cvt3d(mat=cvt3d(rot=[0,0,ABG[0,i]],degrees=Degrees),$ rot=[0,ABG[1,i],ABG[2,i]],degrees=Degrees,vector=Out[*,i]) ENDCASE ENDIF ELSE BEGIN rpm = ToRadians(degrees=Degrees) HalfPi = !dpi/2 TwoPi = 2*!dpi Out = double(Out) ; Switch to double precision Out[0:1,*] = Out[0:1,*]*rpm ; Longitude, latitude (radians) ;nABG = szABG[szABG[0]+2]/szABG[1] ; # Euler angle sets (szABG[1]=3) ;if nABG gt 1 then ABG = reform(ABG, szABG[1], nABG) ; Reform to 2D array leaving 1st dimension the same IF nABG EQ 1 THEN Angle = ABG[0] ELSE Angle = ABG[0,*] Out[0,*] = Out[0,*]-double(Angle)*rpm ; Rotation (1) IF nABG EQ 1 THEN Angle = ABG[1] ELSE Angle = ABG[1,*] SinB = sin(double(Angle)*rpm) CosB = cos(double(Angle)*rpm) SinL = sin(Out[1,*]) CosL = cos(Out[1,*]) SinP = sin(Out[0,*]) CosP = cos(Out[0,*]) ; Rotation (2) CosT = SinL*CosB+CosL*SinB*CosP ; cos(polar angle) CosT = CosT < 1 ; Safety belt ; Note that BadValue(0.0) > (-1) is -1, NOT BadValue(0.0) ; Since we don't want to loose bad values we need to check ; for them explicitly. A = where(finite(CosT)) IF A[0] NE -1 THEN CosT[A] = CosT[A] > (-1) ; Safety belt Out[1,*] = HalfPi-acos(CosT) ; acos returns between 0 and Pi CosT = abs(CosT) ; = 1 on new Z-axis A = where(CosT GE 1.d0) IF A[0] NE -1 THEN Out[0,A] = 0 ; Point on new Z-axis; Lng not defined A = where(CosT LT 1.d0) IF A[0] NE -1 THEN BEGIN CosL = CosL[A] SinP = CosL*SinP[A] CosP = CosL*CosB[A]*CosP[A]-SinL[A]*SinB[A] Out[0,A] = atan(SinP,CosP) IF nABG EQ 1 THEN Angle = ABG[2] ELSE Angle = ABG[2,A] Out[0,A] = Out[0,A]-double(Angle)*rpm ;Rotation (3) ENDIF Out[0,*] = Out[0,*] mod TwoPi ; -360