! _________________________________________________________________________________ !
! _________________CODE SEEKS BEST FIT FOR Gamma MAXIMUM LIKELIHOOD _____________ !
! _________________________________________________________________________________ !
! _________________Correct solution k=19.288624285110700, Thetha=4.4650184____________ !
program Best_Fit
implicit none
integer, parameter ::  ikind=selected_real_kind(p=15)
! ________________Number of specimens_______________________________________________!
integer :: nospec=26
real (kind=ikind):: x1, TOL
real (kind=ikind)::Ak, Bk, Ck, Dk, k1, k2
real (kind=ikind)::AThe, BThe, CThe, DThe, Thet1, Thet2
real, PARAMETER :: R=1.61803399
integer :: i, j, jm, iio, n=50 ! n iteration number
real (kind=ikind) :: gres1
real, dimension(26) :: ulist1
real (kind=ikind), dimension(4) :: klist,Thelist
real (kind=8), dimension(4,4) :: LLL
real (kind=ikind) :: pdfx, mpdf
ulist1 = (/45, 58, 59, 60, 65, 67, 68, 71, 72, 82, 85, 87, 88, 90, 91, 92, 94, 97, 102, 103, 104, 104, 104, 113, 120, 126/)
!______Those above are specimens ulist1___________________!
! __________Thet1, Thet2 ___Range for Theta___
Thet1=1.
Thet2=12.
! __________k1, k2 ___range for k___
k1=1.
k2=60.
!__________________________________________
TOL=1.d-6

!____________________________inital searching range _________________!
!At any given time we will keep trace of 4 k points: Ak, Bk, Ck, Dk
  Ak=MIN(k1,k2)
  Dk=MAX(k1,k2)
  Bk=Dk-(Dk-Ak)/R
  Ck=Ak+(Dk-Ak)/R
klist =(/Ak, Bk, Ck, Dk/)
print *,' klist= ',klist
!!!______________________________________________________________-!
!At any given time we will keep trace of 4 The points: AThe, BThe, CThe, DThe
  AThe=MIN(Thet1,Thet2)
  DThe=MAX(Thet1,Thet2)
  BThe=DThe-(DThe-AThe)/R
  CThe=AThe+(DThe-AThe)/R
Thelist=(/AThe, BThe, CThe, DThe/)
print *,' Thelist= ',Thelist
!!!______________________________________________________________-!
mpdf=0._ikind
do jm=1,4,1 !____Thelist
do j=1, 4,1 !____klist
do i=1, nospec,1
x1=ulist1(i)
 call gammapdf(x1, Thelist(jm), klist(j), gres1)
pdfx=gres1
mpdf=mpdf+log10(pdfx)
print *, 'x=', x1, 'Log10(Gamma_pdf)=', mpdf, 'pdf:',pdfx
 end do

LLL(jm,j)=mpdf
print *, 'Log10(Gamma_pdf)=', mpdf, 'LLL(jm,j)=',LLL(jm,j)
mpdf=0._ikind
 end do
end do
!______________________________________________________________
!____________________________iterations start here _________________!

!____________________________iterations start here _________________!
do iio=1,n,1
!______________________________________________________________

!______________________________________________________________
IF(LLL(2,2)<LLL(2,3)) THEN
print *,' LLL= ',LLL
klist(1)=klist(2)
klist(2)=klist(4)-(klist(4)-klist(1))/R
klist(3)=klist(1)+(klist(4)-klist(1))/R
print *,'klist= ',klist
    ELSE
klist(4)=klist(3)
klist(2)=klist(4)-(klist(4)-klist(1))/R
klist(3)=klist(1)+(klist(4)-klist(1))/R
print *,'klist= ',klist
ENDIF
!______________________________________________________________
IF(LLL(2,2)<LLL(3,2)) THEN
print *,' LLL= ',LLL
Thelist(1)=Thelist(2)
Thelist(2)=Thelist(4)-(Thelist(4)-Thelist(1))/R
Thelist(3)=Thelist(1)+(Thelist(4)-Thelist(1))/R
print *,'Thelist= ',Thelist
    ELSE
Thelist(4)=Thelist(3)
Thelist(2)=Thelist(4)-(Thelist(4)-Thelist(1))/R
Thelist(3)=Thelist(1)+(Thelist(4)-Thelist(1))/R
print *,'Thelist= ',Thelist
ENDIF
!______________________________________________________________

mpdf=0._ikind
do jm=1,4,1 !____Thelist
do j=1, 4,1 !____klist
do i=1, nospec,1
x1=ulist1(i)
 call gammapdf(x1, Thelist(jm), klist(j), gres1)
pdfx=gres1
mpdf=mpdf+log10(pdfx)
print *, 'x=', x1, 'Log10(Gamma_pdf)=', mpdf, 'pdf:',pdfx
 end do
LLL(jm,j)=mpdf
print *, 'Log10(Gamma_pdf)=', mpdf, 'LLL(jm,j)=',LLL(jm,j)
mpdf=0._ikind
end do
end do
!_______   print LLL matrix_____________________________________

  do jm=1,4
             print *,(LLL(jm,j),j=1,4)
  end do
!______LLL matrix printed
!_______   iterations n continues here_________________________
write(*,*) "iteration iio=",iio
print *,' klist= ',klist
print *,' Thelist= ',Thelist
enddo
stop
  !_____________________________________________
  !_____________________________________________
contains
subroutine gammapdf (a, b, c, gam1)
integer, parameter ::  ikind=selected_real_kind(p=15)
real (kind=ikind), intent(in) :: a, b, c
real (kind=ikind), intent(out):: gam1
gam1 = 1/(Gamma_mb(c)*b**c)*a**(c-1._ikind)*exp(-1._ikind*(a/b))
return
end subroutine gammapdf


  !_____________________________________________
real(kind=8) Function Gamma_mb(xx)
  real(kind=8):: xx,cof(6),stp,half,one,fpf,x2,tmp,ser
integer::j2
  DATA cof,stp /76.18009173d0,-86.50532033d0,24.01409822d0,-1.231739516d0,0.120858003d-2,-0.536382d-5,2.50662827465d0/
  DATA half,one,fpf /0.5d0,1.0d0,5.5d0/
  x2=xx-one
  tmp=x2+fpf
  tmp=(x2+half)*DLOG(tmp)-tmp
  ser=one
  do j2=1,6
    x2=x2+one
    ser=ser+cof(j2)/x2
  end do
  Gamma_mb = DEXP(tmp+DLOG(stp*ser))
  return
end Function Gamma_mb
 !_____________________________________________
end program Best_Fit