module gamma
  implicit none
contains
  subroutine calc_gamma(beta,Td,precision,gamma)
    implicit none
    real,intent(in)                     ::beta
    real,intent(in)                     ::Td
    real,intent(in)                     ::precision
    real,intent(inout)                  ::gamma
    !Variables locales
    real                                ::maxfseconde
    integer                             ::n,k

    maxfseconde=0.
    maxfseconde=max(maxfseconde,2.*beta*abs(-1.+2*beta*(0.-Td)**2))
    maxfseconde=max(maxfseconde,2.*beta*abs(-1.+2*beta*(1.-Td)**2))
    
    !evaluation du nombre de points
    n=floor(sqrt(maxfseconde/(12.*precision)))+1


    !calcul de l'integrale
    gamma=0.
    gamma=func(beta,Td,0.)
    gamma=gamma+func(beta,Td,1.)
    gamma=gamma/2.

    do k=1,(n-1)
       gamma=gamma+func(beta,Td,(k*1.)/n)
    enddo

    gamma=gamma/n
    
  end subroutine calc_gamma
  
  function func(beta,Td,x) result(res)
    implicit none
    real,intent(in)              ::beta,Td,x
    real                         ::res

    res=exp(-beta*(x-Td)**2)


  end function func
  
end module gamma
