implicit none
      double precision in(2),out(2), RaPrevious, DecPrevious
      integer flag,i,j

      call pgbegin(0,'?',1,1)
      call pgenv(360.,0.,-90.,90.,0,1)
      call pglabel('Right Ascension', 'Declination',
     $     'Equatorial-Ecliptic-Galactic Conversion')
      call pgslw(3)
      call pgbox('a', 10.0, 1.0, 'a', 10.0, 1.0)
      call pgslw(1)
      call pgbox('g', 10.0, 1.0, 'g', 10.0, 1.0)
C     Draw Galactic Grid with b=const
      do flag = 1, 3, 2
         if(flag.eq.1) call pgsci(2)
         if(flag.eq.3) call pgsci(3)
      do j = 1, 17
         RaPrevious=999.
         if(j.eq.9) then
            call pgslw(3)
         else
            call pgslw(1)
         endif
         do i = 1, 361
            in(1)=1.0*(i-1)
            in(2)=-90.0+j*10.0
            call coord_tran(in,out,flag)
c            write(*,*) i, out(1),out(2)
            if(abs(RaPrevious-out(1)).ge.180.0) then
               call pgmove(REAL(out(1)),REAL(out(2)))
            else
               call pgdraw(REAL(out(1)),REAL(out(2)))
            endif
            RaPrevious=out(1)
         end do
      end do
C     Draw Galactic Grid with l=const
       do j = 1, 37
c       do j = 1, 1
         DecPrevious=999.
         if(j.eq.1) then
            call pgslw(3)
         else
            call pgslw(1)
         endif
         do i = 1, 181
            in(1)=  0.0+10.0*(j-1)
            in(2)= -90.0+(i-1)
            call coord_tran(in,out,flag)
c            write(*,*) i, out(1),out(2)
            if(abs(RaPrevious-out(1)).ge.180.0) then
               call pgmove(REAL(out(1)),REAL(out(2)))
            else
               call pgdraw(REAL(out(1)),REAL(out(2)))
            endif
            RaPrevious=out(1)
         end do
      end do
      end do

      call pgslw(3)
      call pgsci(2)
      call pgtext(216., -68.,'Galactic Plane')
      call pgtext(270., -35.,'GC')
      call pgtext(190.,  30.,'NGP')
      call pgtext( 10., -30.,'SGP')
      call pgptxt(240., -8.,48., 0.5, 'Longitude=0')

      call pgsci(3)
      call pgptxt(216., -20.,26., 0., 'Ecliptic Plane')
      call pgtext(270.,  70.,'NEP')
      call pgtext( 85., -73.,'SEP')
      call pgptxt(344., 20.,70., 0.5, 'Longitude=0')
      call pgptxt(33., -58.,40., 0.5, 'Longitude=0')

      call pgend
      end

CCCCCCCCCCCC
C     Coordinate transformation
C     If flag = 0, equatrial --> Galactic
C     If flag = 1, Galactic  --> equatrial 
C     If flag = 2, equatrial --> ecliptic
C     If flag = 3, ecliptic  --> equatorial
CCCCCCCCCCCC
      subroutine coord_tran(in, out, flag)
      implicit none
CCCCCCCCCCCC
      integer flag
      double precision in(2), out(2)
      double precision matrix1(3,3), matrix2(3,3), matrix3(3,3)
      double precision matrix4(3,3), matrix5(3,3), matrix6(3,3)
      double precision invect(3), outvect(3)
      double precision outvect1(3),outvect2(3)
      double precision dec2rag
      double precision rag2dec

      rag2dec = 45.0D0/atan(1.0D0);
      dec2rag = atan(1.0D0)/45.0D0;

      invect(1)=cos(dec2rag*in(1))*cos(dec2rag*in(2))
      invect(2)=sin(dec2rag*in(1))*cos(dec2rag*in(2))
      invect(3)=sin(dec2rag*in(2))

c      write(*,*) 'invect =', invect

      call check(invect)

      if (flag.eq.0) then
         call set_matrix(matrix1, matrix2, matrix3)
         call multi_matrix(invect,  matrix3, outvect1)
c         write(*,*) 'outvect1 =', outvect1
         call check(outvect1)
         call multi_matrix(outvect1, matrix2, outvect2)
c         write(*,*) 'outvect2 =', outvect2
         call check(outvect2)
         call multi_matrix(outvect2, matrix1, outvect)
         call check(outvect)
c         write(*,*) 'outvect =', outvect
      elseif (flag.eq.1) then
         call set_matrix(matrix1, matrix2, matrix3)
         
         call transpose_matrix(matrix1, matrix4)
         call transpose_matrix(matrix2, matrix5)
         call transpose_matrix(matrix3, matrix6)
         
         call multi_matrix(invect,  matrix4, outvect1)
         call check(outvect1)
         call multi_matrix(outvect1, matrix5, outvect2)
         call check(outvect2)
         call multi_matrix(outvect2, matrix6, outvect)
         call check(outvect)
      elseif(flag.eq.2) then
         call set_matrix2(matrix1)
         call multi_matrix(invect, matrix1, outvect)
         call check(outvect)
      elseif(flag.eq.3) then
         call set_matrix2(matrix1)
         call transpose_matrix(matrix1, matrix2)
         call multi_matrix(invect, matrix2, outvect)
         call check(outvect)
      endif

      if(outvect(1).eq.0.0D0) then
         if(outvect(2).eq.0.0D0) then
            out(1) = 0.0D0
         elseif(outvect(2).gt.0.0) then
            out(1)= 90.0D0
         else
            out(1)= -90.0D0
         endif
      else
         out(1) = atan(outvect(2)/outvect(1))*rag2dec
