! sfit_imsl, sfit using imsl routines
! Copyright (C) 1999,2010 Gauthier Delerce
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 2, or (at your option)
! any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software Foundation,
! Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
!
! Gauthier Delerce <gauthier@delerce.fr>
!* ------------------------------------------------------------------------------------
!* SFIT1 - finds the maximum from a field of grid data using a spline
!* XMAX=COLUMN, YMAX=ROW (FIXED FROM ORIGINAL CIM)
!* ------------------------------------------------------------------------------------
subroutine sfit1_imsl(ncorx,ncory,data,xmax,ymax,max,hor,ver,roi,small,instable)
use peakopt
implicit none
integer,parameter::icn=5,icn2=icn/2 + 1,m2=icn*icn,nit=8
logical::small,instable
integer::ncorx,ncory
real data(ncorx,ncory)
real fdata(icn,icn)
real xd,yd,z,hor,ver,out
real max,xmax,ymax,dxy,max2,fXmax,fYmax
integer::res
real::roi,dd
integer::nx,ny,idex,jdex,igab,jgabi,ii,jj
integer::i,j,loc,jgab,iflag,jflag,it
real::fi,fj,tttt
!* ------------ arrays for SP1,2 interpolation ------------------------------
real(8)::ro
integer loc1(3),loc2(3)
!* ---------------- internal arrays for SMETA -------------------------------
real(8) r(m2*(m2-1)/2),b(m2),aa((m2-3)*(m2-2)/2)
real(8) c1(m2),c2(m2),c3(m2)
!* --------------------------------------------------------------------------
data loc1,loc2/1,2,3,1,icn,3/
data nx,ny/icn,icn/
if(allocated(xdata)) then
deallocate(xdata)
! print*,"xdata was allocated"
endif
if(allocated(ydata)) then
deallocate(ydata)
! print*,"ydata was allocated"
endif
if(allocated(s)) then
deallocate(s)
! print*,"s was allocated"
endif
if(allocated(t)) then
deallocate(t)
! print*,"t was allocated"
endif
allocate(xdata(m2),ydata(m2),s(m2),t(3))
tttt=0.
!* ------------------------------------------------------------------------------------
!* find largest peak in field idex=column (xmax), jdex=row (ymax)
!* ------------------------------------------------------------------------------------
! print*,10
! Print*,2,roi
! print*,ncorx,ncory,data
small=.false.
max=0.
!: idex=ncorx/2
!: jdex=ncory/2
! print*,1110
ro=dble(roi)
! print*,11100
do j=1,ncory
do i=1,ncorx
! print*,1111
z=data(i,j)
! print*,1112
if(z.ge.max) then
max=z
idex=i
jdex=j
endif
end do
end do
! print*,'idex,jdex',idex,jdex,ncorx,ncory
!* ------------------------------------------------------------------------------------
do j=1,icn
fj=float(j)
jj=(j-1)*icn
do i=1,icn
fi=float(i)
loc=jj+i
xdata(loc)=fi
ydata(loc)=fj
end do
end do
!write(*,*) idex,icn2
! ii=234
!write(*,*) ii,ncory,icn
ii=idex-icn2
jj=jdex-icn2
! iflag=0
! jflag=0
do j=1,icn
jgab=jj+j
jflag=0
if((jgab.gt.ncory).or.(jgab.lt.1))then
do i=1,icn
fdata(i,j)=0.
jflag=1
end do
else
do i=1,icn
igab=ii+i
iflag=0
if((igab.gt.ncorx).or.(igab.lt.1))then
fdata(i,j)=0.
iflag=1
else
fdata(i,j)=data(igab,jgab)
endif
end do
endif
end do
! print*,113
! print*,fdata
if((iflag.eq.1).or.(jflag.eq.1)) small=.true.
!* ------------------------- SMETA interpolation ------------------------------------
nxy=nx*ny
loc2(3)=(icn-1)*icn+1
call exchg(xdata,ydata,fdata,nxy,loc1,loc2,3)
! print*,114
!* ------------------ DEBUG ----------------------------
!c do i=1,25
!c write(*,*)xdata(i),ydata(i)
!c end do
! do j=1,5
! write(*,*)(fdata(i,j),i=1,5)
! end do
! write(*,*)'check sp1'
!c pause
!* -----------------------------------------------------
! print*,111,fdata
! print*,''
! print*,xdata;print*,''
! print*,ydata;print*,''
! print*,nxy;print*,''
! print*,ro;print*,''
! print*,s;print*,''
! print*,t;print*,''
! print*,r;print*,''
! print*,b;print*,''
! print*,aa;print*,''
! print*,c1;print*,''
! print*,c2;print*,''
! print*,c3;print*,''
! print*,nxy;print*,''
call sp1(fdata,xdata,ydata,nxy,ro,s,t,r,b,aa,c1,c2,c3,nxy)
!* --------------------------------------------------------------------------------
! print*,222
ymax=float(icn2)
! print*,2221
xmax=icn2
res=2
! print*,222222
max2=0.
! print*,115
!* ------------------------------------------------------------------------------------
!* iterate on smaller and smaller regions
!* ------------------------------------------------------------------------------------
! xmax=3.0
! ymax=3.0
! call find_peak_intel(xmax,ymax,out)
! if( (xmax<1.).or.(ymax<1.).or.(xmax>5.).or.(ymax>5.) ) then
! instable=.true.
! print*,xmax,ymax
! else
! instable=.false.
! endif
! out=out*-1
! goto 888
do it=1,nit
res=res*2
dxy=1./real(res)
dd=2./real(res)
yd=ymax-dd
! print*,'it',it,xmax,ymax,dd,res
fYmax=ymax
do while(yd<=fYmax+dd)
xd=xmax-dd
fXmax=xmax
do while(xd<=fXmax+dd)
call sp2(xd,yd,xdata,ydata,nxy,s,t,out)
! write(*,*) xd,yd,out,xmax,dd
!* ------------------------------------------------------------------------------------
if(out.ge.max2) then
max2=out
xmax=xd
ymax=yd
! print*,idex-icn2+xmax,jdex-icn2+ymax
! write(*,*) max2
endif
xd=xd+dxy
enddo
yd=yd+dxy
enddo
! do yd=(ymax-dd),(ymax+dd),dxy
! do xd=(xmax-dd),(xmax+dd),dxy
!
! call sp2(xd,yd,xdata,ydata,nxy,s,t,out)
! write(*,*) xd,yd,out
!* ------------------------------------------------------------------------------------
!
! if(out.ge.max2) then
! max2=out
! xmax=xd
! ymax=yd
! print*,idex-icn2+xmax,jdex-icn2+ymax
! write(*,*) max2
! endif
! end do
! end do
end do
! print*,116
!* the displacements computed from the spline are relative to the integer
!* peak location at {idex,jdex}, minus one half the surrounding box size
! print*,xmax,ymax,idex,jdex,icn2 !1166666
xmax=idex-icn2+xmax
ymax=jdex-icn2+ymax
max=max2
if( (xmax<1.).or.(ymax<1.).or.(xmax>ncorx).or.(ymax>ncory) ) then
instable=.true.
print*,'instable found at',xmax,ymax
else
instable=.false.
endif
! write(*,*) ' final= ', xmax,ymax,max2
return
end