      implicit none

      integer n
      parameter(n=1950)
      real*8 lambda(n),omega(n),g(n),qext(n)

      character*128 fn

      integer i
      integer i1,i2

      real*8 aveom0,aveg0,aveqe0
      real*8 aveom1,aveg1,aveqe1
      real*8 aveom2,aveg2,aveqe2
      real*8 aveom3,aveg3,aveqe3

      real*8 qe067
      real*8 ratio

c     obtained from vis.f
      qe067=3.04000000000000
      ratio=2.0d0


      fn="gmod.data"
      open(50,file='./'//fn,status='unknown')
      do i=1,n
         read(50,*) lambda(i),g(i)
      enddo
      close(50)
      fn="omegamod.data"
      open(50,file='./'//fn,status='unknown')
      do i=1,n
         read(50,*) lambda(i),omega(i)
      enddo
      close(50)
      fn="qextmod.data"
      open(50,file='./'//fn,status='unknown')
      do i=1,n
         read(50,*) lambda(i),qext(i)
      enddo
      close(50)


      do i=2,n-1
         if(lambda(i).ge.9.0d-6) go to 100
      enddo
 100  continue
      write(*,*) 'check'
      write(*,*) lambda(i)
      aveg0=g(i)
      aveom0=omega(i)
      aveqe0=qext(i)
      write(*,*) 'at 9 micron'
      write(*,*) '    ',aveom0,aveg0,aveqe0
      write(*,*) 

      i1=1
      do i=2,n-1
         if(lambda(i).ge.1.0d0/900.0d2) go to 110
c         if(lambda(i).ge.11.6d-6) go to 110
      enddo
 110  continue
      i2=i
      call average(n,lambda,g,omega,qext,i1,i2,aveom1,aveg1,aveqe1)
      write(*,*) 'average from 5 to 11 micron'
      write(*,*) '    ',aveom1,aveg1,aveqe1
      write(*,*) '            Qe_bar/Qe_9micron',aveqe1/aveqe0
      write(*,*) '            Qe_bar/Qe_0.67micron',
     &     aveqe1*(qe067/ratio/aveqe0)/qe067


      i1=i2
      do i=2,n-1
         if(lambda(i).ge.1.0d0/400.0d2) go to 120
c         if(lambda(i).ge.20.0d-6) go to 120
      enddo
 120  continue
      i2=i
      call average(n,lambda,g,omega,qext,i1,i2,aveom2,aveg2,aveqe2)
      write(*,*) 'average from 11 to 25 micron'
      write(*,*) '    ',aveom2,aveg2,aveqe2
      write(*,*) '            Qe_bar/Qe_9micron',aveqe2/aveqe0
      write(*,*) '            Qe_bar/Qe_0.67micron',
     &     aveqe2*(qe067/ratio/aveqe0)/qe067


      i1=i2
      i2=n
      call average(n,lambda,g,omega,qext,i1,i2,aveom3,aveg3,aveqe3)
      write(*,*) 'average from 25 to 200 micron'
      write(*,*) '    ',aveom3,aveg3,aveqe3
      write(*,*) '            Qe_bar/Qe_9micron',aveqe3/aveqe0
      write(*,*) '            Qe_bar/Qe_0.67micron',
     &     aveqe3*(qe067/ratio/aveqe0)/qe067


      write(*,*)
      write(*,*)
      i1=1
      i2=n
      call average(n,lambda,g,omega,qext,i1,i2,aveom3,aveg3,aveqe3)
      write(*,*) 'average from 5 to 200 micron'
      write(*,*) '    ',aveom3,aveg3,aveqe3
      write(*,*) '            Qe_bar/Qe_9micron',aveqe3/aveqe0
      write(*,*) 'This agrees with that presented in Forget [1998].'

      stop
      end

c******************************************************************************
      subroutine average(n,lambda,g,omega,qext,i1,i2,aveome,aveg,aveqex)

      implicit none

      integer n
      real*8 lambda(n),omega(n),g(n),qext(n)
      integer i1,i2
      real*8 aveome,aveg,aveqex

      real*8 sum1,sum2,sum3
      real*8 temp

      integer i

      real*8 plfunc2

      aveome=0.0d0
      sum1=0.0d0
      aveg=0.0d0
      sum2=0.0d0
      aveqex=0.0d0
      sum3=0.0d0

      temp=200.0d0
      do i=i1+1,i2
         aveome=aveome
     &        +(plfunc2(lambda(i-1),temp)*qext(i-1)*omega(i-1)
     &        +plfunc2(lambda(i),temp)*qext(i)*omega(i))
     &        *(lambda(i)-lambda(i-1))
     &        /2.0d0
         sum1=sum1
     &        +(plfunc2(lambda(i-1),temp)*qext(i-1)
     &        +plfunc2(lambda(i),temp)*qext(i))
     &        *(lambda(i)-lambda(i-1))
     &        /2.0d0

         aveg=aveg
     &        +(plfunc2(lambda(i-1),temp)*qext(i-1)*omega(i-1)*g(i-1)
     &        +plfunc2(lambda(i),temp)*qext(i)*omega(i)*g(i))
     &        *(lambda(i)-lambda(i-1))
     &        /2.0d0
         sum2=sum2
     &        +(plfunc2(lambda(i-1),temp)*qext(i-1)*omega(i-1)
     &        +plfunc2(lambda(i),temp)*qext(i)*omega(i))
     &        *(lambda(i)-lambda(i-1))
     &        /2.0d0

         aveqex=aveqex
     &        +(plfunc2(lambda(i-1),temp)*qext(i-1)
     &        +plfunc2(lambda(i),temp)*qext(i))
     &        *(lambda(i)-lambda(i-1))
     &        /2.0d0
         sum3=sum3
     &        +(plfunc2(lambda(i-1),temp)
     &        +plfunc2(lambda(i),temp))
     &        *(lambda(i)-lambda(i-1))
     &        /2.0d0
      enddo

      aveome=aveome/sum1
      aveg=aveg/sum2
      aveqex=aveqex/sum3



c      sum1=0.0d0
c      do i=1,10000
c         sum2=1000.0d0*1.0d-6/dble(10000)*dble(i)
c         sum3=1000.0d0*1.0d-6/dble(10000)*dble(i+1)
c         sum1=sum1+(plfunc2(sum2,temp)+plfunc2(sum3,temp))
c     &        *(sum3-sum2)/2.0d0
c      enddo
c      sum1=sum1*3.141592
c      write(*,*) sum1,5.6698d-8*temp**4



      return
      end

c******************************************************************************
      real*8 function plfunc(x,t)

      implicit none

      real*8 x,t
      real*8 c,planc,boltz

      c=2.99792458d8
      planc=6.6260755d-34
      boltz=1.380658d-23

      plfunc=2.0d0*planc*c*c*x*x*x
     &     /(dexp(planc*c*x/(boltz*t))-1.0d0)

      return
      end

c******************************************************************************
      real*8 function plfunc2(x,t)

      implicit none

      real*8 x,t
      real*8 c,planc,boltz

      c=2.99792458d8
      planc=6.6260755d-34
      boltz=1.380658d-23

      plfunc2=2.0d0*planc*c*c
     &     /(x*x*x*x*x*(dexp(planc*c/(boltz*x*t))-1.0d0))

      return
      end
