! Programa para calcular los coeficientes de Lagrange y sus
! derivadas primera y segunda. Si definimos los coeficientes como
! b_s^j(\alpha) con j entero, s multiplo de 1/2, y \alpha < 1, 
! tenemos que:
!
! 0.5 b_s^j=[s(s+1)...(s+j-1)]/[1.2.3...j] \alpha^j x
!           x {1+[s(s+j)]/[1.(j+1)] \alpha^2 + [s(s+1)(s+j)(s+j+1)]/
!           / [1.2.(j+1)(j+2)] \alpha^4 + ...}
! 
! Si j=0, el primer factor es 1. Las derivadas se calculan 
! respecto de \alpha.
!
! rgh - Septiembre 2008
! 
program coelag
implicit none
  integer, parameter ::d=selected_real_kind(15,200)
  integer :: i,k,jj,j
  real(d) :: ss,dif,alpha,fac,bb,dbb,ddbb
  real(d),dimension(20) :: coe
 
! ingreso de parametros
! cuidado que jj es entero
!
  write(*,*)'Parametros s, j, y alpha: '
  read(*,*)ss,jj,alpha
  dif=1._d
  coe=0._d
  coe(1)=1._d
  i=1

! calcula los coeficientes de la serie hasta que
! los terminos difieren en menos de 1e-10
!
  bucle: do
    coe(i+1)=coe(i)*(ss+dble(i)-1._d)*(ss+dble(jj)+dble(i)-1._d)/ &
      (dble(i)*(dble(jj)+dble(i)))
    dif=abs(coe(i+1)*alpha**(2*i)-coe(i)*alpha**(2*(i-1)))
    i=i+1
    if((dif < 1.d-10).or.(i == 20))exit bucle
  enddo bucle

  i=i-1

! calcula el primer factor
!
  fac=1._d
  if(jj /= 0)then
    fac=ss
    do k=1,jj-1
      fac=fac*(ss+dble(k))/(dble(k)+1)
    enddo
  endif

  write(*,*)'Numero de terminos: ',i

! calcula el coeficiente de Lagrange y
! sus derivadas primera y segunda
!  
  bb=0._d
  dbb=0._d
  ddbb=0._d
  do j=0,i
    bb=bb+fac*coe(j+1)*alpha**(2*j+jj)
    dbb=dbb+dble(2*j+jj)*fac*coe(j+1)*alpha**(2*j+jj-1)
    ddbb=ddbb+dble(2*j+jj-1)*dble(2*j+jj)*fac*coe(j+1)*alpha**(2*j+jj-2)
  enddo
  write(*,*)'Coef. Lagrange b(',ss,',',jj,'): ',2._d*bb
  write(*,*)'Primera derivada de b(',ss,',',jj,'): ',2._d*dbb
  write(*,*)'Segunda derivada de b(',ss,',',jj,'): ',2._d*ddbb

end program coelag