1!----------------------------------------------------------------------------- 2! 3! Copyright (C) 1997-2013 Krzysztof M. Gorski, Eric Hivon, 4! Benjamin D. Wandelt, Anthony J. Banday, 5! Matthias Bartelmann, Hans K. Eriksen, 6! Frode K. Hansen, Martin Reinecke 7! 8! 9! This file is part of HEALPix. 10! 11! HEALPix is free software; you can redistribute it and/or modify 12! it under the terms of the GNU General Public License as published by 13! the Free Software Foundation; either version 2 of the License, or 14! (at your option) any later version. 15! 16! HEALPix is distributed in the hope that it will be useful, 17! but WITHOUT ANY WARRANTY; without even the implied warranty of 18! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19! GNU General Public License for more details. 20! 21! You should have received a copy of the GNU General Public License 22! along with HEALPix; if not, write to the Free Software 23! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 24! 25! For more information about HEALPix see http://healpix.sourceforge.net 26! 27!----------------------------------------------------------------------------- 28module coord_v_convert 29 30! module provided by C.Burigana & D.Maino (U. Milano) 31! edited by E.Hivon 32 33 use healpix_types 34 use misc_utils, only : strupcase 35 36 Implicit NONE 37 38 private 39 real(DP), parameter, private :: DTOR = PI/180.0_dp 40 41 public :: coordsys2euler_zyz ! produce Euler angles (ZYZ) for coordinate conversion 42 public :: xcc_v_convert ! front end routine 43 public :: xcc_DP_E_TO_Q ! Ecl -> Eq 44 public :: xcc_DP_E_TO_G ! Ecl -> Gal 45 public :: xcc_DP_G_TO_E ! Gal -> Ecl 46 public :: xcc_DP_Q_TO_E ! Eq -> Ecl 47 public :: XCC_DP_PRECESS ! precess ecliptic coord 48 49contains 50 subroutine coordsys2euler_zyz(iepoch, oepoch, isys, osys, psi, theta, phi) 51 ! 52 ! routine producing Euler angles psi, theta, phi 53 ! (for rotation around the unrotated Z, Y, Z axes respectively) 54 ! corresponding to changes of coordinates from ISYS to OSYS 55 ! for respective epochs IEPOCH and OEPOCH 56 ! the angles definition match the one required for rotate_alm. 57 ! EH, March 2005 58 ! April 06 2005, corrected psi calculation 59 ! 60 !use pix_tools, only: angdist 61 use num_rec, only: pythag 62 real(dp), intent(in) :: iepoch, oepoch 63 character(len=*), intent(in) :: isys, osys 64 real(dp), intent(out) :: psi, theta, phi 65 66 real(dp), dimension(1:3) :: v1, v2, v3, v1p, v2p, v3p 67 real(dp) :: st, ct 68 !----------------------- 69 70 ! generate basis vector 71 v1 = (/ 1.0_dp, 0.0_dp, 0.0_dp /) 72 v2 = (/ 0.0_dp, 1.0_dp, 0.0_dp /) 73 v3 = (/ 0.0_dp, 0.0_dp, 1.0_dp /) 74 75 ! rotate them 76 call xcc_V_CONVERT(v1,iepoch,oepoch,isys,osys,v1p) 77 call xcc_V_CONVERT(v2,iepoch,oepoch,isys,osys,v2p) 78 call xcc_V_CONVERT(v3,iepoch,oepoch,isys,osys,v3p) 79 80 ! normalize vectors 81 v1p = v1p / sqrt(dot_product(v1p,v1p)) 82 v2p = v2p / sqrt(dot_product(v2p,v2p)) 83 v3p = v3p / sqrt(dot_product(v3p,v3p)) 84 85 ! figure out Euler angles 86 87! call angdist(v3,v3p, theta) 88! psi = atan2(v2p(3),-v1p(3)) 89! phi = atan2(v3p(2), v3p(1)) 90 91 st = pythag(v3p(1),v3p(2)) ! |sin(theta)| in [0,1] 92 ct = v3p(3) ! cos(theta) in [-1,1] 93 theta = atan2(st, ct) ! theta in [0,Pi] 94 if (st >= 1.d-7) then 95 psi = atan2(v2p(3),-v1p(3)) 96 phi = atan2(v3p(2), v3p(1)) 97 else ! theta=0 or theta=Pi, phi and +/-psi degenerate 98 if (ct >= 0) then ! theta=0 99 ! only phi+psi is known 100 psi = atan2(-v2p(1)+v1p(2), v2p(2)+v1p(1)) ! returns psi+phi for any theta /= 180 101 else ! theta=pi, only psi-phi in known 102 psi = atan2(-v2p(1)-v1p(2), v2p(2)-v1p(1)) ! return psi-phi for any theta /= 0 103 endif 104 phi = 0.0_dp 105 endif 106 107 return 108 end subroutine coordsys2euler_zyz 109 110! 111 Subroutine xcc_V_CONVERT(ivector,iepoch,oepoch,isys,osys,ovector) 112 113! 114! Function to convert between standard coordinate systems, including 115! precession. 116! 117! 118! Q,C: Equatorial coordinates 119! E : Ecliptic coordinates 120! G : Galactic coordinates 121! SG : Super Galactic coordinates (TB implemented) 122! U : Universal Rest-frame coordinates (TB implemented) 123! 124! 03/14/88, a.c.raugh, STX. 125! 11/15/01, C.Burigana & D.Maino, OAT 126! 127! 128! Arguments: 129! 130 real(dp) :: IVECTOR(1:3) ! Input coordinate vector 131 real(dp) :: IEPOCH ! Input epoch 132 real(dp) :: OEPOCH ! Output epoch 133 real(dp) :: OVECTOR(1:3) ! Output coordinate vector 134 Character(len=*) :: ISYS ! Input system 135 Character(len=*) :: OSYS ! Output system 136! 137! Program variables: 138! 139 real(dp) :: XVECTOR(1:3) ! Temporary coordinate holder 140 real(dp) :: YVECTOR(1:3) ! Temporary coordinate holder 141 Character(len=20) :: isysu ! Input system, uppercase 142 Character(len=20) :: osysu ! Output system, uppercase 143! 144 145 isysu = trim(strupcase(isys)) 146 osysu = trim(strupcase(osys)) 147 148! If the input coordinate system was not ecliptic, convert to ecliptic 149! as intermediate form: 150! 151 if (trim(isysu) /= 'E') then 152 if (trim(isysu) == 'Q' .or. trim(isysu) == 'C') then 153 call xcc_dp_q_to_e(ivector,iepoch,xvector) 154 elseif (trim(isysu) == 'G') then 155 call xcc_dp_g_to_e(ivector,iepoch,xvector) 156! elseif (isys == 'SG') then 157! to be implemented 158! elseif (isys == 'U') then 159! to be implemented 160 else 161 print*,' unknown input coordinate system: '//trim(isysu) 162 endif 163 else 164 xvector(1:3)=ivector(1:3) 165 endif 166! 167! If the begin and end epoch are different, precess: 168! 169 if (iepoch /= oepoch) then 170 call xcc_dp_precess(xvector,iepoch,oepoch,yvector) 171 else 172 yvector(1:3)=xvector(1:3) 173 endif 174! 175! If output system is not ecliptic, convert from ecliptic to desired 176! system: 177! 178 if (trim(osysu) /= 'E') then 179 if (trim(osysu) == 'Q' .or. trim(osysu) == 'C') then 180 call xcc_dp_e_to_q(yvector,oepoch,ovector) 181 elseif (trim(osysu) == 'G') then 182 call xcc_dp_e_to_g(yvector,oepoch,ovector) 183 else 184 print*,' unknown output coordinate system: '//trim(osysu) 185 endif 186 else 187 ovector(1:3)=yvector(1:3) 188 endif 189 return 190 end subroutine xcc_v_convert 191!=============================================================== 192 193 Subroutine xcc_DP_E_TO_Q(ivector,epoch,ovector) 194! 195! Routine to convert from ecliptic coordinates to equatorial (celestial) 196! coordinates at the given epoch. Adapted from the COBLIB routine by Dr. 197! E. Wright. 198! 199! 02/01/88, a.c.raugh, STX. 200! 11/16/01, C.Burigana,D.Maino, OAT 201! 202! 203! Arguments: 204! 205 real(dp) :: IVECTOR(3) ! Input coordinate vector 206 real(dp) :: EPOCH ! Equatorial coordinate epoch 207 real(dp) :: OVECTOR(3) ! Output coordinate vector 208! 209! Program variables: 210! 211 real(dp) :: T ! Julian centuries since 1900.0 212 real(dp) :: EPSILON ! Obliquity of the ecliptic 213 real(dp) :: HVECTOR(3) ! Temporary holding place for output vector 214 real(dp) :: dc, ds 215! 216! Set-up: 217! 218 T = (epoch - 1900.d0) / 100.d0 219 epsilon = 23.452294d0 - 0.0130125d0*T - 1.63889d-6*T**2 & 220 & + 5.02778d-7*T**3 221! 222! Conversion: 223! 224 dc = cos(DTOR * epsilon) 225 ds = sin(DTOR * epsilon) 226 hvector(1) = ivector(1) 227 hvector(2) = dc*ivector(2) - ds*ivector(3) 228 hvector(3) = dc*ivector(3) + ds*ivector(2) 229! 230! Move to output vector: 231! 232 ovector(1:3) = hvector(1:3) 233! 234 return 235 end Subroutine xcc_DP_E_TO_Q 236!=============================================================== 237 238 Subroutine xcc_DP_E_TO_G(ivector,epoch,ovector) 239! 240! Routine to convert ecliptic (celestial) coordinates at the given 241! epoch to galactic coordinates. The ecliptic coordinates are first 242! precessed to 2000.0, then converted. 243! 244! -------------- 245! 246! Galactic Coordinate System Definition (from Zombeck, "Handbook of 247! Space Astronomy and Astrophysics", page 71): 248! 249! North Pole: 12:49 hours right ascension (1950.0) 250! +27.4 degrees declination (1950.0) 251! 252! Reference Point: 17:42.6 hours right ascension (1950.0) 253! -28 55' degrees declination (1950.0) 254! 255! -------------- 256! 257! 03/08/88, a.c.raugh, STX. 258! 11/16/01, C.Burigana, D.Maino, OAT 259! 260! 261! Arguments: 262! 263 real(dp) :: IVECTOR(3) ! Input coordinate vector 264 real(dp) :: EPOCH ! Input coordinate epoch 265 real(dp) :: OVECTOR(3) ! Output coordinate vector 266! 267! Program variables: 268! 269 real(dp) :: HVECTOR(3) ! Temporary holding place for coordinates 270 real(dp) :: T(3,3) ! Galactic coordinate transformation matrix 271 integer(i4b) :: I ! Loop variables 272! 273 Data T / -0.054882486d0, -0.993821033d0, -0.096476249d0, & ! 1st column 274 & 0.494116468d0, -0.110993846d0, 0.862281440d0, & ! 2nd column 275 & -0.867661702d0, -0.000346354d0, 0.497154957d0/ ! 3rd column 276! 277! Precess coordinates to 2000.0 if necesssary: 278! 279 if (epoch /= 2000.d0) then 280 call xcc_dp_precess(ivector,epoch,2000.d0,hvector) 281 else 282 hvector(1:3)=ivector(1:3) 283 endif 284! 285! Multiply vector by transpose of transformation matrix: 286! 287 DO I = 1,3 288 ovector(i) = sum(hvector(1:3)*T(1:3,i)) 289 ENDDO 290! 291 RETURN 292 END Subroutine xcc_DP_E_TO_G 293!=============================================================== 294 295 Subroutine xcc_DP_G_TO_E(ivector,epoch,ovector) 296! 297! Routine to convert galactic coordinates to ecliptic (celestial) 298! coordinates at the given epoch. First the conversion to ecliptic 299! 2000 is done, then if necessary the results are precessed. 300! 301! 03/08/88, a.c.raugh, STX. 302! 11/16/01, C.Burigana, D.Maino, OAT 303! 304! 305! Arguments: 306! 307 real(dp) :: IVECTOR(3) ! Input coordinate vector 308 real(dp) :: EPOCH ! Input coordinate epoch 309 real(dp) :: OVECTOR(3) ! Output coordinate vector 310! 311! Program variables: 312! 313 real(dp) :: HVECTOR(3) ! Temporary holding place for coordinates 314 real(dp) :: T(3,3) ! Galactic coordinate transformation matrix 315 integer(i4b) :: I,J ! Loop variables 316! 317 Data T / -0.054882486d0, -0.993821033d0, -0.096476249d0, &! 1st column 318 & 0.494116468d0, -0.110993846d0, 0.862281440d0, &! 2nd column 319 & -0.867661702d0, -0.000346354d0, 0.497154957d0/ ! 3rd column 320! 321! Multiply by transformation matrix: 322! 323 DO I = 1,3 324 hvector(i) = sum(ivector(1:3)*T(i,1:3)) 325 ENDDO 326! 327! Precess coordinates to epoch if necesssary: 328! 329 if (epoch /= 2000.d0) then 330 call xcc_dp_precess(hvector,2000.d0,epoch,ovector) 331 else 332 ovector(1:3)=hvector(1:3) 333 endif 334! 335 return 336 end Subroutine xcc_DP_G_TO_E 337!=============================================================== 338 339 Subroutine xcc_DP_Q_TO_E(ivector,epoch,ovector) 340! 341! Routine to convert equatorial (celestial) coordinates at the given 342! epoch to ecliptic coordinates. Adapted from the COBLIB routine by Dr. 343! E. Wright. 344! 345! 02/01/88, a.c.raugh, STX. 346! 11/16/01, C.Burigana, D.Maino, OAT 347! 348! 349! Arguments: 350! 351 real(dp) :: IVECTOR(3) ! Input coordinate vector 352 real(dp) :: EPOCH ! Epoch of equatorial coordinates 353 real(dp) :: OVECTOR(3) ! Output coordinate vector 354! 355! Program variables: 356! 357 real(dp) :: T ! Centuries since 1900.0 358 real(dp) :: EPSILON ! Obliquity of the ecliptic, in degrees 359 real(dp) :: HVECTOR(3) ! Temporary holding place for output vector 360 real(dp) :: dc, ds 361! 362! Set-up: 363! 364 T = (epoch - 1900.d0) / 100.d0 365 epsilon = 23.452294d0 - 0.0130125d0*T - 1.63889d-6*T**2 & 366 & + 5.02778d-7*T**3 367! 368! Conversion: 369! 370 dc = cos(DTOR * epsilon) 371 ds = sin(DTOR * epsilon) 372 hvector(1) = ivector(1) 373 hvector(2) = dc*ivector(2) + ds*ivector(3) 374 hvector(3) = dc*ivector(3) - ds*ivector(2) 375! 376! Move to output vector: 377! 378 ovector(1:3) = hvector(1:3) 379! 380 return 381 end Subroutine xcc_DP_Q_TO_E 382!=============================================================== 383 384 Subroutine XCC_DP_PRECESS(ivector,iepoch,oepoch,ovector) 385! 386! Subroutine to precess the input ecliptic rectangluar coordinates from 387! the initial epoch to the final epoch. First the precession around the 388! Z-axis is done by a simple rotation. A second rotation about the 389! X-axis is then performed to account for the changing obliquity of the 390! ecliptic. 391! 392! 02/25/88, a.c.raugh, STX. 393! 03/22/88, acr: modified to include changing epsilon 394! 11/16/01, Carlo Burigana & Davide Maino, OAT 395! 396! 397! Parameters: 398! 399 real(dp) :: IVECTOR(3) ! Input coordinates 400 real(dp) :: IEPOCH ! Initial epoch, years (input) 401 real(dp) :: OEPOCH ! Final epoch, years (input) 402 real(dp) :: OVECTOR(3) ! Output coordinates 403! 404! Routine variables: 405! 406 real(dp) :: TEMP ! Holding place for rotation results 407 real(dp) :: GP_LONG ! General precession in longitude 408 real(dp) :: dE ! Change in the obliquity of the ecliptic 409 real(dp) :: OBL_LONG ! -Longitude about which the obliquity is changing 410 real(dp) :: dL ! Longitude difference 411 real(dp) :: TM ! Mid-epoch, in tropical centuries since 1900 412 real(dp) :: TVECTOR(3) ! Temporary holding vector 413 real(dp) :: dco, dso, dce, dse, dcl, dsl 414! 415 Tm = ((oepoch+iepoch)/2.d0 - 1900.d0) / 100.d0 416 gp_long = (oepoch-iepoch) * (50.2564d0+0.0222d0*Tm) / 3600.d0 417 dE = (oepoch-iepoch) * (0.4711d0-0.0007d0*Tm) / 3600.d0 418 obl_long = 180.d0 - (173.d0 + (57.06d0+54.77d0*Tm) / 60.d0) & 419 & + gp_long / 2.d0 420 dL = gp_long - obl_long 421! 422! Z-axis rotation by OBL_LONG: 423! 424 dco = cos(DTOR * obl_long) 425 dso = sin(DTOR * obl_long) 426 tvector(1) = ivector(1)*dco - ivector(2)*dso 427 tvector(2) = ivector(1)*dso + ivector(2)*dco 428 tvector(3) = ivector(3) 429! 430! X-axis rotation by dE: 431! 432 dce = cos(DTOR * dE) 433 dse = sin(DTOR * dE) 434 temp = tvector(2)*dce - tvector(3)*dse 435 tvector(3) = tvector(2)*dse + tvector(3)*dce 436 tvector(2) = temp 437! 438! Z-axis rotation by GP_LONG - OBL_LONG: 439! 440 dcl = cos(DTOR * dL) 441 dsl = sin(DTOR * dL) 442 temp = tvector(1)*dcl - tvector(2)*dsl 443 tvector(2) = tvector(1)*dsl + tvector(2)*dcl 444 tvector(1) = temp 445! 446! Transfer results to output vector: 447! 448 ovector(1:3) = tvector(1:3) 449 450 return 451 end Subroutine XCC_DP_PRECESS 452! 453 454end module coord_v_convert 455 456 457