subroutine fiber_info(mean,width,errorx,errorw,filename,check,turn)
!used for processing the data in the monitor fibers

 use bmad
 use nr

 implicit none

 real(rp) mean                    !mean and width in unit meter
 real(rp) width
 real(rp) meangs,widthgs          !mean and width from gaussian fit
 real(rp) sum1,sum2,sum3,sum4          !used to calculate the error
 real(rp) errorx,errorw            !error in position and width
 real(rp) time/0./,spacing/0.013/ !spacing of the fibers
 real(rp) indexpo                 !index position
 integer  index/0/,total,i,turn
 integer  ios,ios2,unit/10/,lun,lun2
 character(26):: filename
 logical first_g/.true./,check        !whether to record m and v below
 integer,dimension(0:6):: n   !array that counts numbers of muons that hit each fiber
 real(rp),dimension(4)::sums    !quantities used for the gaussian fit 
 real(rp),dimension(3,3)::m1,m2,m3 ! 
 real(rp),dimension(3,1)::v1   !matrix and vector used for the guassian fit  

mean  = 0
width = 0
sum1  = 0
sum2  = 0
sum3  = 0
sum4  = 0
sums(:) = 0.
m1(:,:) = 0.
v1(:,:)   = 0.
meangs = 0.
widthgs = 0.

errorx = 0
errorw = 0

lun = lunget()

open(UNIT=lun,FILE=trim(filename),STATUS='old', IOSTAT=ios)

do i = 1,2      !skip the first two lines
read(lun,*,IOSTAT=ios2)
end do

do i = 0,6
n(i) = 0
end do

do 
read(lun,*,IOSTAT=ios2) time,index
n(index-1) = n(index-1)+1

!print *,index

if (ios2<0) exit

end do

n(0:6) = n(0:6)+1 !avoid problem when using guassian fit

total = sum(n)

do i = 0,6
sums(1) = sums(1)+(i+1-4)*spacing
sums(2) = sums(2)+((i+1-4)*spacing)**2
sums(3) = sums(3)+((i+1-4)*spacing)**3
sums(4) = sums(4)+((i+1-4)*spacing)**4
end do

m1(1,1) = sums(4)
m1(1,2) = sums(3)
m1(1,3) = sums(2)

m1(2,1) = sums(3)
m1(2,2) = sums(2)
m1(2,3) = sums(1) 

m1(3,1) = sums(2)
m1(3,2) = sums(1)
m1(3,3) = 7

do i = 0,6
v1(1,1) = v1(1,1) + ((i+1-4)*spacing)**2*log(real(n(i)))
v1(2,1) = v1(2,1) + ((i+1-4)*spacing)*log(real(n(i)))
v1(3,1) = v1(3,1) + log(real(n(i)))
end do

if (first_g) then
lun2= lunget()

open(UNIT=lun2,FILE="fiber_info_gaussian_fit.dat",STATUS='replace', IOSTAT=ios)

!write(lun2,'(3a6,a14)') 'v1','v2','v3','elements of m'

close(lun2)

first_g = .false.

end if

!if(check) then 
!lun2= lunget()

!open(UNIT=lun2,FILE="fiber_info_gaussian_fit.dat",STATUS='old',access = 'append', IOSTAT=ios)

!write(lun2,'(12es12.4)') v1(:,:),m1(:,:)

!close(lun2)
!end if

m2 = m1

call gaussj(m1,v1) ! solve for the parameters in the gaussian fit 

widthgs = (-1/2/v1(1,1))**0.5
meangs  = -v1(2,1)/2/v1(1,1)
 
do i = 0,6
mean = mean+n(i)*(i+1)                   !(i+1-4)*spacing convert the indices to the actual x positions of the fiber
end do

mean = mean/total


do i = 0,6
width = width+n(i)*(i+1-mean)**2
end do

width = Sqrt(width/total)


do i= 0,6                                 !calculate the error
sum1 = sum1+sqrt(real(n(i)))*(i+1)
sum2 = sum2+n(i)*(i+1)
sum3 = sum3+sqrt(real(n(i)))
end do

errorx = sum1/total - sum2*sum3/(total)**2

sum1  = 0
sum2  = 0

do i= 0,6
sum1= sum1+Sqrt(real(n(i)))*(i+1-mean)**2
sum2= sum2+n(i)*(i+1-mean)**2
sum4= sum4+n(i)*(-2)*(i+1-mean)*errorx
end do

errorw = 1/2/width*(sum1/total+sum4/total-sum2*sum3/total**2)
!print *, errorw


mean  = (mean-4)*spacing                            ! convert to actually position
width = width*spacing                               ! convert to meter
errorx = errorx*spacing
errorw = errorw*spacing

