      implicit none

      integer n
      parameter(n=22+2)
      real*8 lambda(n),omega(n),g(n),qext(n)
      real*8 om1,om2,g1,g2,qe1,qe2
      real*8 aveome,aveg,aveqex
      real*8 sum1,sum2,sum3
      real*8 temp

      character*128 line
      character*128 fn

      integer i

      real*8 plfunc

      fn="visparam.data"
      open(50,file='./'//fn,status='unknown')
      read(50,'(a)') line
      do i=2,n-1
         read(50,*) lambda(i),omega(i),g(i),qext(i)
         lambda(i)=lambda(i)*1.0d-6
      enddo
      close(50)

      i=1
      lambda(i)=0.1d0*1.0d-6
      omega(i)=(omega(3)-omega(2))/(lambda(3)-lambda(2))
     &     *(lambda(i)-lambda(2))+omega(2)
      omega(i)=max(omega(i),0.0d0)
      g(i)=(g(3)-g(2))/(lambda(3)-lambda(2))
     &     *(lambda(i)-lambda(2))+g(2)
      g(i)=max(g(i),0.0d0)
      qext(i)=(qext(3)-qext(2))/(lambda(3)-lambda(2))
     &     *(lambda(i)-lambda(2))+qext(2)
      qext(i)=max(qext(i),0.0d0)

      i=n
      lambda(i)=5.0d0*1.0d-6
      omega(i)=(omega(n-1)-omega(n-2))/(lambda(n-1)-lambda(n-2))
     &     *(lambda(i)-lambda(n-2))+omega(n-2)
      omega(i)=max(omega(i),0.0d0)
      g(i)=(g(n-1)-g(n-2))/(lambda(n-1)-lambda(n-2))
     &     *(lambda(i)-lambda(n-2))+g(n-2)
      g(i)=max(g(i),0.0d0)
      qext(i)=(qext(n-1)-qext(n-2))/(lambda(n-1)-lambda(n-2))
     &     *(lambda(i)-lambda(n-2))+qext(n-2)
      qext(i)=max(qext(i),0.0d0)


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

      temp=6000.0d0
c      do i=2,n
      do i=3,n
         write(*,*) lambda(i)

         aveome=aveome
     &        +(plfunc((1.0d0/lambda(i-1)),temp)*qext(i-1)*omega(i-1)
     &        +plfunc((1.0d0/lambda(i)),temp)*qext(i)*omega(i))
     &        *(1.0d0/lambda(i-1)-1.0d0/lambda(i))
     &        /2.0d0
         sum1=sum1
     &        +(plfunc((1.0d0/lambda(i-1)),temp)*qext(i-1)
     &        +plfunc((1.0d0/lambda(i)),temp)*qext(i))
     &        *(1.0d0/lambda(i-1)-1.0d0/lambda(i))
     &        /2.0d0

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

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

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

      write(*,*) 'at 0.67 micron'
      i=8
      write(*,*) omega(i),g(i),qext(i)
      write(*,*) 'average from 0.2 to 5 micron'
      write(*,*) aveome,aveg,aveqex
      write(*,*) 'ratio of extinction efficiency factor'
      write(*,*) aveqex/qext(i)


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


      stop
      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
