program num_integrate
  implicit none
  double precision, external:: seval, lagint

  integer :: outputunit = 0, nofun, index
  double precision :: startvalue = 2.0, step = 0.1, q8input(11), q8results(11), errest, flag
  double precision :: splineA(11), splineB(11), splineC(11)
  double precision ::  evalX(10), evalq8(10), evalSpline(10), evalLagrange(10)
  double precision :: splineQ8Diff(10), lagrangeQ8Diff(10)
  character(len=255) outfile

  call get_command_argument(1, outfile)

  if(len_trim(outfile) > 0) then

     call process_data()

     open(newunit = outputunit, file = outfile, status="replace", action="write")
     write(outputunit, *) "Значения интеграла, вычисленные через QUANC8"
     write(outputunit, *) "(от 2 до 3 с шагом 0.1)"
     write(outputunit, *) "----------------------------------"
     write(outputunit, *) q8results
     write(outputunit, *) "----------------------------------"

     write(outputunit, *) "Коэффициенты сплайна, построенные в точках (x, a, b, c):"
     write(outputunit, *) "----------------------------------"
     do index = 1, 11
        write(outputunit, *) q8input(index), splineA(index), splineB(index), splineC(index)
     end do
     write(outputunit, *) "----------------------------------"

     write(outputunit, *) "Сравнение точного интеграла и сплайна (x, q8, spline, diff):"
     write(outputunit, *) "----------------------------------"
     do index = 1, 10
        write(outputunit, *) evalX(index), evalq8(index), evalSpline(index), splineQ8Diff(index)
     end do
     write(outputunit, *) "----------------------------------"

     write(outputunit, *) "Сравнение точного интеграла и п. Лагранжа (x, q8, lagrange, diff)"
     write(outputunit, *) "----------------------------------"
     do index = 1, 10
        write(outputunit, *) evalX(index), evalq8(index),  evalLagrange(index), lagrangeQ8Diff(index)
     end do
     close(outputunit)

  else
     print "(a)", "Укажите имя входного и выходного файла"
  end if

contains

  pure double precision function to_integrate(t)
    double precision, intent(in) :: t
    to_integrate = (1 - cos(t)) / t
  end function to_integrate

  subroutine process_data()
    double precision:: current_x

    ! генерируем значения в основных точках
    do index = 0, 10
       current_x = startvalue + index * step
       q8input(index + 1) = current_x
       call quanc8(to_integrate, 0, current_x, 0.88, 0.0, q8results(index + 1), errest, nofun, flag)
       if(flag > 0) then
          WRITE(0,*) "Error"
          write(0,*) current_x
          write(0, *) flag
          write(0, *) errest
          write(0, *) nofun
       end if
    end do

    ! генерируем значения для сплайна
    call spline(11, q8input, q8results, splineA, splineB, splineC)
    ! сравниваем сплайн с точными значениями
    do index = 0, 9
       current_x = startvalue + 0.05 + 0.1 * index
       evalX(index + 1) = current_x

       call quanc8(to_integrate, 0, current_x, 0.88, 0.0, evalq8(index + 1), errest, nofun, flag)
       if(flag > 0) then
          WRITE(0,*) "Error"
          write(0,*) current_x
          write(0, *) flag
          write(0, *) errest
          write(0, *) nofun
       end if

       evalSpline(index + 1) = seval(11, current_x, q8input, q8results, splineA, splineB, splineC)
    end do
    splineQ8Diff = abs(evalq8 - evalSpline)

    ! генерируем значения для полинома Лагранжа
    ! и сравниваем их с точными
    do index = 1, 10
       current_x = evalX(index)
       evalLagrange(index) = lagint(current_x, q8input, q8results, 11, 10)
    end do

    lagrangeQ8Diff = abs(evalq8 - evalLagrange)

  end subroutine process_data

end program num_integrate