!mean  = -v1(2,1)/2./v1(1,1)   !right now use the gaussian fit values
!width = sqrt(-1./2./v1(1,1))    
close(lun)


if(check) then 
lun = lunget()

open(UNIT=lun,FILE="fiber_info_gaussian_fit.dat",STATUS = 'old',ACCESS = 'append', IOSTAT=ios)

!m3 = matmul(m1,m2)

!write(lun,'(12es12.4)') v1(:,:),m1(:,:)
!write(lun,'(a10,9es12.4)') 'product',m3(:,:)
write(lun,'(i5,a8,es12.4,a9,es12.4,a6,es12.4,a)') turn,'N(x) = ',exp(v1(3,1)-v1(2,1)**2/4/v1(1,1)),'*exp(-(x-',mean,')**2/',2*width**2,')' 
!write(lun,'(a5,3es12.4)') 'width',v1(1,1),-1./2./v1(1,1),sqrt(-1./2./v1(1,1))

close(lun)

end if




end subroutine



subroutine fiber_records

!record the fiber information (mean,width)

 use bmad

 implicit none

 real(rp) mean/0./
 real(rp) width/0./
 real(rp) errorx/0./,errorw/0./
 integer  unit1,unit2,unitstart/80/,unitstart2/120/,nturn/128/,uwrite/80/,uwrite2/81/,i,n
 character(26):: filename,filename2
 character(3):: turn
 logical first/.true./

do i=1,nturn    !both 180 and 270 degree x monitors

mean = 0.
width = 0.

write(turn,'(i3)') i

filename =  "x_monitor_180_turn"//trim(adjustl(turn))//".dat"

!print *,filename

call fiber_info(mean,width,errorx,errorw,filename,.true.,i)
write (uwrite,'(es12.4,4es12.4)') i+0.5-1,mean,width,errorx 
write (82,'(es12.4,4es12.4)') i+0.5-1,mean,width,errorx       !separate out the information for x-monitor at 180 degree

print *, "running"     

filename2 ="x_monitor_270_turn"//trim(adjustl(turn))//".dat"

print *, "running2"

call fiber_info(mean,width,errorx,errorw,filename2,.false.,i)

print *, "running3"

write (uwrite,'(es12.4,4es12.4)') i+0.75-1,mean,width,errorx     
write (83,'(es12.4,4es12.4)') i+0.75-1,mean,width,errorx         !separate out the information for x-monitor at 270 degree

filename =  "y_monitor_180_turn"//trim(adjustl(turn))//".dat"  
 
call fiber_info(mean,width,errorx,errorw,filename,.false.,i)               !similarly for y-monitor at 180 degree
write (uwrite2,'(es12.4,4es12.4)') i+0.5-1,mean,width,errorx
write (84,'(es12.4,4es12.4)') i+0.5-1,mean,width,errorx


filename =  "y_monitor_270_turn"//trim(adjustl(turn))//".dat"

call fiber_info(mean,width,errorx,errorw,filename,.false.,i)               !for y monitor at 270 degree
write (uwrite2,'(es12.4,4es12.4)') i+0.75-1,mean,width,errorx
write (79,'(es12.4,4es12.4)') i+0.75-1,mean,width,errorx


end do 

!if (first) then !generate some data for trial

!do n=1,256
!write(92,*)n,(1.59+1.07*cos((1.-0.09375)*real(n-1)*2.*pi+1.4/100)+0.23438*cos((2.-0.1875)*real(n-1)*2.*pi+1.3))/10000
!write(92,*)n,(2.57755+1.19*cos((1.-0.09375)*real(n-1)*2.*pi-1.4)+0.395*cos((2.-0.19531)*real(n-1)*2.*pi-0.57))/10000
!write(92,*)n,(2.67845+14.956/64*cos((1.-0.0859)*real(n-1)*2.*pi-1.1276)+1.2273*cos((1.-0.09375)*real(n-1)*2.*pi-1.175466)+35.6/64*cos((1.-0.10156)*real(n-1)*2.*pi+1.828)+0.39444*cos((2.-0.19531)*real(n-1)*2.*pi-0.891))/10000
!write(94,*)n,(2.83+1.3975*cos((1.-0.09375)*real(n-1)*2.*pi+0.263)+0.348*cos((2.-0.19531)*real(n-1)*2.*pi+2.6223))/10000
!write(94,*)n,(2.7416169+1.2522*cos((1.-0.09375)*real(n-1)*2.*pi-1.1715)+0.42606*cos((2.-0.19531)*real(n-1)*2.*pi-0.62644))/10000

!write(92,*)n,(1.31+0.5714*cos((2.-0.19531)*real(n-1)*2.*pi-1.37))/10000


!end do

!first = .false.

!end if

end subroutine



