VMAT_kurs_15

module params
 real:: L, K , M, g
 real :: a, b, relerr, abserr, result, errest, flag !Starting values for QUANC8
 real :: T, TOUT, Y(4), R_WORK(28) !For RKF45
 integer :: IFLAG, IWORK(5), NEQN !For RKF45
 real :: U_STATIC(7), U_DYNAMIC(7), fm !Vectors for U and their quadratic sum
 DATA U_STATIC/0.0, 0.303, -0.465, 0.592, -0.409, 0.164, 0.180/ !Input data for U
 real :: fm_a , fm_b , tol !For fmin

end module params

program kurs15
use params
implicit none
integer i, nofun;
real:: fmin

 interface

  real function for_l(x)
  real x
  end function for_l

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

  real function sub_fmin(x)
  real x
  end function sub_fmin

end interface

 M = 1.0
 g = 9.81
 fm_a = 36.0
 fm_b = 46.0
 tol = 1.0E-04
 a = 0.0
 b = 1.0
 relerr = 1.E-05
 abserr = 1.E-05
 NEQN = 4

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

 K = FMIN(fm_a, fm_b, sub_fmin, tol)
 write(*,*) "FMIN = ", K

  T = 0.0
  TOUT= 0.0
  IFLAG = 1

  Y(1) = 0.0
  Y(2) = 0.0
  Y(3) = 0.0
  Y(4) = 4.0
  fm = 0.0
  do i=1, 7
    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.4
    fm = fm + (Y(1) - U_STATIC(i))**2
  end do
 write(*,*) "fm = ", sqrt(fm)/7

 !-------------------------------------------------------------------

end program kurs15

real function for_l(x)
implicit none
real x
for_l = exp(x)
end function for_l

subroutine diff(T1, Y1, YP)

use params
implicit none
real T1, Y1(4), YP(4)
real xd, yd, od, zd
xd = Y1(1)
od = Y1(2)
yd = Y1(3)
zd = Y1(4)
YP(1) = yd
YP(2) = zd
YP(3) = -K/M*xd - g*(1 - cos(od)) + (L+xd)*zd**2
YP(4) = -g/(L+xd)*sin(od) - 2/(L+xd)*yd*zd
end subroutine diff

real function sub_fmin(x)
use params

real x
 interface
    subroutine diff(T1, Y1, YP)
    real T1, Y1(3), YP(3)
  end subroutine diff
  end interface
  K = x
  fm = 0.0;
  T = 0.0
  TOUT= 0.0
  IFLAG = 1

  Y(1) = 0.0
  Y(2) = 0.0
  Y(3) = 0.0
  Y(4) = 4.0
  do i=1, 7
    CALL RKF45(diff,NEQN,Y,T,TOUT,RELERR, ABSERR, IFLAG, R_WORK, IWORK)
    TOUT = T + 0.4
    fm = fm + (Y(1) - U_STATIC(i))**2
  end do
 sub_fmin = fm
end function sub_fmin