! Main pt. Halley
Program M_Hal
use work

external SUB, Deriv_Pol

!logical is_pol, not_Pol
!p, is_pol: parametru pt. f(x); f(x) =polinom;

logical exist_p

character(80) text(-2:1), dummy
character(260) fname0, fname
integer width

character(132) f_text(0:2), f_source
logical new_rad, new_p, new_tol

character(1) w1, w
character(20) w_titlu; integer status;
character(8) eps_char

eps_num =1E-12

write(eps_char, '(1pE8.1)') eps_num
text(-2)=' ** Numitor < '//eps_char
text(0) =' * Toleranta EPS s-a atins:'
text(1) =' ** Numarulul limita de iteratii atins:'
text(-1)=' ** Divergenta'

10 format(a, \)

Print '(a, \)', ' Polinom (Y,y/char): '
read(*, '(a1)')  w1
is_pol = (w1 =='Y' .or. w1 =='y')
not_Pol =.NOT. is_Pol

IsPol_IF: IF(is_pol) THEN

w_titlu ='Select Data File'C
call GetFile(fname, status, w_titlu)
OPEN(UNIT=1, FILE=fname, STATUS='OLD')

read(1,'(a80)') dummy
READ(1,*) N

allocate(coef(0:n))
read(1,*) coef(n:0:-1)

! Derivate:
allocate(coef_1(0:n-1))
allocate(coef_2(0:n-2))
call Deriv_Pol()

exist_p =.FALSE.
write(1,'(a)')  ! Before Header

GOTO 150
ENDIF IsPol_IF

f_source ='sub.f90'
call F_Expr(f_source, f_text)
write(*, '(1x, a)') trim(f_text(0))
write(*, '(1x, a)') trim(f_text(1))
write(*, '(1x, a)') trim(f_text(2))

print 10, ' Nume fisier rezultate (Fara extensie): '
read *, fname0
fname =trim(fname0)//'.rez'

open(1, file =fname, form ='formatted', status ='unknown')

print 10,' Functia f(x) include parametrul "p" (Y / char) ?: '
read'(a)',w
exist_p =(w =='Y' .or. w =='y')

150 CONTINUE    ! IF_Pol
! Write header:
write(1, '(a)') ' Metoda HALLEY'

NotPol_IF: IF(not_Pol) THEN
write(1, '(1x, a)') trim(f_text(0))
write(1, '(1x, a)') trim(f_text(1))
write(1, '(1x, a)') trim(f_text(2))
ENDIF NotPol_IF

nrad =1
new_p =.false.
new_tol =.true.

100 if(exist_p) then
print 10,' p=? : '
read*,p
endif

200 continue
if(new_p) then
print '(a, \)',' Alt "eps" si "nlim" ? (Y/ char): '
read '(a)',w
new_tol =w=='Y' .or. w =='y'
endif

if(new_tol) then
print 10,' Introdu eps, nlim : '
read *, eps, nlim0
endif

210 CONTINUE    ! New rad:
print 10,' Introdu x0 si cod tiparire (0/#0): '
read *, x00, kod0

iwn =width(nrad); iwnlim0 =width(nlim0);
write(1, '(/, a, I<iwn>, a1)') ' Radacina nr.', nrad

if(exist_p) then
write(1,*) 'p:', p
endif

write(1,*) 'Aproximatia initiala:', x00
write(1, '(a, 1pG10.3, a, I<iwnlim0>, /)') ' EPS:', EPS, '; NLIM: ', nlim0


x0 =x00; nlim=nlim0; kod=kod0

n_div =0; rad =x0;

call HALLEY(SUB, x0, eps, nlim, rad, kod, n_div)

if(not_Pol) then
call sub(f, f1, f2, rad)

else
f  =polval(n, coef, rad)
endif

!if(kod ==-1) then
!print '(/,a)', trim(text(1))
!write(1,'(/,a,\)'), trim(text(1))
!endif

! Explore divergence:
IF(n_div /=0) THEN
iwn_div =width(n_div)
write(1,'(/, a, I<iwn_div>)') ' * Nr. de divergente locale: ', n_div

ENDIF

if(kod /=2) then
print '(/a)', trim(text(kod))
write(1, '(/, a)'), trim(text(kod))
endif

if(kod ==0) then
write(1,*)
write(1,*) ' Radacina =', rad,'; f(Radacina) =', f
write(*,*) ' Radacina =', rad,'; f(Radacina) =', f

else
write(1,*)
write(1,*) ' Aproximatia curenta: x =', rad,'; f(x) =', f
write(*,*) ' Aproximatia curenta: x =', rad,'; f(x) =', f
endif

iwnlim =width(nlim)
write(1, '(a, I<iwnlim>)') '  Numar de iteratii = ', nlim

! Alta radacina:
print*
print '(/, a, \)',' Alta radacina ? (Y/ char): '
read '(a)',w

new_rad =w=='Y' .or. w =='y'
new_p =.false.

IF(new_rad) then
nrad =nrad+1

	if(exist_p) then
	print '(a, \)',' Alt "p" ? (Y/ char): '
	read '(a)',w
	new_p =w=='Y' .or. w =='y'

	if(new_p) goto 100
	endif


print '(a, \)',' Alt "eps" si "nlim" ? (Y/ char): '
read '(a)',w
new_tol =w=='Y' .or. w =='y'

if(new_tol) then
goto 200
else
goto 210
endif

ENDIF

close(1)

print '(/, a, a)', ' * Fisier de rezultate: ', trim(fname)

END program M_Hal

