How to calculate Legendre polynomials in fortran 90 ?

Legendre polynomials can be calculated using the recurrence relation:

\begin{equation}
P_n = \frac{(2n-1)}{n}.x.P_{n-1}-(n-1).P_{n-2}
\end{equation}

with $P_0=1$ et $P_1=x$

Source code in fortran 90:

program legendre_polynomials

implicit none

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i,j,n

integer :: nb_pl, nb_x

real(pr), dimension(:,:), allocatable :: pl

real(pr) :: pl_n2, pl_n1

real(pr) :: x, x_min, x_max, x_resol

real(pr) :: pi = dacos(-1.0d0)

real(pr) :: cdr

cdr = pi / 180.0

!----------------------------------------------------------------------------------------!
! User Inputs

x_min = -1.0

x_max = 1.0

x_resol = 0.01

nb_pl = 10

nb_x = int((x_max-x_min)/x_resol) + 1

allocate(pl(nb_pl,nb_x))

!----------------------------------------------------------------------------------------!
! Legendre Polynomials

x = x_min

do i = 1, nb_x

    pl(1,i) = 1.0

    pl(2,i) = x

    pl_n2 = pl(1,i)

    pl_n1 = pl(2,i)

    do j = 3 , nb_pl

        n = j - 1

        pl(j,i) = ( 1.0 / n ) * ( ( 2.0 * n - 1.0 ) * x * pl_n1 - ( n - 1.0 ) * pl_n2 )

        pl_n2 = pl_n1

        pl_n1 = pl(j,i)

    end do

x = x + x_resol

end do

!----------------------------------------------------------------------------------------!
! Print results

open(1,file="legendre_polynomials.txt")

x = x_min

do i = 1, nb_x

write(1,*) x,(pl(j,i),j=1,nb_pl)

x = x + x_resol

end do

close(1)

!----------------------------------------------------------------------------------------!

deallocate(pl)

end

Note: output file can be visualized using the python module called matplotlib:

import matplotlib.pyplot as plt
import numpy as np

data = np.loadtxt("../../legendre_polynomials.txt")

for i in [1,2,3,4,5,6]:
        plt.plot(data[:,0],data[:,i],label=str(i))

plt.legend()

plt.savefig('legendre_polynomials.png')

plt.show()

References