subroutine fiber_data_fit(const0,a10,a20,ph10,ph20,tune0,filename,type,nturn)
!used for fitting the simulation data in filename,typically with the input from the FFT results as the initial guess


 use bmad

 implicit none

 real(rp) const0,a10,a20,ph10,ph20,tune0
 real(rp) const,a1,a2,ph1,ph2,tune
 real(rp) constf,a1f,a2f,ph1f,ph2f,tunef
 real(rp) turn,var1,var2,var3
 real(rp) diff,diff0/0./
 real(rp) mean                    !mean and width in unit meter
 real(rp) width                   
 real(rp),dimension(nturn)::data       
 logical  first/.true./
 integer  nturn,n,ix
 integer  ios,ios2,lun
 integer  type                                  !file type:1 for the files generated by fiber_info.f90; 2 for the files from compute_moments
 character(20) ::filename
 


!lun = lunget()

!open(UNIT=lun,FILE=filename,STATUS='old', IOSTAT=ios)

!do i = 1,2      !skip the first two lines
!read(lun,*,IOSTAT=ios2)
!end do

print *, 'fit file=',trim(filename)

lun = lunget()

open(UNIT=lun,FILE=filename,STATUS='old', IOSTAT=ios)

if (type == 1) then
do n = 1,nturn
read(lun,*,IOSTAT=ios2) turn,mean,width
data(n) = width**2
end do

else 

do n = 1,nturn
read(lun,*,IOSTAT=ios2) ix,turn,var1,var2,var3,width
data(n) = width
end do

end if

close(lun)

a1   = a10
a2   = a20
ph1  = ph10
ph2  = ph20
tune = tune0
const = const0


!do tune = tune0-0.1,tune0+0.1,0.0005
tune= tune0-0.1-0.0005
do while(tune <= tune0+0.1)
 tune = tune + 0.0005
diff = 0.

do n = 1,nturn

diff = diff + (data(n)-(const+a1*cos(tune*2*pi*(n-1)+ph1) + a2*cos((2*tune)*2*pi*(n-1)+ph2)))**2
!print *,index

end do


if (first .or. (diff<diff0)) then 
diff0 = diff 
tunef = tune
end if

first = .false.

end do

tune  = tunef
first = .true.

!do const = const0-1.e-4,const0+1.e-4,0.01e-4
const = const0-1.e-4 -0.01e-4
do while(const <= const0+1.e-4)
 const = const+0.01e-4
diff = 0.

do n = 1,nturn

diff = diff + (data(n)-(const+a1*cos(tune*2*pi*(n-1)+ph1) + a2*cos((2*tune)*2*pi*(n-1)+ph2)))**2
!print *,index

end do


if (first .or. (diff<diff0)) then 
diff0 = diff 
constf = const
end if

first = .false.

end do

const  = constf
first = .true.

!do a1 = a10-0.5e-4,a10+0.5e-4,0.005e-4
a1=a10- 0.5e-4 -0.005e-4 
do while (a1 <= a10+0.5e-4)
  a1=a1+0.005e-4
 
diff = 0.

do n = 1,nturn

diff = diff + (data(n)-(const+a1*cos(tune*2*pi*(n-1)+ph1) + a2*cos((2*tune)*2*pi*(n-1)+ph2)))**2
!print *,index

end do


if (first .or. (diff<diff0)) then 
diff0 = diff 
a1f = a1
end if

first = .false.

end do

a1    = a1f
first = .true.

!do ph1 = ph10-pi,ph10+pi,pi/180
ph1=ph10-pi-pi/180
do while(ph1 <= ph10+pi)
 ph1=ph1+pi/180

diff = 0.

do n = 1,nturn

diff = diff + (data(n)-(const+a1*cos(tune*2*pi*(n-1)+ph1) + a2*cos((2*tune)*2*pi*(n-1)+ph2)))**2
!print *,index

end do


if (first .or. (diff<diff0)) then 
diff0 = diff 
ph1f = ph1
end if

first = .false.

end do

ph1    = ph1f
first = .true.

!do a2 = a20-0.5e-4,a20+0.5e-4,0.005e-4
a2=a20-0.5e-4-0.005e-4
do while(a2 <= a20+0.5e-4)
 a2=a2+0.005e-4
diff = 0.

do n = 1,nturn

diff = diff + (data(n)-(const+a1*cos(tune*2*pi*(n-1)+ph1) + a2*cos((2*tune)*2*pi*(n-1)+ph2)))**2
!print *,index

end do


if (first .or. (diff<diff0)) then 
diff0 = diff 
a2f = a2
end if

first = .false.

end do

a2    = a2f
first = .true.

!do ph2 = ph20-pi,ph20+pi,pi/180
ph2=ph20-pi-pi/180
do while(ph2 <=ph20+pi)
  ph2 = ph2+pi/180
 
diff = 0.

do n = 1,nturn

