From b6d8dd5a6f686ddce165472980d1bb4e4c05d909 Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Fri, 13 Jun 2014 09:47:20 +0000
Subject: [PATCH] simpson quadrature replaced with simpler rule for unevenly
 spaced intervals (computation of <rho>_{P,J} and d<rho>_{P,J)

---
 src/gray-externals.f | 48 +++++++++++++++++++-------------------------
 src/grayl.f          | 19 ++++++++++++++++++
 2 files changed, 40 insertions(+), 27 deletions(-)

diff --git a/src/gray-externals.f b/src/gray-externals.f
index e4629bb..d5c3806 100644
--- a/src/gray-externals.f
+++ b/src/gray-externals.f
@@ -6167,36 +6167,30 @@ c
         end do
       end do
  
-      h=1.0d0/dble(nd-1)    !!!!!!!!!! equi-spaced grid expected in
-      rhotpav=0.0d0         !!!!!!!!!! simpson subroutine
-      drhotpav=0.0d0        !!!!!!!!!! wrong integrals until fixed
-      rhotjav=0.0d0         !!!!!!!!!! (d)rhotpav,rhotjav,(d)rhotjava affected
+      rhotpav=0.0d0
+      rhot2pav=0.0d0
+      rhotjav=0.0d0
       rhotjava=0.0d0
       rhot2java=0.0d0
      
-      fi=dpdv/h
-      call simpson (nd,h,fi,spds)
-      fi=rhotv*dpdv/h
-      call simpson (nd,h,fi,rhotpav)
-      fi=rhotv*rhotv*dpdv/h
-      call simpson (nd,h,fi,rhot2pav)
-      rhotpav=rhotpav/spds
-      rhot2pav=rhot2pav/spds
+c     here dpdv is still dP (not divided yet by dV)
+      spds=sum(dpdv(1:nd))
+      if (spds.gt.0.0d0) then
+        rhotpav=sum(rhotv(1:nd)*dpdv(1:nd))/spds
+        rhot2pav=sum(rhotv(1:nd)*rhotv(1:nd)*dpdv(1:nd))/spds
+      end if
 
+c     here ajphiv is still dI (not divided yet by dA)
       if (ieccd.ne.0) then      
-      fi=ajphiv/h
-      call simpson (nd,h,fi,sccs)
-      fi=rhotv*ajphiv/h
-      call simpson (nd,h,fi,rhotjav)    
-      fi=abs(ajphiv)/h
-      call simpson (nd,h,fi,sccsa)
-      fi=rhotv*abs(ajphiv)/h
-      call simpson (nd,h,fi,rhotjava)
-      fi=rhotv*rhotv*abs(ajphiv)/h
-      call simpson (nd,h,fi,rhot2java)     
-      rhotjav=rhotjav/sccs
-      rhotjava=rhotjava/sccsa
-      rhot2java=rhot2java/sccsa
+        sccs=sum(ajphiv(1:nd))
+        if (sccs.gt.0.0d0) then
+          rhotjav=sum(rhotv(1:nd)*ajphiv(1:nd))/sccs
+        end if
+        sccsa=sum(abs(ajphiv(1:nd)))
+        if (sccsa.gt.0.0d0) then
+          rhotjava=sum(rhotv(1:nd)*abs(ajphiv(1:nd)))/sccsa
+          rhot2java=sum(rhotv(1:nd)*rhotv(1:nd)*abs(ajphiv(1:nd)))/sccsa
+        end if
       end if
 
 c     factor sqrt(8)=2 sqrt(2) to match with full width
@@ -6205,7 +6199,7 @@ c     of gaussian profile
       drhotpav=sqrt(8.d0*drhot2pav)
       drhot2java=rhot2java-rhotjava**2
       drhotjava=sqrt(8.d0*drhot2java)
-      
+
       spds=0.0d0
       sccs=0.0d0
       do i=1,nd
@@ -6221,7 +6215,7 @@ c     of gaussian profile
       facjs=1.0d0
       if(spds.gt.0.0d0) facpds=pabs/spds
       if(sccs.ne.0.0d0) facjs=currt/sccs
-
+      
       do i=1,nd
         dpdv(i)=facpds*dpdv(i)
         ajphiv(i)=facjs*ajphiv(i)
diff --git a/src/grayl.f b/src/grayl.f
index 3de8025..761950a 100644
--- a/src/grayl.f
+++ b/src/grayl.f
@@ -257,6 +257,25 @@ c
       end subroutine simpson
 c
 c
+c
+      subroutine trapezoid(n,xi,fi,s)
+c     subroutine for integration with the trapezoidal rule.
+c     fi: integrand f(x); xi: abscissa x;
+c     s: integral Int_{xi(1)}^{xi(n)} f(x)dx
+        implicit none
+        integer n
+        real*8 xi(n),fi(n)
+        real*8 s
+        integer i
+c
+        s  = 0.0d0
+        do i = 1, n-1
+          s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i))
+        end do
+        s = 0.5d0*s
+      end subroutine trapezoid
+c
+c
 c     spline routines: begin
 c
 c