c     Copyright (C) 2008-2014 Vincent Eymet
c
c     KDISTRIBUTION is free software; you can redistribute it and/or modify
c     it under the terms of the GNU General Public License as published by
c     the Free Software Foundation; either version 3, or (at your option)
c     any later version.

c     KDISTRIBUTION is distributed in the hope that it will be useful,
c     but WITHOUT ANY WARRANTY; without even the implied warranty of
c     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c     GNU General Public License for more details.

c     You should have received a copy of the GNU General Public License
c     along with KDISTRIBUTION; if not, see <http://www.gnu.org/licenses/>
c
      subroutine write_quad_results(kpts,kind,Nq,fixed_g,g,w)
      implicit none
      include 'max.inc'
      include 'formats.inc'
c
c     Purpose: to write quadrature results
c     
c     Inputs:
c       + kpts: number of fixed values of "g"
c       + kind: quadrature type
c       + Nq: quadrature order
c       + fixed_g: array of "kpts" fixed values of "g"
c       + g: values of the "g" cumulated function
c       + w: values of quadrature weights
c
      integer kpts,kind,Nq,quad,i,n,kdok
      double precision fixed_g(2),g(Nqmx),w(Nqmx)
      double precision sum_w
      double precision sumwigin(1:2*Nqmx)
      double precision uns1pn(1:2*Nqmx),min
      double precision errkd(1:2*Nqmx),errkd_max
      parameter(errkd_max=1.0D-5)
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine write_quad_results'
      
      sum_w=0.0D+0
      do n=1,Nq
         sum_w=sum_w+w(n)
      enddo

c     validity check for the k-distribution
      kdok=1
      do n=1,2*Nq
         sumwigin(n)=0.0D+0
         do quad=1,Nq
            sumwigin(n)=sumwigin(n)+w(quad)*g(quad)**n
         enddo                  ! quad
         uns1pn(n)=1.0D+0/(1.0D+0+dble(n))
         min=dabs(sumwigin(n))
         if (dabs(uns1pn(n)).lt.min) then
            min=dabs(uns1pn(n))
         endif
         errkd(n)=dabs(sumwigin(n)-uns1pn(n))/min
         if (errkd(n).gt.errkd_max) then
            kdok=0
         endif
      enddo                     ! n

      open(11,file='./results/quadrature_results.txt')
      write(11,*) 'Quadrature type required:',kind
      do i=1,kpts
         write(11,*) 'Using fixed g value:',fixed_g(i)
      enddo
      write(11,*) 'Quadrature order:',Nq
      write(11,*)
      write(11,*) 'Quadrature abscissas and weights: g / w'
      do quad=1,Nq
         write(11,*) g(quad),w(quad)
      enddo
      write(11,*)
      write(11,*) 'Sum of weights=',sum_w
      write(11,*)
      if (kdok.eq.1) then
         write(11,*) 'K-distribution is coherent'
      else
         write(11,*) 'K-distribution error:'
      endif
      write(11,*)
      write(11,*) 'n / sum(wi*gi^n) / 1/(1+n) / err'
      do n=1,2*Nq
         write(11,51) n,sumwigin(n),uns1pn(n),errkd(n)
      enddo
      close(11)
      
      return
      end