diff = diff + (data(n)-(const+a1*cos(tune*2*pi*(n-1)+ph1) + a2*cos((2*tune)*2*pi*(n-1)+ph2)))**2
!print *,index

end do


if (first .or. (diff<diff0)) then 
diff0 = diff 
ph2f = ph2
end if

first = .false.

end do

ph2    = ph2f
first = .true.

a10   = a1f
a20   = a2f
ph10  = ph1f
ph20  = ph2f
tune0 = tunef
const0 = constf

end subroutine


subroutine fft_analysis(filename,nturn,unitw,unitw2,type,fit_index)
!do a Fast Fourier Transform to the data in filename and extract the tunes and amplitudes,assuming the data are real.

 
use bmad
use nr

implicit none
integer nturn   !total number of turns(number of data points) involved,it should be a power of 2
real(rp) mean,width,trialnum
integer lun,n,ios,ios2,ntrial
integer unitw,unitw2,unitinv,unittrial/91/    !units to write the result
integer type                                  !file type:1 for the files generated by fiber_info.f90; 2 for the files from compute_moments;3 and 4 for the files from compute_beam_params
integer,optional:: fit_index                  !the index for the files that contain the fitting results
integer turn2
integer ix,istore/1/,i,j
character(20) fitfile1,fitfile2
character(25) filename,writename,writename2
character(8) filename2
character(3) fortnum,indexnum,writenum,writenum2
real(rp) turn,var1,var2,var3
real(rp) const,a0,ph0,tune0
real(rp) tunefit
real(rp) tuneint,phint,aint                  !used for interchanging variables
real(rp),dimension(2):: a1,ph1,tune1         !used to store the fit parameters
real(rp),dimension(nturn)::data,datainv
real(rp),dimension(nturn)::datatrial
real(rp),dimension(nturn/2)::tune
real(rp) emitx,emity,etax,betax
complex(dp),dimension(nturn/2)::amp 
complex(dp),dimension(nturn/2)::amptrial 
logical first/.true./,fit_param/.true./

!write(99,*) present(fit_index)
!print *,'FFT first=',filename

lun = lunget()

open(UNIT=lun,FILE=trim(filename),STATUS='old', IOSTAT=ios)

if (type == 1) then

do n = 1,nturn
read(lun,*,IOSTAT=ios2) turn,mean,width
data(n) = width**2
end do

else if (type == 2) then

do n = 1,nturn
read(lun,*,IOSTAT=ios2) ix,turn,var1,var2,var3,width
data(n) = width

!write(100,'(2es12.4)') turn,width ! as a check
end do

else

read(lun,*,IOSTAT=ios2) !skip the first line

do n = 1,nturn
read(lun,*,IOSTAT=ios2) turn2,emitx,emity,etax,betax

if(type == 3) data(n) = etax !for the eta function curve
if(type == 4) data(n) = betax   !for the beta function curve
if(type == 4) write(170,*) data(n)
end do

end if !end file type selection

close(lun)

lun = lunget()
open(UNIT=lun,FILE="fort.95",STATUS='old', IOSTAT=ios)
read (lun,*) !
read (lun,*) !skip two lines
do n = 1,nturn
read(lun,*,IOSTAT=ios2) ix,turn,var1,var2,var3,trialnum
!datatrial(n)=(1.59+1.07*cos((1-0.10156)*(n-1)*2*pi)+0.23438*cos((2-0.19531)*(n-1)*2*pi))/10000    !for trial
datatrial(n)=trialnum

end do

close(lun)


call realft_dp(data,1,amp)    !do the FFT using the subroutine in Numerical Recipes
call realft_dp(datatrial,1,amptrial) 

!write(unitw,'(3es12.4)') 0.,sum(data),0.   !for zero frequency
lun = lunget()
write(writenum2,'(i3)') unitw2
writename2 =  "FFT_summary"//trim(adjustl(writenum2))//".dat"
open(unit=lun,file= writename2, status= 'replace')
write(lun,*) "FFT results"
!write(unitw2,*) "constant term=",sum(data)/nturn,"   in meter square" 
write(lun,*) "constant term=",abs(amp(1))/nturn,"   in meter square if fitting for beam width squared" 
if (fit_param) const = abs(amp(1))/nturn

write(lun,'(a10,a12,a9)') "tune","  amplitude","   phase"
close(lun)

lun = lunget()
write(writenum,'(i3)') unitw
writename =  "FFT_results"//trim(adjustl(writenum))//".dat"
open(unit=lun,file= writename, status= 'replace')
close(lun)

do n= 1,nturn/2
tune(n)= (n-1)/1./nturn         

lun = lunget()
open(unit=lun,file= writename, access = 'append')
write(lun,*) tune(n),abs(amp(n)),atan2(aimag(amp(n)),real(amp(n)))
close(lun)

