program exti
c
c      calculate extinction curve using Q values from Bruce Draine's Home Page
c      at: http://www.astro.princeton.edu/~draine/dust/dust.diel.html
c      Please credit the relevant papers from the above site when using this
c      program
c
      parameter (nf=241,nd=81,nff=1205)
      real qg(nd,nf,3),qs(nd,nf,3),qsic(nd,nf,3),w(nf),gmean(nf)
     &,rgr(nd),exts1(nf),exts2(nf),extg1(nf),extg2(nf),extt(nf)
     &,area(nd)
      data exts1,exts2,extg1,extg2,gmean/nff*0./
c
c  Gra_81     - Graphite grains
c  suvSil_81  - UV smoothed Silicate grains
c  SiC_81     - Silicon Carbide grains
c
       open(11,file='Gra_81',status='old')
       open(12,file='suvSil_81',status='old')
       open(13,file='SiC_81',status='old')
c
c  reading grains Q values;
c
      read(11,'(/,/,/)')
      do 1 i=1,81
       read(11,*)
       read(11,*) rgr(i)
       read(11,*)
       do 1 j=1,241
        read(11,*)w(j),qg(i,j,1),qg(i,j,2),qg(i,j,3)
 1    continue
      read(12,'(/,/,/)')
      do 2 i=1,81
       read(12,*)
       read(12,*) rgr(i)
       read(12,*)
       do 2 j=1,241
        read(12,*)w(j),qs(i,j,1),qs(i,j,2),qs(i,j,3)
 2    continue
      read(13,'(/,/,/)')
      do 3 i=1,81
       read(13,*)
       read(13,*) rgr(i)
       read(13,*)
       do 3 j=1,241
        read(13,*)w(j),qsic(i,j,1),qsic(i,j,2),qsic(i,j,3)
 3    continue
      siliden=3.3   ! Silicate grain density (gr cm^-3)
      sicden=3.22   ! Silicon Carbide grain density (gr cm^-3)
      absic=.23
      absil=1.12 !  MRN   Sil/Gr number density ratio
      fsil=1.   ! Sil/Gr/(Sil/Gr)_MRN
      imod=1
      alpha=-3.5  ! grain size distribution power law slope
      rmin=0.005 ! minimum grain radius, 0.005 mic in MRN
      rmax=0.25  ! maximum grain radius, 0.25 mic in MRN
      fmas=0.01  ! Dust/gas mass ratio, 0.01 in MRN model
c
c  Allow modification of MRN parameters
c
      print*,'Type in [1/0], 1 = use MRN model, 0 = specify your model'
      read *,mrn
      If(mrn.eq.1) go to 20
      if(mrn.ne.0) stop 'wrong value for mrn'
      print*,'Type in imod, 1 = Graphite+Silicate, 2 = Graphite+SiC'
      print*,'MRN is imod=1'
      read*,imod
      print*,'Type in fsil = Silicate (or SiC)/Graphite enhancement'
      print*,'MRN is fsil=1'
      read *,fsil
      print*,'Type in alpha = slope of grain size power-law dist'
      print*,'MRN is alpha=-3.5'
      read *, alpha
      print*,'Type in rmin and rmax - min and max grain radii (mic)'
      print*,'values should be in the range 0.001-10'
      print*,'MRN is 0.005 0.25'
      read*, rmin, rmax
      if(rmin.lt..001)stop 'Rmin < 0.001 mic'
      if(rmax.gt.10.) stop 'Rmax > 10 mic'
      print*,'Type in dust/gas mass ratio'
      print*,'normal ISM (MRN) is 0.01'
      read*,fmas
 20   print*,'INPUT: Dust Model Parameters are'
      print*,'imode=',imod,' fsil=',fsil,' alpha=',alpha
      print*,'rmin=',rmin,'  rmax=',rmax,'  fmas=',fmas
      if(imod.eq.2) absil=absic
      if(imod.eq.2) siliden=sicden
      absil=absil*fsil
