FORTRAN

ne.gif (2791 bytes) NE406 Radiation Protection and Shielding
FORTRAN coding to do an S32 Gauss-Legendre integration

      program group
C***********************************************************************
C                                                                      *
C     Find the integral of a function over a given range using         *
C     a Gauss-Legendre procedure of order 32.                          *
C                                                                      *
C     This can be adapted to other functions by changing FUNCTION F.   *
C                                                                      *
C***********************************************************************
      integer nsn
      real*8 snwt(32),snx(32)
      data nsn/32/
      data snx/.001368073,.0071942277,.0176188818,.0325469453,.
     *0518394344,.0753161897,.1027581033,.1339089403,.1684778666,.
     *2061421214,.2465500455,.2893243619,.3340656989,.3803563189,.
     *4277640192,.4758461672,.5241538328,.5722359808,.6196436811,.
     *6659343011,.7106756381,.7534499545,.7938578786,.8315221334,.
     *8660910597,.8972418967,.9246838103,.9481605656,.9674530547,.
     *9823811182,.9928057723,.9986319268/
      data snwt/.003509312,.0081372115,.0126959771,.0171370005,.
     *0214179378,.0254990642,.0293420289,.0329111118,.0361728923,.
     *0390969435,.0416559573,.0438260415,.0455869341,.0469221942,.
     *0478193546,.0482700387,.0482700387,.0478193546,.0469221942,.
     *0455869341,.0438260415,.0416559573,.0390969435,.0361728923,.
     *0329111118,.0293420289,.0254990642,.0214179378,.0171370005,.
     *0126959771,.0081372115,.0035093120/
   10 print*,' Lower limit of integral?'
      read(*,*)xl
      print*,' Upper limit of integral?'
      read(*,*)xu
      val=0.
      do i=1,nsn
        x=xl+(xu-xl)*snx(i)
        val=val+snwt(i)*f(x)
      enddo
      val=val*(xu-xl)
      print*,' Value = ',val
      go to 10
      stop
      end
C***********************************************************************
C***********************************************************************
      function f(x)
      f=0.64*exp(-x/1.175)*sinh(sqrt(1.0401*x))
      return
      end
C***********************************************************************
C***********************************************************************
      function sinh(x)
      sinh=(exp(x)-exp(-x))*0.5
      return
      end