369 lines
11 KiB
Fortran
369 lines
11 KiB
Fortran
program TriCS_CCL
|
|
c -----------------
|
|
c
|
|
c fits data from ccl-files and write *.col file for upals
|
|
c and rafin.dat file for position refining
|
|
c
|
|
|
|
c replaces "none"
|
|
c 8.5.99/21.6.99/15.7.99 Juerg Schefer
|
|
c 29.11.99: Lorentz-Korrektur SJ33
|
|
c 10.4.2000: flat background SJ33
|
|
c 31.10.2000: transfered to UNIX ZO33
|
|
c 22.11.2000: problem with printing is fixed
|
|
c 02.02.01: h,k,l may be real ZM33
|
|
c 17.06.2010: Lorentz correction for tilt Geometry, SJ33
|
|
c
|
|
c link ccl,fit4_shr/opt ! attention: check for last fit_n
|
|
c
|
|
|
|
implicit none
|
|
|
|
integer np_max
|
|
parameter (np_max=1000)
|
|
real output(20),th2,omega_cal,omega_fit,chi,phi,yint,yy1,
|
|
* dint,delta,yoverall,h,k,l,s,temperature,range
|
|
real omega_err
|
|
real output2(20), xx(np_max),lor_cor,fwhm,radian,degrees
|
|
real bgr, wmult
|
|
logical plot,iexp,flat_backgr
|
|
integer lorentz,comens ! SJ33, June 2010
|
|
integer next_flag
|
|
integer iall,iincom,icom,icenter
|
|
character*1 answer,hard
|
|
character*128 input_files
|
|
character*128 output_file1
|
|
character*128 output_file2
|
|
character*128 output_file3 ! SJ33, May 2012
|
|
character*128 output_file4 ! SJ33, May 2012
|
|
character*16 out
|
|
character*2 off,weak
|
|
character*30 hkl
|
|
character*12 hkl_int
|
|
character*18 hkl_float
|
|
character*24 datime
|
|
integer ldt
|
|
|
|
integer ntyp,i,j,numor,ii,iout,np
|
|
radian=3.14159265/180.00
|
|
degrees=1./radian
|
|
|
|
write(*,1000)
|
|
1000 format (/,' FIT CCL Files: Version 2.08, update Jun. 21, 2010',
|
|
* //,' Data processing of files *.ccl from TriCS single'
|
|
* ,' detector mode')
|
|
|
|
call sys_setenv('DAT_DEFSPEC', 'TRICS')
|
|
call dat_ask_filelist(input_files, ' ')
|
|
output_file1=input_files
|
|
i=index(output_file1, '.')
|
|
if (i .eq. 0) i=index(output_file1, '[')
|
|
if (i .eq. 0) i=index(output_file1, ',')
|
|
if (i .eq. 0) i=index(output_file1, ' ')
|
|
do j=i-1,1,-1
|
|
if (output_file1(j:j) .eq. '/') then
|
|
goto 1002
|
|
endif
|
|
enddo
|
|
j=0
|
|
1002 if (i .gt. j+1) then
|
|
output_file2=output_file1(j+1:i-1)//'_rafin.dat'
|
|
output_file1=output_file1(j+1:i-1)//'.col'
|
|
output_file3=output_file1(j+1:i-1)//'.comm'
|
|
output_file4=output_file1(j+1:i-1)//'.incomm'
|
|
else
|
|
output_file2='rafin.dat'
|
|
output_file1='upals.col'
|
|
endif
|
|
c
|
|
write(*,1004)
|
|
1004 format (/,' function type to be used (0=gauss[def],6=strange)',
|
|
* ' typ = ',$)
|
|
read (5,'(i1)') ntyp
|
|
if (ntyp .eq. 6) then
|
|
write(*,'(/,'' width factor [1.0] = '',$)')
|
|
read(5,'(f20.0)') wmult
|
|
if (wmult .eq. 0.0) wmult=1.0
|
|
else
|
|
wmult=1.0
|
|
ntyp=0
|
|
endif
|
|
|
|
write(*,1013)
|
|
1013 format (/,' Monitor to be standardized , monitor=',$)
|
|
read (5,'(f20.0)') yoverall
|
|
if (yoverall.le.1.1) yoverall = 100000.
|
|
|
|
write(*,1015)
|
|
1015 format (/,
|
|
* ' plots (answer=n [no=def], y [yes,terminal], h [hard] ',$)
|
|
read (5,'(a1)') answer
|
|
plot = (answer.eq.'y' .or. answer .eq.'h')
|
|
if (plot) then
|
|
hard='y'
|
|
if (answer.eq.'y') hard='n'
|
|
endif
|
|
|
|
write(*,1019)
|
|
1019 format (' Flat Background (y/n, default=y) : ',$)
|
|
read (5,'(a1)')answer
|
|
flat_backgr = (answer.ne.'n' .and. answer .ne.'N')
|
|
|
|
if (ntyp .eq. 0) then
|
|
write(*,1017)
|
|
1017 format
|
|
* (//,' output intensities: 0=experimental (def), 1=gaussian ',$)
|
|
read (5,'(i1)') Iout
|
|
iexp= (iout.eq.0)
|
|
else
|
|
iexp= .true.
|
|
endif
|
|
if (iexp) then
|
|
out='experimental '
|
|
else
|
|
out='gaussian '
|
|
end if
|
|
|
|
write(*,1018)
|
|
1018 format (' Lorentz-Correction to be applied (0 = 4-circle-geometry(default), 1 = tilt-geometry, 2 = not applied) ?: ',$)
|
|
read (5,'(I1)') lorentz ! SJ33, June 2010, tilt option added (3 options now, def=0)
|
|
if (lorentz.ne.1 .and. lorentz.ne.2) lorentz=0 ! SJ33, June 2010
|
|
write(*,1011) input_files,ntyp,wmult,
|
|
* output_file1,output_file2,out,yoverall,flat_backgr,lorentz
|
|
1011 format (/,' ',78('-'),/,' input from file(s) :',A30,/,
|
|
* ' function type ',i3,' width factor ',f6.3,//,
|
|
* ' output: for UPALS :',A50,/,
|
|
* ' for RAFIN :',A50,/,
|
|
* 1X,A16,' intensities used'/,
|
|
* ' monitor: ',f12.1,/,
|
|
* ' Flat Background (yes/no): ',L33,/,
|
|
* ' Lorentz-Correction (options: 0 = 4-circle-mode,',
|
|
* ' 1 = tilt-mode, 2 = none)',/,
|
|
* ' applied option No. ',i1,/,
|
|
* ' No. 0 Icor=Iobs*sin(2theta)',/,
|
|
* ' No. 1 Icor=Iobs*sin(Gama)*cos(nu)',/,
|
|
* ' ',78('-'),/,' hit return to start ',$)
|
|
c
|
|
read (5,'(i1)') ii
|
|
C if (ii.eq.1) goto 504
|
|
|
|
c
|
|
if (ntyp.ne.0 .and. ntyp.ne.6) ntyp = 0
|
|
c
|
|
c this ends the input for the program, the rest goes automatic
|
|
c ------------------------------------------------------------
|
|
if (plot) then
|
|
call sys_setenv('CHOOSER_PAN','9') ! 9 plots per page
|
|
end if
|
|
c
|
|
c loop over all data following:
|
|
iall =0
|
|
iincom =0
|
|
icom =0
|
|
icenter =0
|
|
open (unit=1,file=output_file1,status='unknown')
|
|
open (unit=2,file=output_file2,status='unknown')
|
|
open (unit=3,file=output_file3,status='unknown')
|
|
open (unit=4,file=output_file4,status='unknown')
|
|
c write headers for fullprof
|
|
c write (3,1027)9.99,0,1,1,0.00,0.00,0.00
|
|
c write (4,1028)9.99,0,1,1,0.00,0.00,0.00
|
|
c1027 format (' title',/,22H(I6,3I10,2F10.2,4F8.2), /,F7.4,i5,i6,/,i1,3f5.2)
|
|
c1028 format (' title',/,24H(I6,3F10.4,2F10.2,4F8.2),/,F7.4,I5,i6,/,i1,3f5.2)
|
|
next_flag=0 ! start flag
|
|
i=0
|
|
call fit_init
|
|
c
|
|
c ------------------------------------------------------ start loop --
|
|
1500 continue
|
|
|
|
call fit_dat_next(input_files, next_flag)
|
|
if (next_flag .eq.0) goto 2222 ! last scan of last file read
|
|
i=i+1
|
|
|
|
call fit_mon(yoverall)
|
|
|
|
C set limits if peak to narrow:
|
|
call fit_get_array('P', output, 20, ii)
|
|
call fit_get_array('X', xx, np_max, np)
|
|
|
|
|
|
c get information out of title:
|
|
call fit_get_real('two_theta',th2)
|
|
call fit_get_real('omega',omega_cal)
|
|
call fit_get_real('chi',chi)
|
|
call fit_get_real('phi',phi)
|
|
call fit_get_real('h',h)
|
|
call fit_get_real('k',k)
|
|
call fit_get_real('l',l)
|
|
call fit_get_real('temp',temperature) ! SJ33, Aug.30,1999
|
|
call fit_get_str('date', ldt, datime)
|
|
|
|
cjs type *,' ***** temp=',temperature
|
|
numor=i
|
|
if (l.eq.88) then
|
|
write(*,1007)
|
|
1007 format (' data numor not found, exit')
|
|
goto 3000
|
|
endif
|
|
c
|
|
if (ntyp .eq. 6) then
|
|
call fit_fun(6,1,wmult,0.0) ! select strange
|
|
else
|
|
call fit_fun(ntyp,0,0.0,0.0) ! select other (guassian)
|
|
endif
|
|
if (flat_backgr) call fit_set(2,0.,0.,0.,0.) ! flat background
|
|
c range = scan with:
|
|
range=xx(np)-xx(1)
|
|
cjs type *,'******',output(5),output(6)
|
|
if (ntyp .ne. 6) then
|
|
if (output(6) .lt. 5.0*range/np) then
|
|
c set position to center:
|
|
call fit_set(3, (xx(np)+xx(1))/2.0, range/np,
|
|
& xx(1)+range/6,xx(np)-range/6)
|
|
c set fwhm
|
|
fwhm = max(0.2,range/5)
|
|
c arbitrary limits: minimal 3 steps, maximal 1 /3 range
|
|
call fit_set(6,fwhm,fwhm/10.,3*range/np, range/3)
|
|
else
|
|
fwhm=max(output(6),5*range/np)
|
|
call fit_set(6,fwhm,range/np
|
|
1 ,3*range/np,range/3)
|
|
endif
|
|
endif
|
|
call fit_fit(0)
|
|
c
|
|
if (plot) call fit_plot(hard) ! hardcopy / terminal
|
|
c
|
|
call fit_get_array('P', output, 20, ii)
|
|
call fit_get_array('E', output2, 20, ii)
|
|
cjs type *,' output ',output
|
|
cjs type *,' output2',output2
|
|
c
|
|
c
|
|
c get results now:
|
|
if (ntyp.eq.0) then
|
|
c function 1 (gauss: use only background, intensity=sum)
|
|
omega_fit = output (3)
|
|
omega_err = output2(3)
|
|
off = ' '
|
|
if (iexp) then
|
|
yint = output(8)
|
|
dint = output2(8)
|
|
else
|
|
c modified 14.7.99 sj33
|
|
yint = output(5)
|
|
dint = output2(5)
|
|
endif
|
|
bgr=output(1)
|
|
fwhm=output(6)
|
|
else
|
|
c function 6 (strange)
|
|
yint = output(7)
|
|
dint = output2(7)
|
|
omega_fit = output(3)
|
|
omega_err = output2(3)
|
|
bgr=output(1)
|
|
fwhm=output(4)
|
|
end if
|
|
c
|
|
c make notes for easy data inspections:
|
|
delta = abs (omega_cal - omega_fit)
|
|
if (delta.ge.0.3) off =' *'
|
|
if (delta.ge.0.5) off ='**'
|
|
weak = ' '
|
|
c
|
|
c write on output files now:
|
|
lor_cor = 1.00000
|
|
c L = 1 / SIN (2-Theta), page 156, Schwarzenbach, EPFL:
|
|
c L applies to the calculated value of the intensities
|
|
c Ical*L = Iobs
|
|
c (measured intensity is always bigger than corrected one)
|
|
if (lorentz.eq.0) lor_cor = sin (abs(th2)*radian)
|
|
c
|
|
c June 2010: SJ33
|
|
c tilt angle nu is read-in as chi from file *.col
|
|
c correction formula: not used (3.59) on page 154, Schwarzenbach/Chapuis Cristallographie, Presses EPFL ! SJ33, June 2010
|
|
c used: Arndt and Wills, L=sin(gamma)*cos(nu), where gamma projection 2-theta, nu tilt
|
|
yy1=lor_cor
|
|
if (lorentz.eq.1) lor_cor =
|
|
* sin(abs(th2)*radian)*cos(abs(chi)*radian)
|
|
c * sqrt( (sin(abs(th2)*radian))**2 -(sin(abs(chi)*radian))**2 ) ! SJ33, June 2010
|
|
c write (6,*)' chi=',chi,th2,radian,sin(chi*radian),sin(th2*radian)
|
|
c write (6,*)' Lorentz=',yy1,lor_cor
|
|
c
|
|
c
|
|
if (abs(nint(h)-h) .gt. abs(h)*1e-4 .or.
|
|
& abs(nint(k)-k) .gt. abs(k)*1e-4 .or.
|
|
& abs(nint(l)-l) .gt. abs(l)*1e-4) then
|
|
write(hkl, '(3(x,f9.4))') h, k, l
|
|
write(hkl_float,'(3f6.2)') h, k, l
|
|
iincom=iincom+1
|
|
iall=iall+1
|
|
comens=1 ! SJ33-May 2012
|
|
else
|
|
write(hkl, '(3(x,i9))') nint(h), nint(k), nint(l)
|
|
write(hkl_int,'(3i4)') nint(h), nint(k), nint(l)
|
|
iall=iall+1
|
|
icom=icom+1 ! SJ33-May 2012
|
|
comens=0
|
|
endif
|
|
c
|
|
write (1, 1005) numor,hkl,lor_cor*yint,lor_cor*dint,
|
|
& Th2/2.,omega_fit,chi,phi,temperature,yoverall/1000.,off,weak
|
|
& ,bgr,fwhm,lor_cor,omega_err,datime(1:ldt)
|
|
c
|
|
if (comens.eq.0) then
|
|
write (3, 1025) numor,hkl_int, lor_cor*yint,lor_cor*dint,
|
|
* Th2/2,omega_fit,chi,phi ! SJ33-May 2012
|
|
else
|
|
write (4, 1026) numor,hkl_float,lor_cor*yint,lor_cor*dint,
|
|
* Th2/2,omega_fit,chi,phi
|
|
endif ! SJ33-May 2012
|
|
1025 format (I6,A12,2F10.2,4f8.2) ! SJ33, nuc
|
|
1026 format (I4,A18,2F10.2,4F8.2) ! SJ33 mag
|
|
1005 format (I6,a,1x,f9.2,1x,f9.3,f7.2,1x,f8.3,1x,f7.2,1x,f7.2,1x,
|
|
& F8.3,1x,F6.0,1x,2A2,1x,f6.1,1x,f5.3,1x,f6.4,1x,f6.3,1x,a)
|
|
|
|
if (dint.ge.1) then
|
|
s = yint/dint
|
|
if (s.le.3) weak='w ' ! weak
|
|
if (s.le.1) weak='vw' ! very weak
|
|
if (s.ge.15 .and. yint.ge.100) then ! write strong hkl to rafin.dat
|
|
write (2,1009) hkl,Th2,omega_fit,chi,phi,.2,temperature
|
|
icenter=icenter+1
|
|
end if
|
|
c
|
|
1009 format (a,f8.2,f8.3,2f8.2,1x,f5.3,f9.3)
|
|
end if
|
|
C
|
|
goto 1500
|
|
c -------------------------------------------------------- end loop --
|
|
C
|
|
2222 write(*,1006)
|
|
1006 format (' Bye bye, normal end of program.')
|
|
C
|
|
C -------------------------------------------
|
|
C error exits here, does not pass label 2222:
|
|
3000 continue
|
|
close (1)
|
|
close (2)
|
|
close (3) ! SJ33-May 2012
|
|
close (4) ! SJ33-May 2012
|
|
c
|
|
c writing the end message for the user
|
|
write (6,2040)
|
|
write (6,2042) input_files(1:60)
|
|
write (6,2041) iall,output_file1(1:30),icenter,output_file2(1:30),
|
|
* icom,output_file3(1:30),iincom,output_file4
|
|
2040 format(///,' trics_ccl prepared the following files for YOU:',//)
|
|
2042 format(' TriCS-input-file used:',A60,//)
|
|
2041 format(i6,' reflections ',a30,' for JANA2006',/
|
|
* ,i6,' positions ',a30,' for rafin (to improve your UB using the dataset)',/
|
|
* ,I6,' reflections ',a30,' for Fullprof - commensurate data',/
|
|
* ,I6,' reflections ',a30,' for Fullprof - incommensurate data',//)
|
|
c
|
|
call fit_exit
|
|
c
|
|
end
|