--------------------------------------------------------
rotm_coo function Generate rotation matrix for
converting J2000.0 coordinates to
galactic l & b coordinates.
Input : - Type of rotation matrix
'e' - equatorial with mean equinox of J2000.0
to ecliptic with mean ecliptic and equinox J2000.0
'E' - ecliptic with mean ecliptic and equinox J2000.0
to equatorial with mean equinox of J2000.0
'g' - galactic to equatorial with mean equinox of J2000.0
'G' - equatorial with mean equinox of J2000.0 to galactic.
'p' - precession matrix from mean equinox
of date to mean equinox of J2000.0.
'P' - precession matrix from mean equinox
J2000.0 to mean equinox of date.
'pd' - precession matrix from true equinox
of date to mean equinox of J2000.0.
'Pd' - precession matrix from mean equinox
J2000.0 to true equinox of date.
- Equinox of coordinates (in Julian Day),
used only in the case of 'p' | 'P' | 'pd' | 'Pd' | 'ed' | 'Ed'
In case of 'E' or 'q' if this parameter is
not given it is taken as 2451545.0 (=J2000.0)
Output : - rotation matrix
Reference : Ex. Supp. to the Astronomical Almanac.
Tested : Matlab 5.3
By : Eran O. Ofek August 1999 / June 2000
URL : http://wise-obs.tau.ac.il/~eran/matlab.html
--------------------------------------------------------0001 function RotM=rotm_coo(Type,EquinoxJD); 0002 %-------------------------------------------------------- 0003 % rotm_coo function Generate rotation matrix for 0004 % converting J2000.0 coordinates to 0005 % galactic l & b coordinates. 0006 % Input : - Type of rotation matrix 0007 % 'e' - equatorial with mean equinox of J2000.0 0008 % to ecliptic with mean ecliptic and equinox J2000.0 0009 % 'E' - ecliptic with mean ecliptic and equinox J2000.0 0010 % to equatorial with mean equinox of J2000.0 0011 % 'g' - galactic to equatorial with mean equinox of J2000.0 0012 % 'G' - equatorial with mean equinox of J2000.0 to galactic. 0013 % 'p' - precession matrix from mean equinox 0014 % of date to mean equinox of J2000.0. 0015 % 'P' - precession matrix from mean equinox 0016 % J2000.0 to mean equinox of date. 0017 % 'pd' - precession matrix from true equinox 0018 % of date to mean equinox of J2000.0. 0019 % 'Pd' - precession matrix from mean equinox 0020 % J2000.0 to true equinox of date. 0021 % - Equinox of coordinates (in Julian Day), 0022 % used only in the case of 'p' | 'P' | 'pd' | 'Pd' | 'ed' | 'Ed' 0023 % In case of 'E' or 'q' if this parameter is 0024 % not given it is taken as 2451545.0 (=J2000.0) 0025 % Output : - rotation matrix 0026 % Reference : Ex. Supp. to the Astronomical Almanac. 0027 % Tested : Matlab 5.3 0028 % By : Eran O. Ofek August 1999 / June 2000 0029 % URL : http://wise-obs.tau.ac.il/~eran/matlab.html 0030 %-------------------------------------------------------- 0031 RADIAN = 180./pi; 0032 J2000 = 2451545.5; 0033 0034 switch Type 0035 case {'e'} 0036 Obl = obliquity(J2000); 0037 RotM = [1 0 0; 0 cos(Obl) sin(Obl); 0 -sin(Obl) cos(Obl)]; 0038 %RotM = [-0.0548755604 +0.4941094279 -0.8676661490; -0.8734370902 -0.4448296300 -0.1980763734; -0.4838350155 +0.7469822445 +0.4559837762]; 0039 case {'E'} 0040 Obl = obliquity(J2000); 0041 RotM = [1 0 0; 0 cos(Obl) sin(Obl); 0 -sin(Obl) cos(Obl)].'; 0042 %RotM = [-0.0548755604 +0.4941094279 -0.8676661490; -0.8734370902 -0.4448296300 -0.1980763734; -0.4838350155 +0.7469822445 +0.4559837762]'; 0043 case {'g'} 0044 RotM = [-0.0548755604 +0.4941094279 -0.8676661490; -0.8734370902 -0.4448296300 -0.1980763734; -0.4838350155 +0.7469822445 +0.4559837762]'; 0045 case {'G'} 0046 RotM = [-0.0548755604 +0.4941094279 -0.8676661490; -0.8734370902 -0.4448296300 -0.1980763734; -0.4838350155 +0.7469822445 +0.4559837762]; 0047 case {'p','P','pd','Pd'} 0048 T = (EquinoxJD - 2451545.0)./36525.0; 0049 0050 ZetaA = 0.6406161.*T + 0.0000839.*T.*T + 0.0000050.*T.*T.*T; 0051 ZA = 0.6406161.*T + 0.0003041.*T.*T + 0.0000051.*T.*T.*T; 0052 ThetaA = 0.5567530.*T - 0.0001185.*T.*T - 0.0000116.*T.*T.*T; 0053 ZetaA = ZetaA./RADIAN; 0054 ZA = ZA./RADIAN; 0055 ThetaA = ThetaA./RADIAN; 0056 0057 RotM = zeros(3,3); 0058 RotM(1,1) = cos(ZetaA).*cos(ThetaA).*cos(ZA) - sin(ZetaA).*sin(ZA); 0059 RotM(2,1) = cos(ZetaA).*cos(ThetaA).*sin(ZA) + sin(ZetaA).*cos(ZA); 0060 RotM(3,1) = cos(ZetaA).*sin(ThetaA); 0061 RotM(1,2) =-sin(ZetaA).*cos(ThetaA).*cos(ZA) - cos(ZetaA).*sin(ZA); 0062 RotM(2,2) =-sin(ZetaA).*cos(ThetaA).*sin(ZA) + cos(ZetaA).*cos(ZA); 0063 RotM(3,2) =-sin(ZetaA).*sin(ThetaA); 0064 RotM(1,3) =-sin(ThetaA).*cos(ZA); 0065 RotM(2,3) =-sin(ThetaA).*sin(ZA); 0066 RotM(3,3) = cos(ThetaA); 0067 0068 switch Type 0069 case {'p'} 0070 RotM = RotM'; 0071 case {'P'} 0072 RotM = RotM; 0073 case {'pd'} 0074 % calculate nutation matrix 0075 [Nut, NutMat]=nutation(EquinoxJD); 0076 RotM = [NutMat']*[RotM']; 0077 case {'Pd'} 0078 % calculate nutation matrix 0079 [Nut, NutMat]=nutation(EquinoxJD); 0080 RotM = NutMat*RotM; 0081 otherwise 0082 error('Unknown rotation matrix type'); 0083 end 0084 0085 0086 otherwise 0087 error('Unknown rotation matrix type'); 0088 end