! 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.
use ifcore

real(4) bess4
real(8) bess8
real(16) bess16

real(4), allocatable :: J7(:), J7R(:)   ! serie; recurenta;
real(8), allocatable :: J16(:), J16R(:)
real(16), allocatable :: J34(:), J34R(:)

integer nr_cifre(3)
character(20) filename(3)

integer numar_termeni(3, 0:100) 
common numar_termeni
character(3) char_nb
integer width
logical pressed

! Precizie:
nr_cifre(1) =7
nr_cifre(2) =16
nr_cifre(3) =34

! Read n:
write(*,'(a)') 'Se calculeaza J(m) pentru m = 0,nb. '
write(*,'(a, \)') 'Introdu "nb" (>=10): '
read(*,*) nb

iwnt =width(nb)
write(char_nb,'(I<iwnt>)') nb

! 
allocate(J7(0:nb), J16(0:nb), J34(0:nb))
allocate(J7R(0:nb), J16R(0:nb), J34R(0:nb))

J7 =0E0; J7R =0E0;
J16 =0D0; J16R =0D0;
J34 =0Q0; J34R =0Q0;

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

do m =0, 1
J7(m) = bess4(m)
J7R(m) =J7(m)
enddo

do m =0, 1
J16(m)= bess8(m)
J16R(m)= J16(m)
enddo

do m =0, 1
J34(m)= bess16(m)
J34R(m)= J34(m)
enddo

do m =2, nb
J7(m) = bess4(m)

J16(m)= bess8(m)

J34(m)= bess16(m)
enddo

! RECURENTA:
! J(m+1)=2mJ(m)-J(m-1)
! J0 si J1, si calculul: cu 3 nivele de precizie

m_do2: DO m =1, nb-1

J7R(m+1)=2.0*m*J7R(m)-J7R(m-1)

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

J34R(m+1)=2Q0*m*J34R(m)-J34R(m-1)
ENDDO m_do2


! Output Files:
filename(1) ='bess7'//'-'//trim(char_nb)//'.rez'
filename(2) ='bess16'//'-'//trim(char_nb)//'.rez'
filename(3) ='bess34'//'-'//trim(char_nb)//'.rez'

7  format(' ',i3, t8, 1pE13.6,  t26, 1pE14.6,  I7, :, 4x, 1pE10.3)
16 format(' ',i3, t8, 1pD22.15, t34, 1pD22.15, I6, :, 5x, 1pD10.3)
34 format(' ',i3, t8, 1pE40.33, t52, 1pE40.33, I6, :, 5x, 1pE10.3)


it1 =9
do io =1,3
open(io, file =filename(io), status ='unknown')

select case(io)
case(1); it2 =28 ; it3 =43; it4 =52
case(2); it2 =35; it3 =58; it4 =68
case(3); it2 =53; it3 =94; it4 =104;
end select

write(io, '(a, I3)') ' Numar cifre semnificative in J0 si J1:', nr_cifre(io)
write(io, '(/, t4,a, t<it1>,a, /)'), 'm', 'Functia Bessel Jm(1)'
    
    wr01: do j =0,1
    select case(io)
    case(1); write(io, '(t4,I1, t8, 1pE13.6, t45,I2)') j, J7(j), numar_termeni(1,j)
    case(2); write(io, '(t4,I1, t8, 1pE22.15, t60,I2)') j, J16(j), numar_termeni(2,j)
    case(3); write(io, '(t4,I1, t8, 1pE40.33, t96, I2)') j, J34(j), numar_termeni(3,j)
    end select
    enddo wr01

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

enddo

io_DO: DO io =1,3
m_DO: DO m =2, nb

SELECT CASE(nr_cifre(io))
CASE(7)
if(J7(m) /=0) then
relerr =(J7(m) -J7R(m))/J7(m)
write(1,7) m,  J7(m), J7R(m), numar_termeni(1,m), relerr
else
write(1,7) m,  J7(m), J7R(m), numar_termeni(1,m)
endif

CASE(16)
if(J16(m) /=0) then
relerr =(J16(m) -J16R(m))/bess8(m)
write(2,16) m, J16(m), J16R(m), numar_termeni(2,m), relerr
else
write(2,16) m, J16(m), J16R(m), numar_termeni(2,m)
endif

CASE(34)
if(J34(m) /=0) then
relerr =(J34(m) -J34R(m))/J34(m)
write(3,34) m, J34(m), J34R(m), numar_termeni(3,m), relerr
else
write(3,34) m, J34(m), J34R(m), numar_termeni(3,m)
endif

END SELECT
ENDDO m_DO

ENDDO io_DO

do j =1,3
write(*,*) 'Fisier de iesire:'//trim(filename(j))
close(j)
enddo

deallocate(J7, J16, J34)
deallocate(J7R, J16R, J34R)

write(*, '(/,a)') 'Apasa orice tasta'
DO WHILE (.NOT. pressed)   
pressed = PEEKCHARQQ ( )
END DO

END


