#
# Procedures for G2S
#
PROCEDURE g2s
PARAMETERS
wg=workspace
RESULT ws
LOCAL g2s res
ws=wg
g2s=fields() ; g2s.lptin=wg.ntc
g2s.Xin=wg.x ; g2s.Yin=wg.y ; g2s.Ein=wg.e
g2s.qmax=inquire("g2s> qmax ")
g2s.npt =inquire("g2s> number of points ")
g2s.ic =inquire("g2s> window function code ")
g2s.rho =inquire("g2s> number density ")
module/load "g2s.so" ;
res=module:execute("g2s", g2s)
ws.ntc=res.lptout
ws.x=res.Xout ; ws.y=res.Yout ; ws.e=res.Eout
ws.xlabel=res.xcaptout ; ws.ylabel=res.ycaptout
ENDPROCEDURE
Now the FORTRAN program that performs the actual calculation:
SUBROUTINE G2S(g2s_get, g2s_put)
EXTERNAL g2s_get, g2s_put
INTEGER mn, lptin, lptout
PARAMETER (mn=33000)
REAL*4 Xin(mn),Yin(mn),Ein(mn),Xout(mn),Yout(mn),Eout(mn)
REAL*4 xnew(mn),yw(mn)
CHARACTER*40 xcaptout, ycaptout
CHARACTER*10 char
LOGICAL LMOD
DANGLE=38.1*1.112/150.0/150.0
PI=4.0*ATAN(1.0D0)
call module_get_int(g2s_get, 'lptin', lptin)
call module_get_real_array(g2s_get, 'Xin', Xin, lptin)
call module_get_real_array(g2s_get, 'Yin', Yin, lptin)
call module_get_real_array(g2s_get, 'Ein', Ein, lptin)
call module_get_real(g2s_get, 'qmax', QMAX)
call module_get_int(g2s_get, 'npt', npt)
call module_get_int(g2s_get, 'ic', ichar)
call module_get_real(g2s_get, 'rho', RHO)
if(lptin.eq.0)then
call module_error(" G2S", 1 "ERROR ** No input data", " ")
endif
if(npt.eq.0)then
call module_error(" G2S>", 1 "ERROR ** Number of output points zero", " ")
RETURN
endif
lptout=npt
if(lptout.GT.mn)then
call module_error(" G2S>", 1 "#points reduced from 33000", " ")
lptout=mn
endif
delq=QMAX/lptout
if(rho.lt.1e-10)then
call module_error(" G2S>", 1 "ERROR ** rho is zero", " ")
RETURN
endif
if(ichar.eq.0)then
LMOD=.false.
else
if(ichar.eq.1)then
LMOD=.true.
else
LMOD=.false.
endif
endif
if(LMOD)then
A=PI/xin(lptin)
do nn=1,lptin
yw(nn)=SIN(xin(nn)*A)/xin(nn)/A
end do
else
do nn=1,lptin
yw(nn)=1.0
end do
endif
C
C FORM VECTOR OF EQUALLY-SPACED R'S AT WHICH THE FOURIER TRANSFORM C IS TO BE COMPUTED (RMAX IN ANGSTROMS), AND THE NUMBER OF R-POINTS.
C
DO NR=1,lptout
xout(NR)=delq*NR
end do
C
C THE NUMBER OF POINTS IN THE RANGE OF DATA TO BE TRANSFORMED.
C
xd=(xin(2)-xin(1))/2. !half x-channel
do n=1,lptin
c xnew(n)=xin(n) +xd !offset - mid channel
xnew(n)=xin(n)
yw(n)=yw(n)*(yin(n)-1.)*xnew(n)
end do
C
C COMPUTE FOURIER TRANSFORM OF THE DATA
C
delr=(xin(lptin)-xin(1))/(lptin-1)
AFACT=delr*2.0/PI
DO 35 NR=1,lptout
FS=0.0
RP=xout(NR)
DO N=2,lptin
SINUS1=SIN(xnew(N-1)*RP)
SINUS=SIN(xnew(N)*RP)
FS=FS+ (SINUS*yw(N)+ 1 SINUS1*yw(N-1))/2.0
end do
yout(NR)=FS*AFACT
35 CONTINUE
C
pir=PI*PI*2.*RHO
do n=1,lptout
yout(n)=yout(n)*pir/xout(n) +1.
end do
xcaptout=' Q (Angstrom-1) ' ycaptout=' Structure Factor S of Q '
call module_put_int(g2s_put, 'lptout', LPTOUT)
call module_put_real_array(g2s_put, 'Xout', Xout, lptout)
call module_put_real_array(g2s_put, 'Yout', Yout, npt)
call module_put_real_array(g2s_put, 'Eout', Eout, npt)
call module_put_string(g2s_put, 'ycaptout', ycaptout)
call module_put_string(g2s_put, 'xcaptout', xcaptout)
call module_information(' G2S> output is in point mode')
RETURN
END