c Program to transform coordinates between equatorial and galactic ones. c 2: equinox 1950 was included c Author: S. Mashchenko (syam@astro.umontreal.ca; Mashchenko@excite.com) program coords implicit NONE c Maximum number of fields: real*8 Pi, RA_ss, DEC_ss, al, de, l, b, Equinox,Numbers(3), # alpha,delta integer RA_hh, RA_mm, DEC_dd, DEC_mm, nn logical Gal character Dir*2, String*256 Pi = dacos(-1.d0) c Reading inputs: write (*,*) ' Program for galactic <-> equatorial coordinates' # //' transformation' write (*,*) ' (Press Ctrl-C to quit)' write (*,*) 2 continue write (*,*) write (*,*) ' Galactic->equatorial (\"ge\"), or equatorial->'// # 'galactic (\"eg\"):' read (*,*) Dir if (Dir.ne.'ge'.and.Dir.ne.'eg'.and.Dir.ne.'GE'.and.Dir.ne.'EG') # then write (*,'(A)') 'Wrong command: '//Dir go to 2 end if write (*,*) ' Equinox for the equatorial coordinates:' read (*,*) Equinox Gal = Dir.eq.'ge'.or.Dir.eq.'GE' 1 continue write (*,*) if (Gal) then c Galactic -> equatorial transformation: write (*,*) ' Galactic coordinates:' write (*,*) 'l (dgr.):' read (*,*) l write (*,*) ' b (dgr.):' read (*,*) b l = l/180.d0*Pi b = b/180.d0*Pi call galactic_equatorial (Gal,Equinox,l,b,al,de) al = al*180.d0/Pi de = de*180.d0/Pi RA_hh = int(al/15) RA_mm = int((al/15-RA_hh)*60) RA_ss = (al/15-RA_hh-RA_mm/60.d0)*3600 DEC_dd = int(dabs(de)) * dsign(1.d0,de) DEC_mm = int((dabs(de)-iabs(DEC_dd))*60) DEC_ss = (dabs(de)-iabs(DEC_dd)-DEC_mm/60.d0)*3600 write (*,'(A,F9.5,A,F9.5,A)') ' RA=',al,' dgr.; Dec=',de, # ' dgr.' write (*,'(2(A,I3,X,I2,X,F6.3),A)') # ' RA=',RA_hh,RA_mm,RA_ss, # '; Dec=',DEC_dd,DEC_mm,DEC_ss else c Equatorial -> galactic transformation: write (*,'(A,F6.1,A)') ' Equatorial coordinates (', # Equinox,'):' 5 continue write (*,*) 'RA:' write (*,*) ' (Accepted formats: \"hh mm ss.s\", \"hh mm.m\",' # //' \"dgr\")' read (*,'(A)') String call Extract_numbers (String, nn,Numbers) if (nn .eq. 0) go to 5 al = Alpha (Numbers,nn) 6 continue write (*,*) 'Dec:' write (*,*) ' (Accepted formats: \"dd mm ss.s\", \"dd mm.m\",' # //' \"dgr\")' read (*,'(A)') String call Extract_numbers (String, nn,Numbers) if (nn .eq. 0) go to 6 de = Delta (Numbers,nn) al = al/180.d0*Pi de = de/180.d0*Pi call galactic_equatorial (Gal,Equinox,l,b,al,de) l = l*180.d0/Pi b = b*180.d0/Pi write (*,'(A,F9.5,A,F9.5,A)') ' l=',l,' dgr.; b=',b, # ' dgr.' end if go to 1 stop end c============================================================================ subroutine galactic_equatorial (Dir,Equinox,l,b, al,de) c If Dir: galactic -> equatorial coordinates transformation, if c .not.Dir: equatorial -> galactic transformation. implicit NONE real*8 pi, e,fa,fl, e2,fa2,fl2, e1,fa1,fl1, Equinox,l,b, al,de, w parameter ( pi = 3.141592653589793115997963468544185161591D+00 ) c Parameters for Equinox=B1950.0: parameter ( e1 = 62.6d0 /180.d0*pi ) parameter ( fa1 = 282.25d0 /180.d0*pi ) parameter ( fl1 = 33.0d0 /180.d0*pi ) c Parameters for Equinox=J2000.0: parameter ( e2 = 62.8717d0 /180.d0*pi ) parameter ( fa2 = 282.8596d0 /180.d0*pi ) parameter ( fl2 = 32.93192d0 /180.d0*pi ) logical Dir if (Equinox .eq. 1950.0d0) then e = e1 fa = fa1 fl = fl1 elseif (Equinox .eq. 2000.0d0) then e = e2 fa = fa2 fl = fl2 else write (*,*) ' In galactic_equatorial only Equinox=1950.0'// # ' or 2000.0 are supported!' stop end if if (Dir) then w = datan2(dsin(l-pi/2.d0-fl),(dcos(l-pi/2.d0-fl)*cos(e)- # dtan(b)*sin(e))) al = (w+fa-1.5d0*pi) de = dasin(sin(b)*cos(e)+cos(b)*sin(e)*cos(l-pi/2.d0-fl)) c al - within interval [0...2*pi] al = al - 2*pi*int(al/(2*pi)) if (al .lt. 0.d0) al = al + 2*pi else l = datan2(dcos(de)*dsin(al-fa)*dcos(e)+dsin(de)*dsin(e), # dcos(de)*dcos(al-fa)) + fl c l - within interval [0...2*pi] l = l - 2*pi*int(l/(2*pi)) if (l .lt. 0.d0) l = l + 2*pi b = dasin(dsin(de)*dcos(e)-dcos(de)*dsin(al-fa)*dsin(e)) end if return end c============================================================================ subroutine Extract_numbers (String, nn,Numbers) c Subroutine to extract 1, 2 or 3 numbers c from the String, and to put them into array Numbers(1..nn). implicit NONE integer nn,i,i1,i2 real*8 Numbers(3) character*(*) String, s*1 logical st st = .false. nn = 0 do i=1,len(string) s = string(i:i) if(s.eq.'0'.or.s.eq.'1'.or.s.eq.'2'.or.s.eq.'3'.or.s.eq.'4'.or. # s.eq.'5'.or.s.eq.'6'.or.s.eq.'7'.or.s.eq.'8'.or.s.eq.'9'.or. # s.eq.'+'.or.s.eq.'-'.or.s.eq.'.'.or.s.eq.'e'.or.s.eq.'E'.or. # s.eq.'d'.or.s.eq.'D') then if (.not.st) i1 = i st = .true. else if (st) then i2 = i-1 if (nn .lt. 3) then nn = nn + 1 read (string(i1:i2),*) Numbers(nn) end if st = .false. end if end if end do if (nn .eq. 0) then write (*,*) 'There are no numbers in the input string!' end if end c============================================================================ function Alpha (Numbers,nn) c The function transforms three different formats of the equatorial coordinate c alpha into degrees: c nn=1: dgr. -> dgr. (no change); {dgr.=Numbers(1)} c nn=2: hh mm.m -> dgr.; {hh=Numbers(1), mm.m=Numbers(2)} c nn=3: hh mm ss.s -> dgr. {hh=Numbers(1), mm=Numbers(2), ss.s=Numbers(3)} c The result is returned in Alpha. implicit NONE integer nn real*8 Numbers(nn), Alpha if (nn.gt.3 .or. nn.lt.1) then write (*,*) ' Wrong nn in Alpha!' stop end if if (nn .eq. 1) then Alpha = Numbers(1) else if (nn .eq. 2) then Alpha = 15.d0 * (Numbers(1) + Numbers(2)/60.d0) else Alpha = 15.d0 * (Numbers(1) + Numbers(2)/60.d0 + # Numbers(3)/3600.d0) end if return end c============================================================================ function Delta (Numbers,nn) c The function transforms three different formats of the equatorial coordinate c delta into degrees: c nn=1: dgr. -> dgr. (no change); {dgr.=Numbers(1)} c nn=2: dd mm.m -> dgr.; {dd=Numbers(1), mm.m=Numbers(2)} c nn=3: dd mm ss.s -> dgr. {dd=Numbers(1), mm=Numbers(2), ss.s=Numbers(3)} c The result is returned in Delta. implicit NONE integer nn real*8 Numbers(nn), Delta if (nn.gt.3 .or. nn.lt.1) then write (*,*) ' Wrong nn in Delta!' stop end if if (nn .eq. 1) then Delta = Numbers(1) else if (nn .eq. 2) then Delta = dsign(1.d0,Numbers(1)) * (dabs(Numbers(1)) + # Numbers(2)/60.d0) else Delta = dsign(1.d0,Numbers(1)) * (dabs(Numbers(1)) + # Numbers(2)/60.d0 + Numbers(3)/3600.d0) end if return end c============================================================================ INTEGER FUNCTION strlen_trail (STRING) C C Return the length of a string less any trailing blanks C C C STRING CH*(*) input String to be measured implicit NONE CHARACTER*(*) STRING INTEGER I, N N = LEN (STRING) DO 1 I = N,1,-1 IF(STRING(I:I).NE.' ') GOTO 2 1 CONTINUE 2 STRLEN_trail = I END c--------------------------------------------------------------------------