Files
fit/pgm/bose.f
2022-08-19 15:22:33 +02:00

181 lines
3.7 KiB
Fortran

program BOSE
! ------------
!
! Simple user function example (straight line).
!
implicit none
external BOSE_FUN
character str*32
integer i,l
!---
! Welcome message
call fit_vers(str)
call str_trim(str, str, l)
type '(X,2A)','Program FIT(BOSE) Version ',str(1:l)
do i=1,l
str(i:i)='-'
enddo
type '(X,2A/)','-----------------------------',str(1:l)
!---
! Function title and parameter names
!
call fit_confun('lorenz/bose', bose_fun) ! function title, function
call fit_userpar('BG:Bgr(0)')
call fit_userpar('DB:dBgr/dX')
call fit_userpar('EF:EF')
call fit_userpar('T:Temp')
do i=1,9
write(str, '(2(a,i1))') 'P',i,':Pos',i
call fit_userpar(str)
write(str, '(2(a,i1))') 'I',i,':Int',i
call fit_userpar(str)
write(str, '(2(a,i1))') 'W',i,':Wid',i
call fit_userpar(str)
enddo
call fit_main
end
real function bose_fun(x,p,n,mode,cinfo)
! -------------------------------------------
implicit none
real x ! x-value
integer n ! number of parameters
real p(n) ! parameters
integer mode ! mode
integer cinfo ! calculation information (see below)
integer i,i0,j,k
parameter (i0=7)
real x0,w0,y0,db,bg,kf,l0
real voigt, bose_fact
real xnew(9),ynew(9),wnew(9)
if (mode .eq. 0) then
bose_fun=0
do i=i0,n-2,3
if (p(i+2) .ne. 0) then ! ignore delta functions (treated later)
bose_fun=bose_fun+p(i+1)*voigt(x-p(i), 0.0, p(i+2))
if (p(i+2) .lt. 0) then ! make a mirror peak for negative width
x0=x+p(i)
if (i .gt. i0 .and. p(i0+2) .eq. 0.0) x0=x0-2*p(i0) ! shift zero
bose_fun=bose_fun+p(i+1)*voigt(x0, 0.0, -p(i+2))
endif
endif
enddo
bose_fun=bose_fun*bose_fact(x/p(6))+p(3)+p(4)*x
elseif (mode .eq. 1) then
! x-independent part
do i=i0,n-2,3
if (p(i+2) .eq. 0) then ! treat delta functions
call fit_delta(p(i),p(i+1))
endif
enddo
call fit_limit_xrange(-p(5),1e6)
elseif (mode .eq. 2) then ! transform x -> t
x0=x+2*p(5)
if (x0 .ge. 0) then
bose_fun=sqrt(x0)
else
bose_fun=0
endif
elseif (mode .eq. 3) then ! transform t -> x
bose_fun=x*x-2*p(5)
else
if (nint(x) .eq. -1 .and. n .ge. 7) then ! convert from multi-voigt
print *
if (n .eq. 7) then
print *,'Convert from voigt'
else
print *,'Convert from multi-voigt'
endif
db=p(2)
bg=p(1)-p(3)*db ! different bg definition
j=3
x0=p(3)
y0=p(5)
w0=max(abs(p(6)), abs(p(7)))
l0=p(7)
do i=8,n,5
if (abs(p(i)) .lt. abs(x0)) then
j=i
x0=p(i)
y0=p(i+2)
w0=max(abs(p(i+3)), abs(p(i+4)))
endif
enddo
k=0
do i=3,n,5
if (p(i) .gt. 0 .and. i .ne. j) then
k=k+1
xnew(k)=p(i)
ynew(k)=p(i+2)
wnew(k)=-max(abs(p(i+3)), abs(p(i+4)))
endif
enddo
p(1)=w0
p(2)=w0*0.05
p(3)=bg
p(4)=db
kf=1.55
call fit_get_real('KFIX', kf)
p(5)=2.0723*kf*kf
p(6)=10 ! default Temp
call fit_get_real('Temp', p(6))
p(7)=x0
p(8)=y0
p(9)=l0
i=10
do j=1,k
p(i)=xnew(k)
p(i+1)=ynew(k)/bose_fact(p(i)/p(6))
p(i+2)=wnew(k)
i=i+3
enddo
x=i-1
else
print *
print *,'Up to 9 Lorenzians multiplied with bose_factor'
1 ,', folded with gaussian'
endif
print *,'Negative Wid makes a mirror peak at -Pos'
print *
endif
end
real function bose_fact(wot)
! argument: omega[meV] over T[K]
real wot
real K_meV, x
parameter (K_meV=11.6048)
x=wot*K_meV
if (abs(x) .lt. 1e-5) then
bose_fact=1.0
else if (x .lt. -30) then
bose_fact=0
else
bose_fact=x/(1-exp(-x))
endif
end