program f08_GPL_FIG
    implicit none
    character(len=100)::linebuf
    character(len=50)::fnameR1,fnameR2,fnameW


    call get_command_argument(1, fnameR1)
    call get_command_argument(2, fnameR2)

    fnameW="gpl_PLOT.txt"
    open(12,file=fnameW,status='replace')
    write(12,'(a)') 'set terminal png'
    write(12,'(a)') 'set output "fig_png_bessel.png"'
    write(12,'(a)') 'set xrange [0:10]'
    write(12,'(a)') 'set samples 1000'
    write(12,'(a)') 'set xlabel "x"'
    write(12,'(a)') 'set ylabel "J0,J1,J2,J3"'
    write(12,'(a)') 'set grid'
    write(12,'(a)') 'plot \'
    write(12,'(a)') 'besj0(x) title "(gnuplot) besj0(x)", \'
    write(12,'(a)') 'besj1(x) title "(gnuplot) besj1(x)", \'
    write(12,'(a)') '2/x*besj1(x)-besj0(x) title "(gnuplot) J2=2/x*J1-J0", \'
    write(12,'(a)') '(8/x/x-1)*besj1(x)-4/x*besj0(x) title "(gnuplot) J3=4/x*J2-J1", \'
    linebuf='"'//trim(adjustl(fnameR1))//'" using 1:2 with points title "(f08) bessel_jn(0.x)", \'
    write(12,'(a)') trim(adjustl(linebuf))
    linebuf='"'//trim(adjustl(fnameR1))//'" using 1:3 with points title "(f08) bessel_jn(1.x)", \'
    write(12,'(a)') trim(adjustl(linebuf))
    linebuf='"'//trim(adjustl(fnameR1))//'" using 1:4 with points title "(f08) bessel_jn(2.x)", \'
    write(12,'(a)') trim(adjustl(linebuf))
    linebuf='"'//trim(adjustl(fnameR1))//'" using 1:5 with points title "(f08) bessel_jn(3.x)"'
    write(12,'(a)') trim(adjustl(linebuf))
    write(12,'()')
    write(12,'(a)') 'reset'
    write(12,'(a)') 'set terminal png'
    write(12,'(a)') 'set output "fig_png_gamma.png"'
    write(12,'(a)') 'set xrange [-5:5]'
    write(12,'(a)') 'set yrange [-30:30]'
    write(12,'(a)') 'set samples 10000'
    write(12,'(a)') 'set xlabel "x"'
    write(12,'(a)') 'set ylabel "Gamma(x)"'
    write(12,'(a)') 'set grid'
    write(12,'(a)') 'set key bottom right'
    write(12,'(a)') 'plot \'
    write(12,'(a)') 'gamma(x) title "(gnuplot) gamma(x)", \'
    linebuf='"'//trim(adjustl(fnameR2))//'" using 1:2 with points title "(f08) gamma(x)", \'
    write(12,'(a)') trim(adjustl(linebuf))
    linebuf='"'//trim(adjustl(fnameR2))//'" using 1:3 with points title "Stirling approximation"'
    write(12,'(a)') trim(adjustl(linebuf))
    write(12,'(a)') '#lngamma singularity error'
    write(12,'(a)') 'quit'
    close(12)

    linebuf='gnuplot "'//trim(adjustl(fnameW))//'"'
    call execute_command_line(trim(adjustl(linebuf)))

end program f08_GPL_FIG