c
c  Calculate total grain volume- volgr,
c  and grain area within each bin- area(id)
c
      volgr=0.    ! initialize total grain volume
      do 4 id=1,nd
       cor=1.
       if(id.eq.1) then
        d1=rgr(id)
        d2=(rgr(id)+rgr(id+1))/2.
        if(rmin.gt.d2) cor=0.
        if(rmin.le.d2) d2=rmin
       elseif(id.eq.nd) then
        d1=(rgr(id-1)+rgr(id))/2.
        d2=rgr(id)
        if(rmax.lt.d1) cor=0.
        if(rmax.ge.d1) d1=rmax
       else
        d1=(rgr(id-1)+rgr(id))/2.
        d2=(rgr(id)+rgr(id+1))/2.
        if(rmin.le.d2.and.rmin.gt.d1) d1=rmin
        if(rmin.gt.d2) cor=0.
        if(rmax.le.d2.and.rmax.gt.d1) d2=rmax
        if(rmax.lt.d1) cor=0.
       endif
       dvol=d2**(4.+alpha)-d1**(4.+alpha)
       volgr=volgr+cor*dvol/(4.+alpha)/rmin**alpha
       area(id)=d2**(3.+alpha)-d1**(3.+alpha)
       area(id)=cor*area(id)/(3.+alpha)/rmin**alpha
 4    continue
       xmas=volgr*(siliden*absil+2.26)  ! total grain mass
c
c              calculating the extinction curve
c
c   extg1 = Graphite absorption
c   extg2 = Graphite scattering
c   exts1 = Silicate absorption
c   exts2 = Silicate scattering
c   extt = Total extinction
c   gmean = Mean g (=cos theta_scattering) 0 = isotropic, 1 = forward
c
      do 5 id=1,nd
        area(id)=area(id)*fmas*1.67e-24*.75/xmas*1.e4
      do 5 if=1,nf
       if(imod.eq.2) qs(id,if,1)=qsic(id,if,1)
       if(imod.eq.2) qs(id,if,2)=qsic(id,if,2)
       if(imod.eq.2) qs(id,if,3)=qsic(id,if,3)
       extg1(if)=extg1(if)+qg(id,if,1)*area(id)
       extg2(if)=extg2(if)+qg(id,if,2)*area(id)
       exts1(if)=exts1(if)+absil*qs(id,if,1)*area(id)
       exts2(if)=exts2(if)+absil*qs(id,if,2)*area(id)
       gmean(if)=gmean(if)+(absil*qs(id,if,2)*qs(id,if,3)+
     &   qg(id,if,2)*qg(id,if,3))*area(id)
       if(id.lt.nd) goto 5
       extt(if)=exts1(if)+exts2(if)+extg1(if)+extg2(if)
       gmean(if)=gmean(if)/(exts2(if)+extg2(if))
       if(w(if-1).ge.0.55.and.w(if).lt.0.55) then
         dvm=(w(if-1)-.55)/(w(if-1)-w(if))
         extv=extt(if-1)*(1.-dvm)+extt(if)*dvm
       endif
 5    continue
c       printing out  A_lambda/A_v   vs.   1/lambda
      write(7,'(2g13.6)') (1./w(if),extt(if)/extv,if=1,nf)
c       printing out lambda, tau_ext, tau_abs(Gra), tau_sca(Gra),
c        tau_abs(Sil), tau_sca(Sil), g (mean cos scattering angle)
      write(8,'(7g13.4)') (w(if),extt(if),extg1(if),extg2(if),
     & exts1(if),exts2(if),gmean(if),if=1,nf)
      print*,'OUTPUT:'
      print*,'1/lambda, A_lambda/A_v written to unit 7'
      print*,'lambda, tau_tot, tau_abs(Gra), tau_sca(Gra), tau_abs(Sil)'
      print*,'tau_sca(Sil), and g=<cos theta_scat> written to unit 8'
      end