! Graficul unei functii f(x) - pe [a,b], cu pasul h;
! Schimbarea de semn - daca exista: marcata cu un *.
! 2015: Cu intervalele radacinilor = max 100
!
subroutine grafic(f)
use work
external f
character(128) fname0, fname*260, w*1

logical exist_p/.FALSE./
character(12) a_str, b_str
character(1) op_br/'['/, cl_br/']'/, comma/','/
! 
real br(100)    !,ar(100)
integer width

iunit =1

write(*,'(a)') fun_expression

print '(a\)', ' Function depends on Parameter? (Y,y/char): '
read '(a1)', w
if(w=='Y'.or.w=='y') then
exist_p =.TRUE.

100 print '(a\)',' Input parameter "p": '
read *, p
endif

print '(a\)', ' Input a, b, h: '
read *, a,b,h
print '(a\)', ' Input file name - No extension: '
read '(a)', fname0


!fname =trim(fname0)//'.dat'
fname =Dat_Directory(:l_Dat)//'\'//trim(fname0)//'.dat'
open(iunit, file =fname)

if(exist_p) then
write(iunit,'(a,\)') trim(fun_expression) 
write(iunit, '(a, 1pE14.6, a1)') ';   p =', p, ';'

else
write(iunit,'(a)') trim(fun_expression) 
endif

write(a_str, '(F12.4)') a
a_str =adjustl(a_str)

write(b_str, '(F12.4)') b
b_str =adjustl(b_str)

write(iunit, '(a, 1pE11.4)') op_br//trim(a_str)//comma//trim(b_str)//cl_br//';   h =', h
write(iunit, *)

f_ant =f(a)

nPas =nint((b-a)/h)
no_roots =0

Pas_Do: do j=0, nPas
x =a +j*h
f_j =f(x)

! zero:
ROOT_IF: IF(f_j ==0.0) THEN
write(iunit, *) x, f_j, '**'
write(*, *) x, f_j, '**'
no_roots =no_roots +1
br(no_roots) =x     !ar(no_roots) =x-h

! schimbare de semn:
ELSE
if(sign(1.,f_ant) ==sign(1., f_j)) then
write(iunit, *) x, f_j  !f(x)

elseif(f_ant /=0.0) then
write(iunit, *) x, f_j, '*'
write(*, *) x, f_j, '*'
no_roots =no_roots +1
br(no_roots) =x     !ar(no_roots) =x-h
endif
ENDIF ROOT_IF

f_ant =f_j
enddo Pas_Do

!write(iunit, *) 'No. of roots: ', no_roots
iwn =width(no_roots)
write(iunit, '(a, I<iwn>)') 'Numar de radacini: ', no_roots
write(iunit, '(a)') 'Intervale cu radacini:'

do ir =1, no_roots
b1 =br(ir)

if(f(b1) ==0.0) then
write(b_str, '(F12.4)') b1
b_str =adjustl(b_str)
write(iunit, '(a,a)') 'Radacina: ', trim(b_str)

else
!b1 =br(ir)
a1 =b1-h
write(a_str, '(F12.4)') a1
write(b_str, '(F12.4)') b1
a_str =adjustl(a_str)
b_str =adjustl(b_str)

write(iunit, '(a)') op_br//trim(a_str)//comma//trim(b_str)//cl_br
endif
enddo

write(*,'(/, a, a)') ' Fisier de iesire: ', trim(fname)
close(iunit)

if(exist_p) then
write(*, '(/, a, \)') ' Another parameter? (Y,y/char): '
read '(a1)', w
if(w=='Y'.or.w=='y') then
iunit =iunit +1
goto 100
endif

endif


END