c         write(*,*) 'atan =', atan(outvect(2)/outvect(1))*rag2dec
      endif
      if(outvect(1).lt.0.0) then
         out(1) = out(1) + 180.0D0
c         write(*,*) 'out(1) =', out(1)
      endif
      if (out(1).lt.0.0) then
         out(1) = out(1) + 360.0D0
      endif
      if(abs(outvect(3)) .gt. 1.0) then
         write(*,*) 'argument of asin below -1 or above 1'
         write(*,*) 'argument of asin',outvect(3)
         stop
      endif
      out(2) = asin(outvect(3))*rag2dec
      
      end 
      
CCCCCCCCCCC
      subroutine multi_matrix(in, matrix, out)
      implicit none
CCCCCCCCC
      double precision in(3),out(3),matrix(3,3)
      integer i, j
      
      do i = 1, 3
         out(i)=0.0D0
         do j = 1, 3
            out(i) = out(i)+matrix(i,j)*in(j)
         end do
      end do

      end

CCCCCCCCCCC
      subroutine set_matrix(matrix1,matrix2,matrix3)
      implicit none
CCCCCCCCCCC
      double precision matrix1(3,3), matrix2(3,3), matrix3(3,3) 
      double precision dec2rag
      double precision alpha, delta, theta
      dec2rag = atan(1.0D0)/45.0D0;

CCCCCCCCCCCCCCCCCC
C     J2000
C     Galactic center alpha and opposite sine of delta
      alpha = 266.404996D0
      delta =  28.936172D0

C     theta is the rotation angle around the axis to the Galactic center
      theta = 58.598666D0;
CCCCCCCCCCCCCCCCCC

      alpha = dec2rag*alpha
      delta = dec2rag*delta
      theta = dec2rag*theta

      matrix1(1,1)= 1.0D0
      matrix1(1,2)= 0.0D0
      matrix1(1,3)= 0.0D0
      matrix1(2,1)= 0.0D0
      matrix1(2,2)= cos(theta)
      matrix1(2,3)= sin(theta)
      matrix1(3,1)= 0.0D0
      matrix1(3,2)= -sin(theta)
      matrix1(3,3)= cos(theta)

      matrix2(1,1)= cos(delta)
      matrix2(1,2)= 0.0D0
      matrix2(1,3)= -sin(delta)
      matrix2(2,1)= 0.0D0
      matrix2(2,2)= 1.0D0
      matrix2(2,3)= 0.0D0
      matrix2(3,1)= sin(delta)
      matrix2(3,2)= 0.0D0
      matrix2(3,3)= cos(delta)

      matrix3(1,1)= cos(alpha)
      matrix3(1,2)= sin(alpha)
      matrix3(1,3)= 0.0D0
      matrix3(2,1)= -sin(alpha)
      matrix3(2,2)= cos(alpha)
      matrix3(2,3)= 0.0D0
      matrix3(3,1)= 0.0D0
      matrix3(3,2)= 0.0D0
      matrix3(3,3)= 1.0D0
      end

CCCCCCCCCCC
      subroutine set_matrix2(matrix1)
      implicit none
CCCCCCCCCCC
      double precision matrix1(3,3)
      double precision e
      double precision dec2rag

      dec2rag = atan(1.0D0)/45.0D0;

CCCCCCCCCCCCCCCCCC
C     obliquity of the ecliptic
c      e = 23.45D0
      e = 23.439291D0
CCCCCCCCCCCCCCCCCC

      e = e * dec2rag

      matrix1(1,1)= 1.0D0
      matrix1(1,2)= 0.0D0
      matrix1(1,3)= 0.0D0
      matrix1(2,1)= 0.0D0
      matrix1(2,2)= cos(e)
      matrix1(2,3)= sin(e)
      matrix1(3,1)= 0.0D0
      matrix1(3,2)= -sin(e)
      matrix1(3,3)= cos(e)

      end
      
CCCCCCCCCCC
      subroutine transpose_matrix(inmatrix, outmatrix)
      implicit none
CCCCCCCCCCC
      integer i, j
      double precision inmatrix(3,3), outmatrix(3,3)
      do i = 1, 3
         do j = 1, 3
            outmatrix(i,j) = inmatrix(j,i)
         end do 
      end do
      end

      subroutine check(vector)
      implicit none
      double precision vector(3),length,tolerance

      tolerance = 1e-5

      length=vector(1)**2+vector(2)**2+vector(3)**2
      
      if(abs(length-1.0D0).gt.tolerance) then
         write(*,*) '*** vector length not 1.0!! ***'
         stop
      endif
      end