if (n>4 .and. abs(amp(n))>abs(amp(n-1)) .and. abs(amp(n))>abs(amp(n+1))&
        .and. abs(amp(n))>abs(amp(n-2)) .and. abs(amp(n))>abs(amp(n+2))&
        .and. abs(amp(n))>abs(amp(n-3)) .and. abs(amp(n))>abs(amp(n+3))) then    !find the peaks (besides the zero frequency)

!write(unitw2,"(a5,es12.4,a12,es12.4,a9,es12.4)") "tune=",tune(n),"  amplitude=",2*abs(amp(n))/nturn,"   phase=",atan2(aimag(amp(n)),real(amp(n))) 
lun = lunget()
open(unit=lun,file = writename2, access= 'append')
write(lun,"(3es12.4)") tune(n),2*abs(amp(n))/nturn,atan2(aimag(amp(n)),real(amp(n))) 
close(lun)

end if
end do

if (present(fit_index)) then   !obtain the fitting curve 

!extract the initial guess of the parameters from ther FFT results
a1(:) = 0              
tune1(:)=0
ph1(:) = 0        


write(fortnum,'(i3)') unitw2
write(indexnum,'(i3)') fit_index


filename2 =  "fort."//trim(adjustl(fortnum))
fitfile1  =  "fit_results"//trim(adjustl(indexnum))//".dat"
fitfile2  =  "fit_curve"//trim(adjustl(indexnum))//".dat"

lun = lunget()
open(UNIT=lun,FILE = writename2,STATUS='old', IOSTAT=ios)

do i = 1,3     !skip the first 3 lines
read(lun,*,IOSTAT=ios2)
end do

do

read(lun,*,IOSTAT=ios2)tune0,a0,ph0

if (a0 > a1(1)) then
tune1(1) = tune0
a1(1)   = a0
ph1(1)  = ph0
end if

if (ios2 < 0) exit
end do

close(lun)

lun = lunget()
open(UNIT=lun,FILE = writename2,STATUS='old', IOSTAT=ios)

do i = 1,3     !skip the first 3 lines
read(lun,*,IOSTAT=ios2)
end do

do 
read(lun,*,IOSTAT=ios2)tune0,a0,ph0

if (a0 > a1(2) .and. a0 < a1(1)) then
tune1(2) = tune0
a1(2)   = a0
ph1(2)  = ph0
end if

if (ios2 < 0) exit
end do

close(lun)

if(tune1(2) < tune1(1)) then !interchange the variables
tuneint  = tune1(2)
phint    = ph1(2)
aint     = a1(2)
tune1(2) = tune1(1)
a1(2)    = a1(1)
ph1(2)   = ph1(1)
tune1(1) = tuneint
a1(1)    = aint
ph1(1)   = phint
end if 

if (type == 3 ) then !for eta fitting, let the dominant amplitude be the first component
if (a1(1)<a1(2)) then
tuneint  = tune1(2)
phint    = ph1(2)
aint     = a1(2)
tune1(2) = tune1(1)
a1(2)    = a1(1)
ph1(2)   = ph1(1)
tune1(1) = tuneint
a1(1)    = aint
ph1(1)   = phint
end if
end if

if (type == 4 ) then !for beta fitting, let the dominant amplitude be the second component
if (a1(1) > a1(2)) then
tuneint  = tune1(2)
phint    = ph1(2)
aint     = a1(2)
tune1(2) = tune1(1)
a1(2)    = a1(1)
ph1(2)   = ph1(1)
tune1(1) = tuneint
a1(1)    = aint
ph1(1)   = phint
end if
end if

!write (94,'(a15,es12.4,a,es12.4,a,es12.4,a,es12.4,a,es12.4,a,es12.4,a,es12.4,a)') "fitting curve: ",const,'+',a1(1),'*cos(',1.-tune1(1),'*2*pi*(n-1)+',ph1(1),') + ',a1(2),'*cos(',2.-tune1(2),&
!      '*2*pi*(n-1)+',ph1(2),')'

!do n = 1,2*nturn
!write (94,*) n, const+a1(1)*cos((1.-tune1(1))*2*pi*n+ph1(1)) + a1(2)*cos((2.-tune1(2))*2*pi*n+ph1(2))
!end do

!do n = 1,2*nturn
!write (98,*) n, const+a1(1)*cos((1.-tune1(1))*2*pi*n+ph1(1)) 
!end do

!improve the fit starting from the FFT results


tunefit = 1.-tune1(1)
if(type == 4) tunefit = 1.- tune1(2)/2
if(type == 3) tunefit = 1.-tune1(1)

lun = lunget()
open(unit=lun,file=fitfile1,status='replace')

