character*30 dis,ns1
	character*12 ac,dc
        pi=4.*atan(1.0)
	cx=pi/12.
	cy=pi/180.

	write(*,*)
	write(*,*)'   Calcultae BH Bmag, EXTINCTION = 4*E(B-V) '
	write(*,*)
	write(*,*)'                                    2000.7'
 	write(*,*)
	ns1(1:2)='  '
	ep_old=2000.

10      if(ns1(1:2).eq.'  ')write(*,9001)ep_old
9001    format('Input: alpha,delta,epoch',
     * ' (Ex.3:46:02.11 23:35:47.8 [',f6.1,'])  "q"~quit')
        if(ns1(1:2).ne.'  ')write(*,9009)ns1
9009    format('Input: alpha,delta,epoch [',a,']  "q"~quit')
        read(*,'(a)')dis
        if(dis(1:2).eq.'  ')dis=ns1(2:)
        if(dis(1:1).eq.'q'.or.dis(1:1).eq.'Q')stop
        if(itohd(dis,alpha,delta,x).eq.0)goto 10        ! char_line  to 3 real
        if(x.ne.0.)ep_old=x
        goto 112                                ! if char_line include 3 item
11      write(*,*)'epoch (ex. 1950.0)'
        read(*,*,err=11)ep_old
112     if(ep_old.lt.1900..or.ep_old.gt.2100.)goto 11
 
        call astprs(alpha,delta,ep_old,x,y,1950.)
        call galactic(x*cx,y*cy,xa,xd,cy)
	xa=redden(xa,xd)
	call toms(alpha,ac,1)
	call toms(delta,dc,0)
	write(*,9005)ac,dc,ep_old,xa
9005	format(a,2x,a,'(',f6.1,')   EXTINCTION = ',f7.4/)
	goto 10
	end

        function redden(al,b)
c this subroutine written by burstein.
      INTEGER*4 IRED(1200),IHI(201)
      XRAD = 57.2958
      ALIN = AL/0.3 +  0.5
      ILIN = ALIN + 0.51
      BB = ABS(B)
      redden= -0.99
      if(bb.lt.10.)  return
      OPEN(UNIT=11,FILE='/EOD/bur/redsouth.dat',
     * STATUS='OLD',ACCESS='DIRECT',recl=4800)
      OPEN(UNIT=12,FILE='/EOD/bur/rednorth.dat',
     * STATUS='OLD',ACCESS='DIRECT',recl=4800)
      OPEN(UNIT=13,FILE='/EOD/bur/hinorth.dat',
     * STATUS='OLD',ACCESS='DIRECT',recl=804)
      OPEN(UNIT=14,FILE='/EOD/bur/hisouth.dat',
     * STATUS='OLD',ACCESS='DIRECT',recl=804)
      IF (B) 11,11,18
   11 IF (B + 62.) 13,12,12
   12 AREC = -(B+10.)/0.6 + 1.
      KZ = AREC + 0.51
      READ (11,rec=KZ) IRED
	call swap4(ired,4800)
      KQ = KZ -1
      IEBV = IRED(ILIN)
      GO TO 19
   13 ALC = AL/XRAD
      AREC = 101. + SIN(ALC)*(90.+B)/0.3
      ALIN = 101. + COS(ALC)*(90.+B)/0.3
      KW = AREC + 0.51
      ILIN = ALIN + 0.51
      READ (14,rec=KW) IHI
	call swap4(ihi,804)
      AHI = IHI(ILIN)
      BMV = -0.0372 + 0.357*AHI/10000.
      IEBV = BMV*1000. + 0.5
      KQ = KW - 1
      GO TO 19
   18 IF (B - 62.) 17,17,16
   17 AREC = (B-10.)/0.6 + 1.
      KY = AREC + 0.51
      READ (12,rec=KY) IRED
	call swap4(ired,4800)
      KQ = KY - 1
      IEBV = IRED(ILIN)
      GO TO 19
   16 ALC = AL/XRAD
      AREC = 101. + SIN(ALC)*(90.-B)/0.3
      ALIN = 101. + COS(ALC)*(90.-B)/0.3
      KX = AREC + 0.51
      ILIN = ALIN + 0.51
      READ (13,rec=KX) IHI
	call swap4(ihi,804)
      AHI = IHI(ILIN)
      BMV = -0.0372 + 0.357*AHI/10000.
      IEBV = BMV*1000. + 0.5
      KQ = KX - 1
   19 ABV = IEBV
      redden= ABV*4./1000. + 0.005
	close(11)
	close(12)
	close(13)
	close(14)
      return
      END

        subroutine galactic(alpha,delta,xl,xb,cy)
        a0=282.25*cy
        theta=62.6*cy
        Xl0=33.0*cy
        x1=cos(delta)*cos(alpha-a0)
        x2=cos(delta)*sin(alpha-a0)*cos(theta)+sin(delta)*sin(theta)
        xl=atan2(x2,x1)+xl0
        if(xl.lt.0.)xl=xl+2.0*pi
        xb=asin(sin(delta)*cos(theta)-
     *          cos(delta)*sin(alpha-a0)*sin(theta))
        xl=xl/cy
        xb=xb/cy
        end