program radec2lb

  !###############################################################
  !
  !     Coordinate transformation from RA(J2000) Dec(J2000) to
  !     Galactic coordinats (l,b). The transformation is taken
  !     from Binney & Merrified (Galactic Astronomy) page 30.
  !
  !     Expects input to be in a file called 'radecs' of the
  !     form: _xx:xx:xx.x_sxx:xx:xx with one entry per line.
  !
  !     Output written to screen.
  !
  !     Bugs: None known
  !
  !     Original version in f77 by Simon Driver (SPD) 23/04/06
  ! 
  !      DATE            MODIFICATION              PROGRAMMER    
  !     24/03/06     Comment statements added  to     SPD
  !                  main programme body.
  !     25/03/06     Converted to f90.                SPD
  !     10/04/06     Converted to strict f90          TR
  !
  !###############################################################

  !     DELCARATIONS

  implicit none

  !     DECLARE AND SET PARAMETERS

  real :: ragc,decgc,lcp,pi

  parameter(ragc=192.85948,decgc=27.12825,lcp=123.932)
  parameter(pi=3.141592654)

  !     VARIABLE DECLARATIONS

  integer :: i,j,k,rah,ram,ras,dd,dm,ds

  real :: rass,dec,ra,l,b,lcpr
  real :: rar,decr,ragcr,decgcr

  character(len=1) :: sgn

  integer :: ioerr

  !     FUNCTION DECLARATIONS

  real :: decimal,rad,deg,getb,getl

  !###############################################################

  !     OPEN DATA FILE 'radecs'

  open(unit=1,file='radecs',status='old')

  !     WRITE TABLE HEADER TO SCREEN

  print *, ' '
  print *, '-----------------------------------------'
  print *, 'RA(J2000)  Dec(J2000)         l      b'
  print *, '-----------------------------------------'

  !     LOOP OVER ALL RAs AND DECs

  do

     !     READ TRADITIONAL RA-DEC FORMAT

     read(1,100,iostat=ioerr) rah,ram,ras,rass,sgn,dd,dm,ds
     if(ioerr.ne.0) exit

     !     CONVERT RA-DEC FORMAT TO DEC DEG

     if(sgn.eq.'+')then
        dec=decimal(real(dd),real(dm),real(ds))
     elseif(sgn.eq.'-')then
        dec=-decimal(real(dd),real(dm),real(ds))         ! remember the minus sign applies to the deg, minutes and seconds and not just the deg.
     else
        dec=decimal(real(dd),real(dm),real(ds))          ! if no sign or garbage assume sign is positive.
     endif

     ra=15.*decimal(real(rah),real(ram),real(ras)+rass) ! remember 24hours=360 deg hence x15. 

     !     CONVERT ALL ANGLES TO RADIANS

     decr=rad(dec)
     rar=rad(ra)
     decgcr=rad(decgc)
     ragcr=rad(ragc)
     lcpr=rad(lcp)

     !     CALCULATE GALACTIC COORDINATES 

     b=getb(decgcr,ragcr,decr,rar)
     l=getl(decgcr,ragcr,lcpr,decr,rar,b)

     !     NOTE: could have written the above as:
     !     b=getb(rad(decgc),rad(ragc),...) etc 
     !     and dispensed with the radian variables altogether
     !     or done the radian to degree conversion inside the
     !     functions. I've left it this way to make the code
     !     more transparent.

     !     CONVERT (l,b) FROM RADIANS TO DEGREES

     b=deg(b)
     l=deg(l)

     !     WRITE INPUT AND OUTPUTS TO THE SCREEN

     write(*,110) rah,':',ram,':',ras,rass,' ',&
          &sgn,dd,':',dm,':',ds,'  ===> ',l,b

     !     END LOOP OVER ALL RAs DECs

  enddo

  !     CLOSE DATA FILE

  close(unit=1)

  !     COMPLETE OUTPUT TABLE

  print *, '-----------------------------------------'
  print *, ' '

  !     FORMAT STATEMENTS

100 format(1x,i2.2,1x,i2.2,1x,i2.2,f2.1,1x,a1,i2.2,1x,i2.2,1x,i2.2)
110 format(1x,i2.2,a1,i2.2,a1,i2.2,f2.1,a1,a1,i2.2,a1,i2.2,a1,i2.2,&
         &a7,f7.2,f7.2)

  !     EXIT PROGRAM

end program radec2lb

!###############################################################
!
!     FUNCTIONS DEFINED BELOW
!
!     decimal - convert dd:mm:ss.s format to dec
!     rad - convert deg to rad
!     deg - convert rad to deg
!     getb - get Galactic Lattitude
!     getl - get Galactic Longitude
!
!###############################################################

function decimal(h,m,s)

  implicit none

  real :: decimal,h,m,s

  decimal=h+m/60.+s/3600.

end function decimal

!###############################################################

function rad(x)

  implicit none

  real :: pi
  parameter(pi=3.141592654)

  real :: rad,x

  rad=x*2.*pi/360.

end function rad

!###############################################################

function deg(x)

  implicit none

  real :: pi
  parameter(pi=3.141592654)

  real :: deg,x

  deg=x*360./(2.*pi)

end function deg

!###############################################################

function getb(dgc,rgc,d,r)

  implicit none

  real :: getb,dgc,rgc,d,r

  getb=asin(sin(dgc)*sin(d)+cos(dgc)*cos(d)*(cos(r-rgc)))

end function getb

!###############################################################

function getl(dgc,rgc,lcp,d,r,b)

  implicit none

  real :: pi
  parameter(pi=3.141592654)

  real :: dgc,rgc,lcp,d,r,b,sin_t,cos_t,t,l
  real :: getl

  !     Calculate sin(theta) and cos(theta)

  sin_t=(cos(d)*sin(r-rgc)/cos(b))
  cos_t=((cos(dgc)*sin(d)-sin(dgc)*cos(d)*cos(r-rgc))/cos(b))

  !     Use signs of sin_t and cos_t to determine which quadrant t lies in

  if((sin_t.ge.0).and.(cos_t.gt.0))then
     t=asin(sin_t)
  elseif((sin_t.ge.0).and.(cos_t.lt.0))then
     t=pi-asin(sin_t)
  elseif((sin_t.lt.0).and.(cos_t.lt.0))then
     t=pi-asin(sin_t)
  elseif((sin_t.lt.0).and.(cos_t.ge.0))then
     t=2.*pi-asin(sin_t)
  endif

  !     Convert t to l

  l=lcp-t

  !     Ensure l lies in the range 0-2pi 

  if(l.lt.0.) l=l+(2.*pi)
  if(l.gt.(2.*pi)) l=l-(2.*pi)

  getl=l

end function getl

!###############################################################