! Main pt. Ordin 15 (Soleymani & Sharifi, 2012)
! Ecuatie: f(x) =0
! - Daca f este polinom: Programul cere un fisier de intrare, cu structura:
!	1) Linie_text: titlul problemei (<= 80 caractere)
!	2) n: gradul polinomului
!	3) a(n),a(n-1), ...,a(0): coeficientii polinomului (în aceasta ordine)
! - Daca f nu e polinom:
!   f si f' se definesc in unitatea de program FUN (f =f(x); f1 =f'(x))
!
! Aproximatia initiala trebuie sa fie suf. de aproape de radacina !
!
Program Main_Ordin_15
use work

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

logical exist_p

character(80) text(-5:3), dummy
character(260) fname0, fname
integer width

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

character(1) w1, w
character(20) w_titlu; integer status;
character(8) eps_num_char, eps_f_char, num_char


ULP(x) =2.0**(exponent(x) -24)

eps_num =1E-12
eps_f =1E-16

write(eps_num_char, '(1pE8.1)') eps_num
write(eps_f_char, '(1pE8.1)') eps_f

text(-5)=' ** w = z !'
text(-4)=' ** z = y !'
text(-3)=' ** y = x !'
text(-2)=' ** Numitor < '//eps_num_char
text(-1)=' ** Divergenta'
text(0) =' * Toleranta EPS s-a atins:'
text(1) =' ** Numarulul limita de iteratii atins:'
text(2) =' ** |f(x)| <= '//eps_f_char
!text(3) =' ** Proces stationar'

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))
call Deriv_Pol()

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

GOTO 150
ENDIF IsPol_IF

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

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

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

! Write header:
write(1, '(/, a)') ' Metoda - Ordinul 15'

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 *, eps0, nlim0
endif


210 CONTINUE    ! New rad:
print 10,' Introdu x0 si Cod tiparire (0/#0; >=2 pt.pas): '
read *, x00, kod0

if(eps0 <ULP(x00)) then
eps =ULP(x00)
write(*,*) ' * EPS ajustat la:', eps

else
eps =eps0
endif

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
if(eps ==eps0) then
write(1, '(a, 1pG10.3, a, I5/)') ' EPS:', EPS, ';   LNIT =', nlim0

else
write(1, '(a, 1pG10.3, a, I5/)') ' EPS (ajustat):', EPS, ';   LNIT =', nlim0
endif


x0 =x00; nlim=nlim0; kod=kod0

n_div =0; rad =x0; num =1.;

call Ordin_15(x0, eps, nlim, rad, kod, n_div)

call compute_f(rad, fa)
! Infinit sau NaN:
IF(kod ==-10) THEN
write(1,*) ' Aproximatia curenta: x =', rad,'; f(x) =', fa
GOTO 800
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 .OR. kod ==2) then
write(1,*) char(13)//char(10)
write(1,*) ' Radacina =', rad,'; f(Radacina) =', fa
write(*,*) ' Radacina =', rad,'; f(Radacina) =', fa

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

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

! Alta radacina:
800 continue
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 Main_Ordin_15

