## 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?'
print*,' Upper limit of integral?'
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```