do j=1,10 
!print *, 'filename_FFT=',filename
call fiber_data_fit_loop_all(const,a1(1),a1(2),ph1(1),ph1(2),tunefit,filename,type,nturn)
write(lun,'(6es12.4)') const,a1(1),a1(2),ph1(1),ph1(2),tunefit

end do

call sleep(10)

!call fiber_data_fit_dfpmin(const,a1(1),a1(2),ph1(1),ph1(2),tunefit,filename,type,nturn)
write(lun,'(6es12.4)') const,a1(1),a1(2),ph1(1),ph1(2),tunefit



close(lun)

lun = lunget()
open(unit=lun,file=fitfile2,status='replace')

write (lun,'(a15,es12.4,a,es12.4,a,es12.4,a,es12.4,a,es12.4,a,es12.4,a,es12.4,a)') "fitting curve: ",const,'+',a1(1),'*cos(',tunefit,'*2*pi*(n-1)+',ph1(1),') + ',a1(2),'*cos(',2.*tunefit,&
       '*2*pi*(n-1)+',ph1(2),')'

do n = 1,2*nturn
write (lun,*) n, const+a1(1)*cos(tunefit*2*pi*(n-1)+ph1(1)) + a1(2)*cos((2*tunefit)*2*pi*(n-1)+ph1(2))
end do

close(lun)

end if
 

if (first) then    !trial data

!write(unittrial,'(3es12.4)') 0.,sum(datatrial),0.   !for zero frequency

do n= 1,nturn/2        
write(unittrial,*) (n-1)/1./(nturn),abs(amptrial(n)),atan2(aimag(amptrial(n)),real(amptrial(n)))
end do

first=.false.
end if

!call realft_dp(datainv,-1,amptrial) !do an inverse FFT to check the results

!do n = 0, nturn-1
!write(unitinv,*) n,datainv(n+1)*2/nturn   !+sum(datatrial)/nturn
!end do

end subroutine


subroutine fiber_data_fit_dfpmin(const0,a10,a20,ph10,ph20,tune0,filename,type,nturn)
!used for fitting the simulation data in filename,typically with the input from the FFT results as the initial guess
!the Fletcher-Reeves-Polak-Ribiere minimization is used


 use bmad
 !use nr,only : frprmn
 use nr
 !use nrtype

 implicit none

 real(rp) const0,a10,a20,ph10,ph20,tune0
 real(rp) const,a1,a2,ph1,ph2,tune
 real(rp) constf,a1f,a2f,ph1f,ph2f,tunef
 real(rp) turn,var1,var2,var3
 real(rp) var
 real(rp) diff,diff0/0./,diff0f
 real(rp) mean                    !mean and width in unit meter
 real(rp) width    
 real(rp) ftol,diffmin,gtol/0.0000001/
 real(rp),dimension(6)    ::fit0,fit1,fitf,fitf_opt,var_range,var_steps               
 real(rp),dimension(nturn)::data       
 integer,dimension(6)     ::orders
 logical  first/.true./,firstf
 integer  nturn,n,ix,iter
 integer  i,i1,i2,i3,i4,i5,i6
 integer  ios,ios2,lun
 integer  type                                  !file type:1 for the files generated by fiber_info.f90; 2 for the files from compute_moments
 character*20 filename
 

lun = lunget()

open(UNIT=lun,FILE=filename,STATUS='old', IOSTAT=ios)

if (type == 1) then
do n = 1,nturn
read(lun,*,IOSTAT=ios2) turn,mean,width
data(n) = width**2
end do

else 

do n = 1,nturn
read(lun,*,IOSTAT=ios2) ix,turn,var1,var2,var3,width
data(n) = width
end do

end if

close(lun)

fit0(1) = const0
fit0(2)   = a10
fit0(3)   = a20
fit0(4)  = ph10
fit0(5) = ph20
fit0(6) = tune0

!call frprmn(fit0,ftol,iter,diffmin)
call dfpmin(fit0,gtol,iter,diffmin,func,dfunc)

print *,"diffmin=", diffmin
a10   = fit0(2)
a20   = fit0(3)
ph10  = fit0(4)
ph20  = fit0(5)
tune0 = fit0(6)
const0 = fit0(1)

contains

function func(p)  
!use nrtype
implicit none
real(rp),dimension(:),intent(in)::p  !dlr 6-> :
real(rp) func
integer n

func = 0.

do n = 1,nturn
func = func + (data(n)-(p(1)+p(2)*cos(p(6)*2*pi*(n-1)+p(4)) + p(3)*cos((2*p(6))*2*pi*(n-1)+p(5))))**2
end do
 
end function func

function dfunc(p)
!use nrtype
implicit none
real(rp),dimension(:),intent(in)::p  !dlr 6 -> :
real(rp),dimension(size(p))::dfunc  !dlr 6 -> :
integer n

dfunc(:) = 0.

