      program main
      dimension dis(1000),x(1000),y(1000),yy(1000),r(1000),
     c     z(1000),zco(1000),zh(1000)
      dimension ah1(1000),aco(1000),zmet(1000),agas(1000),
     c     agastot(1000)


C program computes extinction between Sun and a point with given galactic coordinates
C and distance
      r0=7.5
C  adopted distance of the Galactic center
      conv=3.14159/180.
      step  = 0.05
      dist = 10.0
C distance of the point to which we will calculate the extinction


      print*, 'Interstellar Extinction in the Galaxy',
     c     '(Amores & Lépine - 2004)'
      print*,'Give the galactic longitude and latitude (Degrees)'
      read*, glong,glat
      print*,  'Give distance in kpc'
      read*, dist

      nstep=int(dist/step)+1
C we added +1 to get the same results given by IDL fix

           
C computes  trigonometric functions only once
      yproj=cos(glong*conv)
      xproj=sin(glong*conv)
      bproj=sin(glat*conv)
      dproj=cos(glat*conv)
      av=0.0                  
C   for the integration of the colunar density

C declaring and putting values in the variables. The arrays will contain the
C  quantities like galactic radius or gas density for each step along the line-of sight
C if you work with other language you should probably define these quantities in a loop

      
      do i=1,nstep
         ipas=i
         dis(i)=i*step
         x(i)=dis(i)*xproj*dproj
         y(i)=dis(i)*yproj*dproj
         yy(i)=r0-y(i)
         r(i)=sqrt(x(i)*x(i)+yy(i)*yy(i))
         z(i)=dis(i)*bproj
         
         zco(i)=0.036*exp(0.08*r(i))    
c       H2 scale-height
         zh(i) = zco(i)*1.8
C       H1 scale-height (Guilbert 1978)
         zc = 0.02          
C       shift takes in to account that the sun is not precisely in the galactic plane
	   
         ah1(i)=0.7*exp(-r(i)/7.0-((1.9/r(i))**2))
C function that calculates the HI density
         
         aco(i) = 58.*exp(-r(i)/1.20-((3.50/r(i))**2)) + 
     c        240.*exp(-(r(i)**2/0.095))
C  H2 density; last term is for galactic center
      enddo

      do i=1,nstep
         zmet(i)=1.0
         if (r(i).le.1.2)  zmet(i) = 9.6
         if ((r(i).gt. 1.2).and.(r(i).le.9.0)) zmet(i)=(r0/r(i))**0.5
         if (r(i).gt. 9.0)   zmet(i) = (r0/r(i))**0.1
      enddo
      gam1=1.0
      gam2=2.0

C	See the final tuning (section 4.1) correction factor for interval l=120-200

      tune=1.
      if((glong.ge.120.).and.(glong.le. 200.))  tune=2.
      agas(1)=gam1*(ah1(1)*zmet(1)*exp(-0.5*((z(1)-zc)/zh(1))**2)) +   
     c     gam2*aco(1)*exp(-0.5*((z(1)-zc)/zco(1))**2)
      agastot(1)=agas(1)
      
      do k=2,nstep
         agas(k)=gam1*(ah1(k)*zmet(k)*exp(-0.5*((z(k)-zc)/zh(k))**2)) +   
     c        gam2*aco(k)*exp(-0.5*((z(k)-zc)/zco(k))**2)
         agastot(k)=agastot(k-1) + agas(k)
      enddo

      av=agastot(nstep)*step*3.086*.57*tune

c   agastot is the integral along line-of-sight of the gas density
c   step is in units of kpc= 3.08*10^21 cm and 
c   the conversion factor gamma= .57 10^-21 mag cm2

      rs = 3.05
C	ratio between total to selective extinction
      ebv = av/rs

      print*,'E(B-V)= ',ebv,' mag', 'Av= ', av, 'mag'
      end
