! Functia Bessel prin:
! - SERIE (exact): Bess
!   Insumare pana cand termenul general <= 1D-15.
!
! - RECURENTA: J(m+1)=2.0*m*J(m)-J(m-1)
!   Valori de intrare J0, J1
!   J0, J1: In 3 nivele de precizie: simpla; dubla; cvadrupla.
!
! Restul calculelor: in precizie cvadrupla (maxima disponibila).
!
real(4) bess4, b4(0:1)
real(8) bess8, b8(0:1)
real(16) bess16, J_ex(0:10), J7(0:10), J16(0:10), J34(0:10)

integer nr_cifre(3)/7,16,34/
character(20) filename(3), char_nr_cifre*2
character(12) b0_str, b1_str

integer numar_termeni(0:10) 
common numar_termeni, iprint


! J0 si J1:
iprint =0
J7(0) = bess4(0)
J7(1) = bess4(1)

J16(0)= bess8(0)
J16(1)= bess8(1)

J34(0)= bess16(0)
J34(1)= bess16(1)

! Exact = cvad.
write(*, '(/,a)') ' Calcul prin Serie:'
write(*, '(a, t19, a)') ' m   Nr_Termeni',   'Primul termen neglijat'

iprint =1
do m=0,10
J_ex(m) =bess16(m)
enddo


70  format(' ',i3, t8, 1pE13.6,  t27, 1pE14.7,  I6, :, 4x, 1pE10.3)
7   format(' ',i3, t8, 1pE13.6,  t26, 1pE14.6,  I7, :, 4x, 1pE10.3)

160 format(' ',i3, t8, 1pD22.15, t35, 1pD22.16, I6, :, 3x, 1pD10.3)
16  format(' ',i3, t8, 1pD22.15, t34, 1pD22.15, I7, :, 4x, 1pD10.3)

34 format(' ',i3, t8, 1pE40.33, t52, 1pE40.33, I7, :, 4x, 1pE10.3)

! Citeste nr. de cifre:
!nr_cifre =(/7,16,34/)

! Recurenta:
! J(m+1)=2mJ(m)-J(m-1)
! Precizie calcule: cvad.
! J0 si J1: cu 3 nivele de precizie
DO m =1,9
J7(m+1)=2.0*m*J7(m)-J7(m-1)

J16(m+1)=2D0*m*J16(m)-J16(m-1)

J34(m+1)=2Q0*m*J34(m)-J34(m-1)
ENDDO

DO I =1,3
if(nr_cifre(I) <10) then
write(char_nr_cifre, '(I1)') nr_cifre(I)
else
write(char_nr_cifre, '(I2)') nr_cifre(I)
endif

filename(I) ='bess-'//trim(char_nr_cifre)//'.rez'
open(I, file =filename(I), status ='unknown')

write(I, '(a, I3)') ' Numar cifre semnificative in J0 si J1:', nr_cifre(I)
write(I, '(/, t4,a, t9,a, /)'), 'm', 'Functia Bessel Jm(1)'
ENDDO

iprint =0

do m =0,1
write(1,70) m, J7(m), bess4(m), numar_termeni(m)

write(2,160) m, J16(m), bess8(m), numar_termeni(m)

write(3,34) m, J34(m), bess16(m), numar_termeni(m)
enddo

it1 =9

I_DO2: do I=1,3

if(I ==1) then      !7
it2 =28 ; it3 =43; it4 =52
elseif(I ==2) then  !16
it2 =35; it3 =58; it4 =68
else                !34
it2 =53; it3 =94; it4 =104;
endif

write(I, '(/, t<it1>,a, t<it2>,a, t<it3>,a, t<it4>,a)') &
         'RECURENTA', 'DIRECT (Serie)', 'Termeni', 'Er.Relativa'
ENDDO I_DO2

DO_2: DO m =2, 10

relerr =(bess4(m) -J7(m))/bess4(m)
write(1,7) m, J7(m), bess4(m), numar_termeni(m), relerr

relerr =(bess8(m) -J16(m))/bess8(m)
write(2,16) m, J16(m), bess8(m), numar_termeni(m), relerr

relerr =(bess16(m) -J34(m))/bess16(m)
write(3,34) m, J34(m), bess16(m), numar_termeni(m), relerr

ENDDO DO_2


write(*, '(/,a)') 'Fisiere de iesire:'
do i=1,3
write(*,*)  trim(filename(I))
close(I)
enddo

END


