module params
 real :: L, l1, V_B(3) !Parameters
 DATA V_B/3060.0, 2866.0, 1337.0/
 real :: c != 0.5 * 0.000001
 integer :: nofun !Starting values for QUANC8
 real :: a = 1.0, b = 2.0, relerr = 1.E-06, abserr = 0.0, result, errest, flag, tol !Starting values for QUANC8
 real :: U_STATIC(11), U_DYNAMIC(11), Q_U_STATIC=0.0, Q_U_DYNAMIC=0.0 !Vectors for U and their quadratic sum
 DATA U_STATIC/-1.000, 7.777, 12.017, 10.701, 5.407, -0.843, -5.159, -6.015, -3.668, 0.283, 3.829/ !Input data for U
 real :: T, TOUT, Y(3), R_WORK(21) !For RKF45
 integer :: IFLAG, IWORK(5), NEQN = 3 !For RKF45
 real :: E1, E2, R1, R2, R3,  L3, i1, i3, U, R, fm1 = 0.0, fm_a = 0.5E-6, fm_b = 2.0E-6
end module params

program kurs
 use params
 real :: IPVT(3), M_A(3,3), WORK(3), cond !For DECOMP and SOLVE
 integer i
interface

  real function for_l(x)
  real x
  end function for_l

  real function sub_fmin(x)
  real x
  end function sub_fmin

  real function for_l1(x)
  real x
  end function for_l1

  subroutine diff(T1, Y1, YP)
    real T1, Y1(3), YP(3)
  end subroutine diff

end interface

do i = 1, 11
  Q_U_STATIC = Q_U_STATIC + U_STATIC(i)**2 !Cycle to find squares of the matrix elements
end do

write(*,*) "Q_U_STATIC = ", Q_U_STATIC
  call QUANC8(for_l, a, b, abserr, relerr, result, errest, nofun, flag)
  L = 0.4674158 * result
  write(*,*) "L = ", L

  l1 = 5.718088 * ZEROIN(0.0, 1.0, for_l1, relerr)
  write(*,*) "l1 = ", l1

  DATA M_A/53.0, 46.0, 20.0, 46.0, 50.0, 26.0, 20.0, 26.0, 17.0/

  call DECOMP(3, 3, M_A, cond, IPVT, WORK)
  call SOLVE(3, 3, M_A, V_B, IPVT)
  write(*,*) "V_B = ", V_B
  tol = 1.0E-08
  a = 0.5 * 0.000001
  b = 2.0 * 0.000001

  R = V_B(1)
  R2 = V_B(2)
  E1 = l1
  E2 = V_B(3)
  L1 = L
  L3 = L
  R1 = R
  R3 = R

  C = FMIN(fm_a, fm_b, sub_fmin, tol)
 ! C = 1.0E-6
  write(*,*) "========================="
  write(*,*) "L  = ", L
  write(*,*) "E1 = ", E1
  write(*,*) "R  = ", R
  write(*,*) "R2 = ", R2
  write(*,*) "E2 = ", E2
  write(*,*) "========================"
  write(*,*) "FMIN = ", C

  T = 0.0
  TOUT= 0.0
  IFLAG = 1
  RELERR = 0.1E-2
  ABSERR = 0.0
  Q_U_DYNAMIC = 0.0
  Y(1) = E1/R1
  Y(2) = 0.0
  Y(3) = -(E2)
  do i=1, 11
    CALL RKF45(diff,NEQN,Y,T,TOUT,RELERR, ABSERR, IFLAG, R_WORK, IWORK)
    write(*,*) "T = ", T, "IFLAG = ", IFLAG
    write(*,*) "Y = ", Y
    write(*,*) "----------------------------------"
    TOUT = T + 0.1E-3
    fm1 = fm1 + (y(3) - U_STATIC(i))**2
  end do
  write(*,*) "fm = ", sqrt(fm1)/11
end program kurs

real function for_l(x)
real x
for_l = cos(x) / x
end function for_l

real function for_l1(x)
real x
for_l1 = x - (0.6**x)
end function for_l1

subroutine diff(T1, Y1, YP)
use params
real T1, Y1(3), YP(3)
i1 = Y1(1)
i3 = Y1(2)
U  = Y1(3)
YP(1) = (1 / L1) * (E1 - E2 - U + i3 * R2 - i1 * (R1 + R2))
YP(2) = (1 / L3) * (E2 + U + i1 * R2 - i3 * (R2 + R3))
YP(3) = (1 / C)  * (i1 - i3)
end subroutine diff

real function sub_fmin(x)
use params

real:: x, fm
 interface
    subroutine diff(T1, Y1, YP)
    real T1, Y1(3), YP(3)
  end subroutine diff
  end interface

  T = 0.0
  TOUT= 0.0
  IFLAG = 1
  RELERR = 0.1E-3
  ABSERR = 0.0
  Q_U_DYNAMIC = 0.0
  Y(1) = E1/R1
  Y(2) = 0.0
  Y(3) = -E2
  C = x
  !write(*,*) "FMIN >> C = ", C
  fm = 0.0
  do i=1, 11
    CALL RKF45(diff,NEQN,Y,T,TOUT,RELERR, ABSERR, IFLAG, R_WORK, IWORK)
    TOUT = T + 0.1E-3
    fm = fm + (y(3) - U_STATIC(i))**2
  end do
  sub_fmin = fm
end function sub_fmin