do n = 1,nturn

dfunc(1) = dfunc(1)+ 2*(data(n)-(p(1)+p(2)*cos(p(6)*2*pi*(n-1)+p(4)) + p(3)*cos((2*p(6))*2*pi*(n-1)+p(5))))*(-1)
dfunc(2) = dfunc(2)+ 2*(data(n)-(p(1)+p(2)*cos(p(6)*2*pi*(n-1)+p(4)) + p(3)*cos((2*p(6))*2*pi*(n-1)+p(5))))*(-1)*cos(p(6)*2*pi*(n-1)+p(4))
dfunc(3) = dfunc(3)+ 2*(data(n)-(p(1)+p(2)*cos(p(6)*2*pi*(n-1)+p(4)) + p(3)*cos((2*p(6))*2*pi*(n-1)+p(5))))*(-1)*cos((2*p(6))*2*pi*(n-1)+p(5))
dfunc(4) = dfunc(4)+ 2*(data(n)-(p(1)+p(2)*cos(p(6)*2*pi*(n-1)+p(4)) + p(3)*cos((2*p(6))*2*pi*(n-1)+p(5))))*p(2)*sin(p(6)*2*pi*(n-1)+p(4))
dfunc(5) = dfunc(5)+ 2*(data(n)-(p(1)+p(2)*cos(p(6)*2*pi*(n-1)+p(4)) + p(3)*cos((2*p(6))*2*pi*(n-1)+p(5))))*p(3)*sin((2*p(6))*2*pi*(n-1)+p(5))
dfunc(6) = dfunc(6)+ 2*(data(n)-(p(1)+p(2)*cos(p(6)*2*pi*(n-1)+p(4)) + p(3)*cos((2*p(6))*2*pi*(n-1)+p(5))))*(p(2)*sin(p(6)*2*pi*(n-1)+p(4))*2*pi*(n-1) + p(3)*sin((2*p(6))*2*pi*(n-1)+p(5)) &
           *4*pi*(n-1))

end do

end function dfunc 

end subroutine  fiber_data_fit_dfpmin





subroutine fiber_data_fit_loop_all(const0,a10,a20,ph10,ph20,tune0,filename,type,nturn)
!used for fitting the simulation data in filename,typically with the input from the FFT results as the initial guess
!in the least-square fitting,loop over all the parameters in different orders and find the best match


 use bmad

 implicit none

 real(rp) const0,a10,a20,ph10,ph20,tune0
 real(rp) const,a1,a2,ph1,ph2,tune
 real(rp) constf,a1f,a2f,ph1f,ph2f,tunef
 real(rp) turn,var1,var2,var3
 real(rp) var
 real(rp) diff,diff0/0./,diff0f
 real(rp) mean                    !mean and width in unit meter
 real(rp) width
 real(rp) emitx,emity,etax,betax    
 real(rp),dimension(6)    ::fit0,fit1,fitf,fitf_opt,var_range,var_steps               
 real(rp),dimension(nturn)::data       
 integer,dimension(6)     ::orders
 logical  first/.true./,firstf
 integer  nturn,n,ix
 integer  turn2
 integer  i,i1,i2,i3,i4,i5,i6
 integer  ios,ios2,lun
 integer  type                                  !file type:1 for the files generated by fiber_info.f90; 2 for the files from compute_moments;3 and 4 for the files from compute_beam_params
 character*25 filename
 
firstf = .true.

lun = lunget()

open(UNIT=lun,FILE=filename,STATUS='old', IOSTAT=ios)

if (type == 1) then
do n = 1,nturn
read(lun,*,IOSTAT=ios2) turn,mean,width
data(n) = width**2
end do

else if (type == 2) then

do n = 1,nturn
read(lun,*,IOSTAT=ios2) ix,turn,var1,var2,var3,width
data(n) = width
end do

else 

read(lun,*,IOSTAT=ios2) !skip the first line

do n = 1,nturn
read(lun,*,IOSTAT=ios2) turn2,emitx,emity,etax,betax

if(type == 3) data(n) = etax !for the eta function curve
if(type == 4) data(n) = betax   !for the beta function curve
end do

end if

close(lun)

fit0(1) = const0
fit0(2)   = a10
fit0(3)   = a20
fit0(4)  = ph10
fit0(5) = ph20
fit0(6) = tune0

if (type == 1 .or. type == 2) then
var_range(1) = 1.e-4
var_range(2) = 1e-4
var_range(3) = 0.5e-4
var_range(4) = pi
var_range(5) = pi
var_range(6) = 0.01
end if

if (type ==3 .or. type == 4) then
var_range(1) = 5
var_range(2) = 5
var_range(3) = 5
var_range(4) = pi
var_range(5) = pi
var_range(6) = 0.01
end if

var_steps(:) = 1000

do i1 = 1,6
do i2 = 1,6
do i3 = 1,6
do i4 = 1,6
do i5 = 1,6
do i6 = 1,6

if (i1==i2 .or. i1==i3 .or. i1==i4 .or. i1==i5 .or. i1==i6 .or. i2==i3 .or. i2==i4 .or. i2==i5 .or. i2==i6 .or. i3==i4 .or. i3==i5 .or. i3==i6 .or. i4==i5 .or. i4==i6 .or. i5==i6) cycle

orders(1) = i1
orders(2) = i2
orders(3) = i3
orders(4) = i4
orders(5) = i5
orders(6) = i6



fit1(1) = const0
fit1(2)   = a10
fit1(3)   = a20
fit1(4)  = ph10
fit1(5) = ph20
fit1(6) = tune0

do i = 1,6

!do var = fit0(orders(i))-var_range(orders(i)),fit0(orders(i))+var_range(orders(i)),var_range(orders(i))/var_steps(orders(i))
var = fit0(orders(i))-var_range(orders(i)) - var_range(orders(i))/var_steps(orders(i))
do while(var < fit0(orders(i))+var_range(orders(i)) )
 var = var + var_range(orders(i))/var_steps(orders(i))

fit1(orders(i)) = var
 
diff = 0.

do n = 1,nturn

diff = diff + (data(n)-(fit1(1)+fit1(2)*cos(fit1(6)*2*pi*(n-1)+fit1(4)) + fit1(3)*cos((2*fit1(6))*2*pi*(n-1)+fit1(5))))**2
!print *,index

end do


if (first .or. (diff<diff0)) then 
diff0 = diff 
fitf(orders(i)) = fit1(orders(i))
end if

first = .false.

end do

fit1(orders(i))  = fitf(orders(i))
first = .true.

end do

if (firstf .or. diff0 < diff0f) then
diff0f = diff0
fitf_opt(:) = fitf(:)
write(170+type,'(7es12.4)') diff0f,fitf(:)
end if

firstf = .false.

end do
end do
end do
end do
end do
end do

!print *, "diff0f=",diff0f
a10   = fitf_opt(2)
a20   = fitf_opt(3)
ph10  = fitf_opt(4)
ph20  = fitf_opt(5)
tune0 = fitf_opt(6)
const0 = fitf_opt(1)

end subroutine

subroutine svd_analysis(nturn)
!do a singluar value decomposition analysis to the data from 180 and 270 fiber monitors

use bmad
use nr

integer   nturn,i,j,k
integer   ios,lun
real(rp)  turn,mean,width
real(rp),dimension(nturn,2):: a,a2,a4       !matrix that contains the fiber monitor data and those used for checking
real(rp),dimension(2)  :: w                 !diagonal matrix in the svd decomposition
real(rp),dimension(nturn)  :: a3            !used for checking
real(rp),dimension(2,2):: v                 !matrix in the decomposition
real(rp),dimension(2,2):: w2,w3             



a(:,:)  = 0.
a3(:)   = 0.
a4(:,:) = 0.
w(:)    = 0.
v(:,:)  = 0.
w2(:,:) = 0.
w3(:,:) = 0.
a2(:,:) = 0.

lun = lunget()
open(UNIT=lun,FILE="fort.82",STATUS='old', IOSTAT=ios)

do i = 1,nturn
read(lun,*,IOSTAT=ios) turn,mean,width 
a(i,1) = width**2
a3(i)  = width**2
end do

close(lun)

lun = lunget()
open(UNIT=lun,FILE="fort.83",STATUS='old', IOSTAT=ios)

do i = 1,nturn
read(lun,*,IOSTAT=ios) turn,mean,width 
a(i,2) = width**2
end do

close(lun)

call svdcmp_dp(a,w,v)
w3(1,1) = w(1)
w3(2,2) = w(2)

a4 = matmul(a,matmul(w3,transpose(v)))

lun = lunget()
open(unit=lun,file="svd_results.dat",status='replace')

!write(lun,'(6es12.4)') w(1),w(2),v(1,1),v(1,2),v(2,1),v(2,2)

if (w(1) > w(2)) then
w2(1,1) = w(1)
a2 = matmul(a,matmul(w2,transpose(v)))

do i = 1,nturn
turn = i*1.
write(lun,'(i10,5es12.4)') i,turn,a3(i),a4(i,1),a4(i,2),a2(i,1) !the last a2(i,1) is for the file format in order to use of FFT_analysis
end do

else
w2(2,2) = w(2)
a2 = matmul(a,matmul(w2,transpose(v)))

do i = 1,nturn
turn = i*1.
write(lun,'(i10,5es12.4)') i,turn,a3(i),a4(i,1),a4(i,2),a2(i,1)
end do

end if

close(lun